×
[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。
広告対策
PR
# この記事と同じくoptim関数で
# 井関先生のサイトを参考にした。感謝。
# Bill Prinzmetalのサイトにある論文とExcelファイルを参考
# # Dodson, C. S., Prinzmetal, W., & Shimamura, A. (1998). Using Excel estimate parameters from observed data: An exa mple from source memory data. Behavioral Research Methods, Instruments & Computers, 30, 517-526.
# 一番単純な例でやってみる
# 初期値
values <- rep(0.5, 7)
values
# 観測値
obs <- c(612,151,77,123,643,74,19,18,383)
# 観測回答率
act.prb <- c(
# Sorce=Eric
obs[1]/sum(obs[1:3]), # response=eric
obs[2]/sum(obs[1:3]), # response=julie
obs[3]/sum(obs[1:3]), # response=new
# Sorce=Julie
obs[4]/sum(obs[4:6]), # response=eric
obs[5]/sum(obs[4:6]), # response=julie
obs[6]/sum(obs[4:6]), # response=new
# Sorce=New
obs[7]/sum(obs[7:9]), # response=eric
obs[8]/sum(obs[7:9]), # response=julie
obs[9]/sum(obs[7:9]) # response=new
)
# データの確認
data.frame(obs, act.prb)
# 関数定義
prd.prb <- vector("numeric", 9)
fnc <- function(x, ifm=F){
d1 <- x[1]; d2 <- x[1] # 等値制約のため同じ値を入れる
ld1 <- x[3]; ld2 <- x[3] # 等値制約のため同じ値を入れる
a <- x[5]; g <- x[5] # 等値制約のため同じ値を入れる
b <- x[7]
prd.prb[1] <- (d1*ld1) + (d1*(1-ld1)*a) + ((1-d1)*b*g)
prd.prb[2] <- (d1*(1-ld1)*(1-a)) + ((1-d1)*b*(1-g))
prd.prb[3] <- (1-d1) * (1-b)
prd.prb[4] <- (d2*(1-ld2)*a) + ((1-d2)*b*g)
prd.prb[5] <- (d2*ld2) +(d2*(1-ld2)*(1-a)) + ((1-d2)*b*(1-g))
prd.prb[6] <- (1-d2) * (1-b)
prd.prb[7] <- b * g
prd.prb[8] <- b * (1-g)
prd.prb[9] <- 1 - b
x1 <- 2*obs
x2 <- log(act.prb/prd.prb)
g2v <- x1*x2
estv <- c(d1, d2, ld1, ld2, a, g, b)
reslist <- list(prd.prb, estv, g2v, sum(g2v))
names(reslist) <- c("prd.prb", "estimates", "G2", "sumG2")
if(ifm==T) return(reslist) else return(sum(g2v))
}
(opres <- optim(par=values, fn=fnc, method="BFGS"))
est.v <- opres$par
# 改めて予測値を出す
x <- est.v
d1 <- x[1]; d2 <- x[1] # 等値制約のため同じ値を入れる
ld1 <- x[3]; ld2 <- x[3] # 等値制約のため同じ値を入れる
a <- x[5]; g <- x[5] # 等値制約のため同じ値を入れる
b <- x[7]
est.v <- c(d1,d2, ld1, ld2, a, g, b)
fnc(est.v, ifm=T)
# 適合度の評価
dfv <- 2 # 自由度は2
g2 <- fnc(est.v, ifm=T)$sumG2
pchisq(g2, dfv)
# 前の記事から進歩しておらず、相変わらずconstrOptim関数の使い方がわからない。たぶんこれで制約を色々入れないとダメだろうな…
# Multinomial Processing Treeモデル専用のパッケージmptというのがあるようだ。よくこんなマニアックなものを…。これも使い方がよくわからんのう
http://homepages.uni-tuebingen.de/florian.wickelmaier/wickelmaier_software.html#mpt
# 参考文献
Batchelder, W. H., & Riefer, D. M. (1990). Multinomial processing models of source monitoring. Psychological Review, 97, 548-564.
Batchelder, W. H., & Riefer, D. M. (1999). Theoretical and empirical review of multinomial process tree modeling. Psychonomic Bulletin & Review, 6, 57-86.
Dodson, C. S., Prinzmetal, W., & Shimamura, A. P. (1998). Using Excel to estimate parameters from observed data: An example from source memory data. Behavior Research Methods, Instruments, & Computers, 30, 517-526.
Jacoby, L., Bishara, A., Hessels, S., & Toth, J. (2005). Aging, subjective experience, and cognitive control: Dramatic false remembering by older Adults. Journal of Experimental Psychology: General, 134, 131-148.
Wickelmaier, F. (2010). mpt: Multinomial Processing Tree (MPT) Models. R package version 0.2-0. http://CRAN.R-project.org/package=mpt ワークショップのpdf
Psychoco 2011: International Workshop on Psychometric Computing
# 井関先生のサイトを参考にした。感謝。
# Bill Prinzmetalのサイトにある論文とExcelファイルを参考
# # Dodson, C. S., Prinzmetal, W., & Shimamura, A. (1998). Using Excel estimate parameters from observed data: An exa mple from source memory data. Behavioral Research Methods, Instruments & Computers, 30, 517-526.
# 一番単純な例でやってみる
# 初期値
values <- rep(0.5, 7)
values
# 観測値
obs <- c(612,151,77,123,643,74,19,18,383)
# 観測回答率
act.prb <- c(
# Sorce=Eric
obs[1]/sum(obs[1:3]), # response=eric
obs[2]/sum(obs[1:3]), # response=julie
obs[3]/sum(obs[1:3]), # response=new
# Sorce=Julie
obs[4]/sum(obs[4:6]), # response=eric
obs[5]/sum(obs[4:6]), # response=julie
obs[6]/sum(obs[4:6]), # response=new
# Sorce=New
obs[7]/sum(obs[7:9]), # response=eric
obs[8]/sum(obs[7:9]), # response=julie
obs[9]/sum(obs[7:9]) # response=new
)
# データの確認
data.frame(obs, act.prb)
# 関数定義
prd.prb <- vector("numeric", 9)
fnc <- function(x, ifm=F){
d1 <- x[1]; d2 <- x[1] # 等値制約のため同じ値を入れる
ld1 <- x[3]; ld2 <- x[3] # 等値制約のため同じ値を入れる
a <- x[5]; g <- x[5] # 等値制約のため同じ値を入れる
b <- x[7]
prd.prb[1] <- (d1*ld1) + (d1*(1-ld1)*a) + ((1-d1)*b*g)
prd.prb[2] <- (d1*(1-ld1)*(1-a)) + ((1-d1)*b*(1-g))
prd.prb[3] <- (1-d1) * (1-b)
prd.prb[4] <- (d2*(1-ld2)*a) + ((1-d2)*b*g)
prd.prb[5] <- (d2*ld2) +(d2*(1-ld2)*(1-a)) + ((1-d2)*b*(1-g))
prd.prb[6] <- (1-d2) * (1-b)
prd.prb[7] <- b * g
prd.prb[8] <- b * (1-g)
prd.prb[9] <- 1 - b
x1 <- 2*obs
x2 <- log(act.prb/prd.prb)
g2v <- x1*x2
estv <- c(d1, d2, ld1, ld2, a, g, b)
reslist <- list(prd.prb, estv, g2v, sum(g2v))
names(reslist) <- c("prd.prb", "estimates", "G2", "sumG2")
if(ifm==T) return(reslist) else return(sum(g2v))
}
(opres <- optim(par=values, fn=fnc, method="BFGS"))
est.v <- opres$par
# 改めて予測値を出す
x <- est.v
d1 <- x[1]; d2 <- x[1] # 等値制約のため同じ値を入れる
ld1 <- x[3]; ld2 <- x[3] # 等値制約のため同じ値を入れる
a <- x[5]; g <- x[5] # 等値制約のため同じ値を入れる
b <- x[7]
est.v <- c(d1,d2, ld1, ld2, a, g, b)
fnc(est.v, ifm=T)
# 適合度の評価
dfv <- 2 # 自由度は2
g2 <- fnc(est.v, ifm=T)$sumG2
pchisq(g2, dfv)
# 前の記事から進歩しておらず、相変わらずconstrOptim関数の使い方がわからない。たぶんこれで制約を色々入れないとダメだろうな…
# Multinomial Processing Treeモデル専用のパッケージmptというのがあるようだ。よくこんなマニアックなものを…。これも使い方がよくわからんのう
http://homepages.uni-tuebingen.de/florian.wickelmaier/wickelmaier_software.html#mpt
# 参考文献
Batchelder, W. H., & Riefer, D. M. (1990). Multinomial processing models of source monitoring. Psychological Review, 97, 548-564.
Batchelder, W. H., & Riefer, D. M. (1999). Theoretical and empirical review of multinomial process tree modeling. Psychonomic Bulletin & Review, 6, 57-86.
Dodson, C. S., Prinzmetal, W., & Shimamura, A. P. (1998). Using Excel to estimate parameters from observed data: An example from source memory data. Behavior Research Methods, Instruments, & Computers, 30, 517-526.
Jacoby, L., Bishara, A., Hessels, S., & Toth, J. (2005). Aging, subjective experience, and cognitive control: Dramatic false remembering by older Adults. Journal of Experimental Psychology: General, 134, 131-148.
Wickelmaier, F. (2010). mpt: Multinomial Processing Tree (MPT) Models. R package version 0.2-0. http://CRAN.R-project.org/package=mpt ワークショップのpdf
Psychoco 2011: International Workshop on Psychometric Computing
# 参考は金先生のページ。感謝http://mjin.doshisha.ac.jp/R/
# 足立浩平 (2006). 多変量データ解析法――心理・教育・社会系のための入門―― ナカニシヤ出版の第2章から。感謝
mat <- matrix(c(3.2,3.4,3,3.2,4.2,4,4,3.7,3.6,3.6,3.7,3.4,3.2,3.2,2.7,3.5,3.2,3.2,4.6,4,4.8,4.6,4.3,3.6,3.2,3.7,3.2,3.8,3.7,3.4,3.5,3.5,4.5,3.8,3.9,4.1,3.7,3.5,3.7,3.5,3.6,3.8,2.8,2.5,2.2,2.6,3.1,3.4,3.5,3.4,2.9,3.5,3.9,3.1,2.9,2.8,2.6,2.2,2.1,2.5,3,3.2,3.8,4,3.5,4.2,4.7,2.7,2.2,2.3,2.6,2.6,2.2,2.6,3.2,3.1,3.7,4.1,3.6,4.2,4.7,2.4,2.5,2.4,2.2,3.2,3.3,3.6,2.8,2.4,3.2,4.3,4.7,3.5,4.9,2.3,3.3,3.9,1.4,2.1,3.4,2.9,3.3,1.5,2.1,3.4,4.2,3.5,3.5,1.8,3.3,2.5,1.7,3.6,4.1,4.2,4.1,1.6,2.6,3.5,4.1,3.7,4.2,2.3,3.4,4.7,3.3,4.1,3.4,3.2,4.5,3.7,3.7,4.2,3.9,3.5,3.7,3.3,2.8,3.9,3.8,4.7,1.3,1.5,2.3,3.9,3.6,4.4,3.7,2.5,2.8,2.9,1.8,2.3,1.8,4.2,4.3,4,4.9,3,4.5,4,5,3.5,4.1,3.3,4.3,4.3),
nr=14,
dimnames=list(c("僧侶", "銀行員", "漫画家", "デザイナー", "保母", "大学教授", "医師",
"警察官", "新聞記者", "船乗り", "プロスポーツ選手", "作家", "俳優", "スチュワーデス"),
c("立派な", "役立つ", "よい", "大きい", "力がある", "強い", "速い", "騒がしい", "若い", "誠実な",
"かたい", "忙しい"))
)
distdat <- dist(mat) # ユークリッド距離の算出
res.w <- hclust(distdat, method="ward") # ウォード法による分析
res.a <- hclust(distdat, method="average") # ウォード法による分析
# method引数のメモ
# signle: 最近隣法
# complete: 最遠隣法
# average: 群平均法
# centroid: 重心法
# median: メディアン法
# ward: ウォード法
# mcquitty: McQuitty法
summary(res.w) # リストオブジェクトの要素を確認するのみ
plot(res.w) # ウォード法樹形図の描画 図2.6A
plot(res.a) # 群平均法樹形図の描画 図2.6B
# kmeans法による非階層的クラスター分析はkmeans関数を用いる。そのうちやろう
# 足立浩平 (2006). 多変量データ解析法――心理・教育・社会系のための入門―― ナカニシヤ出版の第2章から。感謝
mat <- matrix(c(3.2,3.4,3,3.2,4.2,4,4,3.7,3.6,3.6,3.7,3.4,3.2,3.2,2.7,3.5,3.2,3.2,4.6,4,4.8,4.6,4.3,3.6,3.2,3.7,3.2,3.8,3.7,3.4,3.5,3.5,4.5,3.8,3.9,4.1,3.7,3.5,3.7,3.5,3.6,3.8,2.8,2.5,2.2,2.6,3.1,3.4,3.5,3.4,2.9,3.5,3.9,3.1,2.9,2.8,2.6,2.2,2.1,2.5,3,3.2,3.8,4,3.5,4.2,4.7,2.7,2.2,2.3,2.6,2.6,2.2,2.6,3.2,3.1,3.7,4.1,3.6,4.2,4.7,2.4,2.5,2.4,2.2,3.2,3.3,3.6,2.8,2.4,3.2,4.3,4.7,3.5,4.9,2.3,3.3,3.9,1.4,2.1,3.4,2.9,3.3,1.5,2.1,3.4,4.2,3.5,3.5,1.8,3.3,2.5,1.7,3.6,4.1,4.2,4.1,1.6,2.6,3.5,4.1,3.7,4.2,2.3,3.4,4.7,3.3,4.1,3.4,3.2,4.5,3.7,3.7,4.2,3.9,3.5,3.7,3.3,2.8,3.9,3.8,4.7,1.3,1.5,2.3,3.9,3.6,4.4,3.7,2.5,2.8,2.9,1.8,2.3,1.8,4.2,4.3,4,4.9,3,4.5,4,5,3.5,4.1,3.3,4.3,4.3),
nr=14,
dimnames=list(c("僧侶", "銀行員", "漫画家", "デザイナー", "保母", "大学教授", "医師",
"警察官", "新聞記者", "船乗り", "プロスポーツ選手", "作家", "俳優", "スチュワーデス"),
c("立派な", "役立つ", "よい", "大きい", "力がある", "強い", "速い", "騒がしい", "若い", "誠実な",
"かたい", "忙しい"))
)
distdat <- dist(mat) # ユークリッド距離の算出
res.w <- hclust(distdat, method="ward") # ウォード法による分析
res.a <- hclust(distdat, method="average") # ウォード法による分析
# method引数のメモ
# signle: 最近隣法
# complete: 最遠隣法
# average: 群平均法
# centroid: 重心法
# median: メディアン法
# ward: ウォード法
# mcquitty: McQuitty法
summary(res.w) # リストオブジェクトの要素を確認するのみ
plot(res.w) # ウォード法樹形図の描画 図2.6A
plot(res.a) # 群平均法樹形図の描画 図2.6B
# kmeans法による非階層的クラスター分析はkmeans関数を用いる。そのうちやろう
