教育と心理のための推計学

```{r}

#p.105

#12.3 連関の強さ

mat=matrix(c(60,2,20,18),ncol=2)

colnames(mat)=c("神経質でない","神経質")

rownames(mat)=c("両親ある","両親なし")

#連関の強さ（<1）

pthi=(mat[1,1]*mat[2,2]-mat[1,2]*mat[2,1])/sqrt((mat[1,1]+mat[2,1])*(mat[1,2]+mat[2,2])*(mat[1,1]+mat[1,2])*(mat[2,1]+mat[2,2]))

t=sum(mat)*log(sum(mat))+sum(mat*log(mat))-sum(apply(mat,2,sum)*log(apply(mat,2,sum)));t=t/sum(mat)

#12.4 順序尺度における相関（ケンドール係数）

#例題1

teacher_A=c(3,2,5,4,7,6,8,1)

teacher_B=c(1,3,4,7,6,5,8,2)

pthi=function(a,b,c,d){

return(ifelse((a-b)*(c-d)>0,1,-1))

}

x=teacher_A;y=teacher_B

n=length(teacher_A)

K=c()

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

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

K=c(K,pthi(x[i],x[j],y[i],y[j]))

}

}

rk=2*sum(K)/(n*(n-1))

#独立性の検定(x,y)

K_star=sum(K)/sqrt(n*(n-1)*(2*n+5)/18)

qnorm(1-0.05/2)


```

```{r}

library(dplyr)

n=100

x=sample(0:1,n,prob=c(0.3,0.7),replace=T)

k=sum(x)

p=k/n

v1=2*(k+1)

v2=2*(n-k)

v1_dash=2*(n-k+1)

v2_dash=2*k

pu=v1*qf(1-0.05,v1,v2)/(v1*qf(1-0.05,v1,v2)+v2)

pl=v2_dash/(v1_dash*qf(1-0.05,v1_dash,v2_dash)+v2_dash)

#中央値（50％）、もしくは～％の値での信頼区間(N>25)

n=100;P=0.5

Ru=n*P+0.5+qnorm(1-0.01/2)*sqrt(n*P*(1-P))

Rl=n*P+0.5-qnorm(1-0.01/2)*sqrt(n*P*(1-P))

#平均値の信頼限界

n=100;x=rpois(n,3)

up=mean(x)+qt(1-0.05/2,df=n-1)*sd(x)/sqrt(n-1)

und=mean(x)-qt(1-0.05/2,df=n-1)*sd(x)/sqrt(n-1)

#母集団が有限の時

N=200

up=mean(x)+qt(1-0.05/2,n-1)*sd(x)*sqrt((N-n)/(N-1))/sqrt(n-1)

und=mean(x)-qt(1-0.05/2,n-1)*sd(x)*sqrt((N-n)/(N-1))/sqrt(n-1)

#第17章　比と中央値

#H:P=Q=0.5(n>25) p.166

n=100

x=sample(0:1,n,prob=c(0.3,0.7),replace=T)

k=sum(x)

p=k/n

CR=(abs(2*k-n)-1)/sqrt(n)

#p.167

mat=matrix(c(111,95,84,122),ncol=2)

colnames(mat)=c("正常成人","精神薄弱者")

rownames(mat)=c("既婚","独身")

#結婚率

P_hat=apply(mat,1,sum)[1]/sum(mat)

p=mat[1,]/apply(mat,2,sum)

sigma_p1=P_hat*(1-P_hat)/apply(mat,2,sum)[1]

sigma_p2=P_hat*(1-P_hat)/apply(mat,2,sum)[2]

CR=(p[1]-p[2])/sqrt(sigma_p1+sigma_p2)






#p.234 共分散分析

library(dplyr)

k=6;n=10

XA=c(108,136,138,159,146,141,175,149,174,176)

YA=c(73,102,118,104,81,107,100,87,117,111)

XB=c(99,117,90,141,106,112,110,117,111,122)

YB=c(98,74,56,111,95,88,82,77,86,92)

XC=c(194,198,196,198,210,196,230,222,220,228)

YC=c(94,79,96,98,102,102,108,91,120,105)

XD=c(165,164,161,159,175,135,132,190,145,142)

YD=c(90,76,90,64,86,51,72,90,95,78)

XE=c(124,95,116,112,123,110,137,105,135,126)

YE=c(107,95,97,80,98,74,74,67,89,58)

XF=c(140,177,189,142,216,200,255,173,153,160)

YF=c(49,82,73,86,81,97,106,70,61,82)

X=as.matrix(cbind(XA,XB,XC,XD,XE,XF))

Y=as.matrix(cbind(YA,YB,YC,YD,YE,YF))


KX=sum(X)^2/(k*n)

KY=sum(Y)^2/(k*n)

KXY=sum(X)*sum(Y)/(k*n)

#全体の平方和と積和を求める

SStX=sum(X^2)-KX

SStY=sum(Y^2)-KY

SPt=sum(X*Y)-KXY

#級間内の平方和と積和を求める

SSbX=sum(apply(X,2,sum)^2)/n-KX

SSbY=sum(apply(Y,2,sum)^2)/n-KY

SPb=sum(apply(X,2,sum)*apply(Y,2,sum))/n-KXY

#級内平方和と積和を求める

SSwX=SStX-SSbX

SSwY=SStY-SSbY

SPw=SPt-SPb

#全体：

SStYt=SStY-(SPt)^2/SStX

#級内（誤差）:

SSwYt=SSwY-(SPw)^2/SSwX

#級間（誤差）：

SSbYt=SStYt-SSwYt

#級内相関

rw=SPw/sqrt(SSwX*SSwY)

#級間相関

rb=SPb/sqrt(n*sum((apply(X,2,mean)-mean(X))^2)*n*sum((apply(Y,2,mean)-mean(Y))^2))

#全体相関

rt=SPt/sqrt(sum((X-mean(X))^2)*sum((Y-mean(Y))^2))

#f-test(H:bA=bB=bC=bD=bE=bF(one-way-anova))

f=(SSbYt/(k-1))/(SSwYt/(n*k-k-1))

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

#bw:

bw=SPw/SSwX

#調整平均値

Y_bar=apply(Y,2,mean)-bw*(apply(X,2,mean)-mean(X))


#各調整平均の組み合わせにおける有意差と予測値

t_mat=array(0,dim=c(k,k))

Y_hat=array(0,dim=c(n,k))

for(i in 1:n){
for(j in 1:k){
  
Y_hat[i,j]=apply(Y,2,mean)[j]+bw*(X[i,j]-apply(X,2,mean)[j])
  
}}

for(i in 1:k){
for(j in 1:k){  
  
t_mat[i,j]=(Y_bar[i]-Y_bar[j])/(sqrt(SSwYt/(n*k-k-1))*sqrt(1/n+1/n+(apply(X,2,mean)[i]-apply(X,2,mean)[j])^2/SSwX))  

}}

#推定誤差平方和

SSwY_hat=sum((Y-Y_hat)^2)

#p.332

trend=c();SSwY2=c()

for(s in 1:k){
  
sswy=sum((Y[,s]-mean(Y[,s]))^2)

spw=sum(X[,s]*Y[,s])-sum(X[,s])*sum(Y[,s])/n

sswx=sum((X[,s]-mean(X[,s]))^2)
  
SSwY2=c(SSwY2,sswy-spw^2/sswx)

sp=sum(X[,s]*Y[,s])-sum(X[,s])*sum(Y[,s])/n

ssx=sum(X[,s]^2)-sum(X[,s])^2/n

trend=c(trend,sp/ssx)

}



#f-test

((SSwYt-sum(SSwY2))/((n*k-k-1)-k*(n-2)))/(sum(SSwY2)/(k*(n-2)))

qf(1-0.05/2,(n*k-k-1)-k*(n-2),k*(n-2))


#H:sigA=sigB=...=sigF

MS=sum(SSwY2)/(k*(n-2))

kai2=(2.3026)*(n-1)*(k*log10(MS)-sum(log10(SSwY2/(n-2))))

qchisq(1-0.05,k-1)



```

```{r}

#要因が2つ以上の場合

mat_x=matrix(c(1,3,3,2,6,1,1,4,0,4,0,0),ncol=3)

mat_y=matrix(c(1,2,4,1,4,2,2,3,2,5,0,1),ncol=3)

tx=apply(mat_x,1,sum)

X2=apply(mat_x^2,1,sum)

ty=apply(mat_y,1,sum)

Y2=apply(mat_y^2,1,sum)

KX=sum(mat_x)^2/prod(dim(mat_x))

KY=sum(mat_y)^2/prod(dim(mat_y))

KXY=sum(mat_x)*sum(mat_y)/prod(dim(mat_x))

SPt=sum(mat_x*mat_y)-KXY

#A要因の平方和・積和：

n=nrow(mat_x)

SSAX=sum(apply(mat_x,2,sum)^2/nrow(mat_x))-KX

SSAY=sum(apply(mat_y,2,sum)^2/nrow(mat_y))-KY

SPA=sum(apply(mat_x,2,sum)*apply(mat_y,2,sum))/n-KXY

#B要因の平方和・積和：

SSBX=sum(apply(mat_x,1,sum)^2)/ncol(mat_x)-KX

SSBY=sum(apply(mat_y,1,sum)^2)/ncol(mat_y)-KY

SPB=sum(apply(mat_x,1,sum)*apply(mat_y,1,sum))/ncol(mat_x)-KXY

#誤差の平方和・積和

SStX=sum(mat_x^2)-KX

SSeX=SStX-SSAX-SSBX

SStY=sum(mat_y^2)-KY

SSeY=SStY-SSAY-SSBY

SPe=SPt-SPA-SPB

#推定の誤差平方和

SSeY_hat=SSeY-(SPe^2)/SSeX

SSAeY_hat=(SSAY+SSeY)-((SPA+SPe)^2/(SSAX+SSeX))

SSBeY_hat=(SSBY+SSeY)-((SPB+SPe)^2/(SSBX+SSeX))

MSeY=SSeY_hat/((dim(mat_y)[1]-1)*(dim(mat_y)[2]-1)-1)

SSAY_hat=SSAeY_hat-SSeY_hat

SSBY_hat=SSBeY_hat-SSeY_hat

fA=(SSAY_hat/(ncol(mat_x)-1))/MSeY

qf(1-0.05,ncol(mat_x)-1,((dim(mat_y)[1]-1)*(dim(mat_y)[2]-1)-1))

fB=(SSBY_hat/(nrow(mat_x)-1))/MSeY

qf(1-0.05,nrow(mat_x)-1,((dim(mat_y)[1]-1)*(dim(mat_y)[2]-1)-1))



```