#### SPAZIO CAMPIONARIO E DISEGNI CAMPIONARI ###

library(sampling)

### SPAZIO CAMPIONARIO ##

P=1:5
N=5
n=3
Snb=writesample(n,N) #costruisce lo spazio campionario binario costituito da tutti i 
                     #possibili campioni di numerosità fissa n=3 (funzioni indicatrici)
                     #da popolazione N=5
Snb                  #combinazioni di 5 elementi presi a gruppi di 3 senza
                     #ripetizione delle unità e prescindendo dall'ordine di 
                     #estrazione(coefficiente binomiale) C(5,3)=10
rownames(Snb)<-c("s1","s2","s3","s4","s5","s6","s7","s8","s9","s10")
colnames(Snb)<-c("u1","u2","u3","u4","u5")
Snb

pik=colSums(Snb)/nrow(Snb)   #somma P(Si)|i appartiene al campione Si. P(Si)=1/10
                             #P(Si)=constante e pik=constante
#sum(pik)                    #sum(pik)=n (somma prob.inclusione = n)
pi12=sum(Snb[,1]==1 & Snb[,2]==1)/nrow(Snb) #prob. inclus. 2° ordine delle unità 1 e 2 
pi13=sum(Snb[,1]==1 & Snb[,3]==1)/nrow(Snb) #pikj=costante
pi14=sum(Snb[,1]==1 & Snb[,4]==1)/nrow(Snb) 
pi15=sum(Snb[,1]==1 & Snb[,5]==1)/nrow(Snb)
pi23=sum(Snb[,2]==1 & Snb[,3]==1)/nrow(Snb)
pi24=sum(Snb[,2]==1 & Snb[,4]==1)/nrow(Snb)
pi25=sum(Snb[,2]==1 & Snb[,5]==1)/nrow(Snb)
pi34=sum(Snb[,3]==1 & Snb[,4]==1)/nrow(Snb)
pi35=sum(Snb[,3]==1 & Snb[,5]==1)/nrow(Snb)
pi45=sum(Snb[,4]==1 & Snb[,5]==1)/nrow(Snb)
#il numero di pikj = combinazioni di 5 elementi presi a gruppi di 2 senza ripetizione = 10
pikj<-c(pik[1],pi12,pi13,pi14,pi15,NA,pik[2],pi23,pi24,pi25,NA,NA,pik[3],pi34,pi35,NA,NA,NA,pik[4],pi45,NA,NA,NA,NA,pik[5])
#inserimento arbitrario di NA per non duplicare le informazioni relative a stesse probabilità (pi12=pi21)
pikjmatrix<-matrix(pikj,5,5,byrow=T)
colnames(pikjmatrix)<-c("u1","u2","u3","u4","u5")
rownames(pikjmatrix)<-c("u1","u2","u3","u4","u5")
pikjmatrix          #vedi esercitazione 5 per il calcolo con il ciclo for                                 

s1<-P[Snb[1,]==1]
s1
Sn<-rbind(P[Snb[1,]==1],P[Snb[2,]==1],P[Snb[3,]==1],P[Snb[4,]==1],P[Snb[5,]==1],
P[Snb[6,]==1],P[Snb[7,]==1],P[Snb[8,]==1],P[Snb[9,]==1],P[Snb[10,]==1])
rownames(Sn)<-c("s1","s2","s3","s4","s5","s6","s7","s8","s9","s10")
colnames(Sn)<-c("uc1","uc2","uc3")
Sn

### SIMPLE RANDOM SAMPLING and PROPORTIONAL PROBABILITY SAMPLING ##

# estrazione di un campione di numerosità fissa n senza reinserimento 
# srswor(n,N)           #Simple Random Sampling without replacement (SRSWOR) 
# esempio 1

P=1:10                  #popolazione
N=10
n=3
srs<-srswor(n,N)        #estrazione campione casuale a probabilità costanti
srs                     #vettore 0,1 per ciascuna unità 
                        #1 se unità estratta
s_srs<-P[srs==1]        #valori campionari srs
s_srs

pik<-inclusionprobabilities(P,n)  #prob. inclus. del 1° ordine delle unità k=1,...N
                                  #probability proportional-to-size!! (PPS)
pik                               #proporzionali alla "dimensione" 
plot(pik,xlab="population size")  #dei valori di P
#sum(pik)                         #sum(pik)=n

pps<-UPtille(pik)                 #Tillé Sampling (Proportional Probability Sampling - PPS)
s_pps<-P[pps == 1]                #valori campionari pps
s_pps

# esempio 2 (dati reali)

data(swissmunicipalities)  #2896 osservazioni comunali svizzere nel 2003
#swissmunicipalities       #data.frame

n=200
SRS<-srswor(n,nrow(swissmunicipalities))
SRS                                        
s_SRS<-swissmunicipalities[SRS == 1,] #campione srs
s_SRS                   
nrow(s_SRS)

pik<-inclusionprobabilities(swissmunicipalities$Pop020,n) #prob. incl. proporzionali a PoP020
pik
plot(pik,xlab="population size",pch=".")

swissmunicipalities$Pop020
plot(swissmunicipalities$Pop020,xlab="population size",pch=".")

PPS<-UPtille(pik)              
s_PPS<-swissmunicipalities[PPS == 1,]  #campione pps
s_PPS
nrow(s_PPS)

### STRATIFIED SAMPLING ##

Nreg<-table(swissmunicipalities[,2])   #2°colonna=codici regionali "REG"
                                       #table() associa una distribuzione 
                                       #di frequanza a ciascun codice
#attach(swissmunicipalities)           #in alternativa!
#table(REG)

sum(Nreg)                              #N=2896
n=200
numpop_s<-round((Nreg/sum(Nreg))*n)    #proporzione della popolazione riportata al campione
                                       #round arrotonda per difetto con n>200
#sum(numpop_s)                         #n=201 (arrotondato)
STRAT<-strata(swissmunicipalities,stratanames="REG",size=as.vector(numpop_s),method="srswor") #stratificato+srswor
s_STRAT<-getdata(swissmunicipalities,STRAT) #campione stratificato+srswor
nrow(s_STRAT)                          #201

save(s_SRS,s_PPS,s_STRAT,file="campioni.rdata")