#Distribució Normal
#Càcul del nombre de mostres per als que encertem que 
#la Mitjana poblacional es troba dins de l'interval

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

set.seed(1000)

#Fiquem les dades

rep<-5000
mu<-0
sigma<-1
alfa<-0.05
nvars=10
matriu.x<-matrix(rnorm(rep*nvars,mu,sigma),rep)
	mu0.hat.10<-rowMeans(matriu.x)
	II_mu0.hat.10<-mu0.hat.10+qnorm(alfa/2)*sigma/sqrt(nvars)
	IS_mu0.hat.10<-mu0.hat.10+qnorm(1-alfa/2)*sigma/sqrt(nvars)
x<-vector()
for (i in 1:rep){
	if (II_mu0.hat.10[i] < mu && mu < IS_mu0.hat.10[i]) {
		x[i]<-1
		} else {
		x[i]<-0
		}
}
	encerts.10<-round(sum(x)/rep,2)*100
nvars=20
matriu.x<-matrix(rnorm(rep*nvars,mu,sigma),rep)
	mu0.hat.20<-rowMeans(matriu.x)
	II_mu0.hat.20<-mu0.hat.20+qnorm(alfa/2)*sigma/sqrt(nvars)
	IS_mu0.hat.20<-mu0.hat.20+qnorm(1-alfa/2)*sigma/sqrt(nvars)
x<-vector()
for (i in 1:rep){
	if (II_mu0.hat.20[i] < mu && mu < IS_mu0.hat.20[i]) {
		x[i]<-1
		} else {
		x[i]<-0
		}
}
	encerts.20<-round(sum(x)/rep,2)*100
nvars=40
matriu.x<-matrix(rnorm(rep*nvars,mu,sigma),rep)
	mu0.hat.40<-rowMeans(matriu.x)
	II_mu0.hat.40<-mu0.hat.40+qnorm(alfa/2)*sigma/sqrt(nvars)
	IS_mu0.hat.40<-mu0.hat.40+qnorm(1-alfa/2)*sigma/sqrt(nvars)
x<-vector()
for (i in 1:rep){
	if (II_mu0.hat.40[i] < mu && mu < IS_mu0.hat.40[i]) {
		x[i]<-1
		} else {
		x[i]<-0
		}
}
	encerts.40<-round(sum(x)/rep,2)*100
nvars=60
matriu.x<-matrix(rnorm(rep*nvars,mu,sigma),rep)
	mu0.hat.60<-rowMeans(matriu.x)
	II_mu0.hat.60<-mu0.hat.60+qnorm(alfa/2)*sigma/sqrt(nvars)
	IS_mu0.hat.60<-mu0.hat.60+qnorm(1-alfa/2)*sigma/sqrt(nvars)
x<-vector()
for (i in 1:rep){
	if (II_mu0.hat.60[i] < mu && mu < IS_mu0.hat.60[i]) {
		x[i]<-1
		} else {
		x[i]<-0
		}
}
	encerts.60<-round(sum(x)/rep,2)*100
nvars=80
matriu.x<-matrix(rnorm(rep*nvars,mu,sigma),rep)
	mu0.hat.80<-rowMeans(matriu.x)
	II_mu0.hat.80<-mu0.hat.80+qnorm(alfa/2)*sigma/sqrt(nvars)
	IS_mu0.hat.80<-mu0.hat.80+qnorm(1-alfa/2)*sigma/sqrt(nvars)
x<-vector()
for (i in 1:rep){
	if (II_mu0.hat.80[i] < mu && mu < IS_mu0.hat.80[i]) {
		x[i]<-1
		} else {
		x[i]<-0
		}
}
	encerts.80<-round(sum(x)/rep,2)*100
nvars=100
matriu.x<-matrix(rnorm(rep*nvars,mu,sigma),rep)
	mu0.hat.100<-rowMeans(matriu.x)
	II_mu0.hat.100<-mu0.hat.100+qnorm(alfa/2)*sigma/sqrt(nvars)
	IS_mu0.hat.100<-mu0.hat.100+qnorm(1-alfa/2)*sigma/sqrt(nvars)
x<-vector()
for (i in 1:rep){
	if (II_mu0.hat.100[i] < mu && mu < IS_mu0.hat.100[i]) {
		x[i]<-1
		} else {
		x[i]<-0
		}
}
	encerts.100<-round(sum(x)/rep,2)*100
#Anàlisi gràfica.
#Generem una matriu gràfica d'una fila i dues columnes
par(mfrow = c(2,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=seq(1,rep,1))
x<-seq(1,rep,1)
plot(x,y, main='n=10',xlab = "Rèpliques", ylab = "", type="n")
abline(h=mu)
lines(x,mu0.hat.10, lty=1, lwd=2, col='red')
lines(x,II_mu0.hat.10, lty=1, lwd=2, col='blue')
lines(x,IS_mu0.hat.10, lty=1, lwd=2, col='blue')
text(0.75*rep,0.95*(mu+2.5*sigma), "Encerts", cex=1.4, col='brown')
text(0.95*rep,0.95*(mu+2.5*sigma), paste(encerts.10,"%"), cex=1.4, col='brown')


#Gràfic a la posició 1,2
y<-seq(mu-2.5*sigma, mu+2.5*sigma, along=seq(1,rep,1))
x<-seq(1,rep,1)
plot(x,y, main='n=20',xlab = "Rèpliques", ylab = "", type="n")
abline(h=mu)
lines(x,mu0.hat.20, lty=1, lwd=2, col='red')
lines(x,II_mu0.hat.20, lty=1, lwd=2, col='blue')
lines(x,IS_mu0.hat.20, lty=1, lwd=2, col='blue')
text(0.75*rep,0.95*(mu+2.5*sigma), "Encerts", cex=1.4, col='brown')
text(0.95*rep,0.95*(mu+2.5*sigma), paste(encerts.20,"%"), cex=1.4, col='brown')

#Gràfic a la posició 1,3
y<-seq(mu-2.5*sigma, mu+2.5*sigma, along=seq(1,rep,1))
x<-seq(1,rep,1)
plot(x,y, main='n=40',xlab = "Rèpliques", ylab = "", type="n")
abline(h=mu)
lines(x,mu0.hat.40, lty=1, lwd=2, col='red')
lines(x,II_mu0.hat.40, lty=1, lwd=2, col='blue')
lines(x,IS_mu0.hat.40, lty=1, lwd=2, col='blue')
text(0.75*rep,0.95*(mu+2.5*sigma), "Encerts", cex=1.4, col='brown')
text(0.95*rep,0.95*(mu+2.5*sigma), paste(encerts.40,"%"), cex=1.4, col='brown')


#Gràfic a la posició 1,1
y<-seq(mu-2.5*sigma, mu+2.5*sigma, along=seq(1,rep,1))
x<-seq(1,rep,1)
plot(x,y, main='n=60',xlab = "Rèpliques", ylab = "", type="n")
abline(h=mu)
lines(x,mu0.hat.60, lty=1, lwd=2, col='red')
lines(x,II_mu0.hat.60, lty=1, lwd=2, col='blue')
lines(x,IS_mu0.hat.60, lty=1, lwd=2, col='blue')
text(0.75*rep,0.95*(mu+2.5*sigma), "Encerts", cex=1.4, col='brown')
text(0.95*rep,0.95*(mu+2.5*sigma), paste(encerts.60,"%"), cex=1.4, col='brown')


#Gràfic a la posició 1,1
y<-seq(mu-2.5*sigma, mu+2.5*sigma, along=seq(1,rep,1))
x<-seq(1,rep,1)
plot(x,y, main='n=80',xlab = "Rèpliques", ylab = "", type="n")
abline(h=mu)
lines(x,mu0.hat.80, lty=1, lwd=2, col='red')
lines(x,II_mu0.hat.80, lty=1, lwd=2, col='blue')
lines(x,IS_mu0.hat.80, lty=1, lwd=2, col='blue')
text(0.75*rep,0.95*(mu+2.5*sigma), "Encerts", cex=1.4, col='brown')
text(0.95*rep,0.95*(mu+2.5*sigma), paste(encerts.80,"%"), cex=1.4, col='brown')


#Gràfic a la posició 1,1
y<-seq(mu-2.5*sigma, mu+2.5*sigma, along=seq(1,rep,1))
x<-seq(1,rep,1)
plot(x,y, main='n=100',xlab = "Rèpliques", ylab = "", type="n")
abline(h=mu)
lines(x,mu0.hat.100, lty=1, lwd=2, col='red')
lines(x,II_mu0.hat.100, lty=1, lwd=2, col='blue')
lines(x,IS_mu0.hat.100, lty=1, lwd=2, col='blue')
text(0.75*rep,0.95*(mu+2.5*sigma), "Encerts", cex=1.4, col='brown')
text(0.95*rep,0.95*(mu+2.5*sigma), paste(encerts.100,"%"), cex=1.4, col='brown')

mtext(side=3, line=0, cex=1.5, outer=T,"Interval de confiança per a la Mitjana Poblacional.Percentatge d'encerts")
mtext(side=1, line=-1, cex=0.75, outer=T, adj=1, "Script creat per Jordi López-Tamayo  ")