差分

このページの2つのバージョン間の差分を表示します。

この比較画面へのリンク

次のリビジョン		前のリビジョン
informatics:2011 [2011/06/16 13:27] – created wataluinformatics:2011 [不明な日付] (現在) – 外部編集 (不明な日付) 127.0.0.1
行 1: 行 1:
 +
 +
 +==== 経営情報学コース分の課題 ====
 +
 +|西先生|これから40年どんな手段で世の中を良くするのか？技術や考え方をビジョン化し、具体的に説明しなさい|
 +|山田先生|講義内での課題を、A4一枚ぐらいまで膨らませて|
 +|水戸先生|[[http://wlgate.inf.uec.ac.jp/j-kiso/|ここにおいてある標示資料]]の一番最後のページ|
 +|由良先生|今日の配付資料のとおり|
 +|山本|コインの表が出る確率は1/2としても、画鋲の針が上になる確率はいくらか、画鋲を30回、投げてみよう ([[http://stat.inf.uec.ac.jp/library/informatics.2011/note-20110616.pdf|配布資料の最後のページ]]と下のコード)|
 +
 +|提出期限|2011年6月30日(木) 0500pmまで|
 +|提出場所|西五号館3階 総合情報学専攻事務室|
 +
 +宜しくお願いします。
 +
 +==== 画鋲を投げて針が上向く確率の推定 ====
 +
 +{{:informatics:529569a.jpg?200|}} 針が上になっている状態を表、下になっている状態を裏とする。
 +
 === モーメント法・最尤法 === === モーメント法・最尤法 ===
-データの与え方。+これはRを用いて計算が可能である。 
 + 
 +== データの与え方 ==
  
 <code> <code>
行 6: 行 27:
 </code> </code>
  
-点推定と区間推定。+== 点推定と区間推定 ==
 <code> <code>
 n <- length(X.raw) n <- length(X.raw)
行 18: 行 39:
  
 === ベイズ法 === === ベイズ法 ===
 +ベイズ法に関する計算、特に事後分布の算出は、WinBUGSというソフトウェアを用いるのが簡単である。
  
-データの与え方。+== データの与え方 ==
 <code> <code>
-DATA list {n=32, x=c(1,0,1,0,0,1,1,1,0,1,0,1,1,0,1,0,0,0,1,1,0,0,1,1,1,0,0,1,0,0,1,0)+DATA list(n=32, x=c(1,0,1,0,0,1,1,1,0,1,0,1,1,0,1,0,0,0,1,1,0,0,1,1,1,0,0,1,0,0,1,0))
 </code> </code>
 +
 +== 点推定と区間推定 ==
  
 モデルと事前分布と初期値の指定。 モデルと事前分布と初期値の指定。
 <code> <code>
-MODEL model {+MODEL model{
   for ( i in 1:n ) {   for ( i in 1:n ) {
     x[i] ~ dbin(p,1)     x[i] ~ dbin(p,1)
   }   }
-  p ~ dgamma(2,3) # objective prior  +  p ~ dbeta(2,3) # objective prior  
-  # p ~ dgamma(1,1)  # flat prior+  # p ~ dbeta(1,1)  # flat prior
 } }
 INIT list {p=0.5} INIT list {p=0.5}
 </code> </code>
  
-ここから先の手順は[[http://web.sfc.keio.ac.jp/~kogure/seminar/07fall/2/WinBUGS_usage.pdf|このマニュアル]]に頼る。+ここでベータ分布はパラメータを変えると、次のように動く。 
 + 
 +{{:informatics:beta-distributions.jpg?400|}} 
 + 
 +ここから先の手順は、多少複雑で[[http://web.sfc.keio.ac.jp/~kogure/seminar/07fall/2/WinBUGS_usage.pdf|この手引き]]に頼るとよい。手順だけ記すと、次の通り。 
 + 
 + 
 +  - 事後分布を得るためのWinBUGSのコードを読み込ませる (〔Model〕→〔Specification〕) 
 +  - モデルを読み込ませる (上で書いたコードのうちmodelをマウスで指定してから、〔Specification Tool〕→〔check model〕、成功すると ``model is syntactically correct''と表示される) 
 +  - データを読み込ませる (上で書いたコードのうち、データの記述に相当するlistをマウスで指定してから、〔Specification Tool〕→〔load data〕、成功すると ``data loaded'') 
 +  - 走らせるマルコフ連鎖の本数を指定する (〔Specification Tool〕→〔num of chains〕) 
 +  - コンパイルする (〔Specification Tool〕→〔compile〕、成功すると ``model compiled'') 
 +  - 初期値を指定する場合、読み込む (〔Specification Tool〕→〔load inits〕、マルコフ連鎖の本数だけ実行する必要あり、不足分があると ``this chain contains uninitialized variables.'') 
 +  - 初期値が不足している分は生成させる (〔Specification Tool〕→〔gen inits〕) 
 +  - 事後分布を得るためにマルコフ連鎖モンテカルロ(MCMC)を走らせる準備をする (〔Model〕→〔Update〕、〔Inference〕→〔Sample〕 
 +  - モニターしたい変数をSample Monitor Toolに登録する (〔Sample Monitor Tool〕→〔node〕にすべての変数を、ひとつひとつ手で入力し、〔set〕) 
 +  - まずは1000歩、走らせてみる (〔Update Tool〕でupdatesに1000とタイプしてから、〔update〕) これはバーンインと呼ぶ準備である。多次元空間を一度に一次元方向ずつだけ移動する点の列なので、初期値の影響がなくなるぐらいまでまず走らせる (変数の数にもよるが、1000程度以上か) 
 +  - 次にまた9000歩、走らせてみてから、historyを眺めて、定常状態になったことを確認する。定常分布に収束できるとは限らないので、複数のマルコフ連鎖を走らせるほうが良い (chains>=2)、 各変数の履歴(history)を見て、定常分布に至っていることを確認する (ホワイトノイズに見えることが望ましい) などの注意を憶えておく。 
 +  - 本格運転 (〔Update Tool〕でマルコフ連鎖の歩を進め、〔Sample Monitor Tool〕で様子を眺める) 
 +  - 十分な長さだけ、マルコフ連鎖モンテカルロを走らせる ($50000$, $100000$, などともかく十分な長さ) 
 +  - 履歴(history)、直近の履歴(trace)、密度関数(density)、分位点の推移(quantiles)、などを眺めて、定常分布となっていることを確認 (〔Sample Monitor Tool〕の〔node〕欄に ``*'' を入力してから、各ボタンをクリックするのが便利) し、自己相関(auto cor)も確認し(妙なパターンがなく、滑らかに減少していれば良い、速く減少していればもっと良いがMCMCなので仕方ない)てから、事後分布を拾う 
 +  - 定常と判断した範囲(beg〜end)までについてstatsを表示させる (〔Sample Monitor Tool〕でbegの値を1ではなく、定常と判断した歩数に変更)
  
 === ブートストラップ法 === === ブートストラップ法 ===
  
-点推定は、モーメント法・最尤法と同じ。データの与え方も同じ。信頼区間の作り方は、+点推定は、モーメント法・最尤法と同じ。データの与え方も同じ。下記のコードで、ブートストラップによる分散の推定と、信頼区間の構成まで行える。
  
 <code> <code>
行 58: 行 103:
 </code> </code>
  
 +これは点推定の方法は正しくても、分布の仮定がそれほど正しくないかもしれない時に、一番最初の信頼区間とは異なった範囲を示す。
 +
 +{{:informatics:bootstrap-sample-distribution-of-p.jpg?400|}}
 +
 +ただし、ブートスラップ反復回数 B の設定には要注意。
 +
 +=== ベータ分布の図 ===
 +<code>
 +# Beta(1,1)
 +plot(c(1:99)/100,dbeta(c(1:99)/100,1,1), xlab="p", ylab="prior probability", type="l", lty=1, xlim=c(0,1), ylim=c(0,4.5))
 +# Beta(3,2)
 +lines(c(1:99)/100,dbeta(c(1:99)/100,3,2), xlab="p", ylab="prior probability", type="l", lty=2)
 +# Beta(2,3)
 +lines(c(1:99)/100,dbeta(c(1:99)/100,2,3), xlab="p", ylab="prior probability", type="l", lty=3)
 +# Beta(10,5)
 +lines(c(1:99)/100,dbeta(c(1:99)/100,10,5), xlab="p", ylab="prior probability", type="l", lty=4)
 +# Beta(5,10)
 +lines(c(1:99)/100,dbeta(c(1:99)/100,5,10), xlab="p", ylab="prior probability", type="l", lty=5)
 +# Beta(15,5)
 +lines(c(1:99)/100,dbeta(c(1:99)/100,15,5), xlab="p", ylab="prior probability", type="l", lty=6)
 +# Beta(5,15)
 +lines(c(1:99)/100,dbeta(c(1:99)/100,5,15), xlab="p", ylab="prior probability", type="l", lty=7)
 +
 +</code>