#Distribuciķ Normal
#Evoluciķ de l'amplada de l'interval de confianįa per a la Variāncia de la poblaciķ
#a mesura que incrementa n i per diferents alfas.

#Fixem la llavor per a que l'experiment sempre surti el mateix

set.seed(1000)

#Fiquem les dades

rep<-500
nvars<-seq(2,200,1)
mu<-10
sigma<-3


#Creem els vectors necessaris de la Mitjana Mostral
e_var0.hat<-vector()
eII_var0.hat<-vector()
eIS_var0.hat<-vector()
e2_var0.hat<-vector()
eII2_var0.hat<-vector()
eIS2_var0.hat<-vector()

e_var1.hat<-vector()
eII_var1.hat<-vector()
eIS_var1.hat<-vector()
e2_var1.hat<-vector()
eII2_var1.hat<-vector()
eIS2_var1.hat<-vector()

e_var2.hat<-vector()
eII_var2.hat<-vector()
eIS_var2.hat<-vector()
e2_var2.hat<-vector()
eII2_var2.hat<-vector()
eIS2_var2.hat<-vector()

for (i in 1:length(nvars)){
	matriu.x<-matrix(rnorm(rep*nvars[i],mu,sigma),rep)
#Interval alfa
	alfa=0.1
	mu0.hat<-rowMeans(matriu.x)
	matriu.xc<-(matriu.x-mu0.hat)^2
	var0.hat<-rowSums(matriu.xc)/(nvars[i]-1)
	II_var0.hat<-(nvars[i]-1)*var0.hat/qchisq(1-alfa/2,nvars[i]-1)
	IS_var0.hat<-(nvars[i]-1)*var0.hat/qchisq(alfa/2,nvars[i]-1)
	e_var0.hat[i]<-mean(var0.hat)
	eII_var0.hat[i]<-mean(II_var0.hat)
	eIS_var0.hat[i]<-mean(IS_var0.hat)
	e2_var0.hat<-rbind(e2_var0.hat,e_var0.hat[i])
	eII2_var0.hat<-rbind(eII2_var0.hat,eII_var0.hat[i])
	eIS2_var0.hat<-rbind(eIS2_var0.hat,eIS_var0.hat[i])
	alfa=0.05
	mu1.hat<-rowMeans(matriu.x)
	matriu.xc<-(matriu.x-mu1.hat)^2
	var1.hat<-rowSums(matriu.xc)/(nvars[i]-1)
	II_var1.hat<-(nvars[i]-1)*var1.hat/qchisq(1-alfa/2,nvars[i]-1)
	IS_var1.hat<-(nvars[i]-1)*var1.hat/qchisq(alfa/2,nvars[i]-1)
	e_var1.hat[i]<-mean(var1.hat)
	eII_var1.hat[i]<-mean(II_var1.hat)
	eIS_var1.hat[i]<-mean(IS_var1.hat)
	e2_var1.hat<-rbind(e2_var1.hat,e_var1.hat[i])
	eII2_var1.hat<-rbind(eII2_var1.hat,eII_var1.hat[i])
	eIS2_var1.hat<-rbind(eIS2_var1.hat,eIS_var1.hat[i])
	alfa=0.01
	mu2.hat<-rowMeans(matriu.x)
	matriu.xc<-(matriu.x-mu1.hat)^2
	var2.hat<-rowSums(matriu.xc)/(nvars[i]-1)
	II_var2.hat<-(nvars[i]-1)*var2.hat/qchisq(1-alfa/2,nvars[i]-1)
	IS_var2.hat<-(nvars[i]-1)*var2.hat/qchisq(alfa/2,nvars[i]-1)
	e_var2.hat[i]<-mean(var2.hat)
	eII_var2.hat[i]<-mean(II_var2.hat)
	eIS_var2.hat[i]<-mean(IS_var2.hat)
	e2_var2.hat<-rbind(e2_var2.hat,e_var2.hat[i])
	eII2_var2.hat<-rbind(eII2_var2.hat,eII_var2.hat[i])
	eIS2_var2.hat<-rbind(eIS2_var2.hat,eIS_var2.hat[i])
}
res0<-cbind(e2_var0.hat,eII2_var0.hat,eIS2_var0.hat)
res1<-cbind(e2_var1.hat,eII2_var1.hat,eIS2_var1.hat)
res2<-cbind(e2_var2.hat,eII2_var2.hat,eIS2_var2.hat)

nomscol<-c("Esperanįa", "I. Inferior", "I. Superior")
colnames(res0)<-nomscol
colnames(res1)<-nomscol
colnames(res2)<-nomscol

res0; res1; res2; 

#Anālisi grāfica.
#Generem una matriu grāfica d'una fila i dues columnes
par(mfrow = c(1,3), oma=c(1, 0, 4, 0))

#Grāfic a la posiciķ 1,1
y<-seq(mu-2.5*sigma, mu+2.5*sigma, along=nvars)
x<-nvars
plot(x,y, main='Alfa=0.1',xlab = "Mida mostral", ylab = "", type="n")
abline(h=sigma^2)
lines(x,e2_var0.hat, lty=1, lwd=3, col='red')
lines(x,eII2_var0.hat, lty=1, lwd=3, col='blue')
lines(x,eIS2_var0.hat, lty=1, lwd=3, col='blue')


#Grāfic a la posiciķ 1,2
y<-seq(mu-2.5*sigma, mu+2.5*sigma, along=nvars)
x<-nvars
plot(x,y, main='Alfa=0.05',xlab = "Mida mostral", ylab = "", type="n")
abline(h=sigma^2)
lines(x,e2_var1.hat, lty=1, lwd=3, col='red')
lines(x,eII2_var1.hat, lty=1, lwd=3, col='blue')
lines(x,eIS2_var1.hat, lty=1, lwd=3, col='blue')

#Grāfic a la posiciķ 1,3
y<-seq(mu-2.5*sigma, mu+2.5*sigma, along=nvars)
x<-nvars
plot(x,y, main='Alfa=0.01',xlab = "Mida mostral", ylab = "", type="n")
abline(h=sigma^2)
lines(x,e2_var2.hat, lty=1, lwd=3, col='red')
lines(x,eII2_var2.hat, lty=1, lwd=3, col='blue')
lines(x,eIS2_var2.hat, lty=1, lwd=3, col='blue')
mtext(side=3, line=0, cex=1.5, outer=T,"Interval de confianįa per a la Variāncia Poblacional.")
mtext(side=1, line=-1, cex=0.75, outer=T, adj=1, "Script creat per Jordi Lķpez-Tamayo  ")