% 95 %-prediksjonintervall for my basert p? 10 observasjoner fra normalfordeling
% med my = 10 og og sigma = 2 (antas ukjent).
% ny observasjon trekkes fra samme populasjon
%
% Kopier alt dette inn i kommandovinduet i MATLAB. Konfidensintervall vises grafisk med
% bl?tt og ny observasjon med r?dt hvis bom og bl?tt ellers. Trykk 'Enter' for ? gjenta. Antall ganger -
% intervallet ikke inneholder ny observasjon angis ?verst.
%
% Trykk 's' 'Enter' for automatisk gjentakelse til 1000 datasett er simulert.

n=10

my=10
sigma=2
alfa=.05
nsim=1000

t=tinv(1-alfa/2,n-1)
kvsd=sqrt(1+1/n)

bom=0
forts=''
clf
plot([my my],[0 nsim],'color','blue')
axis([my-4*sigma my+4*sigma 0 100])
hold on

for i=1:nsim
  x=normrnd(my,sigma,1,n)
  middel=mean(x)
  sd=std(x)
  nedre=middel-t*sd*kvsd
  oevre=middel+t*sd*kvsd
  plot([nedre oevre],[i,i],'Color','blue')
  if i>99
    axis([my-4*sigma my+4*sigma i-99 i+1])
  end
  nyobs=normrnd(my,sigma,1,1)
  
  if nedre>nyobs || oevre<nyobs
    bom=bom+1
    plot(nyobs,i,'*','Color','red')
  else
      plot(nyobs,i,'*','Color','blue')
  end
  title(sprintf('%d bom av %d',bom,i))
  drawnow
  if strcmp(forts,'')
    forts=input('','s')
  end
end