PsyStat.

Just for fun.

Stan Advent Boot Camp 2020 第10日目; おすしでインドスカルしてたら、推しメンが結婚した話。

 

ういっす。この記事は、Stan Advent Calendar 2020 10日目の記事です。

多次元尺度構成法(Multidimensional scaling; MDS)についてStanで実行してみます。

初めに書いておくと、良い/面白い結果が出ていないので、やってみたぐらいのノリです。

 

qiita.com

 

 

~余談~

1ヶ月前のある日、突如某slackでスコココと通知音が鳴った。

 

なにかと思えば、来月に控えたStan Advent Calendarの埋まり具合が悪いらしい。それを危惧した某メンバーが最初の10日間で初心者のためのStan Boot Campをやることを提案した。僕はその様子を静観していた。しかし、それは許されなかった。そして、ここから悲劇は始まる。気付いたら、Stan Boot Campの担当者になっていたのだ。これは道を歩いていたら虎に襲われるぐらいの衝撃である。しかし、驚くべきはこれだけではなかった。最初は自身の研究領域であるIRTについて執筆する予定だったのだが、「専門の人が専門領域の記事を書いても面白くない」などと言った謎の意見により、気付いたらMDSの記事を書くことになっていた。これはまさに道を歩いていたら気が狂って目の前にいる虎を襲ってしまうぐらいの衝撃である。そして、今日に至る。。。

~余談おわり~

 

MDSとは、「距離から地図を作る方法である(小杉, 2018)」と説明がされています。

www.amazon.co.jp

 

距離と言うと、実際に東京から大阪といった物理的な距離をイメージされることが多いと思いますが、ここでは必ずしも物理的な距離のみを指すわけではありません。

 

  • たとえば、心理学にはパーソナルスペースという言葉があります。これは、他人に入ってきて欲しくない(入られると不快になる)距離を指します。電車に乗っていて、そんなに混んでいないのにもかかわらずこの人めっちゃ近いなみたいな場面を想像して貰えば良いと思います。このような心理的距離(測定すれば物理距離として定義もできますが)も、距離として扱うことができます。
  • また、類似度は距離として解釈されることがあります。松屋すき家・𠮷野家の値段の類似度を出したときに、松屋すき家間は同価格帯で類似しており距離が近いとすると、松屋と𠮷野家は前者間よりも価格に開きがあり(類似していなく)距離が遠いと定義することができます。これも物理的な距離ではありませんが距離として扱うことが可能です。
  • そして、対象への選好(好み)も距離として捉えることができます。りんごとれもんのどちらがどれだけ好きかを考え、れもんが好きだとしましょう。このとき、好きな者ほど、距離が近いと定義すれば、ぼくとれもんは、ぼくとりんごよりも近いと解釈することができ、これは様々なフルーツで行うことで、僕の好きなフルーツを地図で表すことができます。このように選好も距離と考えることができます。
  • また、これらの距離は必ずしも対象とは限りません。例えば、仮に僕がある女優Tさんが好きだったととしましょう。あくまでも仮の話です。このとき僕の思うこの気持ち a.k.a 想いは、 距離です。しかし、残念なことに女優Tさんは僕の存在すら認識していませんので、女優Tさんから私に対する気持ち a.k.a 想いは距離ですが、未観測なのでNULLです。つまり、この場合、この私の気持ち a.k.a 恋は、片思いであり、すなわち非対称距離になるのですね。MDSでは、このような非対称距離を扱うこともできるそうです(この例の場合女優から私への距離はNULLなので厳密に非対称距離とはなっていませんかもしれませんが)。

 

以上のように、単純に距離と言っても多種多様です。

では、今回は、対象への選好を距離と捉えて、その関係性を表す地図を作って行きたいと思います。

 

今回用いるデータは、寿司の好みについて評定したデータです。

こちらをお借りしました。

www.kamishima.net

 

これは全国の日本人参加者(N=5000)が各寿司ネタ(全100ネタ)の好みを評定した大変貴重なデータです。

 

今回はここから、「えび」「まぐろ」「とろ」「かっぱ巻き」「はまち」の5ネタを評定した回答者(N=5)に絞り、分析を実施しました。

※本当は、10ネタぐらいでやりたかったのですが、10ネタすべてを評定した回答者がいなかったので(いずれかのネタが欠損になってしまう)、その兼ね合いで5ネタとしました。正直この少なさでMDSを使うのは向いていない気もします。また欠損値を含むMDSができれば良かったのですが、MDSの初心者なのでよくわからなかったです。検索したら僕のボスの博論が出てきました。

 

モデルや分析コードは以下の小杉先生の章を参考にしました。

www.kitaohji.com

 

詳細なモデルの数理等については関連書籍を確認してください(まとめる余裕なかった)。そんなに難しくないです。真のネタ間の距離が、個人ごとの選好の違いと掛け合わされて、誤差とともに出てきたものが今回の評定データであるというモデルになっています。

このようなモデルは、INdivisual Differences multidimensional SCALing(INDSCAL)と呼ぶそうです。なお、INDSCALはインドスカルではなく、インスカルと読むらしい。

 

また、モデルは5ネタでN=5という通常あり得ないサンプルサイズでの解析なので本当にやってみた程度の結果です。

 

これらをもとに作成した地図がこれです(点はEAP,エラーバーは50%信用区間)。

f:id:dastatis:20201211005054p:plain

 

うーん、やはりネタ数が少なすぎて解釈できないですね。

いろいろ寿司ネタを変えたり、試したんですが、無理でした。このデータだと厳しかったかもしれませんね。

別のデータでリベンジしたいですね。

締まらない感じ&オーバータイムになってしまいましたが、とりあえず。

 

 

コード等は以下

https://github.com/dastatis/dastatis.github.io/tree/master/pdf/StanAdvent2020

 

日本視覚学会@2021年冬季大会にてOpen Science Frameworkの解説を行います。

こんにちは。

 

この記事は宣伝です。

以下のAdvent カレンダーの2日目の記事です。

adventar.org

 

オープンサイエンスやっていますか?

やらねばいけないとわかっていながら、結局手を出していないひと。

やりたいけれども、何から始めれば良いかわからないひと。

 

1/20にOSFの使い方について、実際のハンズオンも含め、日本視覚学会冬季大会2021にてトークします。

visionsociety.wixsite.com

 

オンラインですので、ぜひとも聞きに来てくださいませ。

僕の講演の要旨は以下です(後ほど、HP上で公開されると思います)。

近年、心理学を初めとする社会科学の諸領域において、研究の再現性が問題になっている。研究の再現性を高めること、これを示すことは決して受動的にできることではない。これからの時代は、自らの研究が信頼されうる研究であることを積極的に示すことが必要になる。本セッションでは、再現性の問題について簡単に振り返りつつ、自身の研究がReproducible(再現可能)な研究であること、すなわち信頼されうる研究であることを示す方法について解説し、実践する。具体的には、Open Science FrameworkというWebサイトを使用し、研究の目的・方法の再現性と結果の再現性を高める方法について解説し、これらを実際に動かし、その簡便性を体感してもらう。最終的に本セッションの参加者が、この日から再現性を高める研究を実践できることを目指す。

 

また、企画セッションの他の講演では、様々なコンテンツを学ぶことができそうです。

「オンライン実験・調査をする際の諸々」や「計算論的アプローチの最前線研究」、どちらもとても面白そうですね。

 

早期申込は12/7までのようですよ。

 

ということで、 簡単な宣伝で申し訳ないですが、Enjoy!

 

 

余力があれば、データ匿名化の話をこのアドベントカレンダーの何日目かに書きます。

お楽しみに。

 

【Tips】cmdstanr使うときの諸々

こんにちは。

自分用にcmdstanrを使うときの諸々のTipsを随時更新していきます。

 

環境

WindowsからWSLでRstudio Serverからcmdstanrを動かすことを想定。

普段はrstanを使用。

 

cmdstanrのiterationの指定がrstanと異なる

rstanは、warmupを含む反復回数をiterとして指定

cmdstanrは、warmupを含まない(after warmup)をnum_sumplesとして指定(JAGS等と同一)

 

cmdstanrの推定結果(stan_csv)の保存先を変えたい

cmdstanrのsample()のオプションでoutput_dirを指定する

 

WindowsとWSL間でファイルのやりとりがしたい

cmdstanrそのものと関係ありませんが、次の記事を参考にしました。

 

https://qiita.com/quzq/items/1

qiita.com

096c638c0d86795be13

 

(Windows → WSL) のやりとりをしたい

想定例:WSLに保存した分析結果にアクセスしたい

以下のパスからアクセス可能。

//wsl$

なので、WSLのRstudio Server内のファイルへのアクセスは例えば以下から。

//wsl$/Ubuntu-18.04/home/rstudio 

 

(WSL → Windows) のやりとりをしたい

想定例:WSLのRstudio Serverでの分析結果を直接Windowsに保存したい

以下のパス指定で可能(dドライブへのパスの場合)

/mnt/d/

例えば、デスクトップにggplotの図を保存したい。

ggsave(filename = "/mnt/c/Users/dhojo/Desktop/demo.png")

※指定の際のポイントは、ドライブの文字は小文字じゃないとうまくいかない

※先述のoutput_dirと組み合わせれば、推定結果をWindowsに保存可能だと。

 

デフォルトでwarmup情報は保持されない。

save_warmup = FALSEがデフォルトなので、必要に応じてこのオプションをTRUEとする。

 

vbの動作方法

コンパイルしたモデルのオブジェクトにvariationalをつける

model <- cmdstanr::cmdstan_model("demo.stan")

fit_vb <- model$variational

 

== 2020/04/29 追加 ==

 

MCMCサンプルを取り出したい(rstan::extractをしたい)

$draws()で取り出せる。

fit$draws

> fit$draws()
# A draws_array: 250 iterations, 4 chains, and 641 variables
, , variable = lp__

         chain
iteration     1     2     3     4
        1 -3790 -3797 -3837 -3816
        2 -3786 -3786 -3813 -3784
        3 -3794 -3788 -3809 -3781
        4 -3822 -3789 -3803 -3781
        5 -3828 -3811 -3798 -3805

, , variable = theta[1,1]

         chain
iteration      1      2     3    4
        1  0.725  0.448 0.849 1.06
        2  0.261 -0.176 0.625 0.12
...



fit.vb$draws
> fit_vb$draws()
# A draws_matrix: 1000 draws, and 642 variables
    variable
draw  lp__ lp_approx__ theta[1,1] theta[2,1] theta[3,1] theta[4,1]
  1  -3936        -306     -0.101      0.657     -0.095      -1.08      
  2  -3861        -240      0.218     -0.866      0.602       0.55      
  3  -3890        -260      0.390      0.189      0.865      -0.46      
  4  -3896        -268     -0.462     -1.276     -0.112      -0.23      
  5  -3910        -274     -0.179     -1.564      0.194       0.72      
  6  -3906        -281      0.167     -0.087     -0.710      -0.85      
  7  -3909        -271      0.453     -0.735     -0.389      -1.21      

...

 

 

==2020/05/01追記== 

bridgesamplingパッケージを使いたい

現状、cmdstanrのオブジェクトそのままでは使えない?と思うので、

以下の要領で使えるようにする。

1.warmup情報を保存しておく

上記に書いたとおり、warmup保存オプションを書いておく。

> save_warmup = TRUE

2.rstan::stan()でモデルを空回しする

1.の準備を行っても、「モデルが妥当ではない」という旨のエラーをはかれるので、

## 事前に

## fit_cmdstan <- rstan::read_stan_csv(fit$output_files())

## をしておく。

fit_rstan <- rstan::stan("hoge.stan",data,iter=0)

bridgesampling::bridge_sampler(fit_cmdstan,fit_rstan)

これで一応回る。

ただ、これだと結局rstanでモデルのコンパイルをすることになります笑

まぁ、モデル評価の時ぐらいはコンパイル時間を待ちましょうかね!

 

 

 

また、随時追加します。

 

参考:

https://qiita.com/quzq/items/1096c638c0d86795be13

https://mc-stan.org/cmdstanr/reference/cmdstanr-package.html

 

 

マニアックなIRTモデルのパラメータリカバリ性能をひたすらチェックするよ(泣)

こんにちは。

 

この記事は、Stan Advent Calendarの12日目の記事になります。

本当は、今年書く予定はなかったのですが、某先生に名指しで「12日、空いてるよ」と強い圧をかけられたので書くことになりました笑

最初は何言ってんのかわからなくて、デートにでも誘われているのかと思いました。

そんな時期です。

 

今回の記事は書くことは一昨日決まったので、割とやっつけですがすみません。

 

最近、心理学等の社会科学の分野において

(Stanで)ベイズ統計モデリングを行うことが注目されつつあるかと思います。

 

データをとれば自由な統計モデルで推定することができるようになりました。

適当に複数モデル作って、こっちの方がよさそうみたいなことも言えると思います(例えば、情報量規準やBFをもちいて)。

 

ただ、これらの統計モデリングをやる上で

まず、大事なことはそのモデルの挙動を知ることかなと思います。

 

挙動を知るとは、

モデルのあるパラメータの値を動かすと(変化させると)、

どのようなデータが生成されたり、他のパラメータにどのような影響を与えるかを検討することをここでは指しています(他にもあるかもしれませんが)。

 

そして、モデルの挙動を知るために、

統計学の多くのモデル研究でまず検討されることは、

パラメータリカバリ性能の確認かなと思います。

 

心理学で(ベイズ)統計モデリングが今後さらに流行っていく上で、

絶対に必要な研究手続きになると思います。

なので、今回は、IRTモデル(説明変数が潜在変数である特殊なロジスティック回帰モデル)を用いて考えていきます。

 

改めて、

パラメータリカバリの性能の確認とは何でしょうか?

我々が一般的な(応用的な)心理学研究を行う際は、

データを測定し、統計モデルにおけるパラメータ推定を行います。

なので、端的に表すと

実データ → 統計モデル → パラメータ
という流れになると思います。

 

ただ、ここで統計モデルがあまりよくない場合に、

統計モデル → パラメータ

がうまくいかないこともあります。

パラメータリカバリの性能の確認とは、これを確認するシミュレーションになります。

基本的な流れとしては、

パラメータの真値設定 → 統計モデル → 疑似データ → 統計モデル → 推定されたパラメータ

を行い、最初に設定したパラメータ(真値)と推定されたパラメータの差異等を確認することでそのモデル(や推定法)の推定精度をみることです。

 

認知モデルなんかでは、これを行った後に、その真値パラメータを動かしたときのデータの生成挙動や他のパラメータの挙動を確認すべきと思います。

こういうのも所謂、感度分析に近いものだと思います。

 

では、本編に入ります。
IRTの2値のロジスティックモデルにおいて、
最もよく使われるのが1パラメータもしくは2パラメータかなと思います。
当て推量を考慮した3パラメータロジスティックモデルも知っている人はいるかと思いますが、
実はもう少しマニアック?なロジスティックモデルとして、4パラメータ、5パラメータモデルも存在します。
今回は、この二つのモデルのパラメータリカバリを行ってみたいと思います。

ここらへんのモデルは荘島先生のサイトが端的にまとまっています。
http://www.rd.dnc.ac.jp/~shojima/exmk/index.htm

一応、整理すると以下のような感じです。

Parameter 解釈 1PL 2PL 3PL 4PL 5PL
$$\theta_i$$ 回答者iの特性 / 能力 ( trait / ability )
$$\beta_j$$ 項目jにおける困難度 ( difficulty )
$$\alpha_j$$ 項目jにおける識別力 ( discriminaion )
$$\gamma_j$$ 項目jにおける当て推量 / 下限
$$\upsilon_j$$ 項目jにおける当て推量 / 上限
$$\zeta_j$$ 項目jにおける非対称性

モデル

1 Parameter Logistic Model ( Rasch Model )


\begin{align}
P(Y_{ij} | \theta_i, \beta_j)
&= logit^{-1}( \theta_i - \beta_j) \\
&= \frac{1}{1 + \exp( \theta_i - \beta_j)}
\end{align}

項目特性曲線(1PL)

f:id:dastatis:20191212015749p:plain
1PL

2 Parameter Logistic Model

$$
\begin{aligned}
P(Y_{ij} | \theta_i,\alpha_j,\beta_j)
&= logit^{-1}(a_j \theta_i - \beta_j) \\
&= \frac{1}{1 + \exp(a_j \theta_i - \beta_j)}
\end{aligned}
$$

項目特性曲線(2PL)

f:id:dastatis:20191212020032p:plain
2PL

曲線の傾きが項目毎に変わるのが特徴です。

3 Parameter Logistic Model

$$
\begin{aligned}
P(Y_{ij} | \theta_i,\alpha_j,\beta_j, \gamma_j)
&= \gamma_j + ( 1 - \gamma_j ) \times logit^{-1}(a_j \theta_i - \beta_j) \\
&= \gamma_j + ( 1 - \gamma_j ) \times \frac{1}{1 + \exp(a_j \theta_i - \beta_j)} \\
&= \gamma_j + \frac{( 1 - \gamma_j )}{1 + \exp(a_j \theta_i - \beta_j)}
\end{aligned}
$$

項目特性曲線(3PL)

f:id:dastatis:20191212020111p:plain
3PL

曲線の下限が項目毎に変わるのが特徴です。

4 Parameter Logistic Model

$$
\begin{aligned}
P(Y_{ij} | \theta_i,\alpha_j,\beta_j, \gamma_j, \upsilon_j)
&= \gamma_j + ( \upsilon_j - \gamma_j ) \times logit^{-1}(a_j \theta_i - \beta_j) \\
&= \gamma_j + ( \upsilon_j - \gamma_j ) \times \frac{1}{1 + \exp(a_j \theta_i - \beta_j)} \\
&= \gamma_j + \frac{( \upsilon_j - \gamma_j )}{1 + \exp(a_j \theta_i - \beta_j)}
\end{aligned}
$$

項目特性曲線(4PL)

f:id:dastatis:20191212020322p:plain
4PL

曲線の上限が項目毎に変わるのが特徴です。

5 Parameter Logistic Model

$$
\begin{aligned}
P(Y_{ij} | \theta_i,\alpha_j,\beta_j, \gamma_j, \upsilon_j, \zeta_j)
&= \gamma_j + \left[( \upsilon_j - \gamma_j ) \times logit^{-1}(a_j \theta_i - \beta_j)\right]^{\zeta_j} \\
&= \gamma_j + \left[( \upsilon_j - \gamma_j )\times \frac{1}{1 + \exp(a_j \theta_i - \beta_j)}\right]^{\zeta_j} \\
&= \gamma_j + \left[\frac{( \upsilon_j - \gamma_j )}{1 + \exp(a_j \theta_i - \beta_j)}\right]^{\zeta_j}
\end{aligned}
$$

項目特性曲線(5PL)

f:id:dastatis:20191212020454p:plain
5PL

わかりにくいですが、
曲線の対称性が項目毎に変わるのが特徴です。
(反応確率が0.5を境に非対称になっている)


意外とありましたね、IRT2値ロジスティックモデル。
非対称\zeta_jを混ぜた2パラメータロジスティックモデルとか使ったら楽しそう(ここではやらない)。

では、ここから、5PL(パラメータロジスティックモデル)でパラメータリカバリの性能を確認します。
上の各モデルの項目特性曲線含めて、今回は以下のパラメータの真値からデータを生成しました。
詳しいコードは末尾につけるgithubリンクを見てください。

パラメータの真値

I <- 1000 ##  回答者数
J <- 20 ## 項目数

alpha <- rlnorm(J, meanlog = 0, sdlog = 1) ## discrimination
beta <- stats::rnorm(J, mean = 0, sd = 2) ## difficulty
gamma <- stats::runif(J, min = 0, max = 0.2) ## guessing 
upsilon <- stats::runif(J, min = 0.8, max = 1) ## 
zeta <- extraDistr::rhnorm(J, sigma = 1) ## asymmetric
theta <- rnorm(I, mean = 0, sd = 1) ## trait or ability


統計モデル(Stan)

data{
  int<lower = 0> I;
  int<lower = 0> J;
  int<lower = 0, upper = 1> y[I,J];
  int<lower = 0, upper = 1> mode;
}

transformed data {
  real<lower = 0, upper = 1> watanabe_beta;
  if ( mode == 0 ) {
    watanabe_beta = 1.0;
  }
    else { // mode == 1
    watanabe_beta = 1.0/log(I);
  }
}

parameters {
  vector<lower = 0>[J] alpha;
  vector[J] beta;
  vector<lower = 0.0, upper = 0.2>[J] gamma;
  vector<lower = 0.8, upper = 1.0>[J] upsilon;
  vector<lower = 0>[J] zeta;
  vector[I] theta;
}

transformed parameters {
  real<lower = 0, upper = 1> prob[I,J];
  for (i in 1:I) {
    for (j in 1:J) {
      prob[i,j] = gamma[j] + (upsilon[j]-gamma[j]) * inv_logit(alpha[j] * (theta[i] - beta[j]))^zeta[j];
    }
  }
}

model{
  target += lognormal_lpdf(alpha | 0 , 1) - lognormal_lccdf(0 | 0, 1);
  target += normal_lpdf(beta | 0 , 5);
  target += uniform_lpdf(gamma | 0.0 , 0.2) - uniform_lccdf(0 | 0.0 , 0.2) -  uniform_lcdf(0.2 | 0.0 , 0.2);
  target += uniform_lpdf(upsilon | 0.8 , 1.0) - uniform_lccdf(0.8 | 0.8 , 1.0) -  uniform_lcdf(1.0 | 0.8 , 1.0);
  target += normal_lpdf(zeta | 0 , 1) - normal_lccdf(0 | 0, 1);
  target += normal_lpdf(theta | 0 , 1);
  for (i in 1:I) {
    for (j in 1:J) {
      target += watanabe_beta * bernoulli_lpmf(y[i,j] | prob[i,j]);
    }
  }
}

generated quantities {
  int<lower = 0, upper = 1> y_pred[I,J];
  vector[I] log_lik;
  vector[J] log_lik_temp;
  for (i in 1:I) {
    for (j in 1:J) {
      log_lik_temp[j] = bernoulli_lpmf(y[i,j] | prob[i,j]);
      y_pred[i,j] = bernoulli_rng(prob[i,j]);
    }
    log_lik[i] = log_sum_exp(log_lik_temp);
  }
}

細かいところを見てくださったかたは、お気づきかと思いますが、
事前分布の設定と真値の設定を微妙に変えています(項目特性曲線のときとも少し値変わってます)。
本来は、同様の分布が望ましいと思いますが、今回はデモなのでこうしています。

はい、Stanでベイズ推定しました(雑っ!!!)。

5PLのalphaの真値と推定値の相関

はい、あんまりうまく推定できてないね(推定出来ているはずだったw)。
二値データからパラメータが5つあるモデルはやはり難しいね。
ほら、パラメータリカバリの確認は大事でしょ?w

f:id:dastatis:20191212180552p:plain
5PLa

ちなみに、同真値で1PLモデルのβを見てみると下のようになります。

f:id:dastatis:20191212181052p:plain
1PLb


美しい!!!
こんなに綺麗に相関していれば、ちゃんと実践でも使えそうですね。


最後、ほんまにやっつけですみません。。。
忙しくてw


詳細コードは以下に置きました。

github.com






Enjoy!!

【ポエム】どうして僕がOpen Scienceをやるのか。

こんにちは。

 

この記事は、「Open and Reproducible Science Advent Calendar 2019」の3日目の記事です。そして、ポエムです。

気分を害したくないなら、読まない方が身のためです。

 

adventar.org

 

じゃあ、早速ポエムを始めましょうw

 

僕がOpenScienceを行う理由。

それはとても簡単です。僕は統計ができないからです。

だから、私は常にOpenScienceに大きな関心と注意を注いでいます。

 

現在、僕は心理統計学を研究しています。学部3年のゼミから"専攻"していますので、だいたい、学習歴は6年ぐらい?ですかね。

6年といえば、小学校に在籍する期間と一緒です。

ただ、それだけ学習しても(心理)統計はさっぱりわかりません。

僕の人一倍の要領の悪さといい加減な学習方法が、この原因だと思っています。

 

ただ、さすがに6年も専門的にやっているので、

心理学研究を行っている母集団の中では、ほんとにほんとにほんの少しだけ統計がわかるかもしれません。

だけど、やはり統計学を専攻している、所謂、理系の皆様よりはるかに統計ができません。

唯一の救いはわからないことを知っていることです。

 

だから、僕はOpenScienceを行うことを心がけています。

その理由、ちゃんと説明していきます。

 

ここからの話は、特に心理学研究を想定しています。

(量的/定量的/データを用いた)研究を行う上で、統計学及び、それに関連する知識は必須です。絶対にこれから逃れることが出来ません。

これは、理論・応用研究を問わずにそうです。

 

上でいう、統計学に関連する知識とは、

たとえば、

などが該当するのかなと思っています。

 

場面にもよりますが、

僕は、(心理学に問わず)社会科学を研究する上では、これらの知識は絶対に必須だと考えており、完璧に漏らすことなく知っている必要があると考えています。

あたりまえですが、そんなの無理です。全知全能の神でも無い限り。

少し大げさに書きました。*1

 

僕には、これが不可能だからOpenScienceを行います。

OpenScienceは、統計学のできない僕みたいなのこそ、

真っ先に身につけ、実践すべき研究法です。

 

ここで改めて、OpenScienceとは何かを考えます。

研究全体もしくはその過程を可視化することがOpenScienceです(正確な定義ではないかもしれませんが)。

僕がOpenScienceに関して、真っ先に思いつくのが、研究で収集、分析したデータの公開です。みなさんも同じものを想像したかもしれません。

ただ、これだけがOpenScienceではありません。

  • 調査・実験に至るまで
  • 方法(データ収集、手順の詳細)
  • 解析手法(論文に載せていないものも含む)
  • 等々

上記を正確に記述し、公開することもOpenScienceの一つだと思います。

 

これに実施することによって、科学は発展します(と考えています)。

少なくとも、科学の流れを止めません。止めるのがどうかみたいな話はしません。

 

たとえば、

僕は、OpenScienceの作法に則って、データを集め、解析Aを行い、論文を書きました。

しかし、統計学マスターの謙介さん(仮名)*2が、論文を読んだところ、今回の実験に関して、(解析Aでも間違っていないが)解析Bがより適していることに気付きました。

※この状況では、実際、解析Bを行うことがより正しい、わかることが多いとします。

 

ここで、僕がOpenScienceに取り組んでいなければ、 科学の発展は終わりです。少なくとも停滞します。

  • 僕がOpenScienceに取り組んでいることによって、
  1. 謙介さん(仮名)*3が、同様の手順でデータを取り直すことができ、解析Bを行うことができる
  2. もしデータが公開されているならば、データを取り直すことなく解析Bをその場で実行することができる

というように、僕よりも遙かに優秀な人が効率的に科学を発展させる(社会をより良くする)ことにつながるのです。

もちろん、ここで僕もOpenScienceの作法に則り、貴重な研究基礎を提供したので、このようなコミュニティでは評価されることにつながると思います。

一方で、

  • 僕がOpenScienceに取り組んでいなかった場合、
  1. 謙介さん(仮名)*4が、僕の論文同様の手順でデータを取ることが出来ないかもしれない*5。これによって、解析Bが出来たとしても、その結果と元論文の分析Aの結果を比較可能なのか?という問題も出てくる可能性がある。

というように、僕がOpenScienceを取り組んでいないがために、先ほどよりも研究の発展が難しくなるかもしれません。ここで、私は謙介さん(仮名)の先行研究になり得ていますが、謙介さん(仮名)がOpenScienceの作法に則って研究を行い、解析Aと解析Bを行った場合、そちらの研究の方がより研究貢献をしているような気がします*6

 

従って、

僕がOpenScienceの作法に則り、研究を行えば、今後の誰かの研究の大きな一歩につながるかもしれないということです。そして、自分の研究がその研究に大きく貢献する可能性があるということです。

 

 

少し話しが逸れますが、

心理学の分野によってはOpenScienceが行いにくい分野ももちろんあります。

これには様々な理由が考えられますが、

  1. プライバシーの観点からデータ・手順を公開しにくい
  2. 古くからの慣習
  3. 大御所のおっさんが認めてくれない(多分あとで怒られる)

1.は臨床の分野だと多いかもしれませんね。そして、それを理由にOpenScienceを避けているのをたまに見たような気がしますが、あれは夢だったかもしれません。

まだ、心理学領域でほぼ語られていませんが、データを匿名化する方法はたくさん存在し、これらの方法を用いれば、たとえ臨床心理学で使われるような秘匿性の高いデータでも公開できるかもしれません。

また、医療系では既に匿名化を行うことでカルテ等のデータを共有することを考えているという話も聞いています。

おそらく、2-3年の内に(臨床)心理学でも匿名化の話が出てくると思います。*7

 

 

ここまで書いてきたことを踏まえたうえで、 

 

僕は、統計はできないけど(から)、OpenScienceは行います。

 

 

 

*1:千里の道も一歩から、科学の発展は一日にしてならずという教えをボスから教わり、僕もそれを正しいと思っているので、始めから完璧な研究をやる必要は無いと思いますが、あくまで大げさに言っています。

*2:深い意味は全くありませんw

*3:何度も書きますが、深い意味はありません

*4:何度でも書きますが、深い意味はありません

*5:なんか出来てしまいそうな気がしますがw

*6:あくまで1例で、一つの考え方です。これによって、元論文の価値がなくなるとは思いませんが。

*7:今が狙い目です。

CmdStanRのコンパイルが速すぎたお話

こんにちは。

 

今朝、TwitterのStan好きTLを騒がせていた、

 

 

 

こちらについて簡単に検証してみました。

まだ、開発が始まったばかり程度なので、今後どうなるかわかりませんが、とても期待できそうであることがわかりました。

 

結果から書くと、

 

mc-stan.org

 

に記載のモデル(bernoulli.stan)では、

CmdStanRによるコンパイル時間
Time difference of 1.574471 secs
rstanによるコンパイル時間
Time difference of 44.28613 secs

という結果でした。本当にso fucking fast! で笑いました。

 

ということで、もう少し詳しく見ていこうと思います。

基本的には、上のGetting startedを見て、導入、実践してもらえればと思います。

 

まず、私の環境ですが、OSが Windows 10 Proです。

そして、このおかげ?で、Windowsにうまくcmdstanを導入することが出来ませんでした。

というのも、

 

devtools::install_github("stan-dev/cmdstanr")

library(cmdstanr)

 

 ここまではうまくいくのですが、

 install_cmdstan()

 これがうまくいきませんでした。 

Running bash "C:/Users/hogehoge/Documents/R/win-library/3.6/cmdstanr/make_cmdstan.sh" \ -j2 -w
Error in rethrow_call(c_processx_exec, command, c(command, args), stdin, :
Command not found @win/processx.c:977

 という感じのエラーがはかれます。

恐らく、実行に際し、make_cmdstan.shというファイルを実行して、cmdstanをインストールするようなのですが、そこでこけてしまっています。

 

.shファイルは、bashと呼ばれる、linuxとかで開く前提の拡張子らしいので?(あっている?)そもそも、windowsで開けないっぽい。

windowsで.shファイルを開く方法をググったらできそうな感じはありましたが、それを実行するのが面倒だったので window subsystem for linux (ubuntu LTS 18.04) を使って環境を構築しました。なので、環境は、仮想linux(ubuntu)です。

もちろん、通常?のlinuxAWS,AWS,GCP等のlinux環境でも普通に実行できると思いますが。

また、Rからではなく、自力でcmdstanを導入することもできると思いますが、今回は行いませんでした(面倒だった)。

github.com

 

 

まずは、仮想linuxにRやRstudio Serverを構築しました(ここはCmdStanR関係ないので省きます)。 

調べればたくさん出てきますので、そちらを参考にしてください。

 

そしたら、Rで(僕は、Rstudio Serverでやりました。)

devtools::install_github("stan-dev/cmdstanr")
library(cmdstanr)

install_cmdstan()

先ほど同様のコードを走らせます。

すると、cmdstanがインストールされます。

 

その後、cmdstanのpathを指定します(動作環境によって以下でうまくいかないようですが、その際は公式( https://mc-stan.org/cmdstanr/articles/cmdstanr.html )を見てください)。

set_cmdstan_path()

 

続いて、.stanモデルを指定します。

今回は、パッケージインストールと同時に導入される例題モデルをコンパイルしてみます。

 

## .stanファイルのpathを指定する。

file <- file.path(cmdstan_path(), "examples", "bernoulli", "bernoulli.stan")

 

start <- Sys.time()
mod <- cmdstan_model(file) ##モデルのコンパイル
end <- Sys.time()
print(end-start)

 

上記を実行させれば、モデルがコンパイルされます。

僕の環境では、

Time difference of 1.574471 secs

でした。

 

ちなみに、コンパイルされた例題モデルは

mod$print()

で確認することが出来ます。

上記では、

data {
int<lower=0> N;
int<lower=0,upper=1> y[N];
}
parameters {
real<lower=0,upper=1> theta;
}
model {
theta ~ beta(1,1); // uniform prior on interval 0,1
y ~ bernoulli(theta);
}

が、出力されると思います。

 

あとは、

data_list <- list(N = 10, y =c(0,1,0,0,0,0,0,0,0,1))
fit <- mod$sample(data = data_list, seed = 123, num_chains = 2)

を実行すれば、MCMCのサンプリングを行ってくれます。

VBも出来るそうですが、今回は実行しませんでした(内容の本質ではありませんので)。

 

結果オブジェクトは、rstanの結果オブジェクトとは異なります。

class(fit)

[1] "CmdStanMCMC" "R6" 

 という形式で保存されています。

 

そのため、結果の取り出し方も少し違います。

## 要約統計量

fit$summary()

## MCMCの診断
fit$diagnose()

## MCMC samples (rstanのextract関数)

fit$draws()

 

 ここらへん、少し不便かなと思ったのですが、以下のように実行してあげると、rstanの結果オブジェクトと同様に扱ってあげることができます。

stanfit <- rstan::read_stan_csv(fit$output_files())
print(stanfit)

 

やっていることとしては、

cmdstanrの結果オブジェクトは、.csvファイルで保存されている(される?)ようなので、

これらをrstanに標準装備されているread_stan_csv()関数で読み込んであげれば、結果オブジェクトを復元することができます。

 

あとは、いつものようにrstanの結果オブジェクトと同様に扱うことが可能です。

ここらへんは、今後改善されるてきなことが公式にも記載されていたような気がします。

 

それでは、少し話しを戻して、cmdstanrでコンパイル時間を計算できたので、比較してみましょう。

 

先ほどのモデルコードを事前にbernoulli.stanのようにどこかに保存しておいてください。

比較するので、念のため、僕はRを再起動しました。

library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

start <- Sys.time()
mod <- stan_model("Documents/R/model/bernoulli.stan")
end <- Sys.time()
print(end-start)

あとは、こんな感じで書いてあげれば、コンパイル時間を計測することが出来るかと思います。

Time difference of 44.28613 secs

という結果でした。

 

比較してみると、

cmdstanRによるコンパイル時間
Time difference of 1.574471 secs
rstanによるコンパイル時間
Time difference of 44.28613 secs

so fucking fast! でした(言いたいだけ)。

 

経験上、モデルの違いでコンパイル時間は大きく変化しない気がするので(主にCPUのスペックだと思います)、他のモデルでも早いのではないでしょうか。

 

ちなみに、なぜこんなに早いのかについても公式で言及されています。

The RStan interface (rstan package) is an in-memory interface to Stan and relies on R packages like Rcpp and inline call C++ code from R. On the other hand, the CmdStanR interface does not directly call any C++ code from R, instead relying on CmdStan for compilation, running algorithms, and writing results to output files.

https://mc-stan.org/cmdstanr/articles/cmdstanr.html

 

Rcppを一切介さず諸々をやるので速いようです。

実際に計算しているのも、CmdStanくんが頑張っているようです。

 

速いっていいですね。

今回は、サクッと比較したので、記事に誤字、脱字多いかもしれません。すみません。

 

あと、Windowsで簡単に出来る方法(CmdStanをインストールさえしてしまえばという感じですが)をご存じの方は教えていただけると嬉しいです。

 

Enjoy!!

 

 

 





事後最頻値をとる関数(get_posterior_map)を作りました。

どうも。

この記事は、StanAdventも何も関係ありませんw

 

Stanで事後平均値(EAP)をとってくるときにget_posterior_mean()があると思いますが、事後分布の最頻値(MAP)をとってきたい場合もたまにあるかと思います。

 

僕はその都度R codeを書いていたのですが、めんどくさくなったので、MAP推定値をとってくる関数を書きました。
それ以上でもそれ以下でもない記事ですw



依存パッケージは、dplyrだけです(%>%が使いたいだけ)。

構造:get_posterior_map(stanfit, pars)
stanfit : Stan結果オブジェクト
pars : 取りたいパラメータ

注意:抽出したいパラメータの次元が4次元まで対応しています。

使用したい場合、以下のコードをコピーして実行してください。

get_posterior_map <- function(stanfit,pars) {
  if(!require("dplyr",character.only = T)) install.packages("dplyr", dependencies = T); require("dplyr",character.only = T)
  maxpost <- function(par){ density(par)$x[which.max(density(par)$y)]}
  r <- rstan::extract(stanfit,pars=pars)
  Dim <- r[[pars]] %>% dim %>% length - 1
  if (Dim == 1 ){
    map_est <- array(NA,dim=c(dim(r[[pars]])[2]))
    for (l in 1:dim(r[[pars]])[2])
      map_est[l] <- r[[pars]][,l] %>% maxpost(.)
    return(map_est)
  }else if (Dim == 2 ){
    map_est <- array(NA,dim=c(dim(r[[pars]])[2],dim(r[[pars]])[3]))
    for (l in 1:dim(r[[pars]])[2]) 
      for (m in 1:dim(r[[pars]])[3]) 
        map_est[l,m] <- r[[pars]][,l,m] %>% maxpost(.)
    return(map_est)
  }else if (Dim == 3 ){
    map_est <- array(NA,dim=c(dim(r[[pars]])[2],dim(r[[pars]])[3],dim(r[[pars]])[4]))
    for (l in 1:dim(r[[pars]])[2]) 
      for (m in 1:dim(r[[pars]])[3]) 
        for (n in 1:dim(r[[pars]])[4]) 
          map_est[l,m,n] <- r[[pars]][,l,m,n] %>% maxpost(.)
    return(map_est)
  }else if (Dim == 4){
    map_est <- array(NA,dim=c(dim(r[[pars]])[2],dim(r[[pars]])[3],dim(r[[pars]])[4],dim(r[[pars]])[5]))
    for (l in 1:dim(r[[pars]])[2]) 
      for (m in 1:dim(r[[pars]])[3]) 
        for (n in 1:dim(r[[pars]])[4]) 
          for (o in 1:dim(r[[pars]])[5]) 
            map_est[l,m,n,o] <- r[[pars]][,l,m,n,o] %>% maxpost(.)
    return(map_est)
  }else{
    print("Check Parameter Dimensions")
  }
}


って感じです。
もっとうまく書く方法はあると思いますが、僕はこれでうまく取ってこれるので問題無いです:)

Enjoy!