離散データ解析　廣津千尋先生　教育出版　

```{r}

#離散データ解析

library(dplyr)

#ベルヌーイ試行で事象が一度も生起せずに推論に支障をきたす場合、逆に事象が生起する回数rを設定し、それが発生するまでの試行を続けるという考え方が負の二項分布である

#p.4

r=4;p=0.3

negative_binom_values <- rnbinom(10000, size = r, prob = p)

#平均値

mean(negative_binom_values)

r*(1-p)/p

#分散

var(negative_binom_values)

r*(1-p)/(p^2)


```

```{r}

N=c(0:6);n=c(46,86,65,42,21,7,7)

n.=sum(n*N);a=sum(n)

#H0:mu=mu0

mu0=1.8;alpha=0.05

test_statics=abs(n./a-mu0)/sqrt(mu0/a)

qnorm(1-alpha/2)

#confidence interval

under=(sqrt(n./a)-qnorm(1-alpha/2)/sqrt(4*a))^2

upper=(sqrt(n./a)+qnorm(1-alpha/2)/sqrt(4*a))^2


```


```{r}

#p.7

#各確率がポアソン分布(mu_hat)に従う多項分布のあてはめ

N=c(0:9);n=c(46,86,65,42,21,7,5,1,0,1)

mu_hat=sum(N*n)/sum(n)

N=c(0:6);n=c(46,86,65,42,21,7,7)

pi=exp(-mu_hat)*(mu_hat^c(0:6))/factorial(c(0:6))

#chisq_test:H0;

sum((n-sum(n)*pi)^2/(sum(n)*pi))

qchisq(1-0.05/2,df=length(n)-2)



#混合ポアソン分布のあてはめ

N=c(0:9);n=c(46,86,65,42,21,7,5,1,0,1)

mu_hat=sum(N*n)/sum(n)

N=c(0:6);n=c(46,86,65,42,21,7,7)

#mu_hatの密度ｆをN(mu,1)とする

f=function(x){
  
fun=function(mu){  
  
return(dnorm(mu,mu_hat,1)*exp(-mu)*(mu^x)/factorial(x)) 
  
}

return(integrate(fun,0,Inf)$value)
  
}


pi=c()

for(j in 1:length(N)){
  
pi=c(pi,f(N[j]))  
  
}

#chisq_test:H0;

sum((n-sum(n)*pi)^2/(sum(n)*pi))

qchisq(1-0.05/2,df=length(n)-2)



```

```{r}

#p.13

m1=20;n=50;n1=15

#超幾何分布における平均

mean(rhyper(1000,m1,n-m1,n1))

n1*m1/n

#超幾何分布における分散

var(rhyper(1000,m1,n-m1,n1))

(n1*m1/n)*(1-m1/n)*(1-(n1-1)/(n-1))

```

```{r}

#p.16

a=6;b=5

N=array(sample(1:10,a*b,replace = T),dim=c(a,b))

ni=apply(N,1,sum);mj=apply(N,2,sum);n=sum(N)

prob=c(-sum(log(c(1:n))))

for(i in 1:a){
  
prob=c(prob,sum(log(c(1:ni[i]))))  
  
}

for(j in 1:b){
  
prob=c(prob,sum(log(c(1:mj[j]))))  
  
}

E_N=array(0,dim=c(a,b))

for(i in 1:a){
for(j in 1:b){  
  
prob=c(prob,-sum(log(c(1:N[i,j]))))

E_N[i,j]=ni[i]*mj[j]/n
  
}}

exp(sum(prob))



```
```{r}

#p.36

x=c(7,14,27,51)

ni1=c(0,2,7,3)

ni=c(55,157,159,16)

c=1/2

zi=log((ni1+c)/(ni-ni1+c))

v_z=(ni+1)*(ni+2)/(ni*(ni1+1)*(ni-ni1+1))

wi=1/sqrt(v_z)

x.=x-mean(x)

#(イ)

alpha=sum(zi*(x-mean(x)))/sum((x-mean(x))^2)

mu=mean(zi)-alpha*mean(x)

V_alpha=sum(v_z*x.^2)/(sum(x.^2)^2)

cbind(zi,mu+alpha*x)

#(ロ)

wi=1/sqrt((ni-1)/(ni1*(ni-ni1)))

zi_star=wi*log((ni1-c)/(ni-ni1-c))

zi_star=zi_star[2:length(zi_star)]

wi=wi[2:length(wi)]

X=cbind(wi,wi*x[2:length(x)])

params=solve(t(X)%*%X)%*%t(X)%*%zi_star

cov_params=solve(t(X)%*%X)

#(ハ)

t2=sum(ni1*x);t1=sum(ni1)

E_t2=t1*sum(ni*x)/sum(ni)

V_t2=t1*(1-t1/sum(ni))*(sum(ni*x^2)-sum(ni*x)^2/sum(ni))/(sum(ni)-1)

(t2-E_t2)/sqrt(V_t2)

qnorm(1-0.05)

#最尤法

pi=ni1/ni

library(dplyr)

X=as.matrix(cbind(rep(1,length(x)),x))

pis=pi

alpha=0.01;beta=0.01

ite=10

alpha_beta_data=data.frame(ite=1:ite,alpha=0,beta=0)

for(j in 1:ite){

V=array(0,dim=c(length(x),length(x)))

diag(V)=((1/(1+exp(-(X%*%c(alpha,beta)))))*(1-1/(1+exp(-(X%*%c(alpha,beta))))))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%(pis-(1/(1+exp(-(X%*%c(alpha,beta))))))

alpha=alpha+mat[1];beta=beta+mat[2]

print(paste0(alpha,",",beta))

alpha_beta_data$alpha[j]=alpha

alpha_beta_data$beta[j]=beta

}

predict=1/(1+exp(-(X%*%c(alpha,beta))))

```

```{r}

#p.40

#多項超幾何分布の最尤法(修正必要)


a=10;b=12

nij=array(sample(1:10,a*b,replace=T),dim=c(a,b))

Nij=array(sample(15:30,a*b,replace=T),dim=c(a,b))




log_lik=function(lam){
  

value=sum(nij*lam)-sum(log(factorial(nij)))-sum(Nij*lam)+sum(log(factorial(Nij)))

return(value)

}


ite=100000

X=array(0.1,dim=c(a,b))

h=0.01

eta=0.0001

for(l in 1:ite){

vec=X

dX=X

for(i in 1:a){
for(j in 1:b){

vec_sub=vec;vec_sub[i,j]=vec_sub[i,j]+h

dX[i,j]=(log_lik(vec_sub)-log_lik(vec))/h

}}

X=ifelse(X+eta*dX>0,X+eta*dX,X)

print(log_lik(X))

}




```

```{r}

#p.44 

library(dplyr)

data=data.frame(jobs=c(1:10),disease1=c(148,111,645,165,383,96,98,199,59,262),disease2=c(444,352,1911,771,1829,293,330,874,199,1320),disease3=c(86,49,328,119,311,47,58,155,30,236))

mat=data[,-1]

ni.=apply(mat,1,sum)

a=nrow(mat);b=ncol(mat)

log_lik=function(pthis,lam){

p_mat=array(0,dim=c(a,b))

p_mat[a,]=c(exp(pthis),1)/(1+sum(exp(pthis)))

value=cbind(t(exp(t(lam)+c(pthis))),rep(1,a-1))

pij=array(0,dim=c(a-1,b))

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

pij[j,]=value[j,]/sum(value[j,])

}

p_mat[1:(a-1),1:(b)]=pij

vec=0

for(i in 1:a){

val=sum(log(c(1:ni.[i])))

for(j in 1:b){

val=val+mat[i,j]*log(p_mat[i,j])-sum(log(c(1:mat[i,j])))

}

vec=vec+val

}

return(vec)

}


lambda=array(0.1,dim=c(a-1,b-1))

pthi=rep(0.1,b-1)

ite=100000

h=0.01;eta=0.0001

for(l in 1:ite){

vector=pthi

for(j in 1:length(vector)){

vec_sub=vector;vec_sub[j]=vec_sub[j]+h

pthi[j]=pthi[j]+eta*(log_lik(vec_sub,lambda)-log_lik(vector,lambda))/h

}

vector=lambda

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

vec_sub=vector;vec_sub[i,j]=vec_sub[i,j]+h

lambda[i,j]=lambda[i,j]+eta*(log_lik(pthi,vec_sub)-log_lik(pthi,vector))/h

}}

print(log_lik(pthi,lambda))

}

#H:lambda=0

S2=sum((mat-t(t(apply(mat,1,sum)))%*%t(apply(mat,2,sum))/sum(mat))^2/(t(t(apply(mat,1,sum)))%*%t(apply(mat,2,sum))/sum(mat)))

qchisq(1-0.05,df=(a-1)*(b-1))

```

```{r}

#p.48

library(dplyr)

a=5;b=9

nij=matrix(sample(1:10,a*b,replace=T),ncol=b)

f=function(x){
  
k=x[1];sita=x[2:(a+1)];pthi=x[(a+2):(1+a+b)]

value=sum(nij)
  
for(i in 1:a){
for(j in 1:b){  
  
value=value-exp(k+sita[i]+pthi[j])  
  
}}

f1=c()

for(i in 1:a){
  
f1=c(f1,apply(nij,1,sum)[i]-sum(exp(k+sita[i]+pthi)))  
  
}

f2=c()

for(j in 1:b){
  
f2=c(f2,apply(nij,2,sum)[j]-sum(exp(k+sita+pthi[j])))  
  
}

  
return(c(value,f1,f2))

}

ite=10000

X=rep(0.001,a+b+1)

H=diag(f(X))

eta=0.00001

for(l in 1:ite){

X_pre=X  

X=X-eta*H%*%f(X)  

s=X-X_pre

y=f(X)-f(X_pre)

H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H

print(sum(abs(f(X))))


} 

```
```{r}

#p.112

mat=matrix(c(1520,234,177,36,266,1512,362,82,124,432,1772,179,66,78,205,492),ncol=4)

a=nrow(mat)

#(1)対称性の仮説

pij=(t(mat)+mat)/(2*sum(mat))

kai1=sum((mat-sum(mat)*pij)^2/(sum(mat)*pij))

f1=a*(a-1)/2;qchisq(1-0.01,df=f1)

nij=mat;diag(nij)=0;nij2=mat;diag(nij2)=1/2

r1=sum(nij*log(2*nij2/(nij2+t(nij2))))

#(2)周辺対称性の検定

f=function(q){

x=q[1:a];y=q[(a+1):(2*a)]

pij_hat=array(0,dim=c(a,a))

value=array(0,dim=c(a,a))

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

pij_hat[i,j]=mat[i,j]*(1/(1+x[i]-x[j])-(y[i]-y[j])/((1+x[i]-x[j])^2))

value[i,j]=1+x[i]+y[i]-x[j]-y[j]

}}

return(c(c(mat/(sum(mat)*pij)-value),c(apply(pij_hat,1,sum)-apply(pij_hat,2,sum))))

}

alpha=rep(0,a);d=rep(1,a)

X=c(alpha,d)

ite=100000

H=diag(f(X))

eta=0.00000001

for(l in 1:ite){

if(sum(abs(f(X)))>10^(-9)){  

X_pre=X  

X=X-eta*H%*%f(X)  

s=X-X_pre

y=f(X)-f(X_pre)

H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H

print(sum(abs(f(X))))

}

}  


```


```{r}

#Bradley-terry model p.131

library(dplyr)

data=data.frame(力士名=c("北の湖","若乃花","輪島","三重ノ海","貴乃花","旭国"),北の湖=c(0,12,20,9,8,7),貴乃花=c(13,0,15,9,0,13),輪島=c(13,7,0,9,15,3),三重ノ海=c(22,16,27,0,22,15),貴乃花=c(26,0,23,15,0,10),旭国=c(24,13,27,17,24,0))

matrix=data[,2:ncol(data)]

delta=rep(1,nrow(matrix))

r=matrix+t(matrix)

matrix2=matrix;t=apply(matrix2,1,sum)



ite=100;K=300

for(l in 1:ite){

delta2=delta

for(i in 1:length(delta)){

delta[i]=t[i]/sum(r[i,-i]/(delta2[i]+delta2[-i]))

}

delta=K*delta/sum(delta)

print(delta)

}

nij_hat=r*c(delta)/(t(matrix(rep(delta,length(delta)),ncol=length(delta)))+c(delta))

#カイ二乗検定

nij_hat2=nij_hat;nij_hat2[nij_hat2==0]=1

kai2=sum((matrix-nij_hat)^2/nij_hat2)

a=length(t)

f=(a-1)*(a-2)/2-1

qchisq(1-0.05,df=f)

```

```{r}

#各セルの観測値が2項分布に従う場合(交互作用のない場合)

a=5;b=10

nij=array(sample(1:10,a*b,replace=T),dim=c(a,b))/10

rij=nij+array(sample(1:10,a*b,replace=T),dim=c(a,b))/10

alpha=rep(1,a);beta=rep(1,b);

#Broyden Quasi-Newton Algorithm p.213

f=function(x){

mus=x[1];alphas=x[2:(length(alpha)+1)];betas=x[-c(1:(length(alpha)+1))]

f.=-sum(exp(mus+t(t(matrix(rep(alphas,length(beta)),ncol=length(betas)))+c(betas)))*rij/(1+exp(mus+t(t(matrix(rep(alphas,length(betas)),ncol=length(betas)))+c(betas)))))+sum(nij)

f1=c()

for(i in 1:a){

f1=c(f1,-sum(rij[i,]*exp(mus+alphas[i]+betas)/(1+exp(mus+alphas[i]+betas)))+sum(nij[i,]))

}

f2=c()

for(j in 1:b){

f2=c(f2,-sum(rij[,j]*exp(mus+alphas+betas[j])/(1+exp(mus+alphas+betas[j])))+sum(nij[,j]))

}

return(c(f.,f1,f2))

}



X=rep(-0.1,a+b+1)

ite=100000

H=diag(f(X))

eta=0.001

for(l in 1:ite){

if(sum(abs(f(X)))>10^(-9)){  

X_pre=X  

X=X-eta*H%*%f(X)  

s=X-X_pre

y=f(X)-f(X_pre)

H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H

print(sum(abs(f(X))))

}

}  



mu=X[1];alpha=X[2:(length(alpha)+1)];beta=X[-c(1:(length(alpha)+1))]

val=exp(mu+t(t(matrix(rep(alpha,length(beta)),ncol=length(beta))))+c(beta))

mij=rij*val/(1+val)

#H0:rs=0

kai2=2*(sum(nij*log(nij/mij))+sum((rij-nij)*log((rij-nij)/(rij-mij))))

qchisq(1-0.05,df=(a-1)*(b-1))





```

```{r}

#各セルの観測値が2項分布に従う場合(交互作用のある場合)

a=5;b=10

nij=array(sample(1:10,a*b,replace=T),dim=c(a,b))/10

rij=nij+array(sample(1:10,a*b,replace=T),dim=c(a,b))/10

alpha=rep(1,a);beta=rep(1,b);

#Broyden Quasi-Newton Algorithm p.213

f=function(x){

mus=x[1];alphas=x[2:(length(alpha)+1)];betas=x[c((length(alpha)+2):(1+length(alpha)+length(beta)))]

rs=x[-c(1:(1+length(alphas)+length(betas)))];rs=matrix(rs,ncol=b)

f.=-sum(exp(mus+rs+t(t(matrix(rep(alphas,length(betas)),ncol=length(beta)))+c(betas)))*rij/(1+exp(mus+rs+t(t(matrix(rep(alphas,length(betas)),ncol=length(beta)))+c(betas)))))+sum(nij)

f1=c()

for(i in 1:a){

f1=c(f1,-sum(rij[i,]*exp(mus+alphas[i]+betas+rs[i,])/(1+exp(mus+alphas[i]+betas+rs[i,])))+sum(nij[i,]))

}

f2=c()

for(j in 1:b){

f2=c(f2,-sum(rij[,j]*exp(mus+alphas+betas[j]+rs[,j])/(1+exp(mus+alphas+betas[j]+rs[,j])))+sum(nij[,j]))

}

f3=c()

for(i in 1:a){
for(j in 1:b){

f3=c(f3,-as.numeric(rij[i,j]*exp(mus+alphas[i]+betas[j]+rs[i,j])/(1+exp(mus+alphas[i]+betas[j]+rs[i,j])))+nij[i,j])

}}

return(c(f.,f1,f2,f3))

}



X=rep(-0.1,a+b+1+a*b)

ite=100000

H=diag(f(X))

eta=0.00001

for(l in 1:ite){

if(sum(abs(f(X)))>10^(-9)){  

X_pre=X  

X=X-eta*H%*%f(X)  

s=X-X_pre

y=f(X)-f(X_pre)

H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H

print(sum(abs(f(X))))

}

}  







```

```{r}

#p.140

rij=matrix(c(102,64,67,46,49,37),ncol=3)

nij=matrix(c(20,10,26,16,27,14),ncol=3)

a=nrow(rij);b=ncol(rij)

f=function(x){

mu=x[1];x=x[-1]

alpha=x[1:a];x=x[-c(1:a)];

beta=x;

mat_ab=array(0,dim=c(a,b))

for(i in 1:a){
for(j in 1:b){

mat_ab[i,j]=rij[i,j]*exp(mu+alpha[i]+beta[j])/(1+exp(mu+alpha[i]+beta[j]))

}}

mat_a=-apply(mat_ab,1,sum)+apply(nij,1,sum)

mat_a[length(mat_a)]=0

mat_b=-apply(mat_ab,2,sum)+apply(nij,2,sum)

mat_b[length(mat_b)]=0

mat_ab=-sum(mat_ab)+sum(nij)

return(c(mat_ab,mat_a,mat_b))

}



#fiacco,mccormick(1968),murtagh,sargent(1969)

#初期値
X=rep(0,a+b+1)

#ヘッシアンの近似行列の初期値
H=diag(1,length(X))

#反復計算回数
ite=10000

#最小点から非常にかけ離れた点になるのを防ぐためのスカラ
eta=10^(-3)

for(l in 1:ite){

#各代数方程式の値が入ったベクトル
f_val=f(X)

#座標点に対する更新量ベクトル
sigma=t(t(-H%*%f_val))

#etaによって学習幅を設定し座標点を更新
X=X+eta*sigma

#以前の座標点と更新された座標点の関数値の差のベクトル(y)
y=f(X)-f_val

#ヘッシアンの逆行列の近似行列の更新
H=H+(eta*sigma-H%*%y)%*%t(eta*sigma-H%*%y)/as.numeric(t(y)%*%(eta*sigma-H%*%y))

print(f(X))

}

mu=X[1];X=X[-1];alpha=X[1:a];X=X[-c(1:a)];beta=X

pij=array(0,dim=c(a,b))

for(i in 1:a){
for(j in 1:b){

pij[i,j]=exp(mu+alpha[i]+beta[j])/(1+exp(mu+alpha[i]+beta[j]))

}}


mij_hat=rij*pij

2*sum(nij*log(nij/mij_hat)+sum((rij-nij)*log((rij-nij)/(rij-mij_hat))))

df=(a-1)*(b-1)

```
```{r}

#p.140

#母親の総数
rij=matrix(c(102,64,67,46,49,37),ncol=3)

#以前出産事故のあった母親の総数
nij=matrix(c(20,10,26,16,27,14),ncol=3)

a=nrow(rij);b=ncol(rij)

f=function(x){

mu=x[1];x=x[-1]

alpha=x[1:a];x=x[-c(1:a)];

beta=x;

mat_ab=array(0,dim=c(a,b))

for(i in 1:a){
for(j in 1:b){

mat_ab[i,j]=rij[i,j]*exp(mu+alpha[i]+beta[j])/(1+exp(mu+alpha[i]+beta[j]))

}}

mat_a=-apply(mat_ab,1,sum)+apply(nij,1,sum)

mat_a[length(mat_a)]=0

mat_b=-apply(mat_ab,2,sum)+apply(nij,2,sum)

mat_b[length(mat_b)]=0

mat_ab=-sum(mat_ab)+sum(nij)

return(c(mat_ab,mat_a,mat_b))

}



#Broyden Quasi-Newton Algorithm 

#初期値
X=rep(0,a+b+1)

#反復計算回数
ite=10000

#最小点から非常にかけ離れた点になるのを防ぐためのスカラ
eta=10^(-3)

#ヤコビアンの初期値
H=diag(f(X))

for(l in 1:ite){

#以前の座標点を保存
X_pre=X  

#座標点を更新
X=X-eta*H%*%f(X)  

#以前の座標点と更新された座標点の差のベクトル(s)
s=X-X_pre

#以前の座標点と更新された座標点の関数値の差のベクトル(y)
y=f(X)-f(X_pre)

#ヤコビアンの近似行列を更新する
H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H

print((f(X)))

}
 


#計算結果から主効果等のパラメータを取り出す
mu=X[1];X=X[-1];alpha=X[1:a];X=X[-c(1:a)];beta=X

pij=array(0,dim=c(a,b))

for(i in 1:a){
for(j in 1:b){

pij[i,j]=exp(mu+alpha[i]+beta[j])/(1+exp(mu+alpha[i]+beta[j]))

}}


mij_hat=rij*pij

#尤度比検定の統計量
2*sum(nij*log(nij/mij_hat)+sum((rij-nij)*log((rij-nij)/(rij-mij_hat))))

#自由度
df=(a-1)*(b-1)

qchisq(1-0.05,df)

```

```{r}

#p.140

#母親の総数
rij=matrix(c(102,64,67,46,49,37),ncol=3)

#以前出産事故のあった母親の総数
nij=matrix(c(20,10,26,16,27,14),ncol=3)

a=nrow(rij);b=ncol(rij)

f=function(x){

mu=x[1];x=x[-1]

alpha=x[1:a];x=x[-c(1:a)];

beta=x;

mat_ab=array(0,dim=c(a,b))

for(i in 1:a){
for(j in 1:b){

mat_ab[i,j]=rij[i,j]*exp(mu+alpha[i]+beta[j])/(1+exp(mu+alpha[i]+beta[j]))

}}

mat_a=-apply(mat_ab,1,sum)+apply(nij,1,sum)

mat_a[length(mat_a)]=0

mat_b=-apply(mat_ab,2,sum)+apply(nij,2,sum)

mat_b[length(mat_b)]=0

mat_ab=-sum(mat_ab)+sum(nij)

return(c(mat_ab,mat_a,mat_b))

}



#levenberg(1944) 

#初期値
X=rep(0,a+b+1)

#反復計算回数
ite=10000

#最小点から非常にかけ離れた点になるのを防ぐためのスカラ
eta=10^(-2)

#差分メッシュ
h=0.01

#正の小さいスカラ
w=0.1

for(l in 1:ite){

#ヤコビアンを計算する箱を用意
df=array(0,dim=c(length(X),length(X)))

#1次差分を用いて数値微分し、ヤコビアンを計算する
for(j in 1:length(X)){

vec=X;vec[j]=vec[j]+h  

df[j,]=(f(vec)-f(X))/h  

}

#逆行列が特異にならないよう、対角成分をwで調整
A=df%*%t(df)+diag(w,length(X))  

y=df%*%f(X)

#座標点更新の際に必要なベクトルを計算
dx=solve(A)%*%(-y)

#0で固定のパラメータの更新量を0にする
dx[c(3,6)]=0

#学習幅に従って座標点を更新
X=X+eta*dx

print(f(X))

}

 


#計算結果から主効果等のパラメータを取り出す
mu=X[1];X=X[-1];alpha=X[1:a];X=X[-c(1:a)];beta=X

pij=array(0,dim=c(a,b))

for(i in 1:a){
for(j in 1:b){

pij[i,j]=exp(mu+alpha[i]+beta[j])/(1+exp(mu+alpha[i]+beta[j]))

}}


mij_hat=rij*pij

#尤度比検定の統計量
2*sum(nij*log(nij/mij_hat)+sum((rij-nij)*log((rij-nij)/(rij-mij_hat))))

#自由度
df=(a-1)*(b-1)

qchisq(1-0.05,df)

```


```{r}
#p.142


#p.148


library(dplyr)

r_data=data.frame(i=c(1,2),r1=c(102,64),r2=c(67,46),r3=c(49,37))

n_data=data.frame(i=c(1,2),n1=c(20,10),n2=c(26,16),n3=c(27,14))

rij=as.matrix(r_data[,2:ncol(r_data)])

nij=as.matrix(n_data[,2:ncol(n_data)])

a=nrow(r_mat);b=ncol(r_mat)


zij=log(nij/(rij-nij))

wij=sqrt(nij*(rij-nij)/rij)


#対立仮説(r11!=r22!=0)

yij=c(t(zij*wij))

W=c(t(wij))

W_mat=n*W

W_mat2=W_mat[,1:4]

solve(t(W_mat2)%*%W_mat2)%*%t(W_mat2)%*%yij


#帰無仮説(r11=r22=0)

n=matrix(c(1,1,1,0,0,0,1,0,0,1,0,0,0,1,0,0,1,0,0,0,1,0,0,1),ncol=4)

mat=n*c(t(wij))

trends=solve(t(mat)%*%mat)%*%t(mat)%*%yij

alpha=c(trends[1],0);trends=trends[-1]

mu=trends[3];trends=trends[-3]

beta=c(trends-mu,0)

zij_hat=array(0,dim=c(a,b))

pij=array(0,dim=c(a,b))

for(i in 1:a){
for(j in 1:b){

zij_hat[i,j]=mu+alpha[i]+beta[j]

pij[i,j]=exp(zij_hat[i,j])/(1+exp(zij_hat[i,j]))

}}

Sr=sum(((zij_hat-zij)*wij)^2)

qchisq(1-0.05,df=2)

mij_hat=rij*pij

#適合度検定

kai2=sum((nij-mij_hat)^2/mij_hat)+sum(((rij-nij)-(rij-mij_hat))^2/(rij-mij_hat))






```




```{r}

#離散データ解析 p.170

library(dplyr)

#2つの薬をそれぞれn_1,n_2人に任意に投与する

n=50;

n1=sample(1:n,1);n2=n-n1

samples1=c();samples2=c()

print(n1*(n1-1)/(n*(n-1)))

for(j in 1:100){

samples1=c(samples1,(sum(sample(c(0:1),n,prob=c(n2/n,n1/n),replace=T))/(n))*sum(sample(c(0:1),n,prob=c(1-(n1-1)/(n-1),(n1-1)/(n-1)),replace=T))/(n-1))

samples2=c(samples2,(sum(sample(c(0:1),n,prob=c(n2/n,n1/n),replace=T))/(n))*sum(sample(c(0:1),n,prob=c(1-(n2)/(n-1),(n2)/(n-1)),replace=T))/(n-1))

}

#E(d_iv*d_iu)=prob(d_iv=1,d_iu=1)

plot(samples1,type="p",main=paste0("prob:",n1*(n1-1)/(n*(n-1))))

#E(d_1v*d_2u)=prob(d_1v=1,d_2u=1)

plot(samples2,type="p",main=paste0("prob:",n1*(n2)/(n*(n-1))))





#p1をA1が投与されたときに治癒する確率,p2をA2が投与されたときに治癒する確率,が既知とする

Y1_mat=array(0,dim=c(100,n))
Y2_mat=array(0,dim=c(100,n))

for(i in 1:100){

N1=sample(seq(0,1,0.01),n,replace=T)
N2=sample(seq(0,1,0.01),n,replace=T)

d1=sample(seq(0,1,0.01),n,replace=T)
d2=sample(seq(0,1,0.01),n,replace=T)

Y1=c();Y2=c()

for(j in 1:n){
  
Y1=c(Y1,sample(c(0:1),1,prob=c(1-d1[j],d1[j]))*sample(c(0:1),1,prob=c(1-N1[j],N1[j])))  
  
Y2=c(Y2,sample(c(0:1),1,prob=c(1-d2[j],d2[j]))*sample(c(0:1),1,prob=c(1-N2[j],N2[j])))

}

Y1_mat[i,]=Y1;Y2_mat[i,]=Y2

}

p1=apply(Y1_mat,2,mean)
p2=apply(Y2_mat,2,mean)

p_data=data.frame(class=sample(c(1:2),n,replace=T),p1=p1,p2=p2)

#E[Y1/n1]

E_Y1=sum(p1)/n

#E[Y2/n2]

E_Y2=sum(p2)/n

#V[Y1/n1-Y2/n2]

V_Z=sum(p1)/(n*n1)+sum(p2)/(n*n2)-(sum(p1)/n)^2-(sum(p2)/n)^2

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

if(j!=i){  
   
V_Z=V_Z+p1[j]*p1[i]*(n1-1)/(n1*n*(n-1))   
  
V_Z=V_Z+p2[j]*p2[i]*(n2-1)/(n2*n*(n-1))   

V_Z=V_Z-2*p1[j]*p2[i]/(n*(n-1)) 

}  
  
}  
}

V_Z=V_Z+2*(sum(p1)/n)*(sum(p2)/n)

#H0:mean(p1)=mean(p2)

#p1i=p2iならE[V0]-V0(Z)は0になる

E_V0_V0_Z=(((n-2)/(n-1))*sum((p1-p2)^2)+(1/n2-1/n1)*sum(p1^2-p2^2))/(n*(n-1))

plot(c(1:n),p1,type="p",col=2,ylim=c(min(c(p1,p2)),max(c(p1,p2))),ylab="prob")

par(new=T)

plot(c(1:n),p2,type="p",col=3,ylim=c(min(c(p1,p2)),max(c(p1,p2))),ylab="prob")

```

```{r}

#サンプルで漸近正規性をつかって検定してみる

Y1=Y1_mat[1,];Y2=Y2_mat[1,]

value=(n1*n2*(n-1))*(sum(Y1)/n1-sum(Y2)/n2)^2/((sum(Y1)+sum(Y2))*(n-sum(Y1)-sum(Y2)))

qnorm(1-0.05/2)

```
```{r}
#個人による効果を無視した場合

n1=50;n2=60;

n=n1+n2

Y1=30;Y2=30

n1*n2*(n-1)*(Y1/n1-Y2/n2)^2/((Y1+Y2)*(n-Y1-Y2))

qnorm(1-0.05/2)

```

```{r}
#H0:mean(p1)=mean(p2)

#p1i=p2iならE[V0]-V0(Z)は0になる

n1=50;n2=60;

n=n1+n2

p1=rbinom(n,size=50,prob=0.5)

p2=rbinom(n,size=100,prob=0.01)

E_V0_V0_Z=(((n-2)/(n-1))*sum((p1-p2)^2)+(1/n2-1/n1)*sum(p1^2-p2^2))/(n*(n-1))

plot(c(1:n),p1,type="p",col=2,ylim=c(min(c(p1,p2)),max(c(p1,p2))),ylab="prob")

par(new=T)

plot(c(1:n),p2,type="p",col=3,ylim=c(min(c(p1,p2)),max(c(p1,p2))),ylab="prob")


```