分散分析　廣津千尋先生　教育出版

```{r}

#p.10 logistic model

#性別による効果の違いがある場合の線形回帰模型

#混合モデル

I=10;J=5;K=2

yijk=data.frame(i=0,j=0,k=0,y=0)

for(s in 1:I){
for(q in 1:J){
for(l in 1:K){

x_sub=sample(0:1,1);t=j

y_sub=rpois(1,5)*rpois(1,9)*x_sub*t+rpois(1,5)*t+rpois(1,8)*x_sub+1

y_sub=sample(0:1,1)

data_sub=data.frame(i=s,j=q,k=l,y=y_sub)

yijk=rbind(yijk,data_sub)

}}}

yijk=tail(yijk,nrow(yijk)-1)

delta_num1=c();delta_num2=c()

for(j in 1:J){
for(k in 1:K){

delta_num1=c(delta_num1,j)
delta_num2=c(delta_num2,k)

}}

delta_number=data.frame(num=c(1:length(delta_num1)),delta_num1,delta_num2)

#params

mu=0.01;a=rep(0.01,I);b=rep(0.01,J);r=rep(0.01,K);delta=rep(0.01,J*K)

params=c(mu,a,b,r,delta)

f=function(x){

mus=x[1];alpha=x[2:(I+1)];beta=x[(I+2):(1+I+J)];deltas=x[-c(1:(1+I+J))];rs=deltas[1:K];deltas=deltas[-c(1:K)]

delta_num=c()

for(l in 1:nrow(yijk)){

delta_num=c(delta_num,delta_number$num[as.numeric(delta_number$delta_num1)==as.numeric(yijk$j[l]) & as.numeric(delta_number$delta_num2)==as.numeric(yijk$k[l])])

}

values=mu+alpha[yijk$i]+beta[yijk$j]+rs[yijk$k]+deltas[delta_num]

pijk=1/(1+exp(-values))

value=sum(yijk$y*log(ifelse(pijk>0,pijk,1))+(1-yijk$y)*log(ifelse(1-pijk>0,1-pijk,1)))

return(value)

}


X=params

ite=10000;eta=0.001

h=0.01

for(l in 1:ite){

vec=X

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

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

X[j]=X[j]+eta*(f(vec_sub)-f(vec))/h

}

print(f(X))

}
```

```{r}

#医学統計勉強会　p.11

#性別による効果の違いがある場合の線形回帰模型

#混合モデル

I=10;J=5;K=5

yijk=data.frame(i=0,j=0,k=0,y=0)

sigma=1

x_vec=c();t_vec=c()

for(i in 1:I){
for(j in 1:J){
for(k in 1:K){

x_sub=sample(0:1,1);t=j

x_vec=c(x_vec,x_sub);t_vec=c(t_vec,t)

y_sub=x_sub*t+rpois(1,5)*t+rpois(1,8)*x_sub+1+rnorm(1,sigma)

data_sub=data.frame(i=i,j=j,k=k,y=y_sub)

yijk=rbind(yijk,data_sub)

}}}

yijk=tail(yijk,nrow(yijk)-1)

mat=cbind(rep(1,length(x_vec)),x_vec,t_vec,x_vec*t_vec)

mat2=t(mat)%*%mat

params=svd(mat2)$u%*%diag(1/svd(mat2)$d)%*%t(svd(mat2)$v)%*%t(mat)%*%as.numeric(yijk$y)

predict=mat%*%params

plot(c(1:nrow(yijk)),yijk$y,xlim=c(1,nrow(yijk)),ylim=c(min(c(yijk$y,predict)),max(c(yijk$y,predict))),xlab="index",ylab="values",col=2)

par(new=T)

plot(c(1:nrow(yijk)),predict,xlim=c(1,nrow(yijk)),ylim=c(min(c(yijk$y,predict)),max(c(yijk$y,predict))),xlab="index",ylab="values",col=3)

```


```{r}
library(dplyr)

r=10;a=5;b=6;n=r*a

pi=rpois(a,lambda=5);pj=rpois(b,lambda=9)

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

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

mat[i,j]=rpois(1,pi[i])+rpois(1,pj[j]) 

}
}

pthi1=array(0,dim=c(a*b,a))

for(j in 1:a){
  
pthi1[(1+b*(j-1)):(b*j),j]=1  
  
}

pthi2=diag(1,b)

for(j in 1:(a-1)){
  
pthi2=rbind(pthi2,diag(1,b))
  
}


pthi=cbind(rep(1,a*b),pthi1,pthi2)

y=c(mat)

mu=mean(mat);alpha_a=apply(mat,1,mean)-mu;alpha_b=apply(mat,2,mean)-mu

params=c(mu,alpha_a,alpha_b)

predict=pthi%*%params

plot(c(1:prod(dim(mat))),y,type="o",col=2,ylim=c(min(y,predict),max(y,predict)),xlab="samples",ylab="values")

par(new=T)

plot(c(1:prod(dim(mat))),predict,type="o",col=3,ylim=c(min(y,predict),max(y,predict)),xlab="samples",ylab="values")

sigma_hat=sum((y-predict)^2)/prod(dim(mat))

library(mvtnorm)

n=100

for(j in 1:n){

samples=rmvnorm(1,mean=predict,diag(sigma_hat,prod(dim(mat))))

plot(c(1:prod(dim(mat))),y,type="o",col=2,ylim=c(min(y,predict),max(y,predict)),xlab="samples",ylab="values")

par(new=T)

plot(c(1:prod(dim(mat))),samples,type="o",col=3,ylim=c(min(y,predict),max(y,predict)),xlab="samples",ylab="values")


}

J=rep(1,a*b);A=pthi1;B=pthi2

P_pthi=A%*%solve(t(A)%*%A)%*%t(A)+(B)%*%solve(t(B)%*%B)%*%t(B)-J%*%solve(t(J)%*%J)%*%t(J)

P_A=A%*%solve(t(A)%*%A)%*%t(A);P_B=(B)%*%solve(t(B)%*%B)%*%t(B)

P_J=J%*%solve(t(J)%*%J)%*%t(J)

(P_A-P_J)%*%(P_B-P_J)

#cochran theorem

p=a;r=b



#CF

CF=t(y)%*%P_J%*%y

SSA=t(y)%*%P_A%*%y-CF

SSB=t(y)%*%P_B%*%y-CF

SSR=t(y)%*%(diag(1,prod(dim(mat)))-P_pthi)%*%y

#H1:alpha_a

(SSA/(p-1))/((SSR)/((p-1)*(r-1)))

qf(1-0.05,p-1,(p-1)*(r-1))

#H2:alpha_b

(SSB/(r-1))/((SSR)/((p-1)*(r-1)))

qf(1-0.05,r-1,(p-1)*(r-1))

```

```{r}

#p.30

library(mvtnorm)

sigma_vec<- matrix(c(4,0,0,3), ncol=2)

n=500

samples=rmvnorm(n,mean=c(1,2),sigma=sigma_vec)

mean=c(1,2)

dmvnorm(c(1,2),mean=mean,sigma=sigma_vec)

ave=apply(samples,2,mean)

var=diag(0,2)

for(j in 1:n){

var=var+(samples[j,]-mean)%*%t(samples[j,]-mean)

}

var=(1/(n-1))*var

A=matrix(c(1,2,3,4),ncol=2)

E_yAy=sum(diag(A%*%sigma_vec))+ave%*%A%*%t(t(ave))

V_yAy=2*(sum(diag(A%*%sigma_vec))^2+2*ave%*%A%*%sigma_vec%*%A%*%t(t(ave)))


#p.42

mu1=3;mu0=2;sigma0=9;n=100;alpha=0.05

samples=rnorm(n,mu1,sigma0/n)

y_bar=mean(samples)

(y_bar-mu1)/(sqrt(sigma0)/sqrt(n))

#検出力

#k=(mu1-mu0)/sqrt(sigma0)

data=data.frame(k=seq(-1,1,0.01),power=0)

for(j in 1:nrow(data)){
  
data$power[j]=1-pnorm(qnorm(1-alpha/2)-sqrt(n)*data$k[j])+pnorm(-qnorm(1-alpha/2)-sqrt(n)*data$k[j])

}
  
plot(data$k,data$power)  

alpha=0.05;beta=0.10

data=data.frame(k=seq(-1,1,0.01),n=0)

for(j in 1:nrow(data)){
  
data$n[j]=(qnorm(1-alpha)+qnorm(1-beta))^2/(data$k[j]^2)  

}

#p.46 H0:mu0=mu

n=100;mu0=1;alpha=0.05

y=rnorm(n,mean=2,sd=4)

sigma_hat=sum((y-mean(y))^2)/(n-1)

t=(sqrt(n)*(mean(y)-mu0))/sqrt(sigma_hat)

qt(1-alpha,n-1)

#sigma2の検定　p.49 muが未知

n=100;alpha=0.05

y=rnorm(n,mean=2,sd=4)

sigma0=16

sum((y-mean(y))^2)/sigma0

qchisq(1-alpha,n-1)


```


```{r}

#一元配置 p.63

#母平均の推定

library(dplyr)

n=10;row=20

y_mat=array(0,dim=c(n,row))

for(j in 1:n){
  
y_mat[j,]=c(sample(c(1:10),1)+rnorm(row,mean=0,sd=3))  
  
}

data=data.frame(class=1:n,y_mat)

mu_i=apply(y_mat,1,mean)

mu=mean(y_mat)

alpha_i=mu_i-mu

val=0

for(i in 1:n){
for(j in 1:row){
  
val=val+(y_mat[i,j]-sum(y_mat[i,])/row)^2  
  
}
}

sigma_hat=val/(prod(dim(y_mat))-n)

a=n;n=prod(dim(y_mat))

under=array(0,dim=c(a,a));upper=array(0,dim=c(a,a))

for(j in 1:a){
for(i in 1:a){  
 
if(i!=j){
  
under[i,j]=mu_i[i]-mu_i[j]-sqrt((1/row+1/row)*sigma_hat)*qt(1-0.05/2,df=n-a)  
 
upper[i,j]=mu_i[i]-mu_i[j]+sqrt((1/row+1/row)*sigma_hat)*qt(1-0.05/2,df=n-a)  
  
}   
  
}
}

#等分散性の検定

S_vec=c()

for(j in 1:nrow(y_mat)){
  
S_vec=c(S_vec,sum((y_mat[j,]-apply(y_mat,1,mean)[j])^2 )/(row-1)) 

}

s=S_vec*(row-1)

#H0:等分散性のカイ二乗統計量

(n-a)*(log(sum((row-1)*s)/(n-a))-sum((row-1)*log(s))/(n-a))

qchisq(1-0.05/2,df=a-1)




```

```{r}

#2元配置

sigma=1;a=10;b=15

#サンプルサイズ

n=30

#平均値データ作成

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

y=rnorm(array(n,dim=c(a,b)),mu_mat,rep(sigma,prod(dim(mu_mat))))

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

for(j in 1:n){
  
y=rnorm(array(n,dim=c(a,b)),mu_mat,rep(sigma,prod(dim(mu_mat))))

Y=Y+y  
  
}

Y=Y/n

Yi..=apply(Y,2,mean)

Y.j.=apply(Y,2,mean)

Yij.=Y

Y...=mean(Y)

S_ab=c()

for(i in 1:a){
for(j in 1:b){  
  
S_ab=c(S_ab,(Yij.[i,j]-Yi..[i]-Y.j.[j]-Y...)^2)
    
}
}

S_ab=sum(S_ab*n)

S_mu=n*Y...^2

S_a=b*n*sum(Yi..^2)-S_mu

S_b=a*n*sum(Y.j.^2)-S_mu

S_e=sum((Y*n)^2)-S_a-S_b-S_ab-S_mu

#H_ab;(ab)は同一

f_ab=(S_ab/((a-1)*(b-1)))/(S_e/(n-a*b))

qf(1-0.05,(a-1)*(b-1),n-a*b)

#H_a;(a)は同一

f_a=(S_a/(a-1))/(S_e/(n-a*b))

qf(1-0.05,(a-1),n-a*b)

#H_b;(b)は同一

f_b=(S_b/(b-1))/(S_e/(n-a*b))

qf(1-0.05,(b-1),n-a*b)

#交互作用が存在する場合(mu_ijの区間推定)

under=Y-sqrt(1/row)*sqrt(S_e/(n-a*b))*qt(1-0.05/2,(n-a*b))

upper=Y+sqrt(1/row)*sqrt(S_e/(n-a*b))*qt(1-0.05/2,(n-a*b))


```

```{r}

#p.82 分散分析　廣津千尋

library(dplyr)

r=10;a=5;b=6;n=r*a

dir.create("~/variance_samples")

pi=rpois(a,lambda=5);pj=rpois(b,lambda=9)

for(k in 1:r){
  
mat=array(0,dim=c(a,b))  
  
for(j in 1:ncol(mat)){
for(i in 1:nrow(mat)){  

mat[i,j]=rpois(1,pi[i])+rpois(1,pj[j]) 

}
}
  
write.csv(mat,paste0("~/variance_samples/sample_",k,".csv"))  

}   

y1=rep(0,b);y2=rep(0,a);y3=0;ave_mat=0;s_e=0

i_num=1;j_num=1

for(k in 1:r){
  
mat=read.csv(paste0("~/variance_samples/sample_",k,".csv")) 

mat=mat[,2:ncol(mat)]

y1=y1+apply(mat,2,sum)/(r*b)
 
y2=y2+apply(mat,1,sum)/(r*a)

y3=y3+sum(mat)/(r*a*b)

ave_mat=ave_mat+mat/r
  
}

for(k in 1:r){
  
mat=read.csv(paste0("~/variance_samples/sample_",k,".csv")) 

mat=mat[,2:ncol(mat)]

s_e=s_e+sum((mat-ave_mat)^2)

}

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



for(j in 1:b){
for(i in 1:a){  
  
mu[i,j]=y1[j]+y2[i]-y3 
   
}
}

r_e=n/(a+b-1)

sigma_hat=sqrt(s_e/(n-a*b))

alpha=0.1

upper=qt(1-alpha/2,df=n-a*b)*sigma_hat*sqrt(1/r_e)

under=qt(alpha/2,df=n-a*b)*sigma_hat*sqrt(1/r_e)

mu_upper=mu+upper;mu_under=mu+under

```

```{r}

#p.83 4-2-5

library(dplyr)

r=10;a=5;b=6;n=r*a

pi=rpois(a,lambda=5);pj=rpois(b,lambda=9)

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

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

mat[i,j]=rpois(1,pi[i])+rpois(1,pj[j]) 

}
}

alpha=apply(mat,1,mean)-mean(as.matrix(mat))

beta=apply(mat,2,mean)-mean(as.matrix(mat))

g1=0;g2=0;S_e=0;y_ab=mat;g_ab=mat

for(i in 1:nrow(mat)){
for(j in 1:ncol(mat)){  
  
g1=g1+alpha[i]*beta[j]*mat[i,j] 

g2=g2+(alpha[i]*beta[j])^2
   
S_e=S_e+(mat[i,j]-apply(mat,1,mean)[i]-apply(mat,2,mean)[j]+mean(as.matrix(mat)))^2

y_ab[i,j]=y_ab[i,j]-alpha[i]-beta[j]



}
}

g=g1/g2

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

g_ab[i,j]=g*alpha[i]*beta[j]  
  
}}  
  
S_g=g1^2/(sum(alpha^2)*sum(beta^2))

S_e=S_e-S_g

#H:g=0

p=0.05

f=S_g/(S_e/((a-1)*(b-1)-1))

qf(1-p,1,(a-1)*(b-1)-1)

f>qf(1-p,1,(a-1)*(b-1)-1)

#線形模型

mu=mean(as.matrix(y_ab-g_ab))

predict=mat

mat_value=c();predict_value=c()

g_hat=g

for(i in 1:nrow(mat)){
for(j in 1:ncol(mat)){  
  
predict[i,j]=mu+alpha[i]+beta[j]+g*alpha[i]*beta[j] 
mat_value=c(mat_value,mat[i,j])

predict_value=c(predict_value,predict[i,j])
   
}
}

plot(c(1:length(mat_value)),mat_value,type="o",col=2,ylim=c(min(mat_value,predict_value),max(mat_value,predict_value)),xlab="index",ylab="values")

par(new=T)

plot(c(1:length(mat_value)),predict_value,type="o",col=3,ylim=c(min(mat_value,predict_value),max(mat_value,predict_value)),xlab="index",ylab="values")

```

```{r}

#p.86 三元配置

n=20

data=data.frame(i=0,j=0,k=0,l=0,value=0)

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

dat=data.frame(i=i,j=j,k=k,l=1:n,value=(i+j+k)+rnorm(n,sample(1:10,1)))

data=rbind(data,dat)

}  
}  
}

data=tail(data,nrow(data)-1)

sub_i=data %>% group_by(i) %>% summarise(value=mean(value))

sub_j=data %>% group_by(j) %>% summarise(value=mean(value))

sub_k=data %>% group_by(k) %>% summarise(value=mean(value))

sub_all=mean(data$value)

sub_ij=data %>% group_by(i,j) %>% summarise(value=mean(value))

sub_ik=data %>% group_by(i,k) %>% summarise(value=mean(value))

sub_jk=data %>% group_by(j,k) %>% summarise(value=mean(value))

sub_ijk=data %>% group_by(i,j,k) %>% summarise(value=mean(value))



S_a=(n^3)*sum((sub_i$value-sub_all)^2)

S_b=(n^3)*sum((sub_j$value-sub_all)^2)

S_r=(n^3)*sum((sub_k$value-sub_all)^2)

S_ab=c()

for(I in 1:n){
for(J in 1:n){  
 
S_ab=c(S_ab,(sub_ij$value[sub_ij$i==I & sub_ij$j==J]-sub_i$value[sub_i$i==I]-sub_j$value[sub_j$j==J]+sub_all)^2)   
  
}
}

S_ab=sum(S_ab)*n^2

S_ar=c()

for(I in 1:n){
for(K in 1:n){  
  
S_ar=c(S_ar,(sub_ik$value[sub_ik$i==I & sub_ik$k==K]-sub_i$value[sub_i$i==I]-sub_k$value[sub_k$k==K]+sub_all)^2) 
   
}
}

S_ar=sum(S_ar)*n^2

S_br=c()

for(J in 1:n){
for(K in 1:n){
  
S_br=c(S_br,(sub_jk$value[sub_jk$j==J & sub_jk$k==K]-sub_j$value[sub_j$j==J]-sub_k$value[sub_k$k==K]+sub_all)^2)
    
}  
}

S_br=sum(S_br)*n^2

S_e=c();S_T=c();S_ABC=c()

for(I in 1:n){
for(J in 1:n){  
for(K in 1:n){ 
 
vec=data %>% filter(i==I & j==J & k==K)

S_e=c(S_e,sum((vec$value-sub_ijk$value[sub_ijk$i==I & sub_ijk$j==J & sub_ijk$k==K])^2))   

S_T=c(S_T,sum((vec$value-sub_all)^2))

S_ABC=c(S_ABC,n*(sub_ijk$value[sub_ijk$i==I & sub_ijk$j==J & sub_ijk$k==K]-sub_all)^2)
  
}
}
}

S_e=sum(S_e);S_T=sum(S_T);S_ABC=sum(S_ABC)

S_abr=S_ABC-S_a-S_b-S_r-S_ab-S_br-S_ar

pthi_e=(n-1)*n^3

#H:mu_ijk=mu+a+b+r+ab+br+ar

f=(S_abr/((n-1)^3))/(S_e/((n-1)*n^3))

qf(1-0.05,(n-1)^3,(n-1)*n^3)

#H:ab=0

f=(S_ab/((n-1)^2))/(S_e/pthi_e)

qf(1-0.05,(n-1)^2,pthi_e)

#H:ar=0

f=(S_ar/((n-1)^2))/(S_e/pthi_e)

qf(1-0.05,(n-1)^2,pthi_e)

#H:br=0

f=(S_br/((n-1)^2))/(S_e/pthi_e)

qf(1-0.05,(n-1)^2,pthi_e)

#H:r=0

f=(S_r/(n-1))/(S_e/pthi_e)

qf(1-0.05,(n-1),pthi_e)


r_e=n/(n-1+n-1+n-1+3*(n-1)^2+(n-1)^3+1)



```

```{r}

#p.96

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

X=as.matrix(data[,colnames(data) %in% c("x1","x2")])

Y=data$y

c_hat=solve(t(X)%*%X)%*%t(X)%*%Y

intercept=mean(Y)-sum(c_hat*apply(X,2,mean))

s2=sum((Y-X%*%c_hat+intercept)^2)/(nrow(X)-ncol(X))

sigma_c=s2*solve(t(X)%*%X)

sigma_c=sigma_c/s2



#要確認

c(c_hat)%*%solve(sigma_c)%*%t(t(c_hat))/(length(c_hat)*s2)

qf(1-0.05,length(c_hat),nrow(X)-ncol(X))



```

```{r}

# p.101 一元配置における多重比較

a=4;n=100

Y=array(0,dim=c(a,n))

for(j in 1:nrow(Y)){
  
Y[j,]=rnorm(n,sample(1:10,1))  
  
}

alpha=apply(Y,1,mean)

p=(alpha-mean(Y))/sqrt(sum((alpha-mean(Y))^2))

sigma=var(c(Y-alpha))

n*(p%*%t(t(alpha)))^2/((a-1)*sigma)

qf(1-0.05,a-1,a*(n-1))

#信頼区間

alpha-sqrt((a-1)*sigma*qf(1-0.05,a-1,a*(n-1))/n)

alpha+sqrt((a-1)*sigma*qf(1-0.05,a-1,a*(n-1))/n)



```

```{r}

#一元配置での各サンプル数が異なる場合

a=20

r=sample(c(50:100),a)

ave=c();var=c()

for(j in 1:a){
  
y=rnorm(r[j],mean=sample(1:10,1))

ave=c(ave,mean(y));var=c(var,var(y)*r[j]/(r[j]-1))
  
}

#Bartlett test

s2=var

S=s2*(r-1)

sigma_hat=sum(S)/sum(r-1)

N=sum(r)

likelihood_ratio=(N-a)*(log(sum(S)/(N-a))-sum((r-1)*log(s2))/(N-a))

qchisq(1-0.01,a-1)


#c1の場合(標本平均誤差に差があるか)

m=5;n=3

statics=(sum(ave[1:m])/m-sum(ave[(a-n+1):a])/n)*(sum(1/r[1:m])/(m^2)+sum(1/r[(a-n+1):a])/(n^2))^(-1/2)

statics/((a-1)*sigma_hat)>qf(1-0.05,a-1,sum(r)-a)



#c2の場合(重み付き平均誤差に差があるか)

statics=(sum(r[1:m]*ave[1:m])/sum(r[1:m])-sum(r[(a-n+1):a]*ave[(a-n+1):a])/sum(r[(a-n+1):a]))*(1/sum(r[1:m])+1/sum(r[(a-n+1):a]))^(-1)

statics/((a-1)*sigma_hat)>qf(1-0.05,a-1,sum(r)-a)



#p.112 sceffeの方法拡張２


a=20

r=sample(c(50:100),a)

ave=c();var=c()

for(j in 1:a){
  
y=rnorm(r[j],mean=sample(1:10,1))

ave=c(ave,mean(y));var=c(var,var(y)*r[j]/(r[j]-1))
  
}

#Bartlett test

s2=var

S=s2*(r-1)

sigma_hat=sum(S)/sum(r-1)

N=sum(r)

likelihood_ratio=(N-a)*(log(sum(S)/(N-a))-sum((r-1)*log(s2))/(N-a))

qchisq(1-0.01,a-1)

r_plus=t(matrix(rep(r,length(r)),ncol=length(r)))+c(r)

dist=as.matrix(dist(ave)^2)*(t(t(r))%*%t(r))/r_plus

dist/((a-1)*sigma_hat)>qf(1-0.01,a-1,sum(r)-a)


#p.130 tukey method

n=100

y=rnorm(n,mean=sample(1:10,1))

R=max(y)-min(y)

#student化された範囲

q=R/sd(y)







```


```{r}

#p.156

#t-test

a=100;r=10

y=array(0,dim=c(a,r))

for(j in 1:r){
  
y[,j]=rnorm(a,sample(1:10,1),sample(1:10,1))  
  
}

mu=apply(y,2,mean)

sigma=sum((y-apply(y,1,mean))^2)/(a*(r-1))

#非心度

abs(sqrt(r)*sum(-c(1:a)*(mu-mean(y)))/sqrt(sigma*(a*(a^2-1))/12))


alpha=0.01

qt(1-alpha,df=a*(r-1))

```

```{r}

#score method

#H0:mu1=mu2=...=mu

a=10;r=30

y=array(0,dim=c(a,r))

for(j in 1:r){
  
y[,j]=rnorm(a,sample(1:10,1),sample(1:10,1))

}

mu=apply(y,1,mean)

sigma=sum((y-apply(y,1,mean))^2)/(a*(r-1))


C1=-sqrt((c(1:a)-1)*(1-(c(1:a)-1)/a))+sqrt(c(1:a)*(1-c(1:a)/a))

C=-c(1:a)+(a+1)/2

C1=-sqrt((c(1:a)-1)*(1-(c(1:a)-1)/a))+sqrt(c(1:a)*(1-c(1:a)/a))

cor=cor(C,mu);cor1=cor(C1,mu)

#非心度

gamma=r*sum((mu-mean(mu))^2)*(cor^2)/sigma

t=sqrt(r)*sum(C*apply(y,1,mean))/sqrt(sum(C^2)*sigma)

qnorm(1-0.05/2)

2/(2+log(a-1))

#student化された範囲

sqrt(r)*as.matrix(dist(mu))/sigma




#H3:p.166

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

for(j in 1:a){
for(i in 1:a){  
  
if(i<j){  
  
C_ij[i,j]=sqrt((i*(a-j))/(j*(a-i))) 

}
}
}
  
pthi_0=(a-1)^2/(a-1+2*sum(C_ij^2))

statics=c()

val=c();coef=c()

for(i in 1:(a-1)){
  
statics=c(statics,r*i*(a-i)*(sum(mu[1:i])/i-sum(mu[(i+1):a])/(a-i))^2/(a*(a-1)*sigma_hat))

val=c(val,(sum(mu[1:i])/i-sum(mu[(i+1):a])/(a-i))^2)

coef=c(coef,r*i*(a-i)/(a*(a-1)*sigma_hat))

}

sum(statics)>qf(1-0.01,round(pthi_0),a*(r-1))

cumsum(mu)/c(1:a)

```

```{r}


#Bartholomew 尤度比検定 p.168

r=c();n=100

for(j in 1:n){
  
r=c(r,sample(c(1:20),1))  
  
}

y=array(0,dim=c(n,max(r)))

for(j in 1:n){
  
vec=c(rnorm(r[j],sample(1:10,1),sample(1:10,1)))

vec=c(vec,rep(NA,max(r)-r[j]))

y[j,]=vec
  
}

#H2

y_ave=c()

for(j in 1:n){
  
vec=y[j,]
  
vec=vec[is.na(vec)!=T]

y_ave=c(y_ave,mean(vec))

}

mu=y_ave

times=100000

for(k in 1:times){

for(j in length(mu):2){
  
mu_vec=mu[(j-1):(j)];r_vec=r[(j-1):j]  

if(mu_vec[1]<mu_vec[2]){
  
mu_vec=rep((r_vec[1]*mu_vec[1]+r_vec[2]*mu_vec[2])/sum(r_vec),2)

mu[(j-1):j]=mu_vec
  
}

}  
  
}

kai2_1=sum(r*(mu-mean(y[is.na(y)!=T]))^2)

l1=length(unique(floor(mu*10^6)/10^6))

f1=kai2_1/((l1-1)*sigma)

alpha=0.01

qf1=qf(1-alpha,a*max(r)-l,a*(max(r)-1))

mu1=mu



#H3

mu=y_ave

times=100000

for(k in 1:times){

for(j in 2:length(mu)){
  
mu_vec=mu[(j-1):(j)];r_vec=r[(j-1):j]  

if(mu_vec[1]>mu_vec[2]){
  
mu_vec=rep((r_vec[1]*mu_vec[1]+r_vec[2]*mu_vec[2])/sum(r_vec),2)

mu[(j-1):j]=mu_vec
  
}

}  
  
}

kai2_2=sum(r*(mu-mean(y[is.na(y)!=T]))^2)

l2=length(unique(floor(mu*10^6)/10^6))

f2=kai2_2/((l2-1)*sigma)

alpha=0.01

qf2=qf(1-alpha,a*max(r)-l,a*(max(r)-1))

mu2=mu

```

```{r}

#p.172

a=10;r=200

y=array(0,dim=c(a,r))

for(j in 1:a){

y[j,]=rnorm(r,2,sample(1:10,1))

}

sd=apply(y,1,sd)

var=sd^2

y_bar=apply(y/c(var),2,sum)/sum(1/var)

kai2=apply(t((t(y)-y_bar)^2)/var,2,sum)

qchisq(1-0.05,a-1)





```

```{r}

#p.186

b=200

y=3+rnorm(b,2,3)

mu_B=mean(y)

r=100

epsilon=matrix(rnorm(r*b),ncol=r)

X_B=epsilon+c(X_B)

y=epsilon+c(y)

mu_B=mean(y)

E_y=mu_B*rep(1,b)

S_e=sum((y-c(apply(y,1,mean)))^2)

n=b*r

S_b=r*sum((apply(y,1,mean)-mean(y))^2)

f=(S_b/(b-1))/(S_e/(n-b))

qf(1-0.01,b-1,n-b)



sigma_hat=S_e/(n-b)

sigma_B_hat=(S_b/(b-1)-sigma_hat)/r

V_sigma_hat_B=(2/(r^2))*((r*sigma_B_hat+sigma_hat)^2/(b-1)+sigma_hat^2/(n-b))


#confidence interval sigma_B/sigma p.190

under=(((S_b/(b-1))/(S_e/(n-b)))/qf(1-0.05,n-b,b-1)-1)/r

upper=(((S_b/(b-1))/(S_e/(n-b)))*qf(1-0.05,b-1,n-b)-1)/r




```

```{r}

#2元配置 p.190

r=10;a=5;b=6;n=r*a

pi=rpois(a,lambda=5);pj=rpois(b,lambda=9)

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

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

mat[i,j]=rpois(1,pi[i])+rpois(1,pj[j])+rpois(1,pi[i])*rpois(1,pj[j]) 

}
}

alpha=apply(mat,1,mean)-mean(mat)

beta=apply(mat,2,mean)-mean(mat)

alpha_beta=t(t(mat-c(alpha))-c(beta))+mean(mat)

mu

```