-
Notifications
You must be signed in to change notification settings - Fork 5
/
18.html
888 lines (800 loc) · 37.6 KB
/
18.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
<!doctype html>
<html lang="ja">
<head>
<meta charset="utf-8">
<title>プログラマの為の数学勉強会</title>
<!-- For reveal.js -->
<link rel="stylesheet" href="lib/reveal/css/reveal.css">
<link rel="stylesheet" href="lib/reveal/css/theme/night.css">
<link rel="stylesheet" href="lib/reveal/lib/css/ir_black.css">
<!-- For Graphics -->
<link rel="stylesheet" href="css/graphics.css">
<style>
.reveal .chapter-title {
margin-top: 3em;
}
.reveal {
font-size: 32px;
line-height: 1.4em;
}
.reveal .slides {
text-align: left;
}
.reveal section img {
border: none;
background: 0;
margin-left: 1em;
margin-right: 1em;
box-shadow: none;
}
.reveal strong {
color: yellow;
}
.reveal sup {
font-size: 40%;
}
.reveal table {
margin-top: 0.5em;
margin-bottom: 0.5em;
border: 2px solid lightblue;
}
.reveal pre {
font-size: 0.7em;
}
.reveal pre code {
max-height: 600px;
}
.reveal .note {
font-size: 50%;
}
.reveal .controls div.navigate-up,
.reveal .controls div.navigate-down {
display: none;
}
.reveal .block {
border: solid 2px;
position: relative;
border-radius: 8px;
margin-top: 0.8em;
margin-bottom: 0.8em;
padding: 1em 0.8em 1em 0.8em;
}
.reveal .block:after {
content: "";
display: block;
clear: both;
height: 1px;
overflow: hidden;
}
.reveal .answer {
color: #111111;
}
.reveal .block h4 {
position: absolute;
top: -0.5em;
margin: 0 auto;
background: #111111;
font-weight: bold;
}
</style>
<!-- Setup libraries for RequireJS-->
<script src="lib/require.js"></script>
<script>
requirejs.config({
baseUrl: "js",
paths: {
d3: "../lib/d3/d3.v3.min",
numeric: "../lib/numeric-1.2.6",
MathJax: "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS_HTML"
},
shim: {
d3: { exports: "d3" },
numeric: { exports: "numeric" },
MathJax: { exports: "MathJax" }
}
});
</script>
<!-- Initialize MathJax -->
<script type="text/x-mathjax-config">
require(["MathJax"], function (MathJax){
MathJax.Hub.Register.StartupHook("AsciiMath Jax Config", function () {
var AM = MathJax.InputJax.AsciiMath.AM;
AM.symbols.push(
{input:"mathbi",tag:"mstyle",atname:"mathvariant",atval:"bold-italic",
output:"mathbi",tex:null,ttype:AM.TOKEN.UNARY}
);
});
MathJax.Hub.Config({
showProcessingMessages: false,
skipStartupTypeset: false,
tex2jax: {
inlineMath: [ ["\\(","\\)"] ],
displayMath: [ ["\\[","\\]"] ]
}
});
});
</script>
<script>
</script>
</head>
<body>
<div class="reveal">
<div class="slides">
<section style="text-align: center">
<h1> プログラマの為の<br>数学勉強会<br>第18回</h1>
<span>
(於)ワークスアプリケーションズ<br>
中村晃一<br>
2014年1月30日
</span>
</section>
<section>
<h2>謝辞</h2>
<p>
この会の企画・会場設備の提供をして頂きました<br>
㈱ ワークスアプリケーションズ様<br>
にこの場をお借りして御礼申し上げます。
</p>
</section>
<section>
<h2> この資料について </h2>
<p>
<ul>
<li> <a href="http://nineties.github.com/math-seminar">
http://nineties.github.com/math-seminar
</a>に置いてあります。 </li>
<li> SVGに対応したブラウザで見て下さい。主要なブラウザで古いバージョンでなければ大丈夫だと思います。</li>
<li> 内容の誤り、プログラムのバグは<a href="http://twitter.com/9_ties">@9_ties</a>かkoichi.nakamur AT gmail.comまでご連絡下さい。</li>
<li> サンプルプログラムはPython及びMaximaで記述しています。 </li>
</ul>
</p>
</section>
<section>
<h2 class="chapter-title"> ベイズ統計 </h2>
</section>
<section>
<p>
今回はベイズ主義に基づく<strong> ベイズ統計学 </strong>の基礎事項を紹介します.
</p>
<p>
復習すると, ベイズ主義では確率 \(P(A)\) を 事象 \(A\) が起こるという事の確信度の大きさとして解釈するのでした. そして<strong>ベイズの定理</strong>が中心的な役割を果たします.
</p>
</section>
<section>
<p>
ベイズ統計では, 考察の対象 \(A\) に対して予め<strong>事前確率</strong>\(P(A)\)を定めておきます. \(A\) と関係する事象\(B\)が生起したならば <strong>ベイズ改訂</strong> を行って <strong>事後確率</strong> \(P(A|B)\) へと確率を更新する事が出来ます. これを繰り返せば確率の客観性を高めていく事が出来ます.
</p>
<div align="center"> <img width="600px" src="fig/bayesian-revision.png"> </div>
</section>
<section>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> ベイズの定理 </h4>
<p>
\(P(A),P(B)>0\)の時
\[ P(A|B) = \frac{P(A)P(B|A)}{P(B)} \]
が成立する。
</p>
</div>
<p>
ベイズ主義の立場では, \(P(A)\) を<strong>事前確率</strong>, \(P(A|B)\) を \(B\) が生じたという情報(証拠)を得たことによる<strong>事後確率</strong>と解釈します.
</p>
<p style="font-size:60%">
注: この定理は公理的確率論の定理であるので, 頻度主義であってもこの定理は利用されます.
</p>
</section>
<section>
<p> 以下の形も良く使われます. </p>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> ベイズの定理 </h4>
<p>
\(B_1,B_2,\cdots\)が互いに素で,\(\Omega=B_1\cup B_2\cup\cdots\)であるならば
\[ P(B_i|A) = \frac{P(B_i)P(A|B_i)}{\displaystyle\sum_{i=1}^{\infty}P(B_i)P(A|B_i)} \]
である。
</p>
</div>
</section>
<section>
<p>
ベイズの定理
\[ P(B_i|A) = \frac{P(B_i)P(A|B_i)}{\displaystyle\sum_{i=1}^{\infty}P(B_i)P(A|B_i)} \]
の分母は全確率の値が \(1\) という条件を満たす為の規格化定数と考える事が出来るので,
\[ P(B_i|A) \propto P(B_i)P(A|B_i) \]
と表記する場合もあります.
</p>
</section>
<section>
<p> 有名な問題を考えてみましょう. </p>
<div class="block" style="border-color:lightgreen;font-size:90%">
<h4 style="color:lightgreen"> 3囚人問題 </h4>
<p>
\(A,B,C\)の3人の囚人がいます. 3人のうち2人が処刑される事は分かっており, どの1人が助かるのかは囚人達には全く分かりません.
</p>
<p>
ここで囚人 \(A\) が看守に「\(B,C\)の少なくとも1人は処刑になるはずだが, その名前を教えてくれ」と質問した所, 「\(B\)は処刑される」と看守が答えました. それを聞いて \(A\) は自分が処刑されない確率が \(1/3\) から \(1/2\) に上がったと喜びましたが, これは正しいでしょうか?
</p>
</div>
</section>
<section>
<p>
\(A,B,C\)が処刑を免れるという事象をそれぞれ \(A,B,C\) とし, 「\(B\) が処刑される」と看守が証言する事象を \(E\) とします. \(P(A)=P(B)=P(C)=1/3\) として \(P(A|E)\) について考えます.
</p>
<p class="fragment">
(看守が正直だとすれば) 明らかに\(P(E|B)=0, P(E|C)=1\) となります.
また, \(A\) が処刑を免れる時は「\(B\)が処刑される」と「\(C\)が処刑される」両方の証言が可能です. この場合, 看守がどちらの証言を選ぶかは全く分からないので \(P(E|A)=1/2\) としましょう.
</p>
</section>
<section>
<p>
処刑を免れるのは一人と決まっているので \(A,B,C\)は互いに素で \(A\cup B\cup C = \Omega\) ですからベイズの定理より
\[ P(A|E) = \frac{P(A)P(E|A)}{P(A)P(E|A)+P(B)P(E|B)+P(C)P(E|C)} \]
が成立します.
</p>
<p class="fragment" data-fragment-index="1">
従って
\[ P(A|E) = \frac{\frac{1}{3}\cdot\frac{1}{2}}{\frac{1}{3}\cdot\frac{1}{2}+\frac{1}{3}\cdot 0+\frac{1}{3}\cdot 1} = \frac{1}{3}\]
となります. \(P(A)=P(A|E)\) であるので 証言 \(E\) を聞いた所で \(A\) が処刑を免れる確率は変化しません.
</p>
</section>
<section>
<h2> 理由不十分の原則 </h2>
<p>
今の問題では事前確率を
\[ P(A)=P(B)=P(C) \]
としましたが, これは\(A,B,C\)のどれが起こるのか全く分からないからです. この様に「どの事象が生起するかについて全く情報がないときは等確率とする」という原則を <strong>理由不十分の原則 </strong> と呼びます.
</p>
</section>
<section style="font-size:80%">
<h2> ベイズ改訂と事象の独立性 </h2>
<p>
証拠が沢山ある場合の事後確率の計算は
\[ \begin{aligned}
& P(A|B_1,B_2,\ldots,B_n) \\
\propto &P(A)P(B_1,B_2,\ldots,B_n|A) \\
=&P(A)P(B_1|A)P(B_2,\ldots,B_n|A,B_1) \\
=&P(A)P(B_1|A)P(B_2|A,B_1)P(B_3,\ldots,B_n|A,B_1,B_2)\\
=&\cdots \\
=&P(A)P(B_1|A)P(B_2|A,B_1)P(B_3|A,B_1,B_2)\\
&\qquad \cdots P(B_n|A,B_1,B_2,\ldots,B_{n-1})
\end{aligned} \]
となりますが, この計算は非常に複雑です.
</p>
<p>
ここで, 証拠\(B_i,B_j (i < j)\) が独立ならば
\[ P(B_i|\ldots,B_{j-1},B_j,B_{j+1},\ldots)=P(B_i|\ldots,B_{j-1},B_{j+1},\ldots) \]
が成り立ちます(正確には独立性よりも弱い条件付き独立性という性質を持てば十分).
事象の独立性に注目すると計算を簡略化する事が出来ます.
</p>
</section>
<section style="font-size:80%">
<h2> 単純ベイズモデル </h2>
<p>
<strong> 証拠 \(B_1,B_2,\ldots,B_n\) が全て独立である </strong> という強い仮定を置けば,
\[ P(A|B_1,B_2,\ldots,B_n) \propto P(A)P(B_1|A)P(B_2|A)\cdots P(B_n|A) \]
となるので
\[ \color{yellow}{
P(A|B_1,B_2,\ldots,B_n) \propto P(A)\prod_{i=1}^n P(B_i|A)
}\]
によってベイズ改訂の計算が出来ます. この仮定を適用したモデルを <strong> 単純ベイズモデル(ナイーブベイズモデル) </strong> と呼びます.
</p>
</section>
<section style="font-size:80%">
<div class="block" style="border-color:lightgreen;font-size:80%">
<h4 style="color:lightgreen"> 例題 </h4>
<p>
非常に簡単な迷惑メール判定の問題を考えてみましょう. ここでは4つの単語「連絡」「無料」「登録」「経理」のみに注目します.
これらがメールの中に出現する確率は以下のようになっているとします.
\[ \begin{array}{|c|c|c|c|} \hline
& \text{連絡} & \text{無料} & \text{登録} & \text{経理}\\ \hline
\text{通常メール} & 0.8 & 0.2 & 0.3 & \text{0.9} \\ \hline
\text{迷惑メール} & 0.2 & 0.8 & 0.7 & \text{0.1} \\ \hline
\end{array} \]
また, 受け取るメールの6割が通常メールであるとします.
</p>
<p>
あるメールを調べた所「連絡」「無料」「登録」という単語が検出されたとします. このメールが通常メール・迷惑メールのどっちであるか単純ベイズモデルを適用して計算してみましょう.
</p>
</div>
<p class="fragment">
\[ \begin{aligned}
P(\text{通常}|\text{連絡},\text{無料},\text{登録}) &\propto P(\text{通常}) P(\text{連絡}|\text{通常}) P(\text{無料}|\text{通常}) P(\text{登録}|\text{通常})\\
&= 0.6\times 0.8 \times 0.2\times 0.3 = 0.0288 \\
P(\text{迷惑}|\text{連絡},\text{無料},\text{登録}) &\propto P(\text{迷惑}) P(\text{連絡}|\text{迷惑}) P(\text{無料}|\text{迷惑}) P(\text{登録}|\text{迷惑}) \\
&= 0.4\times 0.2 \times 0.8 \times 0.7 = 0.0448
\end{aligned} \]
となるので, この場合は迷惑メールである可能性が高いと考えられます.
</p>
</section>
<section>
<p>
今の問題では, 文章中の単語が全く無関係にランダムに出現するという, 納得しがたい仮定をおいています.
</p>
しかし実際には,単純ベイズモデルに基づく分類器は多くの問題に対して効果的に働くようです.
この事については, 事象の従属関係による確率の変動が実際には上手くキャンセルされるからだという<a href="http://www.cs.unb.ca/profs/hzhang/publications/FLAIRS04ZhangH.pdf"> 理論的な説明 </a> もあります.
</p>
</section>
<section>
<p>
計算上の注意ですが \(P(A)\prod_{i=1}^n P(B_i|A) \) の値は \(n\) が大きいとどんどん小さくなりアンダーフローがおきてしまいます.
そこで実際には,
\[ P(A|B_1,B_2,\ldots,B_n) =\alpha P(A)\prod_{i=1}^n P(B_i|A) \quad (\alpha:\text{規格化定数}) \]
の対数をとって
\[ \ln P(A|B_1,B_2,\ldots,B_n) =\ln \alpha + \ln P(A) + \sum_{i=1}^n \ln P(B_i|A) \]
の形で計算を行います.
</p>
</section>
<section>
<h2> ベイジアンネットワーク </h2>
<p>
今回は詳しい解説をする時間がないのですが, 事象の独立性を<strong>ベイジアンネットワーク</strong>というグラフ構造の上で表す事が出来ます.
グラフ上の幾つかのノードに対して証拠が与えられた時に,それ以外のノードに関する全確率をベイズの定理に基いて計算する事が出来ます.
</p>
<div align="center">
<img width="500px" src="fig/bayesian-network-example.png">
</div>
</section>
<section>
<h2> ベイズ推定 </h2>
<p>
確率分布の母数\(\theta\)に対する推定を行うには以下の定理を利用します.
</p>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> 母数分布に対するベイズの定理 </h4>
<p>
\(X\) が確率密度関数(もしくは確率質量関数) \(f(X|\theta)\) の表す分布に従い, 母数 \(\theta\) の<strong> 事前分布 </strong> が \(p(\theta)\) ならば
\[ p(\theta|X=x) = \frac{f(X=x|\theta)p(\theta)}{f(X=x)} \]
を\(\theta\)の<strong>事後分布</strong>という.
</p>
</div>
<p>
\(p(\theta|X=x)\) を最大にする\(\theta\)の値を <strong> ベイズ推定値 </strong> (もしくはMAP(maximum a posterior estimator)推定値)と呼びます.
</p>
</section>
<section style="font-size:80%">
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen"> 例題 </h4>
<p>
ある新薬を10人の治験者に使用した所, 7人には効果があり3人には効果がなかった.
この薬が効く確率 \(\theta\) は全く分からないとして事後分布とベイズ推定値を計算しましょう.
</p>
</div>
<p>
今回得られた結果を \(A\) としましょう. これは二項分布に従うので
\[ f(A|\theta) = \binom{10}{7}\theta^7(1-\theta)^3 \]
となります. \(\theta\) の値は全く分からないので事前分布は一様分布とします.
\[ p(\theta) = 1\qquad (0\leq \theta \leq 1) \]
すると, 事後分布は
\[ p(\theta|A) \propto \binom{10}{7}\theta^7(1-\theta)^3\times 1 \propto \theta^7(1-\theta)^3 \]
となります.
</p>
</section>
<section style="font-size:80%">
<p>
ここで
\[ \int_0^1 \theta^7(1-\theta)^3\mathrm{d}\theta = \frac{1}{1320} \]
となるので規格化を行って
\[ \color{yellow}{ p(\theta|A) = 1320\theta^7(1-\theta)^3} \]
が求める事後分布となります.(後述しますが, これはベータ分布 \(Be(8,4)\) という分布です.)
</p>
</section>
<section>
<p>
図示してみれば, 以下が \(\theta\) の事前分布です. \(\theta\) の値は全く分からないという様子です.
</p>
<div align="center">
<img width="600px" src="fig/ex18-1-2.png">
</div>
</section>
<section>
<p>
事後分布は以下のようになります. \(\theta\) の値は0.7付近の値である確率が高いという事になります.
</p>
<div align="center">
<img width="600px" src="fig/ex18-1-1.png">
</div>
</section>
<section>
<p>
ベイズ推定値は,
\[ \frac{\mathrm{d}}{\mathrm{d}\theta}\theta^7(1-\theta)^3 = -\theta^7(\theta-1)^2(10\theta-7) \]
より \(\color{yellow}{ \hat{\theta} = 7/10 }\)となります. これは10人中7人に効いたという結果と合致します.
</p>
<p class="fragment">
これは当たり前の事で, 事前分布が一様分布ならば
\[ p(\theta|X=x) \propto f(X=x|\theta) \]
となりますので事後分布は尤度関数に比例します. つまり, <strong> 事前分布が一様分布の時はベイズ推定値と最尤推定値は一致</strong>します.
</p>
</section>
<section>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> ベータ分布 </h4>
確率密度関数が
\[ f(x) = Cx^{m-1}(1-x)^{n-1} \qquad (0\leq x\leq 1,C:\text{規格化定数})\]
と表される確率分布を<strong>ベータ分布</strong>と呼び \(Be(m,n)\) と書く.
<div align="center">
<img width="450px" src="fig/beta-distribution.png"><br>
</div>
</section>
<section>
<p>
\(Be(m,n)\) の平均・分散は以下の様になります.
\[ \text{平均}: \frac{m}{m+n}\quad\text{分散}: \frac{mn}{(m+n)^2(m+n+1)} \]
また, 密度関数が最大値を取る点(最頻値)は
\[ \frac{m-1}{m+n-2} \]
となります. これはベイズ推定値を計算する為に便利です.
</p>
<p>
また, 標準一様分布はベータ分布 \(Be(1,1)\) と一致します.
</p>
</section>
<section>
<p>
最尤推定値とベイズ推定値の違いを見る為に以下の問題を考えましょう.
</p>
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen"> 例題 </h4>
<p>
あるコインを5回投げたら全て表が出た. このコインの表が出る確率 \(p\) を推定しましょう.
</p>
</div>
</section>
<section>
<p>
例題の事象を \(A\) と書きます. 尤度関数は
\[ f(A|p) = p^5 \qquad (0\leq p \leq 1)\]
となるので, <strong> \(p\) の最尤推定値は \(1\)</strong> となります. つまり, 常に表が出るコインであるという事になります.
</p>
</section>
<section>
<p>
一方, ベイズ推定では事前分布として私達の常識的な知識を使う事が出来ます. \(p\) の値は \(p=1/2\) の付近にあるはずなので, 事前分布\(\pi(p)\) として \(Be(2,2)\) を仮定しましょう.
</p>
<div align="center">
<img width="600px" src="fig/ex18-2.png"><br>
</div>
</section>
<section>
<p>
すると, 事後分布 \(\pi(p|A)\) は
\[ \pi(p|A)\propto p(1-p)\times p^5 = p^6(1-p) \]
より \(B(7,2)\) となるので, <strong> \(p\) のベイズ推定量は 6/7 </strong> となります.
</p>
<p class="fragment">
この例で最尤推定法が上手くいかなかったのは, サンプル数が少なかったからだと言えます. ベイズ推定法ではサンプルの不足を事前知識によって補って意味のある結果を得る事が出来ます.
</p>
</section>
<section>
<p> 他の種類の分布の例を見てみましょう. </p>
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen"> 例題 </h4>
<p>
あるパン屋が新しく売りだしたメロンパンは \(100\mathrm{g}\) の重さらしいです. 過去の調査からこのパン屋が作るパンの重さの分散は 10 である事が分かっています.
</p>
<p>
試しに 4つのメロンパンを買って調べた所 \(101\mathrm{g},105\mathrm{g},98\mathrm{g},102\mathrm{g}\) でした. メロンパンの重さの母平均 \(\mu\) を推定しましょう.
</p>
</div>
</section>
<section style="font-size:80%">
<p>
例題の事象を \(A\) とします. パンの重さは正規分布に従うと仮定すると, 尤度関数は
\[\begin{aligned}
f(\mu|A) &\propto
e^{-\frac{(101-\mu)^2}{20}}
e^{-\frac{(105-\mu)^2}{20}}
e^{-\frac{(98-\mu)^2}{20}}
e^{-\frac{(102-\mu)^2}{20}} \\
&\propto \exp\left\{ -\frac{\mu^2-203\mu}{5}\right\}
\end{aligned} \]
となります.
</p>
<p>
\(\mu\) の事前分布 \(p(\mu)\) ですが平均が\(100\)の正規分布としましょう. 分散は大きめに100とします. つまり
\[ p(\mu)\propto \exp\left\{-\frac{(\mu-100)^2}{200}\right\} \]
となります.
</p>
</section>
<section style="font-size:80%">
<p>
よって, \(\mu\) の事後分布は
\[ \begin{aligned}
p(\mu|A) &\propto
\exp\left\{ -\frac{\mu^2-203\mu}{5}\right\}
\exp\left\{-\frac{(\mu-100)^2}{200}\right\} \\
&\propto \exp\left\{-\frac{41\mu^2-8320\mu}{200}\right\} \\
&\propto \exp\left\{-\frac{41(\mu-4160/41)^2}{200}\right\}
\end{aligned} \]
より 平均 \(4160/41\), 分散 \(100/41\) の正規分布となります. よって \(\mu\) のベイズ推定値は \(4160/41=101.5\) となります.
</p>
</section>
<section>
<p>
\(\mu\) の事前分布・事後分布は以下の様になります(それぞれピンクと青). 調査によって分散が減少し \(\mu\) の値の不確実性が減少した事がわかります.
</p>
<div align="center">
<img width="600px" src="fig/ex18-3.png"><br>
</div>
</section>
<section>
<h2> 自然な共役分布 </h2>
<p>
事前分布の取り方によっては事後分布は非常に複雑になってしまいます. すると積分や最大値の計算が面倒になります.
</p>
<p>
<strong> 事前分布と事後分布の種類が一致する</strong> 様に分布を選べばこの問題を回避する事が出来ます. そのような条件を満たす分布を(母集団に対する)<strong> 自然な共役分布 </strong> と呼びます.
</p>
</section>
<section style="font-size:80%">
<p>
これまでの例題で利用したのが
</p>
\[ \begin{array}{|c|c|c|} \hline
\text{母集団の分布} & \text{事前分布} & \text{事後分布} \\ \hline
\text{ベルヌーイ分布・二項分布} & \text{ベータ分布} & \text{ベータ分布} \\ \hline
\text{正規分布} & \text{正規分布} & \text{正規分布} \\ \hline
\end{array} \]
<p>
という関係で, 他には
</p>
\[ \begin{array}{|c|c|c|} \hline
\text{母集団の分布} & \text{事前分布} & \text{事後分布} \\ \hline
\text{ポアソン分布} & \text{ガンマ分布} & \text{ガンマ分布} \\ \hline
\text{正規分布} & \text{逆ガンマ分布} & \text{逆ガンマ分布} \\ \hline
\end{array} \]
<p>
などがあります. 逆ガンマ分布は母分散の推定に使います. 詳しくは省略します.
</p>
</section>
<section>
<p>
公式化してしまいましょう. 証明は練習問題とします.
</p>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> ベイズ改訂の公式:ベータ分布 </h4>
<p>
母集団分布が試行回数 \(n\), パラメータ \(p\) のベルヌーイ分布もしくは二項分布であるとする. 試行の成功回数が \(r\) であるとき, 事前分布に \(Be(p,q)\) を仮定してベイズ改訂を行うと事後分布は
\[ Be(p+r, q+n-r) \]
となる.
</p>
</div>
</section>
<section>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> ベイズ改訂の公式:正規分布 </h4>
<p>
\(N(\mu,\sigma^2)\)に従う独立な確率変数の実現値\(x_1,x_2,\ldots,x_n\)を用いて, 事前分布を \(N(\mu_0,\sigma_0^2)\) としてベイズ改訂を行うと事後分布は
\[ N\left(\frac{\sum x_i/\sigma^2+\mu_0/\sigma_0^2}{n/\sigma^2+1/\sigma_0^2}, \frac{1}{n/\sigma^2+1/\sigma_0^2} \right) \]
となる.
</p>
</div>
</section>
<section>
<p>
自然な共役分布ではなく, 他の種類の確率分布を利用したい場合もあるかもしれません. この場合, 解析的に事後分布を計算する事は難しいので数値計算が必要になります.
</p>
<p>
詳しくは省略しますが, <strong>マルコフ連鎖モンテカルロ法(MCMC法)</strong> という手法が近似計算に利用されます.
</p>
</section>
<section>
<h2> ベイズ線型回帰 </h2>
<p>
線形回帰の問題に事前分布の考え方を導入する事を考えます.
</p>
<p>
\[ y = a_1g_1(\mathbf{x})+\cdots + a_mg_m(\mathbf{x}) \]
という回帰方程式で重回帰を行う事を考えます. 誤差分布の分散を\(\sigma^2\) とします. データセット \((\mathbf{x}_1,y_1),(\mathbf{x}_2,y_2),\ldots,(\mathbf{x}_n,y_n)\) に対して行列 \(G\), ベクトル \(\mathbf{y}\) を前回設定したようにおきます. また, \(\mathbf{a}\) を係数のベクトルとします.
</p>
<p>
まず, このデータセット\(D\)に対する尤度は
\[ f(D|\mathbf{a}) \propto \exp\left\{-\frac{1}{2\sigma^2}(\mathbf{y}-G\mathbf{a})^T(\mathbf{y}-G\mathbf{a})\right\} \]
</p>
</section>
<section style="font-size:80%">
<p>
また, 事前分布として \(\mathbf{a}\) が多変量正規分布 \(N(\mathbf{\mu}_0,\Sigma_0)\) に従う事を仮定します. すなわち,
\[ p(\mathbf{a})\propto \exp\left\{-\frac{1}{2}(\mathbf{a}-\mathbf{\mu}_0)^T\Sigma_0^{-1}(\mathbf{a}-\mathbf{\mu}_0) \right\} \]
とします.
</p>
<p>
すると, 事後分布は
\[ p(\mathbf{a}|D) \propto \exp\left\{-\frac{1}{2}\left(\frac{1}{\sigma^2}(\mathbf{y}-G\mathbf{a})^T(\mathbf{y}-G\mathbf{a})+(\mathbf{a}-\mathbf{\mu}_0)^T\Sigma_0^{-1}(\mathbf{a}-\mathbf{\mu}_0)\right) \right\} \]
となります. これを整理すると 事後分布は多変量正規分布 \(N(\mathbf{\mu},\Sigma)\) に従う事が判ります. 従って \(\mathbf{\mu}\) が \(\mathbf{a}\) のベイズ推定値となります. 但し
\[ \color{yellow}{ \Sigma^{-1} = \frac{1}{\sigma^2}G^TG + \Sigma_0^{-1}, \mathbf{\mu} = \Sigma\left(\frac{1}{\sigma^2}G^T\mathbf{y}+\Sigma_0^{-1}\mathbf{\mu}_0\right) } \]
です.
</p>
</section>
<section style="font-size:70%">
<p>
【証明】<br>
分散共分散行列が対称行列 (\(\Sigma^T=\Sigma\)) である事に注意すると
\[ \begin{aligned}
&(\mathbf{a}-\mathbf{\mu})^T\Sigma^{-1}(\mathbf{a}-\mathbf{\mu}) \\
= &\mathbf{a}^T\Sigma^{-1}\mathbf{a} -2\mathbf{a}^T\Sigma^{-1}\mathbf{\mu}+(\text{定数})\\
= &\mathbf{a}^T\left(\frac{1}{\sigma^2}G^TG+\Sigma_0^{-1}\right)\mathbf{a} - 2\mathbf{a}^T\left(\frac{1}{\sigma^2}G^T\mathbf{y}+\Sigma_0^{-1}\mathbf{\mu}_0\right) + (\text{定数}) \\
=& \frac{1}{\sigma^2}(\mathbf{a}^TG^TG\mathbf{a}-2\mathbf{a}^TG^T\mathbf{y})+(\mathbf{a}^T\Sigma_0^{-1}\mathbf{a}-2\mathbf{a}^T\Sigma_0^{-1}\mathbf{\mu}_0)+(\text{定数}) \\
=& \frac{1}{\sigma^2}(\mathbf{y}-G\mathbf{a})^T(\mathbf{y}-G\mathbf{a}) + (\mathbf{a}-\mathbf{\mu}_0)^T\Sigma_0^{-1}(\mathbf{a}-\mathbf{\mu}_0) + (\text{定数})
\end{aligned} \]
となるから. <span style="float:right">□</span>
</p>
</section>
<section>
<p>
先ほどの式から \(\Sigma\) を消去するとベイズ線型回帰の正規方程式
\[ G^TG\mathbf{\mu}-G^T\mathbf{y}+\sigma^2\Sigma_0^{-1}(\mathbf{\mu}-\mathbf{\mu}_0) = 0 \]
が得られます.
</p>
<p>
もしパラメータ \(\mathbf{a}\) について一切の情報がないならば, 分散が無限大となって(\(\Sigma_0^{-1}=0\) となるから)最小二乗法の正規方程式
\[ G^TG\mathbf{\mu}-G^T\mathbf{y}=0 \]
と一致する事が判ります.
</p>
</section>
<section style="font-size:80%">
<div class="block" style="border-color:lightgreen;font-size:80%">
<h4 style="color:lightgreen"> 例題 </h4>
<p> 以下のデータセットに対して \(9\)次方程式を当てはめてみましょう.
\[ \begin{array}{|c|cccccccccc|} \hline
X & -1.2&1.9&-3.4&-0.4&2.9&2.9&3.5&-0.6&2.0&3.7 \\ \hline
Y & 4.6& -0.6&22.3&-0.0&-0.0&5.0&3.1&5.3&1.9&8.4 \\ \hline
\end{array} \]
</p>
</div>
</section>
<section style="font-size:80%">
<p>
単純な最小二乗法では
\[\small{ y=-0.01x^9+0.05x^8+0.11x^7-0.39x^6-1.0x^5-2.7x^4+16x^3+7.8x^2-32x-11 } \]
となり以下の様にノイズを忠実に拾った回帰方程式となってしまうのでした.
</p>
<div align="center"> <img width="500px" src="fig/ex18-4-1.png"> </div>
</section>
<section style="font-size:70%">
<p>
ベイズ線型回帰でやってみます. ここで次数の高い項の係数ほど \(0\) に近いはずという事を事前分布として導入してみましょう.
</p>
<p>
そこで, \(x^k\) の係数は \(N(0, 10^{-k})\) に従うと仮定します. つまり,
\[ \mathbf{\mu}_0 = \mathbf{0},
\Sigma_0 = \begin{pmatrix}
1 & 0 & \cdots & 0 \\
0 & 10^{-1} & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & \cdots & \cdots & 10^{-9}
\end{pmatrix} \]
とおきましょう. つまり
\[
\Sigma_0^{-1} = \begin{pmatrix}
1 & 0 & \cdots & 0 \\
0 & 10 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & \cdots & \cdots & 10^9
\end{pmatrix} \]
となります. 誤差分布の分散の \(\sigma^2\) の値は分かりませんので適当な値を仮定します.
</p>
</section>
<section>
<p>
正規方程式を解くと
\[\begin{aligned}
&\small{ y=-2.3\cdot 10^{-5}x^9+2.4\cdot 10^{-4}x^8 -3.2\cdot 10^{-4}x^7 + 1.9\cdot 10^{-3}x^6}\\
&\quad\small{ - 3.5\cdot 10^{-3}x^5 + 1.5\cdot 10^{-2}x^4 -3.9\cdot 10^{-2}x^3 +0.13x^2 -0.67x + 2.1}
\end{aligned} \]
となり以下のグラフ(緑の線)のようになります. 高次の急激な変動を除去する事ができている事がわかります.
</p>
<div align="center"> <img width="500px" src="fig/ex18-4-2.png"> </div>
</section>
<section>
<p>
以上でベイズ統計の基本の説明を終わります.
</p>
<p>
事前分布を主観的に定めてしまうという考え方は慣れるまでは(慣れても?)奇妙に感じるかもしれませんが,
データ数が少ない場合の統計的解析には大変強力です. 是非, 復習を行って下さい.
</p>
</section>
<section>
<h2> 以上でプログラマの為の数学勉強会を終わります. </h2>
<p>
全18回お付き合いいただきありがとうございました.
非常に初歩的な所から始めましたが, 結構実用的なレベルの所まで解説する事が出来たと思います.
</p>
<p>
今後ですが, せっかく統計まで講義を行いましたので最近はやりのパターン認識・機械学習の勉強会を開催したいと思います.
そちらではこの勉強会で解説した程度の基礎的な数学知識は前提とし, 積極的に難しいテーマにも取り組んでいこうと思います.
是非お越しください.
</p>
</section>
</div>
</div>
<script src="lib/reveal/lib/js/head.min.js"></script>
<script src="lib/reveal/js/reveal.js"></script>
<script>
Reveal.initialize({
width: 960,
height: 640,
controls: true,
progress: false,
history: true,
overview: false,
touch: true,
center: false,
rollingLinks: false,
transition: "page",
transitionSpeed: "default",
// When scale != 1, positions of mouse events will be incorrect.
minScale: 1.0,
maxScale: 1.0,
dependencies: [
{ src: "lib/reveal/lib/js/classList.js", condition: function() { return !document.body.classList; } },
{ src: "lib/reveal/plugin/highlight/highlight.js", async: true, callback: function() { hljs.initHighlightingOnLoad(); } },
{ src: "lib/reveal/plugin/zoom-js/zoom.js", async: true, condition: function() { return !!document.body.classList; } },
{ src: "lib/reveal/plugin/notes/notes.js", async: true, condition: function() { return !!document.body.classList; } }
]
});
// register event listeners
require(["MathJax"], function (MathJax) {
// Delay typesetting of slides
function typeset (idx) {
for (var i = idx - 2, n = idx + 2; i <= n; i++) {
var slide = Reveal.getSlide(i);
if (!slide) continue;
if (!slide.typeset) {
MathJax.Hub.Typeset(slide);
slide.typeset = true;
}
}
}
function initializeGraphics (idx) {
for (var i = idx - 2, n = idx + 2; i <= n; i++) {
var slide = Reveal.getSlide(i);
if (!slide) continue;
var graphics = slide.getAttribute("graphics");
if (graphics && !slide.initialized) {
slide.initialized = true;
(function () {
var p = slide;
require([graphics], function(g) {
if (g.initialize) g.initialize(p);
});
})();
}
}
}
function start (slide) {
var graphics = slide.getAttribute("graphics");
if (graphics) {
require([graphics], function(g) { if (g.start) g.start(slide); });
}
}
function stop (slide) {
var graphics = slide.getAttribute("graphics");
if (graphics) {
require([graphics], function(g) { if (g.stop) g.stop(slide); });
}
}
function simpleEvent (type) {
var event = document.createEvent("HTMLEvents");
event.initEvent(type, true, true);
return event;
}
Reveal.addEventListener("slidechanged", function (event) {
typeset(event.indexh);
initializeGraphics(event.indexh);
start(event.currentSlide);
stop(event.previousSlide);
});
Reveal.addEventListener("fragmentshown", function (event) {
var slide = Reveal.getCurrentSlide();
var graphics = slide.getAttribute("graphics");
if (graphics) {
require([graphics], function(g) { if (g.proceed) g.proceed(slide); });
}
});
console.log(Reveal.getIndices().h);
typeset(Reveal.getIndices().h);
initializeGraphics(Reveal.getIndices().h);
start(Reveal.getCurrentSlide());
});
</script>
</body>
</htl>