用R語言求概率分布_r語言 概率分布圖

R語言一個很方便的用處是提供了一套完整的統計表集合。函數可以對累積分布函數P(X≤x)，概率密度函數，分位函數（對給定的q，求滿足P(X≤x) > q的最小x）求值，并根據分布進行模擬。

在R中，根據某種分布生成隨機序列的函數如下：

在統計學中，產生隨機數據是很有用的，R可以產生多種不同分布下的隨機數序列。這些分布函數的形式為rfunc(n,p1,p2,…)，其中func指概率分布函數，n為生成數據的個數，p1, p2, . . .是分布的參數數值。上面的表給出了每個分布的詳情和可能的缺省值（如果沒有給出缺省值，則意味著用戶必須指定參數）。數據分析培訓

例：用0~1之間的均勻分布產生10個隨機點

> runif(10)

[1] 0.961465376 0.007521925 0.193619234 0.137027246 0.739370654 0.072907082

[7] 0.674551635 0.650777811 0.984664183 0.796723066

大多數這種統計函數都有相似的形式，只需用d、p或者q去替代r，比如密度函數(dfunc(x, …))，累計概率密度函數（也即分布函數）(pfunc(x,…))和分位數函數(qfunc(p, …)，0<p<1)。最后兩個函數序列可以用來求統計假設檢驗中P值或臨界值。例如，顯著性水平為5%的正態分布的雙側臨界值是：

> qnorm(0.025)

[1] -1.959964

> qnorm(0.975)

[1] 1.959964

對于同一個檢驗的單側臨界值，根據備擇假設的形式使用qnorm(0.05)或1 – qnorm(0.95)。

下面是一些用R語言求解概率問題的例子：

1.       某人進行射擊，每次擊中目標的命中率為0.02，獨立射擊400次，求至少擊中兩次的概率。

解：400重伯努利試驗，用二項分布求解。

P{X = k} = C400k * (0.02)^k * (0.0=98)^(400-k)

P{X≥2} = 1 – P{X = 0} – P{X = 1}

> 1 – sum(pbinom(0:1, 400, 0.02))

[1] 0.9968561

結論：決不能輕視小概率事情，在多次重復試驗的情況下，這一事件的發生幾乎是肯定的。

2.       設X服從平均值為1，標準差為2的正態分布（高斯分布），即X ~ N(1, 4)，求P{0<X≤1.6}

解：這里X是一個連續型隨機變量。求X在某段區間上的概率，用X的分布函數在區間兩端的值的差。

方法一：P{0<X≤1.6} = P{X≤1.6} – P{X≤0} = F(1.6) – F(0)

> pnorm(1.6, 1, 2) – pnorm(0, 1, 2)

[1] 0.3093739

方法二：轉化為標準正態分布。P{x1 < X ≤x2}=P{(x1-μ)/σ < (X-μ)/σ≤(x1-μ)/σ}=φ((x2-μ)/σ) –φ((x1-μ)/σ)

即P{0<X≤1.6}=φ((1.6-1)/2) –φ((0-1)/2)

> pnorm((1.6-1)/2) – pnorm((0-1)/2)   #pnorm函數的缺省參數mean=0,sd=1，即默認標準正態分布

[1] 0.3093739

知識點：設X是一個隨機變量，x是任意實數，函數F(x)=P{X≤x}稱為X的分布函數。

對于任意實數x1,x2(x1＜x2),有P{x1＜X≤x2}=P{X≤x2}-P{X≤x1}=F(x2)-F(x1),

因此，若已知X的分布函數，就可以知道X落在任一區間(x1,x2]上的概率，在這個意義上說，分布函數完整地描述了隨機變量的統計規律性。

分布函數是一個普遍的函數，正是通過它，我們將能用數學分析的方法來研究隨機變量。

如果將X看成是數軸上的隨機點的坐標，那么，分布函數F(x)在x處的函數值就表示X落在區間(-∞,x]上的概率。

3.       求標準正態分布的上α分位點。

知識點：設X~N(0,1)，若Zα滿足條件 P(X>Zα)=α，0<α<1，則稱Zα為標準正態分布的上α分位點.

注意上α分位點和R語言中分位函數（對給定的q，求滿足P(X≤x) > q的最小x）之間的關系。

解：下面給出α=0.001、α=0.005、α=0.01、α=0.025時的上α分位點Zα的值。

> exp <- expression_r(qnorm(1 – alpha))

> alpha = 0.001

> eval_r(exp)

[1] 3.090232

> alpha = 0.005

> eval_r(exp)

[1] 2.575829

> alpha = 0.01

> eval_r(exp)

[1] 2.326348

> alpha = 0.025

> eval_r(exp)

[1] 1.959964