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

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

set.seed(1000)

#Fiquem les dades
p<-0.7
n<-1
rep<-5000
mu<-p
sigma<-sqrt(p*(1-p))
alfa<-0.01

nvars=30
matriu.x<-matrix(rbinom(rep*nvars,n,p),rep)
	mu0.hat.30<-rowMeans(matriu.x)
	II_mu0.hat.30<-mu0.hat.30+qnorm(alfa/2)*sigma/sqrt(nvars)
	IS_mu0.hat.30<-mu0.hat.30+qnorm(1-alfa/2)*sigma/sqrt(nvars)
x<-vector()
for (i in 1:rep){
	if (II_mu0.hat.30[i] < mu && mu < IS_mu0.hat.30[i]) {
		x[i]<-1
		} else {
		x[i]<-0
		}
}
	encerts.30<-round(sum(x)/rep,2)*100
nvars=60
matriu.x<-matrix(rbinom(rep*nvars,n,p),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(rbinom(rep*nvars,n,p),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(rbinom(rep*nvars,n,p),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
nvars=120
matriu.x<-matrix(rbinom(rep*nvars,n,p),rep)
	mu0.hat.120<-rowMeans(matriu.x)
	II_mu0.hat.120<-mu0.hat.120+qnorm(alfa/2)*sigma/sqrt(nvars)
	IS_mu0.hat.120<-mu0.hat.120+qnorm(1-alfa/2)*sigma/sqrt(nvars)
x<-vector()
for (i in 1:rep){
	if (II_mu0.hat.120[i] < mu && mu < IS_mu0.hat.120[i]) {
		x[i]<-1
		} else {
		x[i]<-0
		}
}
	encerts.120<-round(sum(x)/rep,2)*100
nvars=160
matriu.x<-matrix(rbinom(rep*nvars,n,p),rep)
	mu0.hat.160<-rowMeans(matriu.x)
	II_mu0.hat.160<-mu0.hat.160+qnorm(alfa/2)*sigma/sqrt(nvars)
	IS_mu0.hat.160<-mu0.hat.160+qnorm(1-alfa/2)*sigma/sqrt(nvars)
x<-vector()
for (i in 1:rep){
	if (II_mu0.hat.160[i] < mu && mu < IS_mu0.hat.160[i]) {
		x[i]<-1
		} else {
		x[i]<-0
		}
}
	encerts.160<-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
#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=30',xlab = "Rèpliques", ylab = "", type="n")
abline(h=mu)
lines(x,mu0.hat.30, lty=1, lwd=2, col='red')
lines(x,II_mu0.hat.30, lty=1, lwd=2, col='blue')
lines(x,IS_mu0.hat.30, 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.30,"%"), 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')

#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=120',xlab = "Rèpliques", ylab = "", type="n")
abline(h=mu)
lines(x,mu0.hat.120, lty=1, lwd=2, col='red')
lines(x,II_mu0.hat.120, lty=1, lwd=2, col='blue')
lines(x,IS_mu0.hat.120, 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.120,"%"), 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=160',xlab = "Rèpliques", ylab = "", type="n")
abline(h=mu)
lines(x,mu0.hat.160, lty=1, lwd=2, col='red')
lines(x,II_mu0.hat.160, lty=1, lwd=2, col='blue')
lines(x,IS_mu0.hat.160, 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.160,"%"), cex=1.4, col='brown')



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