標本調査法　土屋 隆裕 先生　朝倉書店

```{r}
#標本調査法

#p.30

#非復元抽出の場合

Y=c(576,380,94);pis=c(0.648,0.597,0.569)

tau_y=sum(Y/pis)

pis_ij=matrix(c(0.648,0,0,0.33,0.597,0,0.311,0.278,0.569),ncol=3)

pis_ij=pis_ij+t(pis_ij);diag(pis_ij)=diag(pis_ij)/2

V_HT=array(0,dim=dim(pis_ij));V_YG=array(0,dim=dim(pis_ij))

for(i in 1:nrow(pis_ij)){
for(j in 1:ncol(pis_ij)){

V_HT[i,j]=((pis_ij[i,j]-pis[i]*pis[j])*(Y[i]/pis[i])*(Y[j]/pis[j])*(1/pis_ij[i,j]))

V_YG[i,j]=-(pis_ij[i,j]-pis[i]*pis[j])*((Y[i]/pis[i])^2-2*Y[i]*Y[j]/(pis[i]*pis[j])+(Y[j]/pis[j])^2)*(1/pis_ij[i,j])/2

}}

V_HT=sum(V_HT);V_YG=sum(V_YG)





library(dplyr)

Y=c(576,380,74,292,94);pis=c(0.648,0.597,0.579,0.607,0.569)

X=c(47,31,25,34,22);

tau_y=sum(Y/pis)

data=data.frame(samples=c(1:10),company1=c(1,1,1,1,1,1,2,2,2,3),company2=c(2,2,2,3,3,4,3,3,4,4),company3=c(3,4,5,4,5,5,4,5,5,4),pattern_a=rep(0.1,10),pattern_b=c(0.108,0.117,0.105,0.111,0.099,0.108,0.094,0.082,0.091,0.085))

comp_no=unique(sort(as.matrix(data[,colnames(data) %in% c("company1","company2","company3")])))

pattern_a_mat=array(0.3,dim=c(length(comp_no),length(comp_no)))

diag(pattern_a_mat)=0.6

pattern_b_mat=array(0,dim=c(length(comp_no),length(comp_no)))

pattern_b_mat[1,]=c(0.648,0.33,0.318,0.336,0.311)

pattern_b_mat[2,]=c(0,0.597,0.284,0.303,0.278)

pattern_b_mat[3,]=c(0,0,0.579,0.29,0.265)

pattern_b_mat[4,]=c(0,0,0,0.607,0.284)

pattern_b_mat[5,]=c(0,0,0,0,0.569)

pattern_b_mat=pattern_b_mat+t(pattern_b_mat)

diag(pattern_b_mat)=diag(pattern_b_mat)/2


V_tau_y=c();V_HT=c()

for(i in 1:length(comp_no)){
for(j in 1:length(comp_no)){  
 
V_tau_y=c(V_tau_y,(pattern_b_mat[i,j]-pis[i]*pis[j])*(Y[i]/pis[i])*(Y[j]/pis[j]))   

}
}

V_tau_y=sum(V_tau_y)


#Hansen-Hurwitz estimation 

#復元抽出の場合

p=c(0.296,0.195,0.157,0.214,0.138)

Y=c(576,380,74,292,94)

X=c(47,31,25,34,22);

n=3

pis=1-(1-p)^n

p_sub=c(0.296,0.296,0.157)

Y_sub=c(576,576,74)

V_tau_y=sum(p*(Y/p-sum(Y))^2)/n


tau_y=sum(Y_sub/p_sub)/n

V_tau_y_hat=sum((Y_sub/p_sub-tau_y)^2)/(n*(n-1))

#復元抽出の場合で抽出ウェイトを用いる場合はこうなる(tau_y)

sum(Y_sub/(n*p_sub))

N_hat=sum(1/(n*p_sub))

#(2.25)

V_tau_y_hat_hat=n*sum((Y_sub/(n*p_sub)-tau_y/n)^2)/(n-1)





#p.39 単純無作為抽出法

Y=c(576,380,74,292,94,158,636,479,236,639,465,133,84,565,25,660,65,148,209,62)

n=sample(1:length(Y),1)

randoms=sample(seq(0,1,0.01),length(Y),replace=T)

comps=data.frame(companys=c(1:length(Y)),Y,randoms)

comps=comps[order(comps$randoms,decreasing = T),]

sigma_y=var(Y)

res=comps[1:n,]

N=length(Y)

#非復元単純無作為抽出

Y=res$Y

tau_y=N*mean(Y)

V_tau_y=(N^2)*(1-n/N)*sigma_y/n

V_tau_y_hat=(N^2)*(1-n/N)*var(Y)

#復元単純無作為抽出

tau_y=N*mean(Y)

V_tau_y=N*(N-1)*sigma_y/n

V_tau_y_hat=(N^2)*var(Y)/n



data=data.frame(n=1:N) %>% mutate(非復元単純無作為抽出法=0,復元単純無作為抽出法=0)

for(j in 1:nrow(data)){

data$非復元単純無作為抽出法[j]=sqrt((1-j/N)*j*var(N*Y/j))

data$復元単純無作為抽出法[j]=sqrt(j*var(N*Y/j))

}

plot(data$n,(data$非復元単純無作為抽出法),col=2,xlim=c(min(data$n),max(data$n)),ylim=c(min(data[,!(colnames(data) %in% c("n"))]),max(data[,!(colnames(data) %in% c("n"))])),xlab="n",ylab="values")

par(new=T)

plot(data$n,(data$復元単純無作為抽出法),col=3,xlim=c(min(data$n),max(data$n)),ylim=c(min(data[,!(colnames(data) %in% c("n"))]),max(data[,!(colnames(data) %in% c("n"))])),xlab="n",ylab="values")

#抽出ウェイトを用いた推定量の表現

n=3;N=20;Y=c(380,639,209)

w=N/n

#非復元単純無作為抽出法

tau_y=sum(w*Y)

V_tau_y_hat=(1-n/N)*n*var(w*Y)

#復元単純無作為抽出法

tau_y=sum(w*Y)

V_tau_y_hat=n*var(w*Y)



```

```{r}

#部分母集団総計の推定　p.48

library(dplyr)

data=data.frame(売上高=c(576,380,74,292,94,158,636,479,236,639,465,133,84,565,25,660,65,148,209,62),所在地=c(rep("市",8),rep("郡",12)),delta=c(rep(1,8),rep(0,12)))

tau_y=sum(data$売上高*data$delta)

#非復元単純無作為抽出標本

N=20;n=5;w=N/n

delta=data$delta;y=data$売上高

tau_y_d_hat=sum(w*delta*y)

V_tau_y_d=sum(1-n/N)*n*var(w*delta*y)

#表3.5

n=3;w=N/n;f=n/N

data=data.frame(売上高=c(380,639,639),地域=c("市","郡","郡"),delta=c(1,0,0))

tau_y_city=sum(data$売上高*data$delta)*w

V_tau_y_city=(1-f)*n*sum((w*data$売上高*data$delta-tau_y_city/n)^2)/(n-1)


#確率比例抽出法　p.51

#復元確率比例抽出法　4.2.1

#Ladhiri(1951)

data=data.frame(comp=c(1:20),x=c(47,31,25,34,22,19,57,42,36,60,51,15,15,54,19,52,19,15,28,22),y=c(576,380,74,292,94,158,636,479,236,639,465,133,84,565,25,660,65,148,209,62))

x=data$x;y=data$y;n=3

N=20

X=max(x);tau_x=sum(x);p=x/sum(x);w=1/(n*p)

samp=c();comp=c();ws=c();ys=c()

while(length(samp)<n)({

k=sample(1:N,1)

a=sample(c(0:X),1)

if(x[k]>a){samp=c(samp,x[k]);comp=c(comp,k);ws=c(ws,w[k]);ys=c(ys,y[k])}

})


N_hat=sum(ws)

tau_y_hat=sum(ws*ys)

tau_x_hat=sum(ws*samp)

#(4.6)

V_tau_y_hat=tau_x*sum(x*(y/x-sum(y)/sum(x))^2)/n

#復元確率比例抽出法Vppsと復元単純無作為抽出法Vsirの比較(Vsir-Vpps)

df=(N^2)*sum((x-tau_x/N)*(y^2/x-sum(y^2/x)/N))/(n*(N-1))

df=cov(x,y^2/x)*(N^2)/n

V_hat_tau_y_hat=n*var(ws*ys)



#4.3.1　非復元確率比例抽出法

data=data.frame(comp=c(1:20),x=c(47,31,25,34,22,19,57,42,36,60,51,15,15,54,19,52,19,15,28,22),y=c(576,380,74,292,94,158,636,479,236,639,465,133,84,565,25,660,65,148,209,62))

x=data$x;y=data$y

#poisson sampling

ns=10;N=length(x)

a=runif(N,min=0,max=1)

tau_x=sum(x);pis=ns*x/tau_x;w=1/pis

xs=x[a<pis];ys=y[a<pis];ws=w[a<pis]

tau_y_hat=sum(ws*ys)

V_tau_y_hat=sum(ws*(ws-1)*ys^2)

V_tau_y=sum((1/pis-1)*y^2)

#複数回サンプリング

times=100;samp=c()

n_vec=c()

for(l in 1:times){

ns=10;N=length(x)

a=runif(N,min=0,max=1)

tau_x=sum(x);pis=ns*x/tau_x;w=1/pis

xs=x[a<pis];ys=y[a<pis];ws=w[a<pis]

tau_y_hat=sum(ws*ys)

n_vec=c(n_vec,sum(a<pis))

V_tau_y_hat=sum(ws*(ws-1)*ys^2)

V_tau_y=sum((1/pis-1)*y^2)

samp=c(samp,V_tau_y_hat)

}

hist(samp);hist(n_vec)


#sunter

x=data$x;y=data$y;x=sort(x)

N=length(x);n=15

tau_x=sum(x)

samp=c();sig=c()

while(length(samp)<n)({

samp=c();sig=c()

for(j in 1:N){

if(j==1){

a=runif(1,0,1)

if(a<n*x[j]/tau_x){samp=c(samp,x[j]);sig=c(sig,1)}else{sig=c(sig,0)}

}else{

a=runif(1,0,1)

if(a<(n-length(samp))*x[j]/(tau_x-sum(x[1:(j-1)]))){samp=c(samp,x[j]);sig=c(sig,1)}else{sig=c(sig,0)}

}

}

})

#系統抽出法

x=data$x;y=data$y;

tau_x=sum(x);tau_y=sum(y);N=length(x)

n=3

v=x/tau_x;pis=n*x/tau_x;w=1/pis

cum_v=cumsum(v)

a=runif(1,0,1)

samp=c((x[cum_v>a])[1]);samp_y=c((y[cum_v>a])[1]);samp_pis=c((pis[cum_v>a])[1]);samp_w=c((w[cum_v>a])[1])


d=((1-a)/(n-1))

for(j in 1:(n-1)){

samp=c(samp,(x[cum_v>=(a+j*d)])[1])

samp_y=c(samp_y,(y[cum_v>=(a+j*d)])[1])

samp_pis=c(samp_pis,(pis[cum_v>=(a+j*d)])[1])

samp_w=c(samp_w,(w[cum_v>=(a+j*d)])[1])

}

tau_y=sum(samp_w*samp_y)

V_tau_y=sum(pis*(1-(n-1)*pis/n)*(y/pis-tau_y/n)^2)

V_tau_y_hat=c()

for(i in 1:n){
for(j in 1:n){

V_tau_y_hat=c(V_tau_y_hat,(1-(samp_pis[i]+samp_pis[j])+sum(pis^2)/n)*(samp_y[i]/samp_pis[i]-samp_y[j]/samp_pis[j])^2)

}}

V_tau_y_hat=sum(V_tau_y_hat)/(2*(n-1))


#比推定量 p.71

data=data.frame(comp=c(1:20),x1=c(47,31,25,34,22,19,57,42,36,60,51,15,15,54,19,52,19,15,28,22),x2=c(29,23,73,4,19,32,28,36,11,36,118,12,6,56,11,22,48,52,24,22),y=c(576,380,74,292,94,158,636,479,236,639,465,133,84,565,25,660,65,148,209,62))

y=data$y;x1=data$x1;x2=data$x2

N=nrow(data)

tau_y=sum(y);tau_x1=sum(x1);tau_x2=sum(x2)

R_x1=tau_y/tau_x1;R_x2=tau_y/tau_x2

n=3;f=n/N

V_tau_y_x1=(N^2)*(1-f)*sum((y-R_x1*x1)^2)/(n*(N-1))

V_tau_y_x2=(N^2)*(1-f)*sum((y-R_x2*x2)^2)/(n*(N-1))

p_y_x1=cor(y,x1)

p_y_x2=cor(y,x2)

#samples

x1_samp=x1[sample(1:N,n)]

x2_samp=x2[sample(1:N,n)]

w=N/n

V_tau_y=(N^2)*(1-n/N)*var(y)/n

V_tau_x1=(N^2)*(1-n/N)*var(x1)/n

V_tau_x2=(N^2)*(1-n/N)*var(x2)/n

CV_tau_x1=sqrt(V_tau_x1)/tau_x1

CV_tau_x2=sqrt(V_tau_x2)/tau_x2












```


```{r}

#p.90

library(dplyr)

Y=c(162,170,172,172,161,155,168,171,156,171,168,169,152,154,172,154,152,167,170,154)



sigma_yr=function(w){
  
N_hat=sum(w)

tau_y=sum(w*Y)

mu_y_N=tau_y/N_hat

val1=sum(w*(Y-mu_y_N)^2)

val2=sum((w-1))

return(val1/val2)
  
}



ite=1000000

w=rep(10,length(Y))

h=0.01;eta=0.1

for(i in 1:ite){
 
w_pre=w  
  
for(j in 1:length(w)){  
   
w_sub=w;w_sub[j]=w_sub[j]+h  

if((w[j]-eta*(sigma_yr(w_sub)-sigma_yr(w))/h)>1){
  
w[j]=w[j]-eta*(sigma_yr(w_sub)-sigma_yr(w))/h  
  
}

}
print(sigma_yr(w))  
  
}  

N_hat=sum(w)

tau_y=sum(w*Y)

mu_y_N=tau_y/N_hat

sigma_yr(w)

z=((Y-mu_y_N)^2-sigma_yr(w))/(N_hat-1)

```

```{r}

#p.91

library(dplyr)

Y=c(162,170,172,172,161,155,168,171,156,171,168,169,152,154,172,154,152,167,170,154)

n=sample(3:length(Y),1)

X=sample(Y,n);Y=sample(Y,n)

sigma_yr=function(w){
  
N_hat=sum(w)

tau_y=sum(w*Y)

mu_y_N=tau_y/N_hat

val1=sum(w*(Y-mu_y_N)^2)

val2=sum((w))

return((n/(n-1))*val1/val2)
  
}

ite=1000000

w=rep(10,length(Y))

h=0.01;eta=0.1

for(i in 1:ite){
 
w_pre=w  
  
for(j in 1:length(w)){  
   
w_sub=w;w_sub[j]=w_sub[j]+h  

if((w[j]-eta*(sigma_yr(w_sub)-sigma_yr(w))/h)>1){
  
w[j]=w[j]-eta*(sigma_yr(w_sub)-sigma_yr(w))/h  
  
}

}
#print(sigma_yr(w))  
  
}  

N_hat=sum(w)

tau_y=sum(w*Y);tau_x=sum(w*X)

mu_y_N=tau_y/N_hat;mu_x_N=tau_x/N_hat

sigma_yr(w)

z=(n*(Y-mu_y_N)^2/(n-1)-sigma_yr(w))/N_hat

sigma_yx_hat=(n/(n-1))*sum(w*(Y-mu_y_N)*(X-mu_x_N))/sum(w)

p_xy=sum(w*(Y-mu_y_N)*(X-mu_x_N))/(sqrt(sum(w*(Y-mu_y_N)^2))*sqrt(sum(w*(X-mu_x_N)^2)))

#plot Y

plot(c(1:length(Y)),rep(mu_y_N,length(Y)),type="p",col=2,xlim=c(1,length(Y)),ylim=c(min(Y),max(Y)),xlab="index",ylab="Y")

par(new=T)

plot(c(1:length(Y)),Y,type="p",col=3,xlim=c(1,length(Y)),ylim=c(min(Y),max(Y)),xlab="index",ylab="Y")

#plot X

plot(c(1:length(X)),rep(mu_x_N,length(X)),type="p",col=2,xlim=c(1,length(X)),ylim=c(min(X),max(X)),xlab="index",ylab="X")

par(new=T)

plot(c(1:length(X)),X,type="p",col=3,xlim=c(1,length(X)),ylim=c(min(X),max(X)),xlab="index",ylab="X")



```