#Jordi López-Tamayo, Gener 2011
#Distribució de la relación entre la quasi-variància i l variància poblacional
#a mesura que incrementa la mida mostral. P(0.75<=S2<=1.25)

sigma<-1
IInf<-0.75
ISup<-1.25

# Generació del gràfic
par(mfrow=c(2,2), oma=c(0, 0, 4, 0))
eix <- c(0, 100)

nobs<-10
IInf2<-(nobs-1)*IInf/sigma^2
ISup2<-(nobs-1)*ISup/sigma^2
x <- seq(eix[1], eix[2], length = 1001)
y <- dchisq(x, nobs-1)
plot(x, y, type='l', xlab = expression(italic("n=10")), ylab = expression(italic("Densitat")))
abline(h = 0)
x1 <- x[x >= IInf2 &  x <= ISup2] ; y1 <- dchisq(x1, nobs-1)
x2 <- c(x1[1],x1, x1[length(x1)], x1[1])
y2 <- c(0, y1, 0, 0)
polygon(x2, y2, col = 'yellow')
#Afegim Text
pr<-pchisq(ISup2,nobs-1)-pchisq(IInf2,nobs-1)
Text_I <-paste("[",IInf2,",",ISup2,"]")
text(80, 0.85*dchisq(IInf2, nobs-1), Text_I)
text(80, 0.75*dchisq(IInf2, nobs-1), round(pr, digits=3))

nobs<-20
IInf2<-(nobs-1)*IInf/sigma^2
ISup2<-(nobs-1)*ISup/sigma^2
x <- seq(eix[1], eix[2], length = 1001)
y <- dchisq(x, nobs-1)
plot(x, y, type='l', xlab = expression(italic("n=20")), ylab = expression(italic("Densitat")))
abline(h = 0)
x1 <- x[x >= IInf2 &  x <= ISup2] ; y1 <- dchisq(x1, nobs-1)
x2 <- c(x1[1],x1, x1[length(x1)], x1[1])
y2 <- c(0, y1, 0, 0)
polygon(x2, y2, col = 'yellow')
#Afegim Text
pr<-pchisq(ISup2,nobs-1)-pchisq(IInf2,nobs-1)
Text_I <-paste("[",IInf2,",",ISup2,"]")
text(80, 0.85*dchisq(IInf2, nobs-1), Text_I)
text(80, 0.75*dchisq(IInf2, nobs-1), round(pr, digits=3))

nobs<-40
IInf2<-(nobs-1)*IInf/sigma^2
ISup2<-(nobs-1)*ISup/sigma^2
x <- seq(eix[1], eix[2], length = 1001)
y <- dchisq(x, nobs-1)
plot(x, y, type='l', xlab = expression(italic("n=40")), ylab = expression(italic("Densitat")))
abline(h = 0)
x1 <- x[x >= IInf2 &  x <= ISup2] ; y1 <- dchisq(x1, nobs-1)
x2 <- c(x1[1],x1, x1[length(x1)], x1[1])
y2 <- c(0, y1, 0, 0)
polygon(x2, y2, col = 'yellow')
#Afegim Text
pr<-pchisq(ISup2,nobs-1)-pchisq(IInf2,nobs-1)
Text_I <-paste("[",IInf2,",",ISup2,"]")
text(80, 0.85*dchisq(IInf2, nobs-1), Text_I)
text(80, 0.65*dchisq(IInf2, nobs-1), round(pr, digits=3))

nobs<-60
IInf2<-(nobs-1)*IInf/sigma^2
ISup2<-(nobs-1)*ISup/sigma^2
x <- seq(eix[1], eix[2], length = 1001)
y <- dchisq(x, nobs-1)
plot(x, y, type='l', xlab = expression(italic("n=60")), ylab = expression(italic("Densitat")))
abline(h = 0)
x1 <- x[x >= IInf2 &  x <= ISup2] ; y1 <- dchisq(x1, nobs-1)
x2 <- c(x1[1],x1, x1[length(x1)], x1[1])
y2 <- c(0, y1, 0, 0)
polygon(x2, y2, col = 'yellow')
#Afegim Text
pr<-pchisq(ISup2,nobs-1)-pchisq(IInf2,nobs-1)
Text_I <-paste("[",IInf2,",",ISup2,"]")
text(20, 1.25*dchisq(IInf2, nobs-1), Text_I)
text(20, dchisq(IInf2, nobs-1), round(pr, digits=3))

mtext(side=3, line=0, cex=1.5, outer=T,"Relació entre la quasi-variància i la mida mostral. P(0.75<=S2<=1.25)")
mtext(side=1, line=-2, cex=0.75, outer=T, adj=1, "Script creat per Jordi López-Tamayo")