実験計画法の基礎　サイエンス社　石井吾郎先生

```{r}

# p.15
# 二次曲線のあてはめ

t=c(0:4)
x=c(12,14,23,34,47)
A=cbind(rep(1,length(t)),t,t^2)
param1=solve(t(A)%*%A)%*%t(A)%*%x
e1=x-A%*%param1
Se=sum(e1^2)
sig=Se/(length(t)-ncol(A))


# 帰無仮説：二次曲線でない
A1=A[,1:2]
param2=solve(t(A1)%*%A1)%*%t(A1)%*%x
e2=A1%*%param2
S2=c(x)%*%A%*%param1-c(x)%*%A1%*%param2
S1=c(x)%*%A1%*%param2
f0=S2/sig
f0>qf(1-0.05,1,length(t)-ncol(A))

```

```{r}

#p.18
#多重回帰

x=c(2,5,7,3,4,10)
t=c(5,12,15,8,7,20)
u=c(2,4,5,3,4,6)

A=cbind(rep(1,length(x)),t,u)
param=solve(t(A)%*%A)%*%t(A)%*%x
e=x-A%*%param
sig=sum(e^2)/(length(t)-ncol(A))

```

```{r}
#p.30 例2.4

A=matrix(c(rep(1,8),1,-1,1,-1,1,-1,1,-1,3,2,1,0,0,-1,-2,-3,0,0,2,1,1,2,0,0),nrow=8)

x=c(7,4,9,5,4,5,-4,-6)

param=solve(t(A)%*%A)%*%t(A)%*%x

A1=A[,1:2]

e=x-A%*%param

param1=solve(t(A1)%*%A1)%*%t(A1)%*%x
e1=A1%*%param1
Se=sum(e^2)
S1=c(x)%*%A1%*%param1
S2=c(x)%*%A%*%param-c(x)%*%A1%*%param1

sig=Se/(length(x)-ncol(A))
f0=(S2/(ncol(A)-ncol(A1)))/sig

f0>qf(1-0.01,ncol(A)-ncol(A1),length(x)-ncol(A))

```

```{r}

#p.32 一因子実験(一元配置)

mat=t(matrix(c(10,9,6,8,10,3,6,5,1,4,-3,0,1,-2,0,4,3,0,-1,NA),nrow=5))
rownames(mat)=paste0("A",c(1:4))

r=nrow(mat)
alpha=c()
ns=c()
for(j in 1:r){
vec=mat[j,]  
vec=vec[is.na(vec)!=T]
alpha=c(alpha,mean(vec))
ns=c(ns,length(vec))
}

mat2=mat-alpha
mat2[is.na(mat2)>0]=0
Se=sum(mat2^2)
Sa=sum(ns*alpha^2)-sum(ns)*(sum(alpha*ns)/sum(ns))^2
Smu=sum(ns)*(sum(alpha*ns)/sum(ns))^2

f0=(Sa/(r-1))/(Se/(sum(ns)-r))
#A1,A2,A3,A4には水準による違いがある
f0>qf(1-0.05,r-1,sum(ns)-r)


t_table=array(0,dim=c(r,r))

for(i in 1:r){
for(j in 1:r){
t_table[i,j]=(alpha[i]-alpha[j])/sqrt((1/ns[i]+1/ns[j])*(Se/(sum(ns)-r)))  
}  
}
t_table=abs(t_table)
qt(1-0.05/2,df=sum(ns)-r)

```

```{r}

#乱塊法　p.40

mat=matrix(c(1,5,2,4,3,3,1,2,5,4,2,4,5,1,3,4,2,1,3,5),nrow=5)
rownames(mat)=paste0("A",c(1:5))

n=prod(dim(mat))

mu=mean(mat)
alpha=apply(mat,1,mean)-mu
beta=apply(mat,2,mean)-mu

r=dim(mat)[1]
s=dim(mat)[2]

Smu=n*mu^2
Salpha=s*sum(alpha^2)
Sbeta=r*sum(beta^2)
Se=sum(mat^2)-Salpha-Sbeta-Smu

#H1の検定
(Salpha/(r-1))/(Se/((r-1)*(s-1)))>qf(1-0.05,r-1,(r-1)*(s-1))

#H2の検定
(Sbeta/(s-1))/(Se/((r-1)*(s-1)))>qf(1-0.05,s-1,(r-1)*(s-1))


```

```{r}

#p.48 ラテン方格

#処理水準(A1,A2,A3,A4),作業者(行),実験装置(列)

mat=matrix(paste0("A",c(1,4,3,2,2,3,4,1,4,2,1,3,3,1,2,4)),nrow=4)
x=matrix(c(-10,3,2,9,16,5,5,6,-2,12,0,-6,-1,-2,14,5),nrow=4)
dim=4


mu=mean(x)
beta=apply(x,1,mean)-mu
r=apply(x,2,mean)-mu
alpha=c()
for(j in sort(unique(c(mat)))){
alpha=c(alpha,sum(x[mat==j])/dim-mu)  
}

Smu=prod(dim(x))*mu^2
Sbeta=dim*sum(beta^2)
Sr=dim*sum(r^2)
Salpha=dim*sum(alpha^2)
Se=sum(x^2)-Smu-Sbeta-Salpha-Sr

#H:水準間での違いがない(alpha)
falpha=(Salpha/(dim-1))/(Se/(prod(dim(x))-3*(dim-1)-1))
#H:作業者での違いがない(beta)
fbeta=(Sbeta/(dim-1))/(Se/(prod(dim(x))-3*(dim-1)-1))
#H:実験装置での違いがない(r)
fr=(Sr/(dim-1))/(Se/(prod(dim(x))-3*(dim-1)-1))
#閾値
qf(1-0.05,dim-1,prod(dim(x))-3*(dim-1)-1)

```

```{r}

#2要因配置と交互作用

mat=matrix(c(6,8,8,16,10,12,12,10,14,18,10,8,10,8,18,22,14,12,0,2,0,8,6,8),nrow=6)
rownames(mat)=paste0("A",c(1,1,2,2,3,3))#品種
colnames(mat)=c(paste0("B",c(1:4)))#肥料

r=length(sort(unique(rownames(mat))))#品種の数
s=length(sort(unique(colnames(mat))))#肥料の数
t=2　#繰り返し数
n=prod(dim(mat))#標本数

A=sort(unique(rownames(mat)))#品種のラベル
B=sort(unique(colnames(mat)))#肥料のラベル

mu=sum(mat)/n
alpha=c()
for(j in A){
alpha=c(alpha,sum(mat[rownames(mat)==j,])/(s*t))
}
beta=c()
for(j in B){
beta=c(beta,sum(mat[,colnames(mat)==j])/(r*t))  
}
rij=c()
rnames=c()
for(j in 1:length(A)){
for(i in 1:length(B)){
rij=c(rij,sum(mat[rownames(mat)==A[j],colnames(mat)==B[i]])/t-alpha[j]-beta[i]+mu)
rnames=c(rnames,paste0("r",j,i))
}}
names(rij)=rnames
alpha=alpha-mu
beta=beta-mu

Smu=n*mu^2
Salpha=s*t*sum(alpha^2)
Sbeta=r*t*sum(beta^2)
Sr=t*sum(rij^2)
Se=sum(c(mat)^2)-Salpha-Sbeta-Sr-Smu

p=0.01
#H:品種間での相違がない
falpha=(Salpha/(r-1))/(Se/(s*r*(t-1)))
qf(1-p,r-1,r*s*(t-1))
#H:肥料間での相違がない
fbeta=(Sbeta/(s-1))/(Se/(s*r*(t-1)))
qf(1-p,s-1,r*s*(t-1))
#H:交互作用での相違がない
fr=(Sr/((r-1)*(s-1)))/(Se/(s*r*(t-1)))
qf(1-p,(r-1)*(s-1),r*s*(t-1))


```

```{r}

#p.74 4.直交多項式と反応曲面 
#p.27 参照

#2次曲線に対するあてはめ
u=c(2,3,4,5,6)
x=c(1,1,3,4,8)
A=cbind(rep(1,length(u)),u,u^2)
sita=solve(t(A)%*%A)%*%t(A)%*%x

#1次直線へのあてはめ
B=cbind(rep(1,length(u)),u)
sita1=solve(t(B)%*%B)%*%t(B)%*%x

Sab=x%*%B%*%sita1
Sr=x%*%A%*%sita-x%*%B%*%sita1
Se=sum(x^2)-x%*%A%*%sita

dim_ab=ncol(A)-1
dim_r=ncol(B)-1
dim_e=length(x)-dim_ab-dim_r

f=(Sr/dim_r)/(Se/dim_e)
quantile=qf(1-0.05,dim_r,dim_e)

if(f>quantile){
print("帰無仮説：r=0を棄却する")  
}

```

```{r}

#p.82 反応曲面

t=c(1,1,1,0,0,0,1,-1,-1)
u=rep(c(1,0,-1),3)
x=c(5,3,0,6,8,3,4,3,1)
A=cbind(rep(1,9),t,u,t^2,u^2,t*u)
param=solve(t(A)%*%A)%*%t(A)%*%x

n=nrow(A)

print((sum(x^2)-c(x)%*%A%*%param)/(n-ncol(A)))
sig=sum((x-A%*%param)^2)/(n-ncol(A))

sita1=param[2]
sita2=param[3]
sita11=param[4]
sita22=param[5]
sita12=param[6]

#極大値EXのtvalとuval
val=solve(matrix(c(2*sita11,sita12,sita12,2*sita2),nrow=2))%*%c(-sita1,-sita2)
tval=val[1]
uval=val[2]

vec=c(1,tval,uval,tval^2,uval^2,tval*uval)
Ex=sum(vec*param)
VEx=sig*vec%*%solve(t(A)%*%A)%*%vec



```

```{r}


#p.87 主効果・交互作用と直行多項式の関係
mat1=matrix(c(26,28,30,32,34,30,30,38,34,20,20,26),nrow=3)
mat2=matrix(c(28,36,32,30,38,18,28,42,32,22,28,28),nrow=3)
mat1=(mat1-20)/2
mat2=(mat2-20)/2

n=prod(dim(mat1))+prod(dim(mat2))

mat=mat1+mat2

A=nrow(mat)
B=ncol(mat)


xi=apply(mat,1,mean)/2
xj=apply(mat,2,mean)/2
mu=sum(mat)/n

alpha=xi-mu
beta=xj-mu

r=array(0,dim=c(length(alpha),length(beta)))
for(i in 1:length(alpha)){
for(j in 1:length(beta)){
r[i,j]=mat[i,j]/2-mu-alpha[i]-beta[j]
}}

n1=prod(dim(mat1))
n2=prod(dim(mat2))

delta1=sum(mat1)/n1-mu
delta2=sum(mat2)/n2-mu

Smu=n*mu^2
Salpha=2*A*sum(alpha^2)
Sbeta=2*B*sum(beta^2)
Sdelta=prod(dim(mat))*(delta1^2+delta2^2)
Sr=2*sum(r^2)
Se=sum(mat1^2)+sum(mat2^2)-Smu-Sr-Salpha-Sbeta-Sdelta





```

```{r}

#直交多項式による反応曲面の推定　p.89
pthi_t=function(val,s){
if(s==1){val=rep(1,length(val))}  
if(s==2){val=val}
if(s==3){val=3*val^2-2}
return(val)  
}

pthi_u=function(val,s){
if(s==1){val=rep(1,length(val))}  
if(s==2){val=2*val}
if(s==3){val=val^2-5/4}
return(val) 
}

mat1=matrix(c(26,28,30,32,34,30,30,38,34,20,20,26),nrow=3)
mat2=matrix(c(28,36,32,30,38,18,28,42,32,22,28,28),nrow=3)
mat1=(mat1-20)/2
mat2=(mat2-20)/2
x=c(unlist(t(mat1)),unlist(t(mat2)))

t=c(-1,0,1)
u=c(-1.5,-0.5,0.5,1.5)

pthi_t_mat=array(0,dim=c(length(t),length(t)))
pthi_u_mat=array(0,dim=c(length(u),3))

for(i in 1:ncol(pthi_t_mat)){
pthi_t_mat[,i]=pthi_t(t,i)  
}

for(i in 1:ncol(pthi_u_mat)){
pthi_u_mat[,i]=pthi_u(u,i)  
}

R=kronecker(pthi_t_mat,pthi_u_mat)
A=rbind(cbind(rep(1,nrow(R)),rep(0,nrow(R)),R),cbind(rep(0,nrow(R)),rep(1,nrow(R)),R))



#修正ガウス・ザイデル法
AA=t(A)%*%A　#行列のサンプル
C=t(A)%*%x  #ベクトルのサンプル
rs=C/diag(AA)　#r_iを計算
d=AA/c(diag(AA))　#d_{i,j}を計算
diag(d)=0 #d_{i,j}を計算
mat=cbind(rs,-d) #x_iの更新に必要な行列を作成
ite=30000　#反復計算回数
vec=rep(1,ncol(AA))　#x_iの初期値
eta=0.001
for(l in 1:ite){
vecpre=vec  
z=c(1,vec)
vec[1:2]=vec[1:2]-mean(vec[1:2])
vec=vec*(1-eta)+eta*mat%*%z
print(sum(abs(vecpre-vec)))
}　


alpha_num=c()
for(i in c(0:(ncol(pthi_t_mat)-1))){
alpha_num=c(alpha_num,rep(i,ncol(pthi_u_mat)))
}

delta=vec[1:2]
print(paste0("deltaの値は",sum(mat1)/12-vec[3],"と",sum(mat2)/12-vec[3]))
mu=vec[3]
beta_num=rep(c(0:(ncol(pthi_u_mat)-1)),ncol(pthi_u_mat))
num=data.frame(num1=alpha_num,num2=beta_num,param=vec[-c(1:2)])
alpha=num$param[num$num1!=0 & num$num2==0]
names(alpha)=paste0(num$num1[num$num1!=0 & num$num2==0],num$num2[num$num1!=0 & num$num2==0])
beta=num$param[num$num1==0 & num$num2!=0]
names(beta)=paste0(num$num1[num$num1==0 & num$num2!=0],num$num2[num$num1==0 & num$num2!=0])
r=num$param[num$num1!=0 & num$num2!=0]
names(r)=paste0(num$num1[num$num1!=0 & num$num2!=0],num$num2[num$num1!=0 & num$num2!=0])

df_mu=1
df_delta=length(delta)-1
df_alpha=length(alpha)-1
df_beta=length(beta)-1
df_r=prod(dim(r)-1)
n=length(x)


Sij=(t(A)%*%x*vec)[-c(1:2)]
Se=sum(x^2)-as.numeric(c(x)%*%A%*%vec)
Sdelta=as.numeric(c(x)%*%A%*%vec)-sum(Sij)  
print(paste0("Sdeltaの値は",sum(c(prod(dim(mat1)),prod(dim(mat2)))*vec[1:2]^2)))


```

```{r}

#一対比較　p.95
#p.99のデータの作り方が重要
#N=4人の人に評価してもらった場合

mat=matrix(c(0,0.5,0.25,1,-0.25,0,0,0.75,0.75,-1,0,0,0,0.5,-0.5,0,0,0,0,-1,0,0,0,0,0),nrow=5)

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


```

```{r}

#p.100 一対比較における個人差
#例5.2

t=3
N=4

#k=1:Aさん
mat1=matrix(c(0,-2,-1,2,0,1,1,-1,0),nrow=3)

#k=2:Bさん
mat2=matrix(c(0,-2,0,2,0,2,0,-2,0),nrow=3)

#k=3:Cさん
mat3=matrix(c(0,-2,-2,2,0,1,2,-1,0),nrow=3)

#k=4:Dさん
mat4=matrix(c(0,-1,0,1,0,0,0,0,0),nrow=3)

x=list(A=mat1,B=mat2,C=mat3,D=mat4)

alpha=rep(0,t)
x2=0
for(j in 1:N){
alpha=alpha+apply(x[[j]],1,mean)  
x2=x2+sum(x[[j]]^2)
}
alpha=alpha/N

beta=array(0,dim=c(t,N))
for(j in 1:N){
beta[,j]=apply(x[[j]],1,mean)-alpha
}

alpha_ik=beta+alpha

Salpha=t*N*sum(alpha^2)
Sbeta=(t-1)*(N-1)*sum(beta^2)/2
Se=x2/2-Salpha-Sbeta

#H:alphaが等しい
f=(Salpha/(t-1))/(Se/(N*(t-1)*(t-2)/2))
qf(1-0.05,t-1,N*(t-1)*(t-2)/2)



```

```{r}

#p.105 順序効果による分散分析
t=4

x1=matrix(c(0,1,-1,0,2,0,-1,2,1,-1,0,1,-2,0,1,0),nrow=t)
x2=matrix(c(0,0,0,2,1,0,2,2,-1,1,0,0,1,-2,-1,0),nrow=t)
x3=matrix(c(0,1,1,2,0,0,1,2,-1,1,0,0,-1,-1,0,0),nrow=t)

x=x1
rowx=apply(x,2,sum)
colx=apply(x,1,sum)
sumx2=sum(x^2)
diag(x)=NA
x=c(x)
x=x[is.na(x)!=T]

vec=rep(0,t)

for(i in 1:t){
for(j in 1:t){
vector=rep(0,t)
if(i<j){
vector[i]=1
vector[j]=-1
vec=rbind(vec,vector)
}else{
if(i>j){
vector[j]=-1
vector[i]=1
vec=rbind(vec,vector)
}}
}}

vec=tail(vec,nrow(vec)-1)
A=cbind(vec,rep(1,nrow(vec)))

#修正ガウス・ザイデル法
AA=t(A)%*%A　#行列のサンプル
C=t(A)%*%x  #ベクトルのサンプル
rs=C/diag(AA)　#r_iを計算
d=AA/c(diag(AA))　#d_{i,j}を計算
diag(d)=0 #d_{i,j}を計算
mat=cbind(rs,-d) #x_iの更新に必要な行列を作成
ite=30000　#反復計算回数
vec=rep(1,ncol(AA))　#x_iの初期値
eta=0.001
for(l in 1:ite){
vecpre=vec  
z=c(1,vec)
vec[1:t]=vec[1:t]-mean(vec[1:t])
vec=vec*(1-eta)+eta*mat%*%z
print(sum(abs(vecpre-vec)))
}　

#推定量
alpha=vec[1:t]
delta=vec[length(vec)]

Sa=sum(alpha^2)*2*t
Sdelta=t*(t-1)*sum(delta^2)
Se=sumx2-Sa-Sdelta

#H:alphaはすべて等しい
f=(Sa/(t-1))/(Se/(t*(t-2)))
qf(1-0.01,t-1,t*(t-2))
#H:delta=0
f=(Sdelta)/(Se/(t*(t-2)))
qf(1-0.01,1,t*(t-2))


```

```{r}
#p.108 順序効果における個人差
#例5.4

t=4
N=3

#k=1:Aさん
mat1=matrix(c(0,1,-1,0,2,0,-1,2,1,-1,0,1,-2,0,1,0),nrow=t)

#k=2:Bさん
mat2=matrix(c(0,0,0,-2,1,0,2,2,-1,1,0,0,1,-2,-1,0),nrow=t)

#k=3:Cさん
mat3=matrix(c(0,1,1,2,0,0,1,2,-1,1,0,0,-1,-1,0,0),nrow=t)

x=list(A=mat1,B=mat2,C=mat3)

alpha=rep(0,t)
alpha_ik=array(0,dim=c(t,N))
delta=c()
x2=0
xk=c()
for(j in 1:N){
alpha=alpha+apply(x[[j]],1,sum)-apply(x[[j]],2,sum)  
alpha_ik[,j]=(apply(x[[j]],1,sum)-apply(x[[j]],2,sum))/(2*t)
x2=x2+sum(x[[j]]^2)
xk=c(xk,sum(x[[j]])^2)
delta=c(delta,sum(x[[j]])/(t*(t-1)))
}
alpha=alpha/(2*t*N)
delta_k=delta-mean(delta)
delta=mean(delta)
beta=alpha_ik-alpha

Sa=(2*t*N)*sum(alpha^2)
Sb=(2*t)*sum(beta^2)
Sdelta=t*(t-1)*N*delta^2
Sdelta_k=t*(t-1)*sum(delta_k^2)
print(paste0("Sdelta_kの値は",sum(xk)/(t*(t-1))-Sdelta))
Se=x2-Sa-Sb-Sdelta-Sdelta_k

#H:alphaは全て等しい
f=(Sa/(t-1))/(Se/(t*(t-2)*N))
qf(1-0.05,t-1,t*(t-2)*N)

#H:betaは全て等しい(主効果個人差なし)
f=(Sb/((N-1)*(t-1)))/(Se/(t*(t-2)*N))
qf(1-0.05,(N-1)*(t-1),t*(t-2)*N)

#H:delta=0
f=(Sdelta)/(Se/(t*(t-2)*N))
qf(1-0.05,1,t*(t-2)*N)

#H:delta_kは全て等しい(順序効果による個人差なし)
f=(Sdelta_k/(N-1))/(Se/(t*(t-2)*N))
qf(1-0.05,N-1,t*(t-2)*N)



```

```{r}

#p.110 交互作用、順序効果を含めたパラメータ推定（分散分析）

#p.164 Gauss法

t=4
N=3

#k=1:Aさん
mat1=matrix(c(0,1,-1,0,2,0,-1,2,1,-1,0,1,-2,0,1,0),nrow=t)

#k=2:Bさん
mat2=matrix(c(0,0,0,-2,1,0,2,2,-1,1,0,0,1,-2,-1,0),nrow=t)

#k=3:Cさん
mat3=matrix(c(0,1,1,2,0,0,1,2,-1,1,0,0,-1,-1,0,0),nrow=t)

X=list(A=mat1,B=mat2,C=mat3)

obs=c()
for(j in 1:N){
obs=c(obs,c(X[[j]]))
}

f=function(x){
a=matrix(x[1:(N*t)],nrow=t);x=x[-c(1:(N*t))]
paramr=x[1:((t-1)*(t-2))];x=x[-c(1:((t-1)*(t-2)))]
r=array(0,dim=c(t,t))
r[lower.tri(r)]=-abs(paramr)
r=r-t(r)
delta=array(0,dim=c(t,t))
delta[lower.tri(delta)]=c(x)
delta=delta+t(delta)
pre=c()
for(j in 1:N){
vec=t(matrix(rep(a[,j],t),nrow=t))-c(a[,j])+r+delta  
pre=c(pre,c(vec))  
}
return(pre)  
}

Y=c(obs)
h=0.001
lam=1
ite=10^4
sita=sample(1:9,t*N+((t-1)*(t-2)*2),replace=T)
eta=10^(-3)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(length(Y),length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*beta)

#パラメータに分割
alpha=matrix(sita[1:(N*t)],nrow=t);sita=sita[-c(1:(N*t))]
paramr=sita[1:((t-1)*(t-2))];sita=sita[-c(1:((t-1)*(t-2)))]
r=array(0,dim=c(t,t))
r[lower.tri(r)]=-abs(paramr)
rmat=(r-t(r))
delta=array(0,dim=c(t,t))
delta[lower.tri(delta)]=c(sita)
delta=delta+t(delta)

#平均0の調整
alpha=t(t(alpha)-apply(alpha,2,mean))
as=apply(alpha,1,mean)
r=array(0,dim=c(t,t))
for(i in 1:t){
for(j in 1:t){  
r[i,j]=rmat[i,j]-apply(rmat,1,mean)[i]-apply(rmat,2,mean)[j]-mean(rmat)  
}}
sita=c(c(alpha),c(-r[lower.tri(r)]),c(delta[lower.tri(delta)]))

print(sum((Y-f(sita))^2))
}





```