The One with ...

研究用のメモなどを公開しています

回答のランダム化とStanによる離散パラメータ推定

Stanによる離散パラメータの推定

離散パラメータが含まれている場合のMCMC推定は基本JAGSでやってきたのだが,新しい方法を勉強するためにStanで周辺化消去を試してみる. 以下StanとRコードは,松浦さんのアヒル本11章に依拠しています.

回答のランダム化

被験者が答えにくい質問への回答の母比率を推定するための手法に,「回答のランダム化」という方法がある.これは次のような方法だ.

まず被験者がコインを振る.調査者はその結果を見ない.コインの裏が出たら,被験者は質問に常にYESと答え,表が出たら真実の回答を述べる.自分が真実の回答を述べたか,コインによってYESと答えたのかを,調査者に知られることはないので,被験者は安心して回答できる.よく工夫された設計である.

この調査法は規範的・倫理的に真実の回答を表明しにくい質問に適している.例えばローゼンタールのテキスト『運は数学にまかせなさい』では,マリワナ喫煙者の割合の推定例としてこの手法を紹介している.

人工データの作成

推定がうまくいったかどうかは,そもそも真の分布や真のデータ生成メカニズムが特定できていないと判定のしようがない.そこでまずR上で推定すべき「真のデータ」を生成する.

Rにはベルヌーイ分布の確率変数を生成する関数がどういうわけか用意されていないので,まずは関数bern( )を作成し,それを使ってデータを作る.

人工データを作成する部分の説明

まずbern(0.5)を使い,0 or 1の乱数を発生させる.0がでたら必ず1を出力し,1が出たら確率q0で1を出力(1-q0で0を出力)する.これを1000人分繰り返してデータYとして保存する.

Stanコードはアヒル本記載のコードをそのまま使っています.この例では推定すべき真のパラメータ,つまり母比率はq0=0.85であると仮定している.

Stanの結果は以下の通り.

> fit
Inference for Stan model: model11-1.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

        mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
q       0.86    0.00 0.02    0.83    0.85    0.86    0.87    0.89  2636    1
lp__ -253.67    0.02 0.73 -255.81 -253.80 -253.40 -253.22 -253.17  1439    1

Samples were drawn using NUTS(diag_e) at Fri Jul 07 16:48:00 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

結果qの中央値と平均が0.86で,95%ベイズ信頼区間は(0.83~0.89)となった.

確率論による母比率の(解析的)推定

1000人の回答者に「マリワナ喫煙の経験がありますか」と質問して,コインによってランダム化された回答データがあるとしよう.表が出た場合のみ本当の回答をして,裏が出た場合は常に「はい」と答えてもらう.このとき,MCMCで推定しなくても,母比率の点推定は可能である.

考え方はこうだ.まず集合を4つ定義する.

 A=\{ コインで表が出た人\}

 B=\{ コインで裏が出た人\}

 C=\{ Aの中で実際にマリワナを吸った人\}

 D=\{ Bの中で実際にマリワナを吸った人\}

データを集計したところ,質問に「はい」と回答した人数が930だったとする.つまり,

|B \cup C|=930

である.また大数の法則により,コインは大体半分は裏が出ていると考えられるので|B|=500と仮定する.同様に |A|=500と仮定する.またコインの裏表がランダムに出ているから,

 \displaystyle \frac{|C|}{|A|}=\frac{|D|}{|B|}

と仮定できる.要するに,表が出た人々における「実際のマリワナ喫煙者割合」と,裏が出た人々における「実際のマリワナ喫煙者割合」は等しいはずである.

この値は,全体における「マリワナ喫煙者の割合」に一致する.

 C の要素数が分かれば,\frac{|C|}{|A|}がすぐ分かるので, |C|を計算する.

まず

|C \cup B|-|B| =|C| \quad (B \cap C=\emptyset より)

である.つまり930-500 = |C| だから 430 =  |C|である.

よって, \displaystyle \frac{|C|}{|A|}=\frac{430}{500}=0.86
である.

この数値はベイズ推定による事後平均およびMAP推定値と一致している.