9 第 2 章 R によるべイズ統計解析の基本 2.1 はじめに こでは、べイズ統計学の考え方を、 R の簡単な操作によって学んでいきます。前章でも示したとおり、べイズ統計 学は、最初に事前分布を設定し、それを観測された実態に応じて更新していくというものです。そのような考え方を R で追体験できる方法を、 こで実際に操作しながら学んでいきましよう。 2.2 離散的な事前分布を用いる方法 次のような問題で考えます。 例題 あるべルヌーイ試行 ( 成功するか失敗するかの試行 ) において、その成功率の確度が次の割合であるとする。 成功率 0.25 0.5 0.75 0 1 成功率が上記のものである確度 15 15 15 15 真の成功率を確かめるためこの試行を 30 回行ったところ、 12 回成功して 18 回失敗した。このとき、先の確度の 分布を事前分布として設定し、事後分布を求めよ。 まず、事前分布として用いる成功率の確度の分布を R に入力し ます。成功率を p に、確度の分布を prior に代入します。 > p く一 c ( 0 , 0.25 , 0.5 , 0 .75 , 1 ) > prior く一 c ( 1 / 15 , 4 / 15 , 7 / 15 , 2 / 15 , 1 / 15 ) まず確度の分布を図で見るために、 pl 。 t コマンドを使って中身 を見てみましよう。 > plot(p,prior,type="h") これを見ると、成功率がやや低い方に偏っているかもしれない 試行であると考えることができます。さて、 R におけるべイズ統 計の学習用パッケージに、 LearnBayes があります。この中に収録 されているプログラムを使って、べイズ統計を学んでみましよう。 まずはパッケージを読み込ませます。 15 」 0 こ d 0 第 0-6 02 00
2.3 べータ分布を用いる方法 11 このべータ分布を用いると、ベルヌーイ試行の成功率が〃以下である確率がというふうに表すことができます。 下記のグラフは、左からい , の = ( 1 , 1 ) , ( 2 , 5 ) , ( 5 , 2 ) のべータ分布の確率密度関数です。 こで、成功率の事前分布を、次のような設定に従うべータ分布にしてみましよう。 1. 成功率が 0.3 以下となる確率は 0.5 2. 成功率が 0.5 以下となる確率は 0.9 これを満たすべータ分布を求めるためには、コマンド bet a. select を使います。次のように入力すると、べータ分 布のパラメータ住 , を求めることができます。 > quantilel く一 list ( p = 0.5 , x = 0.3 ) # x = 0.3 以下となる確率が 0.5 > quanti1e2 く一 list ( p = 0.9 , x = 0.6 ) # x = 0.6 以下となる確率が 0.9 > beta. select(quantilel ,quanti1e2) [ 1 ] 1 .60 3 .33 これによると、条件を満たすべータ分布は 0 = 1.60 , = 3.33 であることがわかります。これを用いて推計を行いま す。証明は省略しますが、事後分布は、成功した回数を s 、失敗した回数を / で表すと、次のようになります。この性 質がべータ分布をベイズ推定で使う強みです。 ( 1 ー戸 + / ー 1 0 十 s ー 1 ( 2 ・ 2 ) B(a + s, 0 + 月 これを使うと、事後分布は住 = 1.60 + 12 , 0 = 3.33 十 18 のべータ分布になることがわかります。これをプロットす ると次のようになります ( 左 : 事前分布、右 : 事後分布 ) 。 ( 8 一・の 2 一・ 9 - 区 2 を名 の C ・ 94 ) 里 200 1 0 0 02 0 0 0 なお、。 = 1 , 0 = 1 のべータ分布は、事前情報のない事前分布 ( 無情報事前分布 ) によく使われます。またべータ分
第 2 章 R によるべイズ統計解析の基本 10 > 1ibrary(LearnBayes) 次に、成功した回数と失敗した回数をベクトル形式で読み込ませ、これを事後分布を求めるコマンド pdisc に入れ ます。 > data く一 c ( 12 , 18 ) > post く一 pdisc(), prior, data) これで事後分布の予測が終わります。早速、 post の中身を見てみましよう。 > round(post,5) [ 1 ] 0 .00000 0 .17092 0 .82897 0 .00012 0 .00000 事後分布は次のような感じになりました。実験での成功率が 12 / ( 12 十 18 ) = 0.4000 だったので、 0.5 以下にさらに 偏った形になっています。 1 0.5 0.75 0.25 0 成功率が上記のものである確度 0.0000 0.1709 0.8289 0.00 ( ) 0 0.0001 事前分布と事後分布の違いを見てみるために、事後分布を出力してみます。 > plot(p,post,type="h") 2 つの分布を並べてみると、次のようになります ( 左 : 事前分布、右 : 事後分布 ) 。やはり、事後分布は 0.4 の近くに かなり偏ることが分かるはすです。 成功率 lSOd 」 0 に d 1 -0 0 8 0-6 0 4 0 、 2 0 ℃ 0 第 0 毛 0-4 0-2 2.3 べータ分布を用いる方法 先ほどの作業では事前分布が離散的な場合を使用しましたが、続いて連続的な場合を使います。この例題のような問 題を解く場合、よく使われる分布がべータ分布です。べータ分布は、次の確率密度関数で示される分布です。 0 ー 1 ( 2 ・ 1 ) ただし、 B い , の一 / ( 泳いの = ーソフトで始めるべイズ統計解析 Bayes Analysis Maniax フリ
第 2 章 R によるべイズ統計解析の基本 12 布は、二項分布に対応する尤度関数であるため、このように計算することが可能なのです。 2.4 ヒストグラム事前分布を用いる方法 べータ分布を使うと、任意の離散的な分布から連続的な分布を求めることも可能です。再度、冒頭の例題で用いた事 前分布を持ち出してみます。 成功率 1 0 0.25 0.5 0.75 成功率が上記のものである確度 15 15 本節では予測も行うので、この確度を別の文字 midpt に入れることにします。 > midpt く一 c ( 0 , 0.25 , 0.5 , 0 .75 , 1 ) > prior く一 c ( 1 / 15 , 4 / 15 , 7 / 15 , 2 / 15 , 1 / 15 ) > curve (histprior(x,midpt ,prior) ) 下の行の curve コマンドを使ってこの分布を描くと、次のようなものになります。 15 15 15 0 住 6 0 4 0.2 0.0 これに試行が成功した回数が 12 回、失敗が 18 回である情報を読み込ませ、事後分布を示します。事後分布は、二項 分布に対応する尤度関数であることから、事前分布に住 = 12 + 1 , 0 = 18 + 1 のべータ分布を掛け合わせれば求めるこ とができます。 事後分布の確率密度関数の形だけを知りたいならば、 curve コマンドと histprior コマンドを使って、次のように 入力すれば図として出すことができます。 > curve(histprior(x,midpt,prior)*dbeta(x, 12 + 1 , 18 + 1 ) ) これによって出した図が次の通りです。 Bayes Analysis Maniax フリ ーソフトで始めるべイズ統計解析
2.4 ヒストグラム事前分布を用いる方法 13 96 (L + 9 二一 + 二 x ) 2 qp . ( 」 0 言言翌区」 0 一」 dlS 丘 0.4 02 0 0.8 0 何か奇妙な形になってしまいましたが、形自体は最初に求めた事後分布とだいたい似通っているかと思います。この 事後分布を使ってシミュレーションをしてみましよう。ここでは 1 , 000 個の標本を使ってシミュレーションをしてみま す。 histprior,sample の各コマンドを使って次のように行います。 > p く一 seq(), 1,1ength=1000) > post く一 histprior(p,midpt,prior) * dbeta(), 12 + 1 , 18 + 1 ) > post く - post/sum(post) > ps く一 sample (p ,rep1ace=TRUE,prob=post) このシミュレーション結果をヒストグラムに示してみると次のような形になります。 09 003 な ua コ ba 」」 00 【 0 4 02 0.0 0 名 0.6 こうして見ると、成功率はやはり 0.4 と 0.5 の間にあることが伺えます。
4.4 描画に制限を加える 4.4 描画に制限を加える R によるべイジアンネットワークの描画においては、計算を行わない因果関係を指定することも可能です 37 は、 net2 という文字に、計算を行わない因果関係を含んだべイジアンネットワークの計算を行うネットワークを入れ ることにしましよう。 > net2 く一 network(datasetl) > net2—prior く一 jointprior (net2) lmaginary sample size : 162 本節の分析では、分析の最終的な目的を安倍政権を評価するかどうかに絞るため、「安倍政権」から他の項目への矢 印が出ないようにします。このデータセットでは、「安倍政権」の項目は 2 列目に入力されているため、次のように入 力を行います。項目名ではなく項目が入っている列の番号というのがちょっと厄介ですね・・ ンド banlist で指定を行います。 > net2 > net2 [ 2 , ] [ 3 , ] -ban く一 as . matrix(rbind(c(2,1),c(2,3),c(2,4))) ban [ , 1 ] [ , 2 ] 2 2 2 1 3 4 > ban1ist(net2) く一 net2—ban > p10t(net2) learn,autosearch コマンドを使って次のように行います。 べイジアンネットワークの描画は、計算の制限を行わない場合と同様 を行わない因果関係についてチェックしておきましよう。 たまに成功しないこともあるので、分析を行う前には、このように描画 ているはずです ( 本書は白黒なので分かりづらいかもしれませんが・ plot コマンドで描画してみると、計算しないところが赤い点線になっ 0 0 。行列をつくって、コマ 0
48 第 5 章 KH Coder でのべイズ統計 閑話休題、この分析では 65 冊の本をいくっかの分類に分けていますが、その中の分類として時期があります。これ は、社会の変化に応じて、 65 冊の著作を時期区分に分けたものです。具体的には次のような時期になります ( 「普及版 「劣化言説の時代」のメディアと論客」 pp. 73-74 ) 。 第 1 期 ( ~ 2003 年 ) / 第 2 期 ( 2004 ~ 2006 年 ) 香山的「解離」言説の社会・政治への「適用」・・・ 2004 年に 「就職がこわい」私〉の愛国心」「生きづらいく私〉たち」 ( それぞれ、講談社、ちくま新書、講談社現代新 書 ) が刊行されるが、これらの著書は、従前より香山が「インターネット・マザー』や「多重化するリアル』 などで述べてきた、「解離」概念を中心とする若者論を、それぞれ労働・雇用問題、政治、そして若年層の心 理の掘り下げという各方面に「適用」するようになった。 第 3 期 ( 2007 年 ~ 2011 年 3 月 ) 劣化言説へのコミットと活動範囲の広がり・・・ 2007 年に「なぜ、あの人は仕 事中だけ「うつ」になるのか」の原書である『仕事中だけ「うつ病」になる人たち一一 30 代うつ、甘えと自 己愛の精神分析』 ( 講談社 ) と、「なぜ日本人は劣化したか』 ( 講談社現代新書 ) を上梓し、本格的に若い世 代に対する「攻撃」を開始するようになる。他方で、「しがみつかない生き方一一 - 「ふつうの幸せ」を手に 入れるルール」「人生の法則ー一知るだけでココロがらくになる 10 章」 ( 幻冬舎新書、 2009 年 / ベスト新書、 2010 年 ) などの自己啓発書や、「しがみつかない死に方一一一孤独死時代を豊かに生きるヒント」「いのち問 答 - ーー最後の頼みは医療か、宗教か ? 」 ( 後者は対本宗訓との共著 / 共に角川 One テーマ 21 / それぞれ 2010 年、 2011 年 ) などといった老後に向けた死生観関係の著作も上梓されるようになっている。第 2 期に上梓 された「スピリチュアルにハマる人、ハマらない人」 ( 幻冬舎新書、 2006 年 ) が、「スピリチュアルにハマ る」側としての「若者」を批判的に採り上げる著作だったのが、老後の不安などに「スピリチュアル」的な ものに " ハマる " ことは、これらの著作では肯定的に評価されているのが特徴である。 第 4 期 ( 2011 年 4 月 ~ 2012 年 6 月 ) 東日本大震災の「不安」に応える・・・ 2011 年 3 月 11 日に東北地方太平洋 沖地震が発生、東日本大震災と呼ばれる様々な被害をもたらすが、 2011 年 4 月に講談社現代新書より刊行 された「〈不安な時代〉の精神病理」 ( 講談社現代新書 ) には、 ( 著作の執筆自体は震災以前から用意されてい たものであろうが ) 既に東日本大震災に関する記述がまえがき・あとがきを中心に見られる。また 2011 年 5 月から翌年にかけて、第 2 ・ 3 期に香山が一般書として刊行した自己啓発書が次々と文庫化されているが、 これらにも東日本大震災を意識して加筆されたものが多い。そこで第 4 期に限り、この時期に文庫として再 版された著作もこの時期のカテゴリに入れるものとする。ただし、 2011 年 6 月に刊行された「なせ、あの 人は仕事中だけ「うつ」になるのか」については、原書が第 3 期の入りを特徴付けを象徴する著作と見なし たことを優先し、例外的に第 3 期に入れる。 第 5 期 ( 2012 年 10 月 ~ ) 政治・社会評論への復帰・・・ 2012 年 10 月に香山は「絆ストレス 「つながりた い」という病』 ( 青春新書 intelligence) を刊行するが、そこでは ( 震災後に香山自身もまた「 3 ・ 11 後の心 を立て直す』 ( ベスト新書、 2011 年 ) などで煽っていた ) 「心のつながりの復活」などの言説を批判した本を 上梓し、その後は急速に震災前の言説に回帰していった。政治には「「独裁」入門』 ( 集英社新書、 2012 年 ) 、 若者論には「若者のホンネーーー平成生まれはなにを考えているのか』 ( 朝日新書、 2012 年 ) で復帰する。ま た、 2014 年の「弱者はなぜ救われないのか』「ソーシャルメディアの何が気持ち悪いのか』「劣化する日本 人」 ( それぞれ、幻冬舎新書、朝日新書、ベスト新書 / すべて 2014 年 ) では「リべラルの責任」ということ をしきりに述べるようになるが、特に第 3 期以降に通俗保守的な価値観に立って若者論や劣化言説を産出し てきた香山にそのことを言い立てることはできるのか、と問い直す必要はあるだろう。 この 5 つの時期分類に基づいて、 KHCoder のべイズ学習機能を使ってみようと思います。 KHCoder のべイズ学 習機能は、「ツール > 文書 > べイズ学習による分類 > 外部変数から学習」を使います。ここで学習に使う外部変数を選 択し、「交差妥当化を行う」にチェックを入れて学習を行います。使用する単語は、全体での占有率が 20 % になる、出 現率 264 以上の自立語です ( MeCab の辞書は同書で使ったものをそのまま使用 ) 。 Bayes Analysis Maniax フリ ーソフトで始めるべイズ統計解析
Ca11 : lm(formula = Y Min Residua1s : 20 ー 1 . 61 0 . 00 0 . 01 0 . 13 0 . 18 0 . 21 0 . 42 0 . 05 1 . 32 10.92 15.97 0 .417 seconds (Tota1) 0.185 seconds (Samp1ing) E1apsed Time : 0 . 232 seconds (Warm—up) lteration: 2000 / 2000 [ 100 % ] (Samp1ing) 第 3 章 RStan によるべイズ統計解析 と、このようにべイズ推計のための MCMC ( マルコフ過程モンテカルロ法 ) が行われました ( 実際にはもっと多く の分析が行われている ) 。これによる結果を、通常の回帰分析の結果と比較してみましよう。 > summary(test—lm) data = dataset) IQ Median 0 .04038 3Q 0 . 17107 Max 0 .26980 ー 0 . 91160 ー 0 .08164 Coefficients : Est imat e (lntercept) ー 1 .87193 0 . 15368 Signif . codes: Std . 0 . 001 14 .83 0.01036 ー 15 .08 0 . 12411 Error t value Pr ( > は l) 1 . 18e ー 11 1 .55e ー 11 0 . 05 freedom 1 Residua1 standard error: 0.2672 on 18 degrees Of Mu1tip1e R—squared: 0 .9244 , Adjusted R—squared: F—statistic: 220 0 Ⅱ 1 and 18 DF, p—value : 1 .553e > dataset_fit lnference for Stan model : stan_testl . 4 chains , each with iter=2000; warmup=lOOO; thin=l ; post—warmup draws per mean S e _mean Chain sd 2 . 5 % 0 .9202 ー 11 dra s = 4000. 97 . 5 % n—eff Rhat = 1000 , total post—warmup 25 % 50 % 75 % b ー 1 .87 0 . 15 sigma 0 . 29 14 .46 lp— 0 . 00 0 .04 0 . 00 0 . 13 ー 2 . 15 ー 1 .96 ー 1 .87 ー 1 .79 0 . 15 0 . 25 13.86 0 . 15 0 . 28 14.82 0 . 16 0 . 32 15.41 1664 1600 1537 1209 1 1 1 1 Samp1es were drawn using NUTS(diag-e) at Tue Jul 25 20 : 57 : 42 2017. FO て each parameter, 取 eff crude measure Of effective sample Size , and Rhat is the potential scale reduction factor 0 Ⅱ split chains (at convergence , Rhat=1) . また、 Stan では、 MCMC がちゃんと収束したかを調べることもできます。パッケージ g cmc を使って、次のよう 率は 2.5 % 、という具合に ) 。 モデリングでは、パラメータは平均値及び実際のパラメータがある確率が示されます ( 例えば値が 2.15 以下になる確 通常の回帰分析では、各パラメータの値と標準偏差くらいしか分からないのに対して、 RStan を使ったべイズ推計の に行います。 Bayes Analysis Maniax フリ ーソフトで始めるべイズ統計解析
56 奥付 タイトル Bayes Analysis Maniax ごとうかすとも 著者後藤和智 発行者後藤和智事務所 OffLine 印刷大陽出版 参考文献 とうけいかいせき フリーソフトで始めるべイズ統計解析 発行日 2017 年 8 月 13 日 ( コミックマーケット 92 ) 連絡先 kgot01984@nifty.com 著者ポータルサイト http://www45.atwiki.jp/kazugoto/ サークルプログ http://kazugoto.hatenablog.com/ CircIe. ms プロフィールページ httP://C10000088.circle.ms/oc/CircleProfile.aspx ニコニコ生放送コミュニティ http:/.com/nicovideo.jp/community/c01654257 ニコニコチャンネル http://ch.nicovideo.jp/channel/kazugoto Facebook サークルページ https://www.facebook.com/kazugotooffice twitter https://twitter.com/kazugoto pixiv 5559346
参考文献 5. 松浦健太郎「 stan と R でべイズ統計モデリング (Wonderful R ・ 2 ) 』共立出版、 2016 年 ン、 2012 年 4. J. アルバート : 著、石田基広、石田和枝 : 訳「 R で学ぶべイズ統計学入門」丸善出版 / シュプリンガー 3. 松原望「入門べイズ統言 t ーー意思決定の理論と発展」東京図書、 2008 年 2. 涌井良幸、涌井貞美「 Excel でスッキリわかるべイズ統計入門』日本実業出版社、 2010 年 1. 小島寛之「完全独習べイズ統計学入門』ダイヤモンド社、 2015 年 53 シャノヾ 6. 樋口耕一「社会調査のための計量テキスト分析ーー内容分析の継承と発展をめざして」ナカニシャ出版、 2014 年 7. 「散布図行列を描くには」 http://statmodeling.hatenablog.com/entry/scatter-plot-matrix 8. 「 Rstan でべイズ分析を行うための環境を作る方法メモ」 http://hikaru1122.hatenadiary・jP/entrY/2015/08/04/230000 9. 「 RStan をはじめよう」 https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started-(Japanese) 10. 「 R でべイジアンネットワークメモ」 https.//rpubs.C0m/h0X0ーm/21327 11. 今井俊輔、岡本一志「 R によるべイジアンネットワーク入門」 https : / /www.slideshare.net/0kamot0-laboratory/r-70898327 12. 里洋平「ネットワーク分析 - べイジアン・ネットワーク」 http : / / d. hate Ⅱ a. ne. jp / yokkuns / 20110928 / 1317164851 13. 後藤和智「改訂増補版紅魔館の統計学なティータイム OffLine 、 2013 年 ( コミックマーケット 85 ) 市民のための統計学 Specia12 』後藤和智事務所