Microsoft PowerPoint - 夏の学校2018配布用佐々木.pptx

Transcription

1 最尤推定法 点推定量を求める 般的な 法 L n x n i1 f x i 梨 学佐々 邦明 右上の式を θ の関数とみなしたものが尤度関数 尤度関数を最 化する θ の値を最尤推定量とするのが最尤推定法 平均値の推定を例にすると データ (x:3,5,4) が得られたとき, 平均をいくつとするのがよいか? 平均がいくつの分布だったらデータ (x:3,5,4) が得られやすいか? 2 最 化アルゴリズムの考え 周りが見えない中で, 近傍の情報から頂点を目指す 代表的な繰り返し計算法 対数尤度関数の段階的な最 化 1. 初期値を与える 2. 初期値周りで勾配 (1 次微分 ) 等を いて次の推定値の 向を決める 3. 初期値付近で 1 次微分,2 次微分を いて適切に次の点を決めて推定値を得る 4. 収束基準 ( 時微分ベクトル ) で判定し, 収束していない場合は, 現在の値を初期値として 2 に進む ST L n x ˆML 尤度関数を最 化尤度関数の 階微分 =0 を解く Newon-Raphson 法 テイラー展開の 1 次近似を利 して進める 準 Newon 法 (BFGS,L-BFGS 法 ) ヘッセ 列を, パラメータの差分と 回微分の差分を いて逐次近似する. L-BFGS はヘッセ 列の更新式を展開して, 初期値と差分の関数和で表す. 複雑な場合にメモリの節約などが可能 焼きなまし法 (SANN) 適当にパラメータを遷移しつつ, 関数が きくなる遷移を選ぶ. 悪化する遷移もある確率で許容し, 局所解を避ける g() 1 n1 n H g1 H: 尤度関数の 階微分ヘッセ 列 g: 尤度関数の 階微分 4

2 パラメータ推定がうまくいかない そもそも識別不可 収束するとは θ n+1 と θ n が同じになる g 1 が 0 になる 収束しない 無限に繰り返す θ 2 が計算不能 H -1 が計算できない変数が完全相関変数が効 関数に影響しない式形関数の近似状況初期値の問題 推定プログラムに誤り 最 値において唯 解が求まらない可能性がある ( 最 値となるパラメータベクトルが無数にある ) 意思決定者間では異なるが, 選択肢間では異ならない変数 選択肢間では異なるが, 意思決定者間で異ならない変数 異なるモデル間でのパラメータの直接 較は無意味 例 U 1 U 3 U 2 効 に適当な数を しても選択確率に変化はない 効 に適当な倍率をかけても選択確率に変化がない 局所最適解 かけ上の最 化 g() シミュレーションによる尤度計算 シミュレーションによるパラメータの推定 順 1. 誤差項の密度関数から選択肢の数の次元の ( 準 ) 乱数を発 させる 2. この乱数を誤差の値として, 各代替案の効 値を計算 ( 積分 ) する 3. 代替案 i の効 値とその他の代替案の効 との値を 較し, それらの 関係を 1-0 の変数 G で記述する のステップを繰り返す. その 反復回数を R とする. 5. シミュレーションされた確率は となり, この値は不偏推定量である. 1 Pi R R r1 G r 効 を確定値にする 確定的に選択を決定 率を確率に置き換える これを尤度として最 になるようにパラメータをアップデートする 8

3 準乱数の例 準乱数の例としてHalon 数列がある. 計算 法は素数 pに対して s s, s 1 p, s 2 p, s p 1 1 例えば p=3 ならば, 初期値 0 として 1/3, 2/3, 1/9, 4/9, 7/9, 2/9, 5/9, 8/9 多次元化 数列の異なる素数 p を決めて, それぞれに応じて数列を作り多次元化する. 正規分布化 数列を制約つきの乱数発 と同様の変換をして正規分布化 p シミュレーションベースのパラメータ推定法 シミュレーション尤度最 化 (MSL) シミュレーションによって計算された確率を尤度として, 最 化を う. 特性 サンプル数と乱数発 回数に依存する. 乱数発 回数が 分 きいと 致性や漸近的有効性を持ち解析積分と 緒の特性を持つ. 乱数発 回数がサンプル数に対して さく固定されると 致性もない ベイズの定理と事後確率 データ更新に対応するベイズ推定 A( 原因 ) B( 結果 ) の確率 ( 尤度 ) PAB ( ) B( 結果 ) が得られた時に A が原因である確率 = 事後確率 p x z P( BAPA ) ( ) PB ( ) B の起こる確率 p z x p x p z x p x dx A の起こる確率 = 事前確率 ( 主観確率 ) x: 知りたい量,z: データ センサ普及 尤度ストレージ 容量化 事前情報コンピュータ性能向上

4 ベイズでパラメータ推定 パラメータ推定は可能? そもそもはパラメータの事後分布を推定 どうしても点推定 Expeced a Poseriori (EAP) 推定量 モデルが複雑になると積分で... MCMC 法による推定 前の状態に基づいて (Markov Chain) 新しいパラメータをランダムにサンプリング (Mone Carlo) する 最尤推定と何が違うの? 最尤推定 最適化 = どっかで まる MCMC 分布を再現 = いつまでも動く 乱数を使った単純操作の繰り返しで確率密度の きいところを探しながらサンプルを 成する Gibbs Sampling Meropolis-Hasings Sampling ロジット プロビットモデルの MCMC 推定プログラム ( 兵藤 2009 参照 ) 推定例 (MH ロジット ) library(mcmcpack) ### データファイルの読み込み Da<-read.csv("H:/2014/daa.csv",header=TRUE) hh<-nrow(da) ## データ数 :Daa の 数を数える rime <- Da[, 6]/100; bime <- Da[, 9]/100; cime <- Da[,12]/100 rcos <- Da[, 7]/100; bcos <- Da[,10]/100; ccos <- Da[,13]/100 raged <- marix(0,nrow=hh,ncol=1); baged <- 1*(Da[,3]>=6); caged <- raged rcar <- marix(0,nrow=hh,ncol=1); bcar <- rcar; ccar <- 1*(Da[,4]>=2) ## 選択結果 ch <- marix(0,nrow=hh,ncol=3) colnames(ch) <- c("1", "2", "3") for (i in 1:hh){ if (Da[i, 5]==0) ch[i,1] < if (Da[i, 2]==1) ch[i,1] <- 1 if (Da[i, 8]==0) ch[i,2] < if (Da[i, 2]==2) ch[i,2] <- 1 if (Da[i,11]==0) ch[i,3] < if (Da[i, 2]==3) ch[i,3] <- 1 } pos <- MCMCmnl(ch ~ choicevar(rime, "ime", "1") + choicevar(bime, "ime", "2") + choicevar(cime, "ime", "3") + choicevar(rcos, "cos", "1") + choicevar(bcos, "cos", "2") + choicevar(ccos, "cos", "3") + choicevar(raged, "aged", "1") + choicevar(baged, "aged", "2") + choicevar(caged, "aged", "3") + choicevar(rcar, "car", "1") + choicevar(bcar, "car", "2") + choicevar(ccar, "car", "3"), baseline="3", burnin=1000, mcmc.mehod="rwm", b0=0, B0=0, seed=2348, verbose=1000, mcmc=10000, B=0.001) plo(pos) summary(pos) ロジット library(bayesm) ### データファイルの読み込み Da<-read.csv("h:/2014/daa.csv",header=TRUE) hh<-nrow(da) ## データ数 :Daa の 数を数える al <- 3 rime <- Da[, 6]/100; bime <- Da[, 9]/100; cime <- Da[,12]/100 rcos <- Da[, 7]/100; bcos <- Da[,10]/100; ccos <- Da[,13]/100 raged <- marix(0,nrow=hh,ncol=1); baged <- 1*(Da[,3]>=6); caged <- raged rcar <- marix(0,nrow=hh,ncol=1); bcar <- rcar; ccar <- 1*(Da[,4]>=2) cres <- Da[,2] na <- 4 Xa <- cbind(rime,bime,cime,rcos,bcos,ccos,raged,baged,caged,rcar,bcar,ccar) nd <- 0 X <- creaex(al, na=na, nd=nd, Xa=Xa, Xd=NULL, DIFF=TRUE, base=3) da1 <- lis(p=al, y=cres, X=X) mcmc1 <- lis(r=5000,k=1) res1 <- rmnpgibbs(daa=da1, Mcmc=mcmc1) plo(res1$beadraw) plo(res1$sigmadraw) プロビット 兵藤 2009 を いた

5 MCMC 推定のメリットデメリット メリット 解析的にはパラメータが求まらない複雑なモデルもパラメータ分布が求まる 例 : パネルデータの個 モデルを考慮した階層モデル パラメータの分布がやんわりわかる デメリット 時間がかかる MNL の推定例 MCMC:15 秒最尤推定 :7 秒 パラメータの分布が収束しない場合もある 不完全データの最尤推定 EMアルゴリズム 不完全データの最尤推定? EM アルゴリズム 潜在クラスモデル 集団が異質な複数の集団の混合によって構成されていることを想定した上で 観測されたデータを いて分析対象を複数の潜在クラス ( グループ ) に分割し 各クラスの特徴を把握する 法 マーケティング分野で潜在クラス分析を活 する利点の 1 つは, 消費者セグメントの抽出と各セグメントの特徴把握ができること クラスター分析との違い 消費者の購 状況を確率モデルとして表現できることから他の分析モデルと組合せて拡張しやすい. 分析対象となる各顧客の各クラスへの振り分けをクラスへの所属確率として把握できる 混合ガウス分布の例 不完全データに基づく対数尤度に対して いられた最尤推定のための 法 E ステップ (Expecaion-sep): 測値として扱う変数の期待値の推定 個 i がクラス s に属する場合に 1 となる y is の期待値 y is* をベイズの定理を利 して求める M ステップ (Maximizaion-sep): その期待値を いたパラメータ推定 y is* を いてパラメータ (π s と P i s(j) ) を推定する 推定された π s と P i s(j) をもちいて y is* を更新する Pr(X i1 = x i1,,x im = x im ) =Σ W s=1 π s Π M j=1 P i s (j) xij (1 -P i s (j)) 1-xij 2 つのプロセスを繰り返して最尤推定

6 政策シミュレーション 動モデルを推定し, そのパラメータを いて, 変数の変化による選択の変化を る. 回遊 動の分析 マイクロシミュレーションを いることで, 複雑なモデルの組み合わせを政策評価可能 シミュレーションなので, 複数回実施して平均的な評価を う 政策シミュレーション Sep 1 Sep 2 サンプルのデータ, 個 属性や発ゾーン,LOSデータなどを各段階のモデルに個 nの個 属性やLOS, 各種ダミー等といった説明変数データを代 し, 全ての選択肢ごとの選択確率 P in を求める. 求めた選択確率から確率分布 F in を作成する. 0,1 の 様乱数 γ n を発 させ,γ n の値が,F i1n γ n F in を満たす選択肢 iを選択するものとする. 滞在時間モデル等も同様の 法で可能 参考 献 Discree Choice Mehods wih Simulaion K. Train ベイズ統計学松原望 ベイズモデリングによるマーケティング分析照井伸彦 R による離散選択モデルの推定 法メモ兵藤哲朗 マーケティングのデータ分析岡太彬訓, 守 剛 ベイズ推論と MCMC のフリーソフト岩波データサイエンス