《Modern Statistics for Modern Biology》Chapter 一： 離散數(shù)據(jù)模型的預測（1.1 - 1.3）
1.4、多維正態(tài)分布：DNA
-
當在模擬中有四種情況時候膏燃，例如在圖 1.9的box或者當我們計算
A、T、C聂抢、G四種堿基的數(shù)目時候揽趾，我們需要擴展二項式模型。
回想一下智袭，當使用二項式時奔缠，我們可以考慮兩個結果的不相等的概率，方法是將概率p=P(1)=p1分配給結果1吼野，將1?p=p(0)=p0分配給結果0校哎。當像
ATCG這樣超過兩種可能的輸出結果時候，我們可以考慮把球扔進不同大小的盒子里瞳步，這些盒子的大小對應于不同的概率闷哆。我們可以給這些概率貼上pA，pC单起，pG抱怔，pT的標簽。就像在二項式情形下嘀倒，所有可能結果的概率之和是1：pA+pC+pG+pT=1屈留。嘗試使用隨機數(shù)生成器，該生成器通過名為
runif的函數(shù)生成0到1之間的所有可能的數(shù)字测蘑。用它生成一個4個level(A灌危，C，G碳胳，T)的隨機變量勇蝙，其中pA=，pC=
挨约，pG=
數(shù)學公式. 多項式分布是計數(shù)的最重要的模型，R使用一個通式來計算計數(shù)的多項向量(x1,...,xm)的概率表示具有m個具有概率(p1,...,pm)的box的結果:
? 問題 1.7
- 假設我們有四個同樣可能的盒子翁锡。使用這個公式蔓挖，在第一個盒子中觀察到4，在第二個盒子中觀察到2盗誊，而在另外兩個盒子中觀察不到的概率是多少？
答案：
- 我們經(jīng)常進行模擬實驗隘弊，以檢驗我們所看到的數(shù)據(jù)是否與最簡單的四箱模型一致哈踱，其中每個箱的概率是1/4。在某種意義上梨熙，它是稻草人(沒有發(fā)生什么有趣的事情)开镣。我們將在第2章中看到更多的例子。在這里咽扇，我們使用一些R命令來生成這樣的計數(shù)矢量邪财。首先，假設我們有四種不同的质欲、同樣可能的類型的8個字符：
> pvec = rep(1/4, 4)
> pvec
[1] 0.25 0.25 0.25 0.25
> t(rmultinom(1, prob = pvec, size = 8))
[,1] [,2] [,3] [,4]
[1,] 1 0 4 3
rmultinom
## Usage（用法)
rmultinom(n, size, prob) #【產(chǎn)生隨機樣本】
dmultinom(x, size = NULL, prob, log = FALSE) # 【密度函數(shù)】
## 參數(shù)說明
x：長度為k的整數(shù)向量树埠，從0到size。
n：要抽取的隨機變量的個數(shù)
size：整數(shù)嘶伟，指定多維正態(tài)試驗中被放入K個盒子對象的總數(shù)量怎憋，對于dmultnom，它默認為sum(x)
prob：長度為K的非負數(shù)值向量九昧，指定K類的概率绊袋，內(nèi)部規(guī)范為1，無限和缺失值是不允許的
log：邏輯值铸鹰，如果為真癌别，計算的是對數(shù)概率
## 比如
# https://blog.csdn.net/zhaozhn5/article/details/77933670
rmultinom(n, size, prob)
#拋 10 次骰子為一次實驗，做 1000 次實驗蹋笼。則 n=1000展姐，size=10。
#prob為每個獨立結果出現(xiàn)的概率剖毯，其總和為1诞仓。
#結果為 k×n 的矩陣，k即length(prob)
? 問題 1.8：t函數(shù)表示轉置
? 問題 1.9：rmultinom(n = 8, prob = pvec, size = 1)和 rmultinom(n = 1, prob = pvec, size = 8)的區(qū)別? 提示：1.3.1 和 1.3.2.
> rmultinom(n = 8, prob = pvec, size = 1)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0 0 0 0 0 0 0 0
[2,] 1 0 0 0 0 0 0 0
[3,] 0 1 0 1 1 1 0 1
[4,] 0 0 1 0 0 0 1 0
> rmultinom(n = 1, prob = pvec, size = 8)
[,1]
[1,] 3
[2,] 2
[3,] 1
[4,] 2
n = 8, size = 1表示丟1一個球為一次實驗速兔，做8次實驗墅拭，得到的矩陣為length(prob) * n(即4 * 8)的矩陣，每一列代表四種情況中哪一種會成功涣狗，標記為1。四行對應的是四種情況太示，八列代表的是八次試驗香拉。n = 1, size = 8表示丟八個球為，做一次實驗镀迂，得到一個 4 * 1的矩陣。由于只做一次實驗唤蔗，所以只有一列探遵，每一行表示每種可能出現(xiàn)的次數(shù)，最后總次數(shù)等于size妓柜。
1.4.1 Simulating for power
讓我們來看一個例子箱季，用蒙特卡羅來表示多項式，這與科學家們在計劃他們的實驗時經(jīng)常要解決的問題有關：我需要多大的樣本量呢棍掐？藏雏。
權重一詞在統(tǒng)計學中有特殊的含義。它是發(fā)現(xiàn)某事件是否發(fā)生的概率作煌，也被稱為真正的陽性率掘殴。
按照慣例，實驗者在計劃實驗時的目標是80%(或更多)的權重粟誓。這意味著奏寨，如果相同的實驗運行多次，即使本應該發(fā)生鹰服，大約20%的試驗次數(shù)將不能產(chǎn)生顯著的結果服爷。
零假設H0，我們收集的DNA數(shù)據(jù)來自一個公平的過程获诈，在這個過程中仍源，4個堿基中的每一個都有相同的可能(pA，pC舔涎，pG笼踩，pT)=(0.25，0.25亡嫌，0.25嚎于，0.25)。這意味著：沒有什么有意思的事情發(fā)生挟冠。這就是我們試圖反駁(用統(tǒng)計學家的術語來說于购，是“拒絕”)的稻草人，所以零假設應該是這樣的知染，偏離它是很有趣的肋僧。。
正如您看到的，通過對8個字符和4個等概率結果(由大小相等的方框表示)運行R命令嫌吠，我們并不總是在每個box中得到2個止潘。這是不可能的，從只看8個字符辫诅，核苷酸是否來自一個公平的過程凭戴。
讓我們來確定，通過觀察一個長度為n=20的序列炕矮，我們是否能夠檢測出最初的核苷酸分布是否公平么夫，或者它是否來自其他的(“替代”)過程。
我們使用
rmulom函數(shù)從零假設中生成了1000個模擬肤视。我們只顯示前11列以節(jié)省空間档痪。
> obsunder0 = rmultinom(1000, prob = pvec, size = 20)
> dim(obsunder0)
[1] 4 1000
> obsunder0[, 1:11]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[1,] 4 7 7 8 7 5 7 9 1 6 3
[2,] 3 4 7 4 8 7 4 3 6 3 6
[3,] 9 2 2 4 3 2 5 2 6 3 3
[4,] 4 7 4 4 2 6 4 6 7 8 8
- 矩陣
obsunder0中的每一列都是一個模擬實例。您可以看到钢颂，框中的數(shù)字差異很大：有些大到8钞它，而期望值是5 =拜银。
創(chuàng)建測試
-
記资獗蕖：我們知道這些值來自一個公平的過程。顯然尼桶，僅僅知道流程的期望值是不夠的操灿。我們還需要一種可變性的度量，使我們能夠描述預期的可變性有多大泵督，以及有多大就是太大了趾盐。我們使用以下統(tǒng)計量作為度量，該統(tǒng)計量是觀測值與期望值相對于期望值的差值的平方和小腊。因此救鲤，對于每個實例，
- 生成的數(shù)據(jù)的前三列與我們預期的有多大不同秩冈？我們得到：
> expected0 = pvec * 20
> sum((obsunder0[, 1] - expected0)^2 / expected0)
[1] 4.4
> sum((obsunder0[, 2] - expected0)^2 / expected0)
[1] 3.6
> sum((obsunder0[, 3] - expected0)^2 / expected0)
[1] 3.6
- 度量的值可能會不同-您可以查看3個以上的列本缠，我們將了解如何研究所有1000個列。為了避免重復鍵入入问，我們將
stat公式(表達式(1.1))封裝在一個函數(shù)中：
> stat = function(obsvd, exptd = 20 * pvec) {
+ sum((obsvd - exptd)^2 / exptd)
+ }
> stat(obsunder0[, 1])
[1] 4.4
- 為了更全面地了解這種變化丹锹，我們計算了所有1000個實例的度量，并將這些值存儲在一個我們稱為 S0 的向量中：它包含在
H0下生成的值芬失。我們可以考慮圖1.10所示的S0值的直方圖楣黍，這是對空分布的估計。
> S0 = apply(obsunder0, 2, stat)
> summary(S0
+ )
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 1.200 2.800 3.115 4.400 17.200
> hist(S0, breaks = 25, col = "lavender", main = "")
-
summary函數(shù)向我們展示了S0具有不同值的分布棱烂。例如租漂，從模擬數(shù)據(jù)中，我們可以近似地得到95%的分位數(shù)(這個值將較小的95%的值與5%的最大值分開)。
> q95 = quantile(S0, probs = 0.95)
> q95
95%
7.6
- 因此窜锯，我們看到5%的
S0值大于7.6张肾。我們將提出這作為我們測試數(shù)據(jù)的臨界值，并將拒絕這樣的假設锚扎，即如果加權平方和統(tǒng)計數(shù)據(jù)大于7.6吞瞪，數(shù)據(jù)來自一個公平的過程，具有同樣可能的核苷酸驾孔。
Determining our test’s power
- 我們必須計算我們的檢驗-基于加權平方和差-將檢測到數(shù)據(jù)實際上不是來自零假設的概率芍秆。通過仿真計算了拒絕率。我們從一個由
pvecA參數(shù)化的可選進程生成1000個模擬實例翠勉。
> pvecA = c(3/8, 1/4, 3/12, 1/8) ## 概率
> observed = rmultinom(1000, prob = pvecA, size = 20) ## 實驗一千次妖啥，每次20個球
> dim(observed)
[1] 4 1000
> observed[, 1:7]
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 10 10 8 5 9 12 7
[2,] 2 3 3 5 3 3 6
[3,] 7 4 6 7 4 3 4
[4,] 1 3 3 3 4 2 3
> apply(observed, 1, mean)
[1] 7.388 5.078 5.009 2.525
> expectedA = pvecA * 20
> expectedA
[1] 7.5 5.0 5.0 2.5
與零假設的模擬一樣，觀測值也有很大的差異对碌。問題是：我們的測試多久(在1000個實例中)檢測數(shù)據(jù)是否偏離NULL荆虱？
測試不會拒絕第一個觀察結果(10、3朽们、5怀读、2)，因為統(tǒng)計數(shù)據(jù)的值在95%以內(nèi)骑脱。
> stat(observed[, 1])
[1] 10.8
> S1 = apply(observed, 2, stat)
> q95
95%
7.6
> sum(S1 > q95)
[1] 187
> power = mean(S1 > q95)
> power
[1] 0.187
在1000個模擬中菜枷，測試確定了187個來自
alternative distribution。因此叁丧，我們計算出概率P(拒絕H0|HA)為0.187啤誊。With a sequence length of n=20 we have a power of about 20% to detect the difference between the fair generating process and our alternative在實踐中，正如我們提到的拥娄，一個可接受的加權值（power）是0.8或更多蚊锹。重復模擬實驗，并建議一個新的序列長度n稚瘾，以確保加權值（power）是可接受的牡昆。
經(jīng)典數(shù)據(jù)的經(jīng)典統(tǒng)計
我們不需要用蒙特卡羅模擬數(shù)據(jù)來計算第95個百分位數(shù)，有足夠的理論可以幫助我們進行計算孟抗。
我們的統(tǒng)計數(shù)據(jù)實際上有一個眾所周知的分布迁杨，稱為卡方分布(具有3個自由度)，寫作
凄硼。
這里用markdwon表示為 $χ_2 ^3$
- 在第2章中铅协，我們將看到如何使用Q-Q圖比較分布(參見圖2.8)。我們本可以使用更標準的測試摊沉，而不是運行手工模擬狐史。然而，我們所學到的過程擴展到許多情況下，卡方分布并不適用骏全。例如苍柏，當某些盒子的
概率極低，并且它們的計數(shù)大多為零時姜贡。
1.5试吁、本章節(jié)總結
我們使用數(shù)學公式和R來計算各種離散事件的概率，我們可以用一些基本的分布來建模：
Bernoulli分布是我們最基本的建模-它用于表示單個二進制試驗(即結果只有無兩種情況楼咳，用1和0表示)熄捍，如硬幣翻轉。我們可以將結果編碼為0和1母怜。我們稱p為成功概率(成功輸出1)余耽。在n個二元試驗中，對
數(shù)字1s的個數(shù)采用二項分布(The binomial distribution)苹熏，并利用R函數(shù)dbinom生成k次成功的概率碟贾。我們還看到，我們可以使用函數(shù)rbinom來模擬n次試驗二項式轨域。泊松分布(The Poisson distribution)最適合于p較小時(數(shù)字1s很少)的情形袱耽。它只有一個參數(shù)
λ，當p較小時疙挺，λ=np的泊松分布近似于(n扛邑，p)的二項分布怜浅。我們使用泊松分布來模擬在沿序列檢測表位的試驗中隨機出現(xiàn)的假陽性數(shù)铐然，假設每個位置的假陽性率很小。只要我們知道所有的參數(shù), 我們就能看到了這樣一個參數(shù)模型如何使我們能夠計算極端事件的概率恶座。多維正態(tài)分布（The multinomial distribution）用于具有
兩個以上可能結果或級別的離散事件搀暑。POWER示例向我們展示了如何使用蒙特卡洛(MonteCarlo)模擬來決定我們需要收集多少數(shù)據(jù)，如果我們想要檢驗具有等概率的多項模型是否與數(shù)據(jù)一致的話跨琳。我們使用概率分布和概率模型來評估關于我們的數(shù)據(jù)是如何生成的假設自点，通過對生成模型的假設。在給定假設的情況下脉让，我們將看到數(shù)據(jù)的概率稱為p值桂敛。這和假設是真的概率是不一樣的！
課后閱讀
The elementary book by Freedman, Pisani, and Purves (1997) provides the best introduction to probability through the type of box models we mention here.
The book by Durbin et al. (1998) covers many useful probability distributions and provides in its appendices a more complete view of the theoretical background in probability theory and its applications to sequences in biology.
Monte Carlo methods are used extensively in modern statistics. Robert and Casella (2009) provides an introduction to these methods using R.
Chapter 6 will cover the subject of hypothesis testing. We also suggest Rice (2006) for more advanced material useful for the type of more advanced probability distributions, beta, gamma, exponentials we often use in data analyses.
參考資源
Freedman, David, Robert Pisani, and Roger Purves. 1997. Statistics. New York, NY: WW Norton.
Durbin, Richard, Sean Eddy, Anders Krogh, and Graeme Mitchison. 1998. Biological Sequence Analysis. Cambridge University Press.
Robert, Christian, and George Casella. 2009. Introducing Monte Carlo Methods with R. Springer Science & Business Media.
Rice, John. 2006. Mathematical Statistics and Data Analysis. Cengage Learning.
