May 20 2017

心理統計法('17) その10

放送大学「心理統計法(‘17)」第12章

授業ではrstanパッケージを使ってますが、MCMCpackパッケージを使ってみます。初学者ですので間違いが多々あると思います。

（参考）
ベイジアンMCMCによる統計モデル
 Using the ggmcmc package
Bayesian Inference With Stan ~番外編~
Exercise 2: Bayesian A/B testing using MCMC Metropolis-Hastings
政治学方法論 I

1 2	library(MCMCpack) library(knitr)

カテゴリ数がkの比率の推測

ペット問題

 x <- c(32,29,18,15,10)
#
post <- MCmultinomdirichlet(x, c(1,1,1,1,1), mc=100000)
df<-data.frame(post[,1],post[,2],post[,3],post[,4],p5=post[,5])
#
ryou<-rbind(EAP=apply(df, 2, mean),post.sd=apply(df, 2, sd),
	apply(df, 2, quantile, probs=c(0.025,0.05, 0.5,0.95, 0.975)))
colnames(ryou)<-c("p1(犬)","p2(猫)","p3(魚)","p4(鳥)","p5(その他)")
kable(t(round(ryou,3)))

	EAP	post.sd	2.5%	5%	50%	95%	97.5%
p1(犬)	0.303	0.044	0.220	0.232	0.302	0.377	0.392
p2(猫)	0.275	0.043	0.195	0.207	0.274	0.348	0.363
p3(魚)	0.174	0.036	0.109	0.118	0.172	0.237	0.251
p4(鳥)	0.147	0.034	0.087	0.095	0.145	0.206	0.219
p5(その他)	0.101	0.029	0.052	0.058	0.098	0.152	0.164

mat<-matrix(rep(NA,25),ncol=5)
colnames(mat)<-rownames(mat)<-c("p1(犬)","p2(猫)","p3(魚)","p4(鳥)","p5(その他)")
for ( i in 1:5){
	for ( j in 1:5){
	  mat[i,j]<-round(sum((df[,i]-df[,j] > 0 )=="TRUE")/nrow(df),3)
	}
}
kable(mat)

	p1(犬)	p2(猫)	p3(魚)	p4(鳥)	p5(その他)
p1(犬)	0.000	0.649	0.975	0.993	1.000
p2(猫)	0.351	0.000	0.945	0.982	0.999
p3(魚)	0.025	0.055	0.000	0.695	0.932
p4(鳥)	0.007	0.018	0.305	0.000	0.838
p5(その他)	0.000	0.001	0.068	0.162	0.000

連言命題が正しい確率

u_p1>p2 　u_p2>p3 u_p3>p4

1	sum((df[,1]-df[,2] > 0 )=="TRUE" & (df[,2]-df[,3] > 0 )=="TRUE" & (df[,3]-df[,4] > 0 )=="TRUE")/nrow(df)

0.40009

u_p1>p3 　u_p2>p3 u_p1>p4 * u_p2>p4

1	sum((df[,1]-df[,3] > 0 )=="TRUE" & (df[,2]-df[,3] > 0 )=="TRUE" & (df[,1]-df[,4] > 0 )=="TRUE" & (df[,2]-df[,4] > 0 )=="TRUE")/nrow(df)

0.91187

対応ある2×2のクロス表の分析

ブランド認知問題2

x <- matrix(c(
       70,30,
       28,72),2,2,T)                         #反応数 
colnames(x)<-c("B認知","B非認知")
rownames(x)<-c("A認知","A非認知")
kable(x)

	B認知	B非認知
A認知	70	30
A非認知	28	72

#png("sinri12_1.png")
par(family="TakaoMincho")
mosaicplot(x,col=T,cex=1.5,main="2つのブランドの認知")
#dev.off()

mc<-1002000
burnin<-2000
thin <-10
y<-c(70,30,28,72)
alpha0<-c(1,1,1,1)
set.seed(1)
l<-MCmultinomdirichlet(y ,alpha0 ,mc=mc)
df<-data.frame(l[,1],l[,2],l[,3],l[,4])
names(df)<-c("p11","p12","p21","p22")
df<-df[(burnin+thin*(1:((mc-burnin)/thin))-(thin-1)),]
#
ryou<-rbind(EAP=apply(df, 2, mean),post.sd=apply(df, 2, sd),
	apply(df, 2, quantile, probs=c(0.025,0.05, 0.5,0.95, 0.975)))
kable(t(round(ryou,3)))

	EAP	post.sd	2.5%	5%	50%	95%	97.5%
p11	0.348	0.033	0.284	0.294	0.347	0.404	0.415
p12	0.152	0.025	0.106	0.113	0.151	0.195	0.205
p21	0.142	0.024	0.098	0.104	0.141	0.184	0.193
p22	0.358	0.034	0.293	0.303	0.357	0.414	0.425

独立と連関

#### ピアソン残差
#
p1n<-df$p11+df$p12
p2n<-df$p21+df$p22
pn1<-df$p11+df$p21
pn2<-df$p12+df$p22
#
e11<-(df$p11-p1n*pn1)/sqrt(p1n*pn1)
e12<-(df$p12-p1n*pn2)/sqrt(p1n*pn2)
e21<-(df$p21-p2n*pn1)/sqrt(p2n*pn1)
e22<-(df$p22-p2n*pn2)/sqrt(p2n*pn2)
#
#### クラメルの連関係数
#
V<-sqrt(e11^2+e12^2+e21^2+e22^2)
#
hyouka<-data.frame(V,e11,e12,e21,e22)
#
ryou<-rbind(EAP=apply(hyouka, 2, mean),post.sd=apply(hyouka, 2, sd),
	apply(hyouka, 2, quantile, probs=c(0.025,0.05, 0.5,0.95, 0.975)))
kable(t(round(ryou,3)))

	EAP	post.sd	2.5%	5%	50%	95%	97.5%
V	0.412	0.063	0.284	0.305	0.413	0.514	0.532
e11	0.208	0.034	0.141	0.152	0.208	0.264	0.275
e12	-0.203	0.032	-0.266	-0.256	-0.204	-0.150	-0.139
e21	-0.207	0.033	-0.270	-0.261	-0.208	-0.153	-0.142
e22	0.204	0.034	0.138	0.149	0.204	0.259	0.270

対応あるa×bのクロス表の分析

パスタ問題

x <- matrix(c(
       19,  9,  6,
       10, 19,  5,
       15, 14, 18),3,3,T)                    #反応数 
rownames(x)<-c("トマトの冷製","カルボナーラ","ペペロンチーノ")
colnames(x)<-c("バジル","トリュフ","なし")
kable(x)

	バジル	トリュフ	なし
トマトの冷製	19	9	6
カルボナーラ	10	19	5
ペペロンチーノ	15	14	18

#png("sinri12_2.png")
par(family="TakaoMincho")
mosaicplot(x,col=T,cex=1.5,main="パスタに選ばれたトッピング")
#dev.off()

mc<-1002000
burnin<-2000
thin <-10
y<-c(19,9,6,10,19,5,15,14,18)
alpha0<-c(1,1,1,1,1,1,1,1,1)
set.seed(1)
l<-MCmultinomdirichlet(y ,alpha0 ,mc=mc)
df<-data.frame(l[,1],l[,2],l[,3],l[,4],l[,5],l[,6],l[,7],l[,8],l[,9])
names(df)<-c("p11","p12","p13","p21","p22","p23","p31","p32","p33")
df<-df[(burnin+thin*(1:((mc-burnin)/thin))-(thin-1)),]
#
ryou<-rbind(EAP=apply(df, 2, mean),post.sd=apply(df, 2, sd),
	apply(df, 2, quantile, probs=c(0.025,0.05, 0.5,0.95, 0.975)))
kable(t(round(ryou,3)))

	EAP	post.sd	2.5%	5%	50%	95%	97.5%
p11	0.161	0.033	0.102	0.110	0.159	0.218	0.231
p12	0.081	0.024	0.040	0.045	0.078	0.124	0.134
p13	0.057	0.021	0.023	0.027	0.054	0.094	0.103
p21	0.088	0.025	0.045	0.051	0.086	0.134	0.144
p22	0.161	0.033	0.103	0.110	0.159	0.219	0.231
p23	0.048	0.019	0.018	0.022	0.046	0.084	0.092
p31	0.129	0.030	0.076	0.083	0.127	0.181	0.193
p32	0.121	0.029	0.070	0.077	0.119	0.172	0.184
p33	0.153	0.032	0.095	0.104	0.151	0.209	0.221

独立と連関

#### ピアソン残差
#
p1n<-df$p11+df$p12+df$p13
p2n<-df$p21+df$p22+df$p23
p3n<-df$p31+df$p32+df$p33
pn1<-df$p11+df$p21+df$p31
pn2<-df$p12+df$p22+df$p32
pn3<-df$p13+df$p23+df$p33
#
e11<-(df$p11-p1n*pn1)/sqrt(p1n*pn1)
e12<-(df$p12-p1n*pn2)/sqrt(p1n*pn2)
e13<-(df$p13-p1n*pn3)/sqrt(p1n*pn3)
e21<-(df$p21-p2n*pn1)/sqrt(p2n*pn1)
e22<-(df$p22-p2n*pn2)/sqrt(p2n*pn2)
e23<-(df$p23-p2n*pn3)/sqrt(p2n*pn3)
e31<-(df$p31-p3n*pn1)/sqrt(p3n*pn1)
e32<-(df$p32-p3n*pn2)/sqrt(p3n*pn2)
e33<-(df$p33-p3n*pn3)/sqrt(p3n*pn3)
#
#### クラメルの連関係数
#
# 3*3 だからmin(3,3)-1 = 2
V<-sqrt((e11^2+e12^2+e13^2+e21^2+e22^2+e23^2+e31^2+e32^2+e33^2)/2)
#
hyouka<-data.frame(V,p1n,p2n,p3n,pn1,pn2,pn3)
#
ryou<-rbind(EAP=apply(hyouka, 2, mean),post.sd=apply(hyouka, 2, sd),
	apply(hyouka, 2, quantile, probs=c(0.025,0.05, 0.5,0.95, 0.975)))
colnames(ryou)<-c("V","p1.","p2.","p3.","p.1","p.2","p.3")
kable(t(round(ryou,3)))

	EAP	post.sd	2.5%	5%	50%	95%	97.5%
V	0.253	0.059	0.138	0.155	0.252	0.351	0.370
p1.	0.298	0.041	0.222	0.233	0.297	0.368	0.381
p2.	0.298	0.041	0.222	0.233	0.297	0.367	0.382
p3.	0.403	0.044	0.319	0.332	0.403	0.476	0.490
p.1	0.379	0.043	0.296	0.308	0.378	0.451	0.465
p.2	0.363	0.043	0.281	0.293	0.362	0.435	0.450
p.3	0.258	0.039	0.185	0.196	0.257	0.325	0.339

ピアソン残差が0より大きい確率

mat1<-matrix(rep(NA,9),ncol=3)
mat1[1,1]<-round(sum((e11>0)=="TRUE")/length(e11),3)
mat1[1,2]<-round(sum((e12>0)=="TRUE")/length(e11),3)
mat1[1,3]<-round(sum((e13>0)=="TRUE")/length(e11),3)
mat1[2,1]<-round(sum((e21>0)=="TRUE")/length(e11),3)
mat1[2,2]<-round(sum((e22>0)=="TRUE")/length(e11),3)
mat1[2,3]<-round(sum((e23>0)=="TRUE")/length(e11),3)
mat1[3,1]<-round(sum((e31>0)=="TRUE")/length(e11),3)
mat1[3,2]<-round(sum((e32>0)=="TRUE")/length(e11),3)
mat1[3,3]<-round(sum((e33>0)=="TRUE")/length(e11),3)
rownames(mat1)<-c("トマトの冷製","カルボナーラ","ペペロンチーノ")
colnames(mat1)<-c("バジル","トリュフ","なし")	
kable(mat1)

	バジル	トリュフ	なし
トマトの冷製	0.992	0.072	0.113
カルボナーラ	0.101	0.997	0.045
ペペロンチーノ	0.128	0.111	0.994

ピアソン残差が0より小さい確率

mat2<-matrix(rep(NA,9),ncol=3)
mat2[1,1]<-round(sum((e11<0)=="TRUE")/length(e11),3)
mat2[1,2]<-round(sum((e12<0)=="TRUE")/length(e11),3)
mat2[1,3]<-round(sum((e13<0)=="TRUE")/length(e11),3)
mat2[2,1]<-round(sum((e21<0)=="TRUE")/length(e11),3)
mat2[2,2]<-round(sum((e22<0)=="TRUE")/length(e11),3)
mat2[2,3]<-round(sum((e23<0)=="TRUE")/length(e11),3)
mat2[3,1]<-round(sum((e31<0)=="TRUE")/length(e11),3)
mat2[3,2]<-round(sum((e32<0)=="TRUE")/length(e11),3)
mat2[3,3]<-round(sum((e33<0)=="TRUE")/length(e11),3)
rownames(mat2)<-c("トマトの冷製","カルボナーラ","ペペロンチーノ")
colnames(mat2)<-c("バジル","トリュフ","なし")	
kable(mat2)

	バジル	トリュフ	なし
トマトの冷製	0.008	0.928	0.887
カルボナーラ	0.899	0.003	0.955
ペペロンチーノ	0.872	0.889	0.006

連言命題が正しい確率

u_e11>0 u_e22>0 u_e33>0 u_e12<0 u_e23<0

1	sum((e11>0)=="TRUE" & (e22>0)=="TRUE"& (e33>0)=="TRUE" & (e12<0)=="TRUE" & (e23<0)=="TRUE" )/length(e11)

0.8814

u_e11>0 u_e22>0 u_e33>0

1	sum((e11>0)=="TRUE" & (e22>0)=="TRUE"& (e33>0)=="TRUE" )/length(e11)

0.98397

放送大学「心理統計法(‘17)」第12章

授業ではrstanパッケージを使ってますが、MCMCpackパッケージを使ってみます。初学者ですので間違いが多々あると思います。

カテゴリ数がkの比率の推測

ペット問題

連言命題が正しい確率

u_p1>p2 u_p2>p3 u_p3>p4

u_p1>p3 u_p2>p3 u_p1>p4 * u_p2>p4

対応ある2×2のクロス表の分析

ブランド認知問題2

独立と連関

対応あるa×bのクロス表の分析

パスタ問題

独立と連関

ピアソン残差が0より大きい確率

ピアソン残差が0より小さい確率

連言命題が正しい確率

u_e11>0 u_e22>0 u_e33>0 u_e12<0 u_e23<0

u_e11>0 u_e22>0 u_e33>0

u_p1>p2 　u_p2>p3 u_p3>p4

u_p1>p3 　u_p2>p3 u_p1>p4 * u_p2>p4