## Simulazione: proprietà e distribuzione dello stimatore media campionaria


## creo popolazione di N=30.000, un terzo maschi e 2/3 femmine,
## estraendo da una uniforme. Le femmine sono nella prima parte.
N <- 30000
pop <- c(runif(N*2/3, min=150, max=180), runif(N/3, min=160, max=190))
hist(pop)

## I momenti della popolazione sono:

## media
Y.bar <- sum(pop)/length(pop)
mean(pop)

## varianza
sigma.Y <- sum((pop-mean(pop))^2)/length(pop)
var(pop) ## NB var() calcola la varianza campionaria!
sum((pop-mean(pop))^2)/(length(pop)-1)

## Definiamo il piano di campionamento: CCS senza reimbussolamento

## numero di estrazioni
K <- 1000

## numero di estratti per ciascun campione (prendiamo apposta un campione
## molto piccolo)
n <- 10

## organizzo i campioni in una matrice 10x1000, ogni colonna un campione
campioni <- matrix(ncol=K, nrow=n)
for(k in 1:K) { 
  campioni[ ,k] <- sample(pop, n)
}

## Stimatore "media campionaria"

## come si distribuisce lo stimatore "media campionaria" se n=10?
medie <- apply(campioni, 2, mean)
plot(density(medie), col="red")

## (illustrazione: "peschiamo" tutti i campioni e situiamoli nel grafico)
library(scales)
plot(density(medie), col="white")
for(i in 1:K) abline(v=medie[i], col=alpha("red", 0.1))
lines(density(medie), col="red", lwd=3)

## Valore atteso dello stimatore e correttezza: 
## confronto con la "vera" media della popolazione (in nero)
plot(density(medie), col="red")
abline(v=mean(pop))
abline(v=mean(medie), col="red", lty=3)

## Distribuzione dello stimatore: confronto con la Normale
lines(density(rnorm(1000000, mean=mean(medie), sd=sd(medie))), lty=2,
      col="blue")
qqnorm(medie)


## ES della media campionaria:
## ho bisogno di stimare la varianza campionaria

## Stimatore "varianza campionaria corretta"

varianze <- apply(campioni, 2, var)
plot(density(varianze), col="green3")

## confronto con la "vera" varianza della popolazione (in nero)
abline(v=sum((pop-mean(pop))^2)/length(pop))
abline(v=mean(varianze), col="green3", lty=3)

## ...e la varianza campionaria "non corretta"?
varianze.nc <- apply(campioni, 2,
                     FUN=function(x) sum((x-mean(x))^2)/length(x))
lines(density(varianze.nc), col="orange")
abline(v=mean(varianze.nc), col="orange", lty=3)


## Calcoliamo gli ES della media campionaria

## Dalla teoria, il "vero" ES è
veroES <- sqrt((1-n/N)*sigma.Y/n)

## l'equivalente stimatore campionario:
ES.y <- function(x, f) sqrt((1-f)*var(x)/length(x)) 

errori.standard <- apply(campioni, 2, ES.y, f=n/N)
es <- mean(errori.standard)

## calcolo il valore critico
z.alfa <- qnorm(0.975)

## test: quanti intervalli di confidenza "contengono" il vero parametro?
test <- 0
for(k in 1:K) {
  media.k <- mean(campioni[ ,k])
  ES.y.k <- ES.y(campioni[ ,k], f=n/N)
  ci.l <- media.k - z.alfa*ES.y.k
  ci.u <- media.k + z.alfa*ES.y.k
  if((Y.bar > ci.l)&(Y.bar < ci.u)) {
    test <- test + 1
  }
}
## percentuale di successo:
test/K*100

## Problema: i nostri campioni sono piccoli (n=10): meglio usare
## i valori critici della t di Student! (e anche così ci manca la
## normalità della variabile di partenza). Vediamo che succede...

t.alfa <- qt(0.975, df=n-1)

test <- 0
for(k in 1:K) {
  media.k <- mean(campioni[ ,k])
  ES.y.k <- ES.y(campioni[ ,k], f=n/N)
  ci.l <- media.k - t.alfa*ES.y.k
  ci.u <- media.k + t.alfa*ES.y.k
  if((Y.bar > ci.l)&(Y.bar < ci.u)) {
    test <- test + 1
  }
}
## percentuale di successo:
test/K*100


## Cosa succede all'aumentare del campione?

## numero di estratti per ciascun campione ora è 100
## organizzo i campioni in una matrice 100x1000, ogni colonna un campione
campioni100 <- matrix(ncol=K, nrow=100)
for(k in 1:K) { 
  campioni100[ ,k] <- sample(pop, 100)
}

## come si distribuisce lo stimatore "media campionaria" se n=10?
## e se n=100?
medie100 <- apply(campioni100, 2, mean)
plot(density(medie), col="red", ylim=c(0, 0.8))
lines(density(medie100), col="blue")

## numero di estratti per ciascun campione ora è 300
campioni300 <- matrix(ncol=K, nrow=300)
for(k in 1:K) { 
  campioni300[ ,k] <- sample(pop, 300)
}
medie300 <- apply(campioni300, 2, mean)
lines(density(medie300), col="green3")

abline(v=mean(pop), lty=3)

#########################################

## Campionamento stratificato:

## dalla conoscenza della popolazione, è
N.f <- 20000
N.m <- 10000

W.f <- N.f/N
W.m <- N.m/N

## Definiamo il piano di campionamento: CCS senza reimbussolamento,
## ma adesso stratificato per femmine e maschi

## numero di estrazioni
K <- 1000

## numero di estratti per ciascun campione e strato: la stratificazione
## è invertita (i maschi sono sovrarappresentati) rispetto alle reali
## proporzioni nella popolazione
n.f <- 4
n.m <- 8

## organizzo i campioni in una matrice 10x1000, ogni colonna un campione
campioni.w <- matrix(ncol=K, nrow=(n.f + n.m))
for(k in 1:K) { 
  campioni.w[ ,k] <- c(sample(pop[1:N.f], n.f), sample(pop[(N.f+1):N], n.m))
}


## La media campionaria pesata (correttamente) per gli strati *della
## popolazione*, NB non per le frazioni di campionamento
mean.w <- function(x, W.f, W.m) {
   mean.w <- W.f * mean(x[1:n.f]) + W.m * mean(x[(n.f+1):(n.f+n.m)])
   return(mean.w)
}

## come si distribuisce lo stimatore "media campionaria pesata"?

medie.w <- apply(campioni.w, 2, FUN=mean.w, W.f=W.f, W.m=W.m)
plot(density(medie.w), col="red")

## Valore atteso dello stimatore e correttezza: 
## confronto con la "vera" media della popolazione (in nero)
abline(v=mean(pop))
abline(v=mean(medie.w), col="red", lty=3)


## Calcoliamo gli ES della media campionaria:

## (Radice della) varianza campionaria corretta pesata per gli strati:
ES.y.w <- function(x, W.f, W.m) {
   ES.y.w <- sqrt(W.f^2 * (1-n.f/N.f) * var(x[1:n.f])/n.f +
                  W.m^2 * (1-n.m/N.m) * var(x[(n.f+1):(n.f+n.m)])/n.m)
   return(ES.y.w)
}

es.w <- apply(campioni.w, 2, FUN=ES.y.w, W.f=W.f, W.m=W.m)
plot(density(es.w), col="green3")
## confronto con la "vera" varianza della popolazione (in nero)
abline(v=mean(es.w), col="green3", lty=3)

## Dalla teoria, il "vero" ES è (servono le varianze di ogni strato
## della popolazione):
pop.f <- pop[1:N.f]
pop.m <- pop[(N.f+1):(N.f+N.m)]
var.f <- sum((pop.f-mean(pop.f))^2)/length(pop.f)
var.m <- sum((pop.m-mean(pop.m))^2)/length(pop.f)
sigma2.w <- W.f^2 * (1-n.f/N.f) * var.f/n.f +
            W.m^2 * (1-n.m/N.m) * var.m/n.m
veroES <- sqrt(sigma2.w)
abline(v=veroES)

## NB se uno avesse usato le proporzioni degli strati *nel campione*:
badES.y.w <- function(x, W.f, W.m) {
   badES.y.w <- sqrt((4/12)^2 * (1-n.f/N.f) * var(x[1:n.f])/n.f +
                  (8/12)^2 * (1-n.m/N.m) * var(x[(n.f+1):(n.f+n.m)])/n.m)
   return(badES.y.w)
}

## confronto con il procedimento corretto:
plot(density(es.w), col="green3", ylim=c(0, 1))
abline(v=mean(es.w), col="green3", lty=3)
abline(v=veroES)

bades.w <- apply(campioni.w, 2, FUN=badES.y.w, W.f=W.f, W.m=W.m)
lines(density(bades.w), col="orange")
abline(v=mean(bades.w), col="orange", lty=3)



#####################################################

## Proporzioni:

## Sappiamo che la proporzione femmine/maschi=2, femmine/totale=2/3
## La variabile da campionare è "femmina", ovvero nella popolazione:
prop <- c(rep(1, N.f), rep(0, N.m))

## media della popolazione:
Tp <- mean(prop)

## varianza della popolazione:
Vp <- Tp*(1-Tp)



## Definiamo il piano di campionamento: CCS senza reimbussolamento

## numero di estrazioni
K <- 1000

## numero di estratti per ciascun campione
n <- 10

## organizzo i campioni in una matrice 10x1000, ogni colonna un campione
campionip <- matrix(ncol=K, nrow=n)
for(k in 1:K) { 
  campionip[ ,k] <- sample(prop, n)
}

## Stimatore "media campionaria della proporzione"

## come si distribuisce lo stimatore "media campionaria" se n=10?
mediep <- apply(campionip, 2, mean)
plot(density(mediep), col="red")

## (illustrazione: "peschiamo" tutti i campioni e situiamoli nel grafico)
library(scales)
plot(density(mediep), col="white")
for(i in 1:K) abline(v=mediep[i], col=alpha("red", 0.1))
lines(density(mediep), col="red", lwd=3)

## Valore atteso dello stimatore e correttezza: 
## confronto con la "vera" media della popolazione (in nero)
plot(density(mediep), col="red")
abline(v=mean(prop))
abline(v=mean(mediep), col="red", lty=3)

## Distribuzione dello stimatore: confronto con la Normale
lines(density(rnorm(1000000, mean=mean(mediep), sd=sd(mediep))), lty=2,
      col="blue")
qqnorm(mediep)

## ES della proporzione campionaria:
## ho bisogno di stimare la varianza della proporzione campionaria

## Stimatore "ES della proporzione campionaria (corretto)"
ES.y.p <- function(x, f) sqrt((1-f)*mean(x)*(1-mean(x))/(length(x)-1)) 

errori.standard.p <- apply(campionip, 2, ES.y.p, f=n/N)
es.p <- mean(errori.standard.p)

## test: quanti intervalli di confidenza "contengono" il vero parametro?
test <- 0
for(k in 1:K) {
  prop.k <- mean(campionip[ ,k])
  ESp.y.k <- ES.y(campionip[ ,k], f=n/N)
  ci.l <- prop.k - z.alfa*ESp.y.k
  ci.u <- prop.k + z.alfa*ESp.y.k
  if((Tp > ci.l)&(Tp < ci.u)) {
    test <- test + 1
  }
}
## percentuale di successo:
test/K*100

## In piccoli campioni, i valori critici della Normale sono inappropriati
## Inoltre, nemmeno la t di Student va bene perché la popolazione di
## selezione non è normale.




#################################################

## allarghiamo il campione
## numero di estratti per ciascun campione
n <- 100

## organizzo i campioni in una matrice 100x1000, ogni colonna un campione
campionip100 <- matrix(ncol=K, nrow=n)
for(k in 1:K) { 
  campionip100[ ,k] <- sample(prop, n)
}

## Stimatore "media campionaria della proporzione"

## come si distribuisce lo stimatore "media campionaria" se n=100?
mediep100 <- apply(campionip100, 2, mean)
plot(density(mediep100), col="red")

## (illustrazione: "peschiamo" tutti i campioni e situiamoli nel grafico)
library(scales)
plot(density(mediep100), col="white")
for(i in 1:K) abline(v=mediep100[i], col=alpha("red", 0.1))
lines(density(mediep100), col="red", lwd=3)

## Valore atteso dello stimatore e correttezza: 
## confronto con la "vera" media della popolazione (in nero)
plot(density(mediep100), col="red")
abline(v=mean(prop))
abline(v=mean(mediep100), col="red", lty=3)

## Distribuzione dello stimatore: confronto con la Normale
lines(density(rnorm(1000000, mean=mean(mediep100),
                    sd=sd(mediep100))), lty=2, col="blue")
qqnorm(mediep)

errori.standard.p100 <- apply(campionip100, 2, ES.y.p, f=n/N)
es.p100 <- mean(errori.standard.p100)

## test: quanti intervalli di confidenza "contengono" il vero parametro?
test <- 0
for(k in 1:K) {
  prop.k <- mean(campionip100[ ,k])
  ESp.y.k <- ES.y(campionip100[ ,k], f=n/N)
  ci.l <- prop.k - z.alfa*ESp.y.k
  ci.u <- prop.k + z.alfa*ESp.y.k
  if((Tp > ci.l)&(Tp < ci.u)) {
    test <- test + 1
  }
}
## percentuale di successo:
test/K*100

## In campioni più grandi, i valori critici della Normale sono appropriati
## (a prescindere dalla distribuzione della popolazione di selezione).