#############################################################################
###
###  labIV.txt - Mauro Gasparini, Dicembre 2002
###              quarto laboratorio di R per studenti di ingegneria
###
###              Simulazione di un campione casuale.
###              QQplots.
###              Simulazione da una normale bivariata.
###              Leggi dei grandi numeri e stima.
###              Approssimazione di distribuzione campionaria
###                        di una statistica tramite simulazione.
###
#############################################################################


##### Simulazione di campioni casuali
### Ricordiamo che in R la simulazione di un campione casuale
### di dimensione n si ottiene con il comando r*
### dove al posto di * occorre mettere il nome di una distribuzione

rnorm(10) # simulazione di 10 osservazioni normali standard
rnorm(10,-3,2) # simulazione di 10 osservazioni N(-3,4)

rbinom(1, 1, 0.5)  # una Bernoulli(p) è una binomiale di parametri n=1 e p
rbinom(10, 1, 0.5) # simulazione di 10 Bernoulli p=0.5

##### QQplot 
### Il plot quantile-quantile di un vettore di osservazioni x si fa con 
### la funzione qqnorm(x). 
### Più rettilineo è il qqplot, più crediamo alla normalità dei dati.

# Esempio
qqnorm(rnorm(1000))
# Controesempio
qqnorm(rexp(1000))

##### Simulazione da una normale bivariata
### R non ha una funzione nativa per simulare da una normale bivariata;
### per tale fine si può usare la seguente funzione

rbivnorm <- function(n,       # dimensione del campione voluta
                     mux,     # valore atteso di X
                     muy,     # valore atteso di Y
                     sigmax,  # deviazione std di X
                     sigmay,  # deviazione std di Y
                     rho){    # coefficiente di correlazione lineare
            # genera x dalla marginale di X
            x <- rnorm(n,mux,sigmax) 
            # ora genera Y dalla condizionata di Y dato X=x
            y <- rnorm(n,muy+rho*sigmay*(x-mux)/sigmax, sigmay*sqrt(1-rho^2))
            # ritorna x,y
            cbind(x,y) 
            }

# parametro grafico per fare 3x2 grafici nella stessa pagina
# con un margine esterno (Outer MArgin) sul lato 3 (in alto)
par(mfrow=c(3,2),oma=c(0,0,2,0)) 

plot(rbivnorm(1000,0,0,1,1,rho=-0.8),xlim=c(-4,4),ylim=c(-5,4))
legend(-2,-3,expression(paste(mu[x]==0,", ",mu[y]==0,", ",sigma[x]==1,", ",sigma[y]==1,", ", rho==-0.8)), text.width=5)
plot(rbivnorm(1000,0,0,1,1,rho=-0.5),xlim=c(-4,4),ylim=c(-5,4))
legend(-2,-3,expression(paste(mu[x]==0,", ",mu[y]==0,", ",sigma[x]==1,", ",sigma[y]==1,", ", rho==-0.50)), text.width=5)
plot(rbivnorm(1000,0,0,1,1,rho=0),xlim=c(-4,4),ylim=c(-5,4))
legend(-2,-3,expression(paste(mu[x]==0,", ",mu[y]==0,", ",sigma[x]==1,", ",sigma[y]==1,", ", rho==0.00)), text.width=5)
plot(rbivnorm(1000,0,0,1,1,rho=0.99),xlim=c(-4,4),ylim=c(-5,4))
legend(-2,-3,expression(paste(mu[x]==0,", ",mu[y]==0,", ",sigma[x]==1,", ",sigma[y]==1,", ", rho==0.99)), text.width=5)
plot(rbivnorm(1000,0,1,2,1,rho=0.00),xlim=c(-4,4),ylim=c(-5,4))
legend(-2,-3,expression(paste(mu[x]==0,", ",mu[y]==1,", ",sigma[x]==2,", ",sigma[y]==1,", ", rho==0.00)),text.width=5)
plot(rbivnorm(1000,0,1,2,1,rho=0.99),xlim=c(-4,4),ylim=c(-5,4))
legend(-2,-3,expression(paste(mu[x]==0,", ",mu[y]==1,", ",sigma[x]==2,", ",sigma[y]==1,", ", rho==0.99)), text.width=5)
mtext("Sei simulazioni da normali bivariate",3,outer=T)
par(mfrow=c(1,1)) #rimetti un grafico per pagina

##### Leggi dei grandi numeri e stima in pratica

mean(rnorm(10,-3,2))
mean(rnorm(100,-3,2))
mean(rnorm(1000,-3,2))
mean(rnorm(10000,-3,2))      # cosa si nota in questa progressione?

sd(rnorm(10,-3,2))
sd(rnorm(100,-3,2))
sd(rnorm(1000,-3,2))
sd(rnorm(10000,-3,2))        # e in questa?

par(mfrow=c(2,2)) # 2x2 grafici nella stessa pagina
hist(rnorm(10,-3,2),xlim=c(-9,3),main="Istogramma da 10 osservazioni")
hist(rnorm(100,-3,2),xlim=c(-9,3),main="Istogramma da 100 osservazioni")
hist(rnorm(1000,-3,2),xlim=c(-9,3),main="Istogramma da 1000 osservazioni")
hist(rnorm(10000,-3,2),xlim=c(-9,3),main="Istogramma da 10000 osservazioni")      

##### Simulazione di distribuzioni campionarie in R

xbar <- matrix(0,1000,4) ### inizializzazione di xbar

### simulazione di quattro diverse (n=10,100,1000,10000) 
### distribuzioni campionarie della media campionaria
### nelle quattro colonne di xbar
for (i in 1:1000) {       
xbar[i,1] <- mean( rnorm(10,-3,2) )           
xbar[i,2] <- mean( rnorm(100,-3,2) )
xbar[i,3] <- mean( rnorm(1000,-3,2) )
xbar[i,4] <- mean( rnorm(10000,-3,2) )
}

### stima grafica della densità tramite istogramma
hist(xbar[,1],xlim=c(-4,-2),main=expression(paste("Densità stimata di ",bar(X)," con 10 osservazioni")))
hist(xbar[,2],xlim=c(-4,-2),main=expression(paste("Densità stimata di ",bar(X)," con 100 osservazioni"))) 
hist(xbar[,3],xlim=c(-4,-2),main=expression(paste("Densità stimata di ",bar(X)," con 1000 osservazioni"))) 
hist(xbar[,4],xlim=c(-4,-2),main=expression(paste("Densità stimata di ",bar(X)," con 10000 osservazioni"))) 

par(mfrow=c(1,1)) #rimetti un grafico per pagina

### Esercizio da fare in classe:
### 1) generare un campione di 10, 100, 1000 e 10000 osservazioni normali bivariate
###    di medie 3 e -2, varianze 2 e 1 e covarianza .7.
###    Calcolare media e varianza di ciascuna. 
### 2) Nel contesto dell'esercizio precedente, suggerire una formula per la stima di rho
###    e implementarla in R.
### 3) Scrivere una funzione che simuli la distribuzione campionaria di S^2
###    basato su un campione casuale di dimensione 10 da una esponenziale standard (lambda=1)