集団遺伝学の数学的理論　木村資夫先生

#量的形質の分散分析

P=matrix(c(0.5,0.36,0.44,0.12,0.23,0.32,0.56,0.98,0.36,0.54,0.75,0.48,0.15,0.48,0.84,0.52),ncol=4)

Y=matrix(c(50,42,40,15,28,40,50,85,45,42,60,40,20,50,80,60),ncol=4)

P=P+t(P);Y=Y+t(Y);P=P/sum(P)

alpha=c(1:ncol(P))

Y_bar=mean(Y)

#Broyden Quasi-Newton Algorithm p.213

f=function(x){

alpha=x  

alpha_mat=t(matrix(rep(alpha,length(alpha)),ncol=length(alpha)))  

mat=P*(Y-mean(Y))

f=c()

for(k in 1:length(alpha)){

f=c(f,sum(mat[k,]-P[k,]*(alpha[k]+alpha))+sum(mat[,k]-P[,k]*(alpha[k]+alpha)))  

}
return(c(f))

}

X=alpha

ite=100

H=diag(f(X))

for(l in 1:ite){

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

X_pre=X  

X=X-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))))

}

}  


predicts=t(matrix(rep(X,length(X)),ncol=length(X)))+c(X)+Y_bar

d=Y-predicts

a=apply(P*(Y-Y_bar),1,sum)/apply(P,1,sum)

#遺伝的分散

Vg=2*alpha*a*apply(P,1,sum)

#優性分散

Vd=sum(P*d^2)

sita=P/(t(t(c(apply(P,1,sum))))%*%t(c(apply(P,1,sum))))