JAGSで対数尤度を得る方法(Stanも)
WAICなどを計算するには，対数尤度のMCMCサンプルが必要。
個々のMCMCサンプルからがんばって計算する方法もあるが，より簡単にはMCMCサンプラーのなかで対数尤度を得ると楽。
Stanではgenerated quantityが容易されているが，JAGSでは少し工夫が必要。
参考
https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/a3bef062/
https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/87902622/
その他参考
ito-hi.blog.so-net.ne.jp
hamada.hatenablog.jp
JAGSのモデルの中でのようなものでも対数尤度が得られる
一変量の分布に関しては，密度関数などがRのように用意されている。
model{ for(i in 1:I) logLik[i] <- log(dnorm(x[i], mu, 1/sigma)) }
そのほか, logdensity.**系の関数が用意されているが，これはマニュアルに書かれていないので必要に応じて調べる。
多分，すべての分布に対してこれらの関数が利用できる？
（未確認）
# # 相関係数については，以下の頁のものを参考にした。 # https://kyougokumakoto.blogspot.jp/2016/09/jags_91.html # # logdensity.mnormで多変量正規分布の対数密度を得られる。 # https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/a3bef062/ # library(MASS) set.seed(1234) mu <- c(0,0) Sig <- matrix(c(1,0.75,0.75,1),2,2) x <- mvrnorm(100,mu,Sig) cor.test(x[,1],x[,2]) data = list(x = x, N = nrow(x)) # JAGS code cat("model{ # likelihood for(i in 1:N){ x[i,1:2] ~ dmnorm(mu[], prec[ , ]) logLik[i] <- logdensity.mnorm(x[i,1:2], mu, prec) } # corresponding the precision matrix prec[1:2, 1:2] <- inverse(cov[,]) # constructing the covariance matrix cov[1,1] <- sigma[1]^2 cov[1,2] <- sigma[1]*sigma[2]*rho cov[2,1] <- sigma[1]*sigma[2]*rho cov[2,2] <- sigma[2]^2 #priors rho ~ dunif(-1,1) mu[1] ~ dnorm(0, 1.0e-4) mu[2] ~ dnorm(0, 1.0e-4) sigma[1] ~ dt(0, 1, 1)T(0,) sigma[2] ~ dt(0, 1, 1)T(0,) }", file="cor.jag") fit_jag <- jags(data=data, parameters.to.save=c("mu", "sigma","rho","logLik"), model.file="cor.jag", n.chains = 3, n.adapt = 100, n.iter = 5000, n.burnin = 2500, n.thin = 2) library("loo") # チェインがたくさんある場合には，長い行列をとしてみる？ log_lik1 <- sapply(1:nrow(x), function(i) do.call("rbind",fit_jag$samples)[,paste0("logLik[",i,"]")]) waic(log_lik1) log_lik2 <- sapply(1:nrow(x), function(i) fit_jag$samples[[1]][,paste0("logLik[",i,"]")]) waic(log_lik2) library(rstan) # # Stanのコード参考：https://qiita.com/hoxo_m/items/0f1b05681f5d6c4b560a # ベクトルの配列について注意。y[i]で得られるのはi番目の長さが2のベクトル。 # y[i][1]でi番目の配列にあるベクトルの1番目の要素を取り出すということになる。 # cat("data{ int<lower=0> N; vector[2] y[N]; } parameters{ vector[2] mu; real<lower=0> var_x; real<lower=0> var_y; real cov; } transformed parameters{ matrix[2,2] sigma; sigma[1,1] = var_x; sigma[2,1] = cov; sigma[1,2] = cov; sigma[2,2] = var_y; } model{ for(i in 1:N){ target += multi_normal_lpdf(y[i]| mu, sigma); //target += normal_lpdf(y[i][2]|0, 1); //target += normal_lpdf(y[i][1]|0, 1); //y[i] ~ multi_normal(mu, sigma); } } generated quantities { real cor; real log_lik[N]; cor = cov/sqrt(var_x * var_y); for(i in 1:N) log_lik[i] = multi_normal_lpdf(y[i] | mu, sigma); } ", file="cor.stan") datastan <- list(N=nrow(x), y = x) fit <- stan(file = "cor.stan", data=datastan, iter = 5000, chains = 2) log_lik_stan <- extract_log_lik(fit) str(log_lik_stan) waic(log_lik_stan) waic(log_lik1) waic(log_lik2) # 値は大体一致している。
MplusのMCMCサンプルをGetする方法(+Mplusの結果をRでプロットする方法)その1
気が向いたらStanやJAGSなどとの対応関係をきちんと取りならが検算を行っていきたい。
MCMCサンプルは「.gh5」ファイルに入っているので，これを読み込む。
以下の資料のp.65あたりの「mplus.get.bayesian.」を利用すると，ギブスサンプラーのMCMCサンプルが得られるはず。
http://www.statmodel.com/mplus-R/Mplus%20R%20tutorial.pdf
以下のinpファイルを作成して，Mplusで推定を行う。
データは
https://www.statmodel.com/usersguide/chap4/ex4.1a.dat
にあるものを，ダウンロードしてくる。
TITLE: this is an example of an exploratory factor analysis with continuous factor indicators DATA: FILE IS ex4.1a.dat; VARIABLE: NAMES ARE y1-y12; ANALYSIS: TYPE = general; ESTIMATOR = BAYES; PROCESSORS = 2; BITERATIONS = (2000); model: y1-y12 with y1-y12; y1-y12; [y1-y12]; OUTPUT: CINTERVAL; plot: type = plot1 plot2 plot3;
library("MplusAutomation")
# MacではRを使ってプロットをする必要がある。
# 以下のsourceを読み込むと必要な関数がいろいろと読み込まれる。
source("https://www.statmodel.com/mplus-R/mplus.R")
# アウトプットの読み込み
hoge1 <- readModels("ex4.1bayes.out")
# 結果の表示
hoge1$parameters
# Rstudioのコンソールでmplusと入力して，タブを押すと色々な関数が候補として出てくる。
# プロットできるもの一覧を表示する。
mplus.view.plots("ex4.1bayes.gh5")
# 変数名を見る
(v_list <- mplus.list.variables("ex4.1bayes.gh5"))
# ヒストグラムの作成
mplus.plot.histogram("ex4.1bayes.gh5","Y1")
# 密度曲線を追加
mplus.plot.densityplot("ex4.1bayes.gh5","Y1")
# qqplotの追加
mplus.plot.qqnorm("ex4.1bayes.gh5","Y1")
# 散布図の作成
mplus.plot.scatterplot("ex4.1bayes.gh5","Y1","Y2")
# データの取得：1変数だけ。
mplus.get.data("ex4.1bayes.gh5","Y1")
ベイズ推定の結果をプロットするためのもの。
まだ利用しきれていない関数もあるので，今後追記する。
#
# ベイズ推定の結果をいろいろと得るためのもの
# mplus.get.bayesian.●●
#
# パラメタ名を取得する
# ついでに，パラメタ名が必要な関数も表示してくれる。
(par_list <- mplus.list.bayesian.parameters("ex4.1bayes.gh5"))
#
# 特定のパラメタのMCMCサンプルを取得する。
# 列がチェインに対応している。
#
par_data <- mplus.get.bayesian.parameter.data("ex4.1bayes.gh5", par_list[1])
# トレースプロット
matplot(par_data, type="l")
mplus.plot.bayesian.traceplot("ex4.1bayes.gh5",paramstr = par_list[1])
# 自己相関
mplus.plot.bayesian.autocorrelation("ex4.1bayes.gh5",paramstr = par_list[1])
# 事後分布
mplus.plot.bayesian.distribution("ex4.1bayes.gh5",paramstr = par_list[1])
Mplusは便利なので，もう少し使いこなしていきたいところ。
よく必要になるけど，よく忘れるR(とExcel)の小技メモ
文字（数字）の桁数を揃えたい。
参考URL：
http://d.hatena.ne.jp/teramonagi/20100902/1283429021
formatC("hoge", widh = 5, flag="0")flagは '0 + - #' だけが指定可能とのこと。
不要な空白を消したい。
# いらない空白を全て削除する
# 全角半角の指定に注意。
gsub(" ", "", x, fixed = TRUE)
# 文字列の前後の空白を削除する
stringr::str_trim(x, side="both")全角半角を変換したい。
Nipponパッケージを利用するといいらしいが，上手く動作しない場合がある。
そんなときは，Excelを使う。
参考URL：全角を半角に、半角を全角に変換する エクセルの関数技
文字を半角にするにはASC関数を使い，半角文字を全角にするにはJIS関数を使う。
不要な全角半角スペースを削除するには，=SUBSTITUTE(SUBSTITUTE(A1," ","")," ","")という関数を利用するとよい。
(trimは文字列前後の空白を削除するための関数。）
Rで図を出力する場合のあれこれ。
参考URL
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/53.html
R のグラフで日本語を使う - 知に至る病
https://oku.edu.mie-u.ac.jp/~okumura/stat/first.html
Rのグラフで軸の目盛りの刻み幅を変更する方法 - Rプログラミングの小ネタ
# 文字のフォントを指定する。 # これを指定しておけば，とりあえず文字化けしない。 pdf(file="hoge.pdf", family="Japan1GothicBBB") # これはMac向き par(family="HiraKakuProN-W3") # ヒラギノ角ゴシックProN W3 # 空のプロットを作成する。 plot.new() # tclで軸を内側に向ける。 # 正の値が内側で，枠の大きさに対する割合で指定で指定する。0.3ぐらいがけっこうちょうどいい。 plot(0,1,tcl = 0.3) # メモリを独自に振り直したいときに使う。 #メモリを消す curve(dnorm, xlim=c(-3, 3), xaxt="n") # atに書き込みたい目盛りのベクトルを指定する axis(side=1, at=-3:3) # sideは下から時計回りに1,2,3,4になっている。 # ラベルを各軸に並行して描く las = 0 # ラベルをすべて水平に描く las = 1 # ラベルを軸に対して垂直に描く las = 2 # ラベルをすべて垂直に描く las = 3
Mplusでターゲット行列を指定した回転を行う方法
Analysisコマンドで，rotation = targetを指定し，
modelコマンド内で目標とする値を~で指定する。
また，各因子の最後の(*t)などとつけて，各因子負荷が1つの因子負荷行列で表すことができることを示さないといけない。
Mplusのマニュアルは，versin 8のp681などに記載されている。正直，いろいろと回転がありすぎて何が何だかわからない。
本文をtargetで検索するのがはやい。
これを使えば，任意の値を指定したターゲット行列に向けた回転を行ってくれる。
（つまり，プロクラステス(procrustes)回転ができる。）
確認的因子分析よりもゆるく，探索的因子分析よりも仮説を反映した分析ができるので，もう少し利用されてもいい方法のように感じる。
使用したデータ
https://www.statmodel.com/usersguide/chap4/ex4.1a.dat
TITLE: DATA: FILE IS ex4.1a.dat; VARIABLE: NAMES ARE y1-y12; ANALYSIS: ESTIMATOR IS ML; ROTATION=TARGET; MODEL: f1 BY y1-y12 y1~0 y2~0 y3~0 y4~0 y5~1 y6~1(*t); f2 BY y1-y12 y7~0 y8~0 y9~0 y10~0 y11~1 y12~1(*t); OUTPUT: MODINDICES CINTERVAL;
Wordの小技（大文字小文字の変換，テキストボックス内の文字列の検索）
Wordで便利な小技を見つけた。
どちらも，「高度な検索と置換」を使う。
特に，ワイルドカードと，検索対象を指定するのがポイント。
「[０-９]」で検索を行うと，大文字の数字が選択される。その上で，「A->a」で「半角」を選択すればOK。
テキストボックス内の文字列の選択では，ワイルドカードの「*」を使って，検索対象を「メイン文書」ではなく「メイン文書のテキストボックス」にして検索をする。その後，検索ウィンドを閉じて「F9」をつかって，フィールドの更新をすればよい。
Wordで数字を大文字から小文字に変換する。
office-qa.com
Wordで複数のテキストボックスの文字列を選択する。
detail.chiebukuro.yahoo.co.jp
そして，フィールドについて，何処かでまとめて勉強したほうがいい気がする。
あとで，まとめて読んでみること。
平均値の信頼区間のシミュレーションプロット
Rで，信頼区間を書いていいプロットするコードのメモ。
layout関数で，上手くプット画面を分割して，さも1枚のグラフかのように見せかけるのがポイント。
けっこう可愛い感じのグラフを描画できる。
# # 縦向きの，信頼区間の図。 # ポイントは，レイアウトとparの使い方。 # layout(matrix(c(1,2),ncol=1), heights=c(1/5,4/5)) par(mar=c(0,3,1,1)) mu <- 0 sigma <- 1 n <- 5 size <- 100 curve(dnorm(x,mu,sigma/sqrt(n)),from=-3, to = 3,axes=F, xlab = "",ylab="",lty =1, col = rgb(1,0,1,alpha = 0.7),lwd = 2) arrows(mu, 0, mu, dnorm(mu,mu,sigma/sqrt(n)), lwd=2, angle = 0,lty =1, col = rgb(1,0,1,alpha = 0.7)) x <- matrix(rnorm(size*n,mean = mu,sd = sigma), ncol=n) x_mean <- rowMeans(x) upp <- x_mean + 1.96*sigma/sqrt(n) low <- x_mean - 1.96*sigma/sqrt(n) is_sig <- ((upp < 0) + (low < 0)) !=1 par(mar=c(3,3,0,1)) plot(x_mean,1:size, xlim = c(-3,3),col =ifelse(is_sig, rgb(1,0,0,alpha=0.7), rgb(0,0,1,alpha=0.7)), pch = 1- is_sig ,xlab="x平均") arrows(x_mean, 1:size, upp, 1:size, lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7), rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig) arrows(x_mean, 1:size, low, 1:size, lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7), rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig) abline(v = mu , lwd = 2,lty=2,col = rgb(1,0,1,alpha = 0.7))
出来上がる図は次のようなもの。
# # 横向きの，信頼区間の図。 # layout(matrix(c(2,1),nrow=1), widths =c(4/5,1/5)) par(mar=c(3,0,1,1)) # 下，左，上，右の順で，余白を設定する。上が0なので，いい感じでくっつけられる。 mu <- 0 sigma <- 1 n <- 10 size <- 100 # 右に90度，図を転置させる。 t_dnorm <- function(x, mu, sigma){ cbind(x, dnorm(x,mu,sigma)) %*% matrix(c(0,-1,1,0),ncol=2,byrow=T) } plot(t_dnorm(seq(from=-3, to =3, by=0.1), mu, sigma/sqrt(n)),type="l",axes=F, xlab = "",ylab="",lty =1, col = rgb(1,0,1,alpha = 0.7),lwd = 2) arrows(0, mu, dnorm(mu,mu,sigma/sqrt(n)), mu, lwd=2, angle = 0,lty =3, col = rgb(1,0,1,alpha = 0.7)) x <- matrix(rnorm(size*n,mean = mu,sd = sigma), ncol=n) x_mean <- rowMeans(x) upp <- x_mean + 1.96*sigma/sqrt(n) low <- x_mean - 1.96*sigma/sqrt(n) is_sig <- ((upp < 0) + (low < 0)) !=1 par(mar=c(3,3,1,0)) # この余白の設定で，絶妙に図をくっつけるのがポイント。 plot(x_mean, ylim = c(-3,3),col =ifelse(is_sig, rgb(1,0,0,alpha=0.7), rgb(0,0,1,alpha=0.7)), pch = 1- is_sig ,xlab="x平均") arrows(1:size,x_mean, 1:size, upp, lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7), rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig) arrows(1:size,x_mean, 1:size, low,lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7), rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig) abline(h = mu , lwd = 2,lty=3,col = rgb(1,0,1,alpha = 0.7))
出来上がる図は次のようなもの。
