#Distribuciķ Normal

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

set.seed(1000)

#Fiquem les dades

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


#Creem els vectors necessaris de la Mitjana Mostral
e_mu0.hat<-vector()
b_mu0.hat<-vector()
v_mu0.hat<-vector()
eqm_mu0.hat<-vector()
e2_mu0.hat<-vector()
b2_mu0.hat<-vector()
v2_mu0.hat<-vector()
eqm2_mu0.hat<-vector()


#Creem els vectors necessaris del primer estimador
e_mu1.hat<-vector()
b_mu1.hat<-vector()
v_mu1.hat<-vector()
eqm_mu1.hat<-vector()
e2_mu1.hat<-vector()
b2_mu1.hat<-vector()
v2_mu1.hat<-vector()
eqm2_mu1.hat<-vector()

#Creem els vectors necessaris del segon estimador
e_mu2.hat<-vector()
b_mu2.hat<-vector()
v_mu2.hat<-vector()
eqm_mu2.hat<-vector()
e2_mu2.hat<-vector()
b2_mu2.hat<-vector()
v2_mu2.hat<-vector()
eqm2_mu2.hat<-vector()


for (i in 1:length(nvars)){
	matriu.x<-matrix(rnorm(rep*nvars[i],mu,sigma),rep)
#Mitjana mostral
	mu0.hat<-rowMeans(matriu.x)
	e_mu0.hat[i]<-mean(mu0.hat)
	v_mu0.hat[i]<-var(mu0.hat)
	b_mu0.hat[i]<-mean(mu0.hat)-mu
	eqm_mu0.hat[i]<-b_mu0.hat[i]^2+v_mu0.hat[i]
	e2_mu0.hat<-rbind(e2_mu0.hat,e_mu0.hat[i])
	b2_mu0.hat<-rbind(b2_mu0.hat,b_mu0.hat[i])
	v2_mu0.hat<-rbind(v2_mu0.hat,v_mu0.hat[i])
	eqm2_mu0.hat<-rbind(eqm2_mu0.hat,eqm_mu0.hat[i])
#Primer estadístic
	mu1.hat<-rowSums(matriu.x)/(nvars[i]-1)
	e_mu1.hat[i]<-mean(mu1.hat)
	v_mu1.hat[i]<-var(mu1.hat)
	b_mu1.hat[i]<-mean(mu1.hat)-mu
	eqm_mu1.hat[i]<-b_mu1.hat[i]^2+v_mu1.hat[i]
	e2_mu1.hat<-rbind(e2_mu1.hat,e_mu1.hat[i])
	b2_mu1.hat<-rbind(b2_mu1.hat,b_mu1.hat[i])
	v2_mu1.hat<-rbind(v2_mu1.hat,v_mu1.hat[i])
	eqm2_mu1.hat<-rbind(eqm2_mu1.hat,eqm_mu1.hat[i])
#Segon estadístic
	mu2.hat<-rowSums(matriu.x)/(nvars[i]+1)
	e_mu2.hat[i]<-mean(mu2.hat)
	v_mu2.hat[i]<-var(mu2.hat)
	b_mu2.hat[i]<-mean(mu2.hat)-mu
	eqm_mu2.hat[i]<-b_mu2.hat[i]^2+v_mu2.hat[i]
	e2_mu2.hat<-rbind(e2_mu2.hat,e_mu2.hat[i])
	b2_mu2.hat<-rbind(b2_mu2.hat,b_mu2.hat[i])
	v2_mu2.hat<-rbind(v2_mu2.hat,v_mu2.hat[i])
	eqm2_mu2.hat<-rbind(eqm2_mu2.hat,eqm_mu2.hat[i])
}
res0<-cbind(e2_mu0.hat,b2_mu0.hat,v2_mu0.hat,eqm2_mu0.hat)
res1<-cbind(e2_mu1.hat,b2_mu1.hat,v2_mu1.hat,eqm2_mu1.hat)
res2<-cbind(e2_mu2.hat,b2_mu2.hat,v2_mu2.hat,eqm2_mu2.hat)
nomscol<-c("Esperanįa", "Biaix", "Variāncia", "EQM")
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(2,2))
mu0<-expression(bar(italic(X)))
mu1<-expression(hat(italic(mu[1])))
mu2<-expression(hat(italic(mu[2])))
mmu1<-expression(hat(italic(mu[2]))/bar('X'))
mmu2<-expression(hat(italic(mu[2]))/bar('X'))
mmu3<-expression(hat(italic(mu[1]))/hat(italic(mu[2])))

#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='Esperanįa Matemātica',xlab = "Mida mostral", ylab = "", type="n")
abline(h=mu)
lines(x,e2_mu0.hat, lty=1, lwd=3, col='red')
lines(x,e2_mu1.hat, lty=1, lwd=3, col='blue')
lines(x,e2_mu2.hat, lty=1, lwd=3, col='violet')
legend("topright", c(mu0,mu1,mu2), lty=1, col=c("red","blue","violet"), lwd=3, bty="n")

#Grāfic a la posiciķ 1,2
y<-seq(0,max(v2_mu1.hat,v2_mu2.hat), along=nvars)
plot(x,y, main='Variāncia',xlab = "Mida mostral", ylab = "", type="n")
lines(x,v2_mu0.hat, lty=1, lwd=3, col='red')
lines(x,v2_mu1.hat, lty=1, lwd=3, col='blue')
lines(x,v2_mu2.hat, lty=1, lwd=3, col='violet')
legend("topright", c(mu0,mu1,mu2), lty=1, col=c("red","blue","violet"), lwd=3, bty="n")

#Grāfic a la posiciķ 2,1
y<-seq(0,max(eqm2_mu1.hat, eqm2_mu2.hat), along=nvars)
plot(x,y, main='Error Quadrātic Mitjā',xlab = "Mida mostral", ylab = "", type="n")
lines(x,eqm2_mu0.hat, lty=1, lwd=3, col='red')
lines(x,eqm2_mu1.hat, lty=1, lwd=3, col='blue')
lines(x,eqm2_mu2.hat, lty=1, lwd=3, col='violet')
legend("topright", c(mu0,mu1,mu2), lty=1, col=c("red","blue","violet"), lwd=3, bty="n")

#Grāfic a la posiciķ 2,2
q_eqm1<-eqm2_mu1.hat/eqm2_mu0.hat
q_eqm2<-eqm2_mu2.hat/eqm2_mu0.hat
q_eqm3<-eqm2_mu1.hat/eqm2_mu2.hat
y<-seq(0, max(q_eqm1,q_eqm2), along=nvars)
plot(x,y, main="Quocient d'EQM",xlab = "Mida mostral", ylab = "", type="n")
abline(h=1)
lines(x,q_eqm1, lty=1, lwd=3, col='red')
lines(x,q_eqm2, lty=1, lwd=3, col='blue')
lines(x,q_eqm3, lty=1, lwd=3, col='violet')
legend("topright", c(mmu1,mmu2,mmu3), lty=1, col=c("red","blue","violet"), lwd=3, bty="n")
mtext(side=1, line=-2, cex=0.75, outer=T, adj=1, "Script creat per Jordi Lķpez-Tamayo")