forked from winterwang/LSHTMlearningnote
-
Notifications
You must be signed in to change notification settings - Fork 0
/
08-Intro-to-Bayes.Rmd
1258 lines (887 loc) · 64.9 KB
/
08-Intro-to-Bayes.Rmd
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
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# (PART) 貝葉斯統計 Introduction to Bayesian Statistics {-}
# 貝葉斯統計入門 {#intro-Bayes}
> A Bayesian statistician is one who, vaguely expecting a horse and catching a glimpse of a donkey, strongly concludes he has seen a mule.
> ~ Guernsey McPearson's Drug Development Dictionary^[http://www.senns.demon.co.uk/wdict.html]
```{block2, note-thankJavi, type='rmdnote'}
The Introduction to Bayesian Statistics lectures were orgainised and taught by Dr. [F. Javier Rubio](https://sites.google.com/site/fjavierrubio67/) and Professor [Karla DiazOrdaz](https://scholar.google.com/citations?user=GSGcT-EAAAAJ&hl=en).
```
本章節之目的：
- 介紹 (啓發) 貝葉斯推斷 Bayesian inference 的基本概念。並且與概率論 frequentist inference 推斷實例作比較。
- 介紹共軛分佈的概念 conjugate distributions。用單一參數家族 (single parameter family) ，特別是二項分佈的圖形來描述共軛分佈；用方差已知的正態分佈均值來描述共軛分佈。
- 介紹貝葉斯預測分佈 Bayesian prediction distribution。
推薦書目：
1. ["Principles of Statistical Inference"](https://books.google.co.jp/books?id=nRgtGZXi2KkC&dq=principles+of+statistical+inference&lr=&source=gbs_navlinks_s) by D.R. Cox [@cox2006principles]
2. ["Bayesian Data Analysis"](https://books.google.co.jp/books?id=nRgtGZXi2KkC&dq=principles+of+statistical+inference&lr=&source=gbs_navlinks_s) by Gelman, Carlin, Stern, Dunson, Vehtari, and Rubin [@gelman2013bayesian], [website for the book](www.stat.columbia.edu/~gelman/book/)
3. ["Bayesian Biostatistics"](https://books.google.co.uk/books?id=WV7KVjEQnJMC&printsec=frontcover&dq=Bayesian+biostatistics&hl=zh-CN&sa=X&ved=0ahUKEwjct_3vqcbXAhXFJ8AKHar3BbQQ6AEIJjAA#v=onepage&q=Bayesian biostatistics&f=false) by Vehtari and Rubin [@lesaffre2012bayesian]
貝葉斯統計推斷，提供了不同於概率論推斷的另一種考察和解決問題的思路。所有的思考，都源於貝葉斯定理 Bayes' Theorem (Section \@ref(Bayes-Definition))。起源於英國統計學家[托馬斯貝葉斯 (Thomas Bayes)](https://en.wikipedia.org/wiki/Thomas_Bayes) 死後被好友 [Richard Price](https://en.wikipedia.org/wiki/Richard_Price) 整理發表的論文: ["An essay towards solving a problem in the doctrine of chances."](www.stat.ucla.edu/history/essay.pdf)
概率論推斷與貝葉斯推斷的中心都圍繞似然 likelihood (Section \@ref(likelihood-definition)) 的概念。然而二者對似然提供的信息之理解和解釋完全不同。即在對於觀察數據提供的信息的理解，和如何應用已有信息來影響未來決策（或提供預測）的問題上常常被認爲是統計學中形成鮮明對比的兩種哲學理念。過去幾個世紀二者之間孰優孰劣的爭論相當激烈。但是，從實際應用的角度來看，我們目前更關心哪種思維能更加實用地描述和模擬真實世界。幸運地是，多數情況下，二者的差距不大。所以無法簡單地從一個實驗或者一次爭論中得出誰更出色的結論。現在的統計學家們通常不再如同信仰之爭那樣的互相水火不容，而是從實用性角度來判斷一些實際情況下，採用哪種思想能使計算過程更加簡便或者計算結果更加接近真實情況。
請思考如下的問題：
什麼是概率？ What is probability?
1. 概率論思想下的定義：某事件在**多次重複觀察**實驗結果中發生次數所佔的比例。<br> The probability of an event is the limit of its relative frequency in a large number of trials."
2. 貝葉斯思想下的定義：概率是你相信某事件會發生的可能性。 <br> Probability is a measure of the degree of belief
about an event.
## 概率論推斷的複習
思考不同場景：
- 場景 A：假如我們在監測一個製造鐵絲的工廠，需要測量該工廠生產的鐵絲的強度。
- 場景 B：假如我們正在進行一個大型隊列研究，該研究是關於心臟病和與之相關的某個危險因子的評價。數據來源是家庭醫生的診療數據庫。
- 場景 C：假如一名警察凌晨三點在空無一人的街頭巡邏時，突然聽見防盜自動警鈴的報警聲。他立刻循聲望去，對面街上的珠寶店玻璃碎了一地。一個戴着巴拉克拉瓦頭套的人正揹着一個大包從破碎玻璃窗中爬出。該警察毫不猶豫地判定該人就是劫匪，立刻將其逮捕。
在這些場景下，請用概率論思想思考如下幾個問題：
1. 事件是什麼？
2. 如何解讀總體參數？
3. 如何使用參數進行概率推斷？
4. 用經典概率論時，有什麼缺點嗎？
- 場景 A：
1. 事件：該工廠製造的鐵絲，長期以來的強度大小是多少。
2. 總體參數：鐵絲的真實強度，或者與鐵絲強度相關的特性。
3. 概率推斷：我們進行鐵絲的強度實驗，即從該工廠已經生產的鐵絲中大量抽取樣本逐一進行強度檢測。用相應的概率模型來模擬抽取的樣本數據，並且使用極大似然估計找到最能體現抽樣數據的參數估計，然後對獲得的極大似然估計進行95%信賴區間的計算。然後如果我們重複這樣相同的實驗無數次，那麼我們計算的所有的信賴區間中，有95%包含了真實的鐵絲強度大小。
4. 在鐵絲強度測量的場景中，經典概率論顯得十分自然，因爲我們真的可以重複這樣的實驗很多很多次以獲得想要的參數的精確估計。
- 場景 B：
1. 事件：由於我們用的是整個隊列研究的數據。所以從概率論的角度來看，本事件就是假定我們可以在人數無限多的人羣中重複同樣的隊列研究。
2. 總體參數：我們感興趣的心臟病相關危險因子，在抽取該隊列作爲樣本的人羣中的真實值大小。
3. 概率推斷：我們用泊松分佈的概率模型來模擬人羣中從開始觀察時起，至心臟病發病這段時間內和該危險因子之間的關係大小。然後用傳統的極大似然估計法計算獲得 HR, OR 等值來表示危險因素和心臟病的關係。
4. 缺點：實際情況是，經費時間和人力資源的限制下，我們無法“重複相同的隊列研究”。而且該對列本身可能就是十分獨特的，比如只有男性，或者有年齡限制，或者其他的特性使隊列本身在理論上就是不可能被重複的。所以，在這樣的場景下，用經典的概率論思想作統計推斷常常會被認爲是不自然不妥當的。
- 場景 C：
1. 事件：警察無數次在同樣的時間同樣的地點巡邏時，聽見防盜自動警鈴的報警聲，他看見頭戴巴拉克拉瓦頭套的人從破碎的玻璃窗中爬出......
2. 總體參數：在無數次上面描述的場景時，發生盜竊案的真實概率。
3. 概率推斷：使用某種可以描述該事件（巡邏時。。。發生盜竊案的概率）的數學模型，我們用極大似然估計來計算發生盜竊案概率的估計和95%信賴區間，**然後警察同志再來決定是否要去抓眼前這個頭戴巴拉克拉瓦頭套的人**。
4. 缺點：經典概率論在如此場景下很明顯是完全不適用的。<br>1) 這裏經典概率論思維下的概率實際上無法準確定義，充其量是一種發生盜竊案可能性的估計。<br>2) 在如此場景下，警察會根據已經觀察到的現象（已知信息），來判斷一場盜竊案發生的概率是多少。
通過上面不同場景的下的思考，應該能看到傳統概率論中始終假設我們可以**重複相同的實驗**多次，然後從長遠來估測相關事件發生的概率。許多場景下，即使事件概率能被準確定義，我們是很難知道我們關心的參數的分佈的，從而導致我們常常要用到漸進法估算 (asymptotic approximation)。
## 貝葉斯概率推理/逆概率 Bayesian reasoning/inverse probability
首先，不得不承認的一個事實是，**所有的概率都是條件概率**。
- 要麼是根據已知的信息。
- 要麼是一般性大家都接受的某種假設條件。
其次，概率，並不是“長遠”地重複觀察獲得的事件發生頻率。相反地，概率的大小取決與**你自己**和**你感興趣的話題（事件）**。思考下列例子：
1. 明天會下雨嗎？
2. 阿森納下一場比賽會贏還是會輸？
3. 你的期末考試能不能過？
### 演繹推理 deductive reasoning 和 三段論 weak syllogisms
數學要用到邏輯，假設我們用 $A,B,C$ 標記不同的事件。
- 如果 $A\Rightarrow B$ (事件 A 可以推導出事件 B)
- 那麼當我們知道“事件 B 爲真”時，雖然B不一定能倒推回 A，但是我們會相信**事件 A 很可能發生了**。
例如，A 表示“正在下雨”這件事，B 表示 “天上有烏雲”。那麼從邏輯學上來說，$A\Rightarrow B$ 。然而有烏雲本身不一定會下雨。但是會讓我們覺得下雨的可能性增加了。
再來思考警察巡邏的例子。A 表示 “在珠寶店正在發生盜竊案”，B 表示 “一個頭戴巴拉克拉瓦頭套的人正在從玻璃窗中爬出”。也是一樣的道理。
所以警察薯熟在做判斷的時候，需要判斷 Pr(A|B)。他需要如下的信息：
1. 珠寶店發生盜竊案的前提下，有個人從碎玻璃窗中爬出來的概率。
2. 該警察薯熟正處於的環境（半夜三點無人的街頭，等場景）
所以，看到這裏是不是覺得貝葉斯使用的是我們的“常識”在思考決斷問題？因爲我們的先驗概率 (prior) 至關重要。這是我們的背景知識和解釋參數似然（推斷）的依據。
### 如何給可能性定量 Quantifying plausibility
進行可能性定量之前，R.T. Cox 制定了如下的規則[@Cox1946]：
1. $\text{plausibility}(A)$ 是一個有邊界的實數；
2. 傳遞性，transitivity：如果
- $\text{plaus}(C)>\text{plaus}(B)$ and
- $\text{plaus}(B)>\text{plaus}(A)$ then
- $\text{plaus}(C)>\text{plaus}(A)$
3. 一致性，consistency：事件 $A$ 發生的可能性只取決於所有與 $A$ 直接相關的信息，而不包括那些推理到與 $A$ 相關信息之前的信息。<br> The plausibility of proposition $A$ depends only on the relevant information on $A$ and not on the path of reasoning followed to arrive at $A$.
R.T. Cox 證明了他提出的這些規則可以完全適用於所有的可能性計算，而且可能性 (plausibility) 的這些規則和概率 (probability) 的微積分計算完全一致。
所以利用上面的可能性規則，我們可以對條件概率進行更深層次的定義：
$$\text{Pr}(A|B)=\frac{\text{Pr}(B|A)\text{Pr}(A)}{\text{Pr}(B)}\propto \text{Pr}(B|A)\text{Pr}(A)$$
用文字表述爲：
<center>
事後概率 $\propto$ 似然 $\times$ 先驗概率
</center>
其中：
- **事後概率，posterior probability**：$B$ 發生的條件下, $A$ 發生的概率；
- $\propto$ ：與...成正比；
- **似然，likelihood**：$A$ 發生的條件下，$B$ 發生的概率；
- **先驗概率，prior probability**：事件 $A$ 發生的概率。
這就是**貝葉斯定理**。這個定理也告訴我們爲什麼貝葉斯論證在18，19世紀時被叫做“逆概率推理, inverse probability reasoning”。因爲似然 ($A$ 發生的條件下，$B$ 發生的概率) 在與先驗概率相乘以後，概率發生了逆轉--事後概率 ($B$ 發生的條件下, $A$ 發生的概率)。
回頭再來看之前的珠寶店盜竊案：
- 事件 $A$：珠寶店正在發生盜竊案；
- 事件 $B$：一個頭戴巴拉克拉瓦頭套的人正在從玻璃窗中爬出。
所以：
- $\text{Pr}(A)=$ 珠寶店發生盜竊案的概率 -- 先驗概率 (prior probability);
- $\text{Pr}(B|A)=$ 當珠寶店發生盜竊案時，觀察到“一個頭戴巴拉克拉瓦頭套的人正在從玻璃窗中爬出”事件的可能性 -- 似然 (likelihood);
- $\text{Pr}(A|B)$ 當觀察到“一個頭戴巴拉克拉瓦頭套的人正在從玻璃窗中爬出”事件時，倒推珠寶店發生了盜竊案的概率 -- 事後概率 (posterior probability)。
用例子來解釋貝葉斯推理之後你會發現，其實貝葉斯思想也是純粹的概率理論。與經典概率論不同的是，**我們沒有必要認爲某些事件發生的概率需要被重複實驗驗證**。貝葉斯對整個世界的理解源於我們每個人自己認爲的事件發生概率 (personalisitic probability)，或者叫信念度（degree of belief）。
## 貝葉斯推理的統計學實現
在經典概率論中，概率分佈的標記 $f_X(x;\theta)$ 的涵義爲：
對於一個隨機變量 $X$，它在我們假設的某種固定的真實（上帝才知道是多少的）參數 $\theta$ 的分佈框架下，不斷重複相同的實驗之後獲得的概率分佈。
在貝葉斯統計推理中，一切都被看作是一個服從概率分佈的隨機變量。利用貝葉斯定理，我們將先驗隨機概率分佈 (prior probability distribution)，和觀察數據作條件概率 (condition on the observed data)，從而獲得事後概率分佈 (posterior probability distribution)。
### 醫學診斷測試 diagnostic testing
貝葉斯推理最常用的實例是在診斷測試中，即當一個人拿着陽性的檢驗報告結果來找你，你如何判斷這個人有多大的概率真的患有該疾病。
用 $D$ 標記患病， $\bar{D}$ 標記不患病；$T$ 標記檢查結果爲陽性，$\bar{T}$ 標記檢查結果爲陰性。那麼，陽性檢查結果時，真的患病的概率 $\text{Pr}(D|T)$：
$$
\begin{aligned}
\text{Pr}(D|T) &= \frac{\text{Pr}(T|D)\text{Pr}(D)}{\text{Pr}(T)}\\
&=\frac{\text{Pr}(T|D)\text{Pr}(D)}{\text{Pr}(T|D)\text{Pr}(D)+\text{Pr}(T|\bar{D})\text{Pr}(\bar{D})}
\end{aligned}
$$
其中分母的轉換用到了 Law of Total Probability (L.T.P):
$$
\begin{aligned}
\text{Pr}(T) &= \text{Pr}(T \cap D) + \text{Pr}(T \cap \bar{D}) \\
&= \text{Pr}(T|D)\text{Pr}(D)+\text{Pr}(T|\bar{D})\text{Pr}(\bar{D})
\end{aligned}
$$
所以說，貝葉斯定理在這裏告訴我們，要計算 $\text{Pr}(D|T)$ 我們只需要下列幾個信息：
1. 患病率： $\text{Pr}(D)$
2. 檢測手段的敏感度 (sensitivity)： $\text{Pr}(T|D)$
3. 檢測手段的 1 - 特異度 (specificity)： $\text{Pr}(T|\bar{D})=1-\text{Pr}(\bar{T}|\bar{D})$
### HIV 檢查時的應用
假設人羣中患病率爲 $1/1000$，所用的 HIV 檢測手段的敏感度爲 $0.99$， 特異度爲 $0.98$。試計算該檢測HIV手段的事後概率（即拿到陽性結果時，患病的概率 $\text{Pr}(D|T)$）。
**解**
令 $D=\text{HIV positive}, \bar{D}=\text{HIV negative}\\
T=\text{test postive}, \bar{T}=\text{test negative}$
$$
\begin{aligned}
\text{Pr}(D|T) &= \frac{\text{Pr}(T|D)\text{Pr}(D)}{\text{Pr}(T|D)\text{Pr}(D)+\text{Pr}(T|\bar{D})\text{Pr}(\bar{D})} \\
&= \frac{0.99\times0.001}{0.99\times0.001+(1-0.98)\times0.999} \\
&= 0.0472
\end{aligned}
$$
如果 特異度能達到 $0.99$
$$
\begin{aligned}
\text{Pr}(D|T) &= \frac{\text{Pr}(T|D)\text{Pr}(D)}{\text{Pr}(T|D)\text{Pr}(D)+\text{Pr}(T|\bar{D})\text{Pr}(\bar{D})} \\
&= \frac{0.99\times0.001}{0.99\times0.001+(1-0.99)\times0.999} \\
&= 0.0901
\end{aligned}
$$
如果特異度能達到 $0.999$
$$
\begin{aligned}
\text{Pr}(D|T) &= \frac{\text{Pr}(T|D)\text{Pr}(D)}{\text{Pr}(T|D)\text{Pr}(D)+\text{Pr}(T|\bar{D})\text{Pr}(\bar{D})} \\
&= \frac{0.99\times0.001}{0.99\times0.001+(1-0.999)\times0.999} \\
&= 0.497
\end{aligned}
$$
可見，對於像 HIV 這樣人羣中患病率較爲罕見的疾病，其檢驗手段的敏感度，特異度都要達到極高才能讓檢驗結果可靠，即拿到陽性結果的人的確患有該疾病。其中當敏感度爲 $0.99$，特異度爲 $0.999$ 時，才能讓這樣的檢驗手段達到接近一半的可靠程度 (即只有接近一半的陽性結果是真陽性)。
注意本例爲貝葉斯理論的特例，即我們使用的是一個固定的先驗概率 (prior) 和似然 (likelihood)。一般情況下，先驗概率和似然會有自己的概率分佈 (probability distribution)，而很少會是一個固定的值， 其相應的事後概率 (posterior) 也擁有概率分佈，並且使用它本身的均值和方差來描述。
### 說點小歷史
```{r Bayes00, cache=TRUE, echo=FALSE, fig.asp=.7, fig.width=4, fig.cap='Sir Ronald Fisher', fig.align='center', out.width='50%'}
knitr::include_graphics("img/Fisher.jpg")
```
[Ronald Aylmer Fisher (1890-1962)](https://en.wikipedia.org/wiki/Ronald_Fisher) 推動了統計學在20世紀前半頁的重大發展。他鞏固了概率論統計學堅實的基礎，並且積極提倡這一套理論[@Fisher1922]。但是 Fisher 本人對於統計學的“統計學意義, level of significance” 的認識卻是隨着時間和他年齡的變化而變化的：
```{r Bayes01, cache=TRUE, echo=FALSE, warning=FALSE, eval=FALSE}
dt <- read.csv("backupfiles/Fisher.csv", header = T)
kable(dt, "html",align = "c",caption = "Fisher's interpretation of 'level of significance' and the Neyman-Pearson interpretation") %>%
kable_styling(bootstrap_options = c("striped", "bordered"))
```
<table class="table table-striped table-bordered" style="margin-left: auto; margin-right: auto;">
<caption>表 40.1: Fisher's interpretation of 'level of significance' and the Neyman-Pearson interpretation</caption>
<thead><tr>
<th style="text-align:center;"> 早期 Fisher (1935) </th>
<th style="text-align:center;"> 晚期 Fisher (1956) </th>
<th style="text-align:center;"> Neyman and Pearson </th>
</tr></thead>
<tbody><tr>
<td style="text-align:left;"> 統計學有意義的水平（傳統上使用 $\alpha=5\%$），必須在實施統計檢驗之前就被決定。因此，統計學意義的水平是相應統計學檢驗本身的性質之一。<br>
Thus, the level of significance is a property of the _test_. </td>
<td style="text-align:left;"> 統計學意義的水平，應該被精確計算並且在報告中明確 $p$ 值的大小，故統計學意義的水平本身是在實施了統計檢驗之後計算的。它應該是屬於觀察數據的固有性質。 <br>
Here the level of significance is a property of the _data_. </td>
<td style="text-align:left;"> $\alpha$ 和 $\beta$ 作爲統計檢驗的第一類錯誤和第二類錯誤指標，應該在實施統計檢驗之前被決定。所以 $\alpha, \beta$ 是屬於統計檢驗的性質。<br>
Yet, to determine $\alpha, \beta$ no convention is required, but rather a cost-benefit estimation of the severity of the two kinds of error. </td>
</tr></tbody>
</table>
隨着马尔科夫蒙特卡洛 (Markov-Chain Monte Carlo, MCMC) 法的廣泛應用，貝葉斯統計學在事後概率計算上（計算量超大的）棘手問題，得到了解決。
## 練習題
1. 從經典概率論的角度，準確定義 $95\%$ 信賴區間。思考，在貝葉斯統計理論中，它會如何被定義。
**解**
<u>概率論：</u>
對於一個總體參數 $\theta$ 來說，$95\%$ 信賴區間是一個從觀察數據中計算得到的數值區間。如果重複相同的實驗無數次，我們從無數個觀察數據中計算這個區間，那麼這些無數多的信賴區間 (confidence interval, CI) 裏有 $95\%$ 包含了總體參數 $\theta$。
<u>貝葉斯：</u>
對於一組觀察數據，它可以計算獲得可信區間 (credible interval, CI)。如果使用 $L, U$ 分別表示下限和上限的值，$\theta$ 表示參數，$x$ 表示觀察數據，$\pi(\theta|x)$ 表示事後概率分佈的密度方程， posterior distribution。那麼有：
$$\text{Pr}(\theta \in (L,U)) = \int_L^U\pi(\theta|x)\text{d}\theta = 95\%$$
**即，在貝葉斯理論下，95% 可信區間就是這一個區間包含了參數的概率是95%。**
2. 證明貝葉斯定理。
a. 並且用二項分佈隨機變量的例子來證明：<br>$\text{posterior odds} = \text{prior odds}\times\text{likelihood ratio}$
b. 用前面提到的 HIV 的案例來說明這個公式的實際應用。
**解**
參照上面的標記法：
- $\theta$ 表示參數
- $x$ 表示觀察數據
- $\pi(\theta|x)$ 表示事後概率分佈的密度方程， posterior distribution
- $f(\theta,x)$ 表示參數和數據的聯合分佈， joint distribution
- $f(x)$ 表示先驗概率分佈的密度方程， prior distribution
$$
\begin{aligned}
\pi(\theta|x) &= \frac{f(\theta, x)}{f(x)} \\
&=\frac{f(\theta, x)}{f(x)}\cdot\frac{1/\pi(\theta)}{1/\pi(\theta)} \\
&=\frac{\frac{f(\theta,x)}{\pi(\theta)}}{\frac{f(x)}{\pi(\theta)}}
\end{aligned}
$$
其中分子部分 $\frac{f(\theta,x)}{\pi(\theta)}$ 就是條件概率 $f(x|\theta)$。
分母的 $f(x)$ 部分
$$
\begin{aligned}
f(x) &= \int f(x,\theta) \text{d}\theta \\
&= \int \frac{f(x,\theta)}{\pi(\theta)} \cdot \pi(\theta) \text{d}\theta \\
&= \int f(x|\theta) \cdot \pi(\theta) \text{d}\theta
\end{aligned}
$$
所以，
$$\pi(\theta|x)=\frac{f(x|\theta)\pi(\theta)}{\int f(x|\theta) \cdot \pi(\theta) \text{d}\theta}$$
a. 用二項分佈隨機變量 ($\theta=1, 0$) 來證明：**$\text{posterior odds} = \text{prior odds}\times\text{likelihood ratio}$**
**解**
假設 $\theta$ 是一個二項分佈的隨機變量，那麼 $f(\theta|x)=\text{Pr}(\theta |x)$。
$$
\begin{aligned}
\text{posterior odds} &= \frac{\text{Pr}(\theta=1|x)}{\text{Pr}(\theta=0|x)} \\
&= \frac{\frac{\text{Pr}(x|\theta=1)\text{Pr}(\theta=1)}{\text{Pr}(x)}}{\frac{\text{Pr}(x|\theta=0)\text{Pr}(\theta=0)}{\text{Pr}(x)}}\\
&=\frac{\text{Pr}(\theta=1)}{\text{Pr}(\theta=0)}\cdot\frac{\text{Pr}(x|\theta=1)}{\text{Pr}(x|\theta=0)} \\
&=\text{prior odds}\times\text{likelihood ratio}
\end{aligned}
$$
b. 用前面提到的 HIV 案例來驗證：
HIV的患病率爲 $1/1000$，所以 $\text{prior odds}=1:999$，似然比 $\text{likelihood ratio}=0.99:(1-0.98)$。所以就有：
$$
\begin{aligned}
\text{posterior odds} &=\text{prior odds}\times\text{likelihood ratio} \\
&= \frac{1}{999}\times\frac{0.99}{1-0.98} \\
&= \frac{0.99}{19.98} \\
&= \frac{1}{20.18182}
\end{aligned}
$$
所以事後概率（陽性結果患病的概率）爲 $1/(1+20.18182)=0.0472$。
3. 史密斯先生有2個孩子，其中之一是男孩。另一個孩子是女孩的概率是多少？ 如下前提默認成立：
1. 男女比例爲: 50-50。
2. 這個家庭中沒有對男孩或者女孩的偏好。
3. 這兩個孩子不是同胞雙胞胎。
一個家庭有兩個孩子的性別組合的所有可能性：
<table class="table table-striped table-bordered" style="width: auto !important; margin-left: auto; margin-right: auto;">
第一個孩子性別 | 第二個孩子性別
:------------:|:--------------:
男孩 | 男孩
男孩 | 女孩
女孩 | 男孩
女孩 | 女孩
</table>
所以根據已知條件，其中之一是男孩，所以最後一種情況：“兩個女孩” 是不可能的。故另一孩子是女孩的概率就是 $\frac{2}{3}$。
如果用貝葉斯理論來正式計算的話：
$$
\begin{aligned}
& \text{Pr (1 girl in family of 2 | family does not have 2 girls)} \\
&= \frac{\text{Pr(family doesn't have 2 girls|1 girl in a family of 2)}\times \\ \text{Pr(1 girl in a family of 2 )}}{\sum_{j=0,1,2}\text{Pr(family doesn't have 2 girls|j girl in a family of 2)}\times\\\text{Pr(j girl in a family of 2)}} \\
&= \frac{1\times\frac{1}{2}}{1\times\frac{1}{4}+1\times\frac{1}{2}+0\times\frac{1}{4}} \\
&= \frac{\frac{1}{2}}{\frac{3}{4}}=\frac{2}{3}
\end{aligned}
$$
也是一樣的結論。
4. 下表是全國普查以後得出的家庭有兩個孩子，**且至少一個是男孩的數據分佈**：
<table class="table table-striped table-bordered" style="margin-left: auto; margin-right: auto;">
第一個孩子性別 | 第二個孩子性別 | 家庭數量
:------------:|:--------------:|:----:
男孩 | 男孩 | 657
男孩 | 女孩 | 591
女孩 | 男孩 | 610
女孩 | 女孩 | 0
</table>
求同樣的概率問題：
**解**
另一個孩子是女孩的概率是：$\frac{610+591}{610+591+657}=0.646$
# 貝葉斯定理的應用：單一參數模型 {#single-para-model}
從前一章節我們可以深切體會到，貝葉斯統計是如何讓我們的先驗概率，在觀察到數據之後，更新信息，獲得事後概率 (這是一個通過數據自我學習，進化的過程)。(How a prior belief about an event can be updated, given data, to a posterior belief.)
所以說，在貝葉斯模型中，我們期待使用觀察數據來學習，以增加現有的對相關參數的知識和信息。本章我們把重點放在二項分佈，用二項分佈作爲單一參數模型來瞭解怎樣推導事後分佈。
## 貝葉斯理論下的事後二項分佈概率密度方程 notation for probability density functions
- $R$ 用來表示服從一個二項分佈的隨機變量， $R\sim Bin(n, \theta)$。
- $r$ 表示觀察到 $r$ 次成功實驗，實驗次數爲 $n$。
- 先驗概率分佈： $\pi_\Theta(\theta)$
- 應用貝葉斯定理：
$$
\begin{aligned}
\pi_{\Theta|R}(\theta|r) &= \frac{f_R(r|\theta)\pi_\Theta(\theta)}{\int_0^1f_R(r|\theta)\pi_\Theta(\theta)\text{ d}\theta}\\
&= \frac{f_R(r|\theta)\pi_\Theta(\theta)}{f_R(r)}
\end{aligned}
$$
如果我們的先驗概率分佈：
$$\begin{equation}
\pi_\Theta(\theta)=\begin{cases}
1 \text{ if } \theta=0.2\\
0 \text{ otherwise}
\end{cases}
\end{equation}$$
意思就是，我們 100% 相信 $\theta$ 絕對就等於 0.2，不相信 $\theta$ 竟然還能取任何其他值（霸道自大又狂妄的我們）。
如果先驗概率分佈：
$$\begin{equation}
\pi_\Theta(\theta)=\begin{cases}
0.4 \text{ if } \theta=0.2\\
0.6 \text{ if } \theta=0.7
\end{cases}
\end{equation}$$
意思就是，我們有 60% 的把握相信 $\theta=0.7$，有 40% 的把握相信 $\theta=0.2$，稍微傾向於 $\theta=0.7$。
假設進行10次實驗，觀察到3次成功。當 $\theta=0.2$ 時，觀察數據的似然 (likelihood) 爲：
$$f_R(3|\theta=0.2)=\binom{10}{3}0.2^3(1-0.2)^7$$
當 $\theta=0.7$ 時，觀察數據的似然爲：
$$f_R(3|\theta=0.7)=\binom{10}{3}0.7^3(1-0.7)^7$$
應用貝葉斯定理計算事後概率分佈：
$$\begin{equation}
\pi_{\Theta|R}(\theta|3)=\begin{cases}
\frac{\binom{10}{3}0.2^3(1-0.2)^7\times0.4}{\binom{10}{3}0.2^3(1-0.2)^7\times0.4+\binom{10}{3}0.7^3(1-0.7)^7\times0.6}=0.937 \text{ if } \theta=0.2\\
\frac{\binom{10}{3}0.7^3(1-0.7)^7\times0.6}{\binom{10}{3}0.7^3(1-0.7)^7\times0.6+\binom{10}{3}0.2^3(1-0.2)^7\times0.4}=0.063 \text{ if }\theta=0.7
\end{cases}
\end{equation}$$
所以，我們從一開始認爲只有40%的把握相信 $\theta=0.2$，觀察數據告訴我們 10 次實驗，3次獲得了成功。所以我們現在有 93.7% 的把握相信 $\theta=0.2$。也就是說，觀察數據讓我們對參數 $\theta$ 的取值可能性發生了質的變化，從原先的傾向於 $\theta=0.7$ 到現在幾乎接近 100% 的認爲 $\theta=0.2$。也就是，觀察數據獲得的信息改變了我們的立場。
上面的例子很直觀，但是有下面幾個問題：
1. 如果我們無法對參數 $\theta$ 賦予先驗概率的點估計時，該怎麼辦？
2. 如果事後概率不是一個離散的分佈時，該如何才能表達事後概率？
## $\theta$ 的先驗概率
一種選擇是，我們用均一分佈 (uniform distribution)，即我們對數據一無所知，認爲所有的 $\theta$ 的可能性都一樣，概率密度方程爲 $1$。在這一情況下，先驗概率爲 1： $\pi_\Theta(\theta)=1$，其事後概率分佈爲：
$$
\begin{equation}
\pi_{\Theta|R}(\theta|r)=\frac{\binom{n}{r}\theta^r(1-\theta)^{n-r}}{\int_0^1\binom{n}{r}\theta^r(1-\theta)^{n-r} \text{ d}\theta}
(\#eq:uniformBayes)
\end{equation}
$$
看到即使在如此簡單的先驗概率下，我們還是要使用複雜的微積分進行計算。幸運的是，像 \@ref(eq:uniformBayes) 的分母這樣的積分公式其實是有跡可循的。這就是 beta ($\beta$) 分佈。
### beta 分佈 the beta distribution {#beta-distribution-intro}
```{r beta-distr, echo=FALSE, fig.height=7, fig.width=6, fig.cap='Beta distribution functions for various values of a, b', fig.align='center', out.width='90%', cache=TRUE}
par(mfrow=c(3,2))
pi <- Vectorize(function(theta) dbeta(theta, 1,1))
curve(pi, xlab=~ theta, ylab="Density", main="Beta prior: a=1, b=1", frame=FALSE, lwd=2)
pi <- Vectorize(function(theta) dbeta(theta, 0.3,1))
curve(pi, xlab=~ theta, ylab="Density", main="Beta prior: a=0.3, b=1",frame=FALSE, lwd=2)
pi <- Vectorize(function(theta) dbeta(theta, 0.3,0.3))
curve(pi, xlab=~ theta, ylab="Density", main="Beta prior: a=0.3, b=0.3",frame=FALSE, lwd=2)
pi <- Vectorize(function(theta) dbeta(theta, 2,2))
curve(pi, xlab=~ theta, ylab="Density", main="Beta prior: a=2, b=2",frame=FALSE, lwd=2)
pi <- Vectorize(function(theta) dbeta(theta, 8,2))
curve(pi, xlab=~ theta, ylab="Density", main="Beta prior: a=8, b=2",frame=FALSE, lwd=2)
pi <- Vectorize(function(theta) dbeta(theta, 2,8))
curve(pi, xlab=~ theta, ylab="Density", main="Beta prior: a=2, b=8",frame=FALSE, lwd=2)
```
你會發現，beta分佈的圖形特徵由它的兩個超參數 (hyperparameter) `a, b` 決定。當相對地 `a` 比較大的時候，beta分佈的概率多傾向於較靠近橫軸右邊，也就是1的位置（較高概率），相對地 `b` 比較大的時候，beta分佈的曲線多傾向於靠近橫軸左邊，也就是0的位置（較低概率）。如果 `a` 和 `b`同時都在變大的話，beta分佈的曲線就變得比較“瘦”。這樣決定概率分佈形狀的參數，又被叫做**形狀參數 (shape parameters)**
我們定義 $a>0$ 時[伽馬方程](https://zh.wikipedia.org/wiki/%CE%93%E5%87%BD%E6%95%B0)爲
$$\Gamma(a)=\int_0^\infty x^{a-1}e^{-ax}\text{ d}x$$
當 $a$ 取正整數時， $\Gamma(a)$ 是 $(a-1)!$。例如，當 $a=4, \Gamma(a)=3\times2\times1=6$。
對於 $\theta\in[0,1]$ 時，beta 方程 $Beta(a,b)$ 被定義爲：
$$
\begin{aligned}
\pi\Theta(\theta|a,b) &= \theta^{a-1}(1-\theta)^{b-1}\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\\
&= \frac{\theta^{a-1}(1-\theta)^{b-1}}{B(a,b)}
\end{aligned}
$$
其中
$$B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)} = \int_0^1\theta^{a-1}(1-\theta)^{b-1}\text{d}\theta$$
**莫要混淆 B 函數和 Beta 方程。B 函數的性質和詳細定義可以參照[其維基百科頁面](https://zh.wikipedia.org/wiki/%CE%92%E5%87%BD%E6%95%B0)。**
利用 Beta 方程作爲先驗概率顯得十分便捷且靈活。圖 \@ref(fig:beta-distr) 展示的是 6 種不同的 $(a,b)$ 取值下的先驗概率分佈示意圖。其實我們可以看到，包括均一分佈在內的各種可能性都可以通過 Beta 分佈實現。其中 $(a,b)$ 被叫做超參數 (hyperparameter)。$(a,b)$ 取值越大，先驗概率分佈的方差越小。
關於 Beta 分佈的幾個性質：
1. 均值：$\text{mean}=\frac{a}{a+b}$；
2. 衆數：$\text{mode}=\frac{a-1}{a+b-2}$；
3. 方差：$\text{variance}=\frac{ab}{(a+b)^2(a+b+1)}$。
回到均一分佈的簡單例子 \@ref(eq:uniformBayes) 上：
$\pi_\Theta(\theta)=Beta(1,1)$ 是 $\theta\in[0,1]$ 上的均一分佈。所以事後概率 posterior 和下面的式子成正比：
$$\theta^r(1-\theta)^{n-r}$$
換句話說，事後概率分佈服從 $Beta(r+1,n-r+1)$ 分布，均值爲 $\frac{r+1}{n+2}$，方差爲 $\frac{(1+r)(n-r+1)}{(n+2)^2(n+3)}$。
由此可見，在貝葉斯統計思維下，先驗概率爲均一分佈的二項分佈數據，其事後概率分佈的均值和方差，和經典概率論下的極大似然估計 $r/n$ 不同，和它的漸進樣本方差 $r(n-r)/n^3$ 也不同。但是，當 $n$ 越來越大，獲得的觀察數據越多提供的信息越來越多以後，我們會發現事後概率分佈的均值和方差也會越來越趨近於經典概率論下的極大似然估計和它的方差。
於是這裏可以總結以下兩點：
1. 即使先驗概率對參數毫無用處（不能提供有效信息，或者我們對所觀察的數據一無所知），也可能會對事後概率分佈結果提供一些意外的信息。
2. 當樣本量增加，似然就主導了整個貝葉斯方程，在數學計算上，經典概率論和貝葉斯推理的估計結果將會十分接近。當然，其各自的意義還是截然不同的。
### 二項分佈數據事後概率分佈的一般化：共軛性 {#conjugate}
當 $r\sim \text{Binomial}(n,\theta)$ 時，如果先驗概率 $\pi_\Theta(\theta)=\text{Beta}(a,b)$。那麼參數 $\theta$ 的事後概率分佈的密度方程滿足：
$$\pi_{\Theta|r}(\theta|r)=\text{Beta}(a+r, b+n-r)$$
它的事後概率分佈均值爲：
$$E[\theta|r]=\frac{a+r}{a+b+n}$$
事後概率分佈的衆數爲：
$$\text{Mode}[\theta|r]=\frac{a+r-1}{n+a+b-2}$$
事後概率分佈方差爲：
$$\text{Var}[\theta|r]=\frac{(a+r)(b+n-r)}{(a+b+n)^2(a+b+n+1)}$$
因此，我們看到先驗概率服從 $\text{Beta}(a,b)$ 分佈，觀察數據爲二項分佈時，事後概率分佈還是服從 $\text{Beta}$ 分佈，僅僅只是超參數發生了轉變（更新）。這就是**共軛分佈**的實例。$\text{Beta}$ 分佈是二項分佈的共軛先驗概率分佈 (the $\text{Beta}(a,b)$ is the conjugate prior for the binomial likelihood)。
在經典概率論的框架下，參數 $\theta$ 的估計就是極大似然估計 (MLE)。在二項分佈的例子中， $\text{MLE}=\hat\theta=r/n$，當樣本量 $n\rightarrow\infty$ 時，事後概率分佈均值：
$$E[\theta|r]=\frac{a+r}{a+b+n}=\frac{\frac{r}{n}+\frac{a}{n}}{1+\frac{a+b}{n}}\approx\frac{r}{n}=\text{MLE}$$
事後概率分佈的衆數爲：
$$
\begin{aligned}
\text{Mode}[\theta|r] &=\frac{a+r-1}{n+a+b-2} \\
&= \frac{\frac{r}{n}+\frac{a-1}{n}}{1+\frac{a+b-2}{n}}\\
&\approx \frac{r}{n}
\end{aligned}
$$
事後概率分佈的方差爲：
$$\frac{(a+r)(b+n+r)}{(a+b+n)^2(a+b+n+1)}\approx0$$
當 $n$，樣本越來越大時，我們獲得更多的來自數據的信息，所以來自數據的信息逐漸主導 (dominate) 了整個貝葉斯推斷的過程，事後均值等衆多統計結果都越來越趨近於概率論統計思想下的極大似然估計等結論。
我們也可以注意到，當 $a\rightarrow0, b\rightarrow0$ 時，事後概率分佈的均值 $E[\theta|r] = \frac{a+r}{a+b+n} \rightarrow \frac{r}{n}$，方差也趨向於樣本漸進方差 (asymptotic sample variance)。但是當 $a\rightarrow 0, b\rightarrow0$ 時，先驗概率是沒有被定義的，可是此例下事後概率卻可以正常被定義。所以當先驗概率分佈無法被定義，或者被定義的不恰當時，事後概率分佈依然不受太大影響。所以特別是對於均值（或迴歸係數，regression coefficients）等參數，我們常常會使用均一分佈這樣的無信息先驗概率。
## 附贈--加量不加價
**當樣本數據服從二項分布時使用 $\text{Beta}(a,b)$ 作爲先驗概率分布，試推導事後概率分布。並且證明 Beta 分布是二項分布的共軛分布 (conjugate distribution)。**
- 當數據服從二項分布時，其似然方程 (或概率方程) 爲: <br> $$f(n,r|\theta) = \binom{n}{r}\theta^r(1-\theta)^{n-r}$$
- 用 $\text{Beta}(a,b)$ 做先驗概率分布 ($a,b$ 是 $\theta$ 的超參數)，那麼 <br> $$\begin{aligned}\pi(\theta) &= \text{Beta}(\theta|a,b) \\
& = \theta^{a-1}(1-\theta)^{b-1}\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \\
& = \frac{\theta^{a-1}(1-\theta)^{b-1}}{B(a,b)} \\
\text{Where } B(a,b) &=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)} = \int_0^1\theta^{a-1}(1-\theta)^{b-1}\text{d}\theta
\end{aligned}$$
- 利用貝葉斯定理，後概率分布就可以推導爲 <br>
$$
\begin{aligned}
\pi(\theta|n,r) & = \frac{f(n,r|\theta)\pi(\theta)}{\int f(n,r |\theta)\pi(\theta)\text{d}\theta} \\
& = \frac{\binom{n}{r}\theta^r(1-\theta)^{n-r}\pi(\theta)}{\int_0^1\binom{n}{r}\theta^r(1-\theta)^{n-r} \pi(\theta) \text{ d}\theta} \\
& = \frac{\binom{n}{r}\theta^r(1-\theta)^{n-r}\frac{\theta^{a-1}(1-\theta)^{b-1}}{B(a,b)}}{\int_0^1\binom{n}{r}\theta^r(1-\theta)^{n-r} \frac{\theta^{a-1}(1-\theta)^{b-1}}{B(a,b)} \text{ d}\theta} \\
& = \frac{\binom{n}{r}\theta^{r+a-1}(1-\theta)^{n+b-r-1}}{\int_0^1\binom{n}{r}\theta^{r+a-1}(1-\theta)^{n+b-r-1} \text{ d}\theta} \\
& = \frac{\theta^{r+a-1}(1-\theta)^{n+b-r-1}}{B(r+a, n+b-r)} \\
& \sim \text{Beta}(r+a, n+b-r)
\end{aligned}
$$
這個推導的結果告訴我們，先驗概率 $\pi(\theta) = \frac{\theta^{a-1}(1-\theta)^{b-1}}{B(a,b)}\sim \text{Beta}(a,b)$，通過和二項分布的似然方程相乘，獲得後驗概率 $\pi(\theta|n,r) = \frac{\theta^{r+a-1}(1-\theta)^{n+b-r-1}}{B(r+a, n+b-r)}\sim\text{Beta}(r+a, n+b-r)$。這兩個概率都服從 $\text{Beta}$ 分布，也就是證明了 $\text{Beta}$ 分布，是二項分布的共軛分布。新產生的事後概率的性質有: <br>
1. $\theta$ 的均值 mean $\mu= \frac{a+r}{a+b+n}$;
2. $\theta$ 的方差 variance $\sigma^2 = \frac{(a+r)(n+b-r)}{(a+b+n)^2(a+b+n+1)}$;
## 練習題
### Q1
當先驗概率分佈服從 $\text{Beta}(0.5,0.5)$，觀察數據記錄到 $5$ 個患者中 $3$ 人死亡的事件。
試求：
死亡發生概率 $\theta$ 的 95% 可信區間 (credible intervals)。
**解**
根據 Section \@ref(conjugate) 的公式，當先驗概率爲 $\pi_{\Theta|r}(\theta|r)=\text{Beta}(a=0.5,b=0.5)$ ，數據 $n=5, r=3$。參數 $\theta$ 的事後概率分佈 $\pi_{\Theta|r}(\theta|r)=\text{Beta}(a+r,b+n-r)=\text{Beta}(3.5,2.5)$。
在 [R](https://www.r-project.org/) 裏進行貝葉斯計算十分簡便：
```{r Bayes02, cache=TRUE}
# 95% Credible Intervals
L <- qbeta(0.025, 3.5, 2.5)
U <- qbeta(0.975, 3.5, 2.5)
print(c(L,U))
```
事後分佈 $\pi_{\Theta|r}(\theta|r)=\text{Beta}(3.5,2.5)$ 的分佈圖形如下：
```{r Bayes03, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Posterior distribution of Beta(3.5,2.5)', out.width='80%', cache=TRUE}
post <- Vectorize(function(theta) dbeta(theta, 3.5, 2.5))
# Illustration
x <- seq(0,1,length=10000)
y <- post(x)
plot(x,y, type = "l", xlab=~theta, ylab="Density", lwd=2, frame.plot = FALSE)
polygon(c(L, x[x>=L & x<= U], U), c(0, y[x>=L & x<=U], 0), col="grey")
```
我們可以自己寫一個求可信區間的公式來計算：
```{r Bayes04, cache=TRUE}
# Credible Interval function:
# a, b : shape / super parameters
# level: probability level (0,1)
cred.int <- function(a,b,level){
L <- qbeta((1-level)/2, a, b) # Lower limit
U <- qbeta((1+level)/2, a, b) # Upper limit
return(c(L,U))
}
cred.int(3.5,2.5,0.95)
```
95%可信區間 $(0.2094, 0.9056)$ 告訴我們，參數 $\theta\in(0.2094, 0.9056)$ 的概率是 $0.95$。
下面我們嘗試寫 Beta 分佈的其他統計量：均值，衆數，方差等。
```{r Bayes05, cache=TRUE}
# a, b: shape / super parameters
MeanBeta <- function(a,b) a/(a+b)
ModeBeta <- function(a,b) {
m <- ifelse(a>1 & b>1, (a-1)/(a+b-2), NA)
return(m)
}
VarianceBeta <- function(a,b) (a*b)/((a+b)^2*(a+b+1))
# mean
MeanBeta(3.5,2.5)
# mode
ModeBeta(3.5,2.5)
# Variance
VarianceBeta(3.5, 2.5)
# SD
sqrt(VarianceBeta(3.5, 2.5))
```
### Q2
假如數據還是 Q1 的數據，然而先驗概率讓我們認爲可能在 5 名受試對象中觀察到 1 次事件。
1. 試求超參數 $(a,b)$ 滿足先驗概率的 Beta 分佈。(不止一組)
**解**
我們認爲最有可能發生 “5 名受試對象中觀察到 1 次事件” 的情況，那麼先驗概率的均值爲 $\frac{a}{a+b}=0.2$。所以，在實數中有無數組超參數都可以用來模擬先驗概率分佈。例如 $a=1, b=4; a=10, b=40; a=100, b=400; a=0.317, b=1.268, \cdots$。
2. 假如觀察數據是 $n=5, r=1$，計算事後概率分佈及其均值，標準差。
**解**
先來嘗試寫一個計算貝葉斯二項分佈的方程：
```{r 08-Intro-to-Bayes-1}
# binbayes function in R
#------------------------------
# a, b: shape / super parameters
# r : number of successes
# n : number of trials
binbayes <- function(a, b, r, n) {
prior <- c(a, b, NA, MeanBeta(a,b), sqrt(VarianceBeta(a, b)), qbeta(0.025, a, b), qbeta(0.5, a, b), qbeta(0.975, a, b))
posterior <- c(a+r, b+n-r, r/n, MeanBeta(a+r, b+n-r), sqrt(VarianceBeta(a+r, b+n-r)), qbeta(0.025, a+r, b+n-r), qbeta(0.5, a+r, b+n-r), qbeta(0.975, a+r, b+n-r))
out <- rbind(prior, posterior)
out <- round(out, 4)
colnames(out) <- c("a","b","r/n","Mean", "SD", "2.5%", "50%", "97.5%")
return(out)
}
# a=1, b=4, r=1, n=5
binbayes(1,4,1,5)
# a=10, b=40, r=1, n=5
binbayes(10,40,1,5)
```
通過繪製先驗概率分佈圖和事後概率分佈圖來比較二者的變化：
```{r 08-Intro-to-Bayes-2, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Prior (dashed) Beta(1,4) vs. Posterior (cont.) Beta(2,8)', out.width='80%', cache=TRUE}
# Prior vs posterior graphs
# a,b : shape / super parameters
# r : number of successes
# n : number of trials
graph.binbayes <- function(a,b,r,n) {
prior <- Vectorize(function(theta) dbeta(theta, a,b))
posterior <- Vectorize(function(theta) dbeta(theta, a+r, b+n-r))
YL <- max(prior(seq(0.001,0.999,by=0.001)),posterior(seq(0.001,0.999,by=0.001)))
curve(prior, xlab=~theta, ylab="Density", lwd=1, lty=2,n=10000,ylim=c(0,YL), frame.plot = FALSE)
curve(posterior, xlab=~theta,lwd=2,lty = 1, add=T,n=10000)
}
graph.binbayes(1,4,1,5)
```
```{r 08-Intro-to-Bayes-3, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Prior (dashed) Beta(10,40) vs. Posterior (cont.) Beta(11, 44)', out.width='80%'}
graph.binbayes(10,40,1,5)
```
我們可以很清楚的看見，先驗概率相同時，$\text{Beta} (a,b)$ 的超參數如果越大，先驗概率的分佈就越趨近與對稱圖形，且極大值也就越出現在均值的地方 (本例中是 $0.2$)。而且也會使事後概率的 HPD (highest posterior density) 的區間更狹窄 (意爲對事後概率的預測越準確)，同時事後概率分佈也更加接近左右對稱。
### Q3
我們事先估計某個事件在 $n=20$ 名患者中發生的概率爲 $15\%$。當實際觀察數據爲 $n=15,r=3$ 時，計算相應的事後概率。
**解**
```{r 08-Intro-to-Bayes-4}
# because 15% events happened in 20 subjects, assuming prior Beta(a=3, b=17)
# observed n=15, r=3
binbayes(3,17,3,15)
```
```{r 08-Intro-to-Bayes-5, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Prior (dashed) Beta(3,17) vs. Posterior (cont.) Beta(6,29)', out.width='80%'}
graph.binbayes(3,17,3,15)
```
試着繪製先驗概率服從 $\text{Beta} (1,1)$，回憶之前本章開頭的圖 \@ref(fig:beta-distr)，這個先驗概率的含義就是我們沒有任何背景知識，對數據完全陌生的情況：
```{r 08-Intro-to-Bayes-6, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Prior (dashed) Beta(1,1) vs. Posterior (cont.) Beta(4,13)', out.width='80%'}
graph.binbayes(1,1,3,15)
```
### Q4
試給出上面各題中參數 $\theta$ 落在 $(0.1,0.25)$ 之間的概率。
```{r 08-Intro-to-Bayes-7}
# function to calculate probabilities in a interval
# a, b: super parameters
# r : number of successes
# n : number of trials
# L : Lower limit of the probability interval
# U : Upper limit of the probability interval
prob.int <- function(a,b,r,n,L,U){
prior0 <- pbeta(U,a,b) - pbeta(L,a,b)
posterior0 <- pbeta(U,a+r,n-r+b) - pbeta(L,a+r,n-r+b)
prob <- as.matrix(c(prior0, posterior0))
prob <- round(prob,4)
colnames(prob) <- paste("Probability of theta lies between the Interval", L, U)
rownames(prob) <- c("Prior","Posterior")
return(prob)
}
# Prior Beta(0.317,1.286) n = 5, r=1
binbayes(0.317, 1.286, 1, 5)
prob.int(0.317, 1.286, 1, 5,0.1,0.25)
```
```{r 08-Intro-to-Bayes-8, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Prior (dashed) Beta(0.317,1.286) vs. Posterior (cont.) Beta(1.317, 5.286)', out.width='80%'}
graph.binbayes(0.317, 1.286, 1, 5)
```
```{r 08-Intro-to-Bayes-9}
# Prior Beta(10,40) n = 5, r=1
binbayes(10, 40, 1, 5)
prob.int(10, 40, 1, 5,0.1,0.25)
```
所以在範圍固定的時候，事後概率分佈總是能夠比先驗概率分佈給出更高的累計概率。
### Q5
一個臨牀試驗要進行兩個階段 (two phases)，第一階段我們觀察到 $10$ 個患者中 $1$ 個事件。第二階段，觀察到 $n=50, r=5$。
1. 兩個階段都使用 $\text{Beta}(1,1)$ 作先驗概率。求兩個實驗階段參數 $\theta<0.1$ 的概率。
```{r 08-Intro-to-Bayes-10}
# Phase I
binbayes(1, 1, 1, 10)
prob.int(1,1,1,10,0,0.1)
```
```{r 08-Intro-to-Bayes-11, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Prior (dashed) Beta(1,1) vs. Posterior (cont.) Beta(2, 10)', out.width='80%'}
graph.binbayes(1, 1, 1, 10)
```
```{r 08-Intro-to-Bayes-12}
# Phase II
binbayes(1, 1, 5, 50)
prob.int(1,1,5,50,0,0.1)
```
```{r 08-Intro-to-Bayes-13, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Prior (dashed) Beta(1,1) vs. Posterior (cont.) Beta(6, 46)', out.width='80%'}
graph.binbayes(1, 1, 5, 50)
```
2. 繼續使用先驗概率分佈 $\text{Beta}(1,1)$，合併兩個實驗階段，求此時的事後概率分佈，以及參數 $\theta<0.1$ 的概率。
```{r 08-Intro-to-Bayes-14}
# Combining both phases
binbayes(1, 1, 6, 60)
prob.int(1,1,6,60,0,0.1)
```
```{r 08-Intro-to-Bayes-15, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Prior (dashed) Beta(1,1) vs. Posterior (cont.) Beta(7, 55)', out.width='80%'}
graph.binbayes(1, 1, 6, 60)
```
3. 用第一階段的實驗結果做第二階段實驗的先驗概率分佈，再計算事後概率分佈，以及 $\theta<0.1$ 的概率。
```{r 08-Intro-to-Bayes-16}
# Using Phase I results as a prior for Phase II
binbayes(2, 10, 5, 50)
prob.int(2,10,5,50,0,0.1)
```
```{r 08-Intro-to-Bayes-17, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Prior (dashed) Beta(2,10) vs. Posterior (cont.) Beta(7, 55)', out.width='80%'}
graph.binbayes(2, 10, 5, 50)
```
第2，3兩個小問題提示我們，無論是將第一階段實驗結果作爲第二階段實驗的先驗假設還是將兩次實驗合併，最終的結果是不會改變的。Both approaches are equivalent.
### Q6
藥物 A 和藥物 B 都被批准用於治療某種疾病。在 5000 例病例中使用藥物 A，發現有 3 人發生了不良副作用。在另外 7000 例病例中使用藥物 B，發現只有 1 例發生了副作用。
1. 先使用單一分佈作爲先驗概率 (uniform prior: $\text{Beta}(1,1)$)。求藥物 A 和藥物 B 各自發生不良反應的事後概率。
**藥物 A**
```{r 08-Intro-to-Bayes-18, cache=TRUE}
# Drug A
binbayes(1,1, 3, 5000)
```
```{r 08-Intro-to-Bayes-19, echo=FALSE, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Prior (dashed) Beta(1,1) vs. Posterior (cont.) Beta(4, 4998)', out.width='80%', cache=TRUE}
priorA <- Vectorize(function(theta) dbeta(theta,1,1))
posteriorA <- Vectorize(function(theta) dbeta(theta,1+3,5000-3+1))
curve(priorA,0,0.0025, xlab=~theta, ylab="Density",
lwd=1, lty=2,n=10000,ylim=c(0,1200),frame=FALSE)
curve(posteriorA, xlab=~theta,lwd=2,lty = 1, add=T,n=10000)
```
**藥物 B**
```{r 08-Intro-to-Bayes-20}
# Drug B
binbayes(1,1, 1, 7000)
```
```{r 08-Intro-to-Bayes-21, echo=FALSE, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Prior (dashed) Beta(1,1) vs. Posterior (cont.) Beta(2, 7000)', out.width='80%', cache=TRUE}
priorB <- Vectorize(function(theta) dbeta(theta,1,1))
posteriorB <- Vectorize(function(theta) dbeta(theta,1+1,7000-1+1))
curve(priorB,0,0.0025, xlab=~theta, ylab="Density",
lwd=1, lty=2,n=10000,ylim=c(0,2700),frame=FALSE)
curve(posteriorB, xlab=~theta,lwd=2,lty = 1, add=T,n=10000)
```
2. 使用 $\text{Beta}(0.00001,0.00001)$ 作爲先驗概率，重複上面的計算
**藥物 A**
```{r 08-Intro-to-Bayes-22}
# Drug A
binbayes(0.00001, 0.00001, 3, 5000)
```
```{r 08-Intro-to-Bayes-23, echo=FALSE, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Prior (dashed) Beta(0.00001,0.00001) vs. Posterior (cont.) Beta(3, 4997)', out.width='80%', cache=TRUE}
priorA <- Vectorize(function(theta) dbeta(theta,0.00001,0.00001))
posteriorA <- Vectorize(function(theta) dbeta(theta,0.00001+3,5000-3+0.00001))
curve(priorA,0,0.0025, xlab=~theta, ylab="Density",
lwd=1, lty=2,n=10000,ylim=c(0,1500),frame=FALSE)
curve(posteriorA, xlab=~theta,lwd=2,lty = 1, add=T,n=10000)
```
**藥物 B**
```{r 08-Intro-to-Bayes-24}
# Drug B
binbayes(0.00001, 0.00001, 1, 7000)
```
```{r 08-Intro-to-Bayes-25, echo=FALSE, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Prior (dashed) Beta(0.00001,0.00001) vs. Posterior (cont.) Beta(1, 6999)', out.width='80%', cache=TRUE}
priorB <- Vectorize(function(theta) dbeta(theta,0.00001,0.00001))
posteriorB <- Vectorize(function(theta) dbeta(theta,1+0.00001,7000-1+0.00001))
curve(priorB,0,0.0025, xlab=~theta, ylab="Density",
lwd=1, lty=2,n=10000,ylim=c(0,7000),frame=FALSE)
curve(posteriorB, xlab=~theta,lwd=2,lty = 1, add=T,n=10000)
```
3. 現在使用概率論的計算信賴區間 (confidence intervals) 的方法，求上面數據的精確二項分佈 95% 信賴區間。之前兩問中使用的哪個先驗概率更加接近概率論算法？
```{r 08-Intro-to-Bayes-26}
#------------------------------------------------
# Binomial confidence intervals
#------------------------------------------------
# r : number of successes
# n : number of trials
# level: confidence level
binom.confint <- function(r,n,level){
p <- r/n
conf <- as.vector(binom.test(r,n,conf.level = 0.95)$conf.int)
out <- c(p,conf)
out <- as.matrix(t(round(out,8)))
colnames(out) <- c("MLE", "L", "U")
return(out)
}
# Drug A
binom.confint(3,5000,0.95)
# Drug B
binom.confint(1,7000,0.95)
```
明顯可以看到，先驗概率使用 $\text{Beta}(0.00001,0.00001)$ 時，事後概率的均值和可信區間的下限值更接近概率論算法。使用先驗概率 $\text{Beta}(1,1)$ 時，事後概率的可信區間的上限值更接近概率論算法。
4. 如果需要你來下結論說，藥物 B 和藥物 A 哪個更加安全？ 求 $\text{Pr}(\theta_B < \theta_A|data)$。
**解**
**貝葉斯**
在計算機的輔助下，這是一個十分簡單的計算。我們從各自的事後分佈中採集大量隨機樣本，然後求 $\theta_B-\theta_A$ 然後看有多少比例這個數值是小於零的就可以得出結論：
```{r 08-Intro-to-Bayes-27, cache=TRUE}
# Simulating from each posterior
set.seed(1001)
post.thetaA <- rbeta(1000000, 3, 4997)
post.thetaB <- rbeta(1000000, 1, 6999)
# Taking the differences
theta.diff0 <- post.thetaB - post.thetaA
```
```{r 08-Intro-to-Bayes-28, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Histogram of Drug B - Drug A', out.width='80%'}
# Histogram of the differences
hist(theta.diff0,probability = TRUE, breaks = 50, xlab=expression(theta[B] - theta[A]),main = "")
abline(v=0, col="red", lwd=2)
box()
```
也可以不採用直方圖而是使用連續曲線：
```{r 08-Intro-to-Bayes-29, fig.asp=.7, fig.width=5, fig.align='center',fig.cap='Density of Drug B - Drug A', out.width='80%'}
# Continuous version