摘要: 本文由digging4發表於：http://www.cnblogs.com/digging4/p/5054594.html
統計建模與R軟件-第四章 參數估計
4.1設總體的分布密度為
$$f(x;\alpha)=\begin{cases}(\alpha+1)x^\alpha&,x\in(0,1)\0&,other\end{cases}$$
\(X_1,X_2,\cdots,X_n\)為其樣本，求參數\(\alpha\)的矩估計量\(\widehat{\alpha_1}\)和極大似然估計量\(\widehat{\alpha_2}\).現測得樣本觀測值為
$$0.1,0.2,0.9,0.8,0.7,0.7$$
求參數\(\alpha\)的估計值。
矩估計：用樣本的一階原點矩估計總體均值；用樣本的二階中心矩估計總體方差。
如果總體的分布已知，那么總體的均值和方差就可以用分布中的參數表示，再等於樣本的一階原點矩和二階中心矩，可以計算出總體分布中的參數。
矩估計優點：在其能用的情況下，計算往往簡單
矩估計缺點：相對其他估計方法，如極大似然法，其效率往往較低。
已知分布密度，求隨機變量的期望（均值）如下，該期望值等於樣本的均值：
\[\int_{-\infty}^{\infty}xf(x;\alpha )dx=\int_{0}^{1}x(\alpha+1)x^{\alpha}dx=\frac{\alpha+1}{\alpha+2}x^{\alpha+2}\mid_0^1=\frac{\alpha+1}{\alpha+2}=mean(x) \]
x <- c(0.1, 0.2, 0.9, 0.8, 0.7, 0.7)
(2 * mean(x) - 1)/(1 - mean(x))
## [1] 0.3077
極大似然估計：我們所估計的模型參數，要使得產生這個給定樣本的可能性最大。似然函數如下：
\[L(x;\alpha )=\prod_{i=1}^{n}f(x;\alpha)=(\alpha+1)^n\prod_{i=1}^{n}x_i^\alpha \]
取對數得：
\[lnL(x;\alpha )=nln(\alpha+1)+\alpha \sum_{i=1}^{n}lnx_i \]
在對對數似然函數求倒，可得到
\[\frac{\partial lnL(x;\alpha )}{\partial \alpha}=\frac{n}{\alpha+1}+\sum_{i=1}^{n}lnx_i \]
即為求上述偏導數等於0的\(\alpha\)值，此例中\(n=3\)
x <- c(0.1, 0.2, 0.9, 0.8, 0.7, 0.7)
f <- function(a) 6/(a + 1) + sum(log(x))
uniroot(f, c(0, 1))
## $root
## [1] 0.2112
##
## $f.root
## [1] -3.845e-05
##
## $iter
## [1] 5
##
## $estim.prec
## [1] 6.104e-05
root為估計值，iter為迭代次數，optimize函數也可以用來求解方程。
4.2設元件無故障工作時間\(X\)具有指數分布，取1000個元件工作時間的記錄數據，經分組后得到它的頻數分布為
組中值\(x_i\) 5 15 25 35 45 55 65
頻數\(v_i\) 365 245 150 100 70 45 25
如果各組中數據都取為組中值，試用極大似然估計求\(\lambda\)的點估計。
指數分布\(\lambda\)的估計量為
\[\hat{\lambda}=\frac{n}{\sum_{i=1}^{n}x_i} \]
x <- c(rep(5, 365), rep(15, 245), rep(25, 150), rep(35, 100), rep(45, 70), rep(55,
45), rep(65, 25))
1000/sum(x)
## [1] 0.05
4.3為檢驗某自來水消毒設備的效果，現從消毒后的水中隨機抽取50升，化驗每升水中大腸桿菌的個數（假設一升水中大腸桿菌個數服從\(Poisson\)分布），其化驗結果如下：
大腸桿菌數/升 0 1 2 3 4 5 6
升數 17 20 10 2 1 0 0
試問平均每升水中大腸桿菌個數為多少時，才能使上述情況的概率為最大？
本題實際上是求泊松分布的\(\lambda\)估計，泊松分布的均值和方差均為\(\lambda\)，使用點估計即可
x <- c(rep(0, 17), rep(1, 20), rep(2, 10), rep(3, 2), rep(4, 1), rep(5, 0),
rep(6, 0))
## 得到$\lambda$的估計值
mean(x)
## [1] 1
4.4利用R軟件中的nlm()函數求解無約束優化問題
$$minf(x)=(-13+x_1+((5-x_2)x_2-2)x_2)^2+(-29+x_1+((x_2+1)x_2-14)x_2)^2$$
取初始點\(x^{(0)}=(0.5,-2)^T\).
# 目標函數
obj <- function(x) {
f <- c(-13 + x[1] + ((5 - x[2]) * x[2] - 2) * x[2], -29 + x[1] + ((x[2] +
1) * x[2] - 14) * x[2])
sum(f^2)
}
x0 <- c(0.5, -2)
nlm(obj, x0)
## $minimum
## [1] 48.98
##
## $estimate
## [1] 11.4128 -0.8968
##
## $gradient
## [1] 1.415e-08 -1.435e-07
##
## $code
## [1] 1
##
## $iterations
## [1] 16
# 最優目標值為 $minimum 48.98
4.5正常人的脈搏平均每分鍾72次，某醫生測得10例四乙基鉛中毒患者的脈搏數（次/分）如下：
54 67 68 78 70 66 67 70 65 69
已知人的脈搏次數服從正態分布，試計算這10名患者平均脈搏次數的點估計和\(95%\)的區間估計。並做單側區間估計，試分析這10名患者的平均脈搏次數是否低於正常人的平均脈搏次數。
x <- c(54, 67, 68, 78, 70, 66, 67, 70, 65, 69)
# 由於人的脈搏次數服從正態分布，因此平均次數的點估計為樣本均值
mean(x)
## [1] 67.4
# t.test 可以進行區間估計
t.test(x)
##
## One Sample t-test
##
## data: x
## t = 35.95, df = 9, p-value = 4.938e-11
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 63.16 71.64
## sample estimates:
## mean of x
## 67.4
#
# 得到95%置信度下的區間估計為[63.16,71.64]，可認為患者低於常人每分鍾72次的正常水平。
4.6甲、乙兩種稻種分布播種在10塊試驗田中，每塊試驗田甲、乙稻種各種一半，假設兩稻種產量\(X,Y\)均服從正態分布，且方差相等，收獲后10塊試驗田的產量如下所示（單位：千克）
甲種 140 137 136 140 145 148 140 135 144 141
乙種 135 118 115 140 128 131 130 115 131 125
求出兩稻種產量的期望差\(\mu_1-\mu_2\)的置信區間（\(\alpha=0.05\)）.
# 兩個正態總體，方差未知但相等情況下的均值差估計
x <- c(140, 137, 136, 140, 145, 148, 140, 135, 144, 141)
y <- c(135, 118, 115, 140, 128, 131, 130, 115, 131, 125)
t.test(x, y, var.equal = TRUE)
##
## Two Sample t-test
##
## data: x and y
## t = 4.629, df = 18, p-value = 0.0002087
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 7.536 20.064
## sample estimates:
## mean of x mean of y
## 140.6 126.8
# 得到95%置信度下的區間估計為[7.36,20.24]
4.7甲、乙兩組生產同種導線，現從甲組生產的導線中隨機抽取4根，從乙組生產的導線中隨機抽取5根，它們的電阻值（單位：\(\Omega\)）
甲組 0.143 0.142 0.143 0.137
乙組 0.140 0.142 0.136 0.138 0.140
假設兩組電阻值分別服從正態分布\(N(\mu_1,\alpha^2)\)和\(N(\mu_2,\alpha^2)\)，\(\alpha^2\)未知，試求\(\mu_1-\mu_2\)的置信系數為0.95的區間估計。
# 兩個正態總體，方差未知但相等情況下的均值差估計
x <- c(0.143, 0.142, 0.143, 0.137)
y <- c(0.14, 0.142, 0.136, 0.138, 0.14)
t.test(x, y)
##
## Welch Two Sample t-test
##
## data: x and y
## t = 1.164, df = 5.701, p-value = 0.2909
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.002315 0.006415
## sample estimates:
## mean of x mean of y
## 0.1412 0.1392
4.8對習題4.6中甲乙兩種稻種的數據作方差比的區間估計，並用其估計值來判定兩數據是否等方差，若兩數據方差不相等，試重新計算兩稻種產量的期望差\(\mu_1-\mu_2\)的置信區間（\(\alpha=0.05\)）。
# var.test() 雙樣本方差比的區間估計
x <- c(140, 137, 136, 140, 145, 148, 140, 135, 144, 141)
y <- c(135, 118, 115, 140, 128, 131, 130, 115, 131, 125)
var.test(x, y)
##
## F test to compare two variances
##
## data: x and y
## F = 0.2353, num df = 9, denom df = 9, p-value = 0.04229
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.05845 0.94744
## sample estimates:
## ratio of variances
## 0.2353
#
# 得到方差比的區間估計為[0.05845,0.94744]，區間中不包含1，所以可認為兩個樣本的方差不等
# 重新計算兩樣本的期望差
t.test(x, y, var.equal = FALSE)
##
## Welch Two Sample t-test
##
## data: x and y
## t = 4.629, df = 13.01, p-value = 0.0004712
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 7.36 20.24
## sample estimates:
## mean of x mean of y
## 140.6 126.8
4.9設電話總機在某段時間內接到的呼喚的次數服從參數未知的\(Poisson\)分布\(P(\lambda)\),現搜集了42個數據
接到呼喚的次數 0 1 2 3 4 5 6
出現的頻數 7 10 12 8 3 2 0
試求出平均呼喚次數\(\lambda\)的估計值和它的置信系數為0.95的置信區間。
# 由於分布不是正態分布，需要樣本量比較大，利用中心極限定理進行分析
# 按書P190的公式進行計算
x <- c(rep(0, 7), rep(1, 10), rep(2, 12), rep(3, 8), rep(4, 3), rep(5, 2))
n <- length(x)
tmp <- sd(x)/sqrt(n) * qnorm(1 - 0.05/2)
mean(x) - tmp #估計區間下限
## [1] 1.494
mean(x) + tmp #估計區間上限
## [1] 2.315
4.10已知某種燈泡壽命服從正態分布，在某星期所生產的該燈泡中隨機抽取10只，測得其壽命（單位：小時）為：1067,919,1196,785,1126,936,918,1156,920,948，求燈泡壽命平均值的置信度為0.95的單側置信下限。
# t.test()可完成單側區間估計
x <- c(1067, 919, 1196, 785, 1126, 936, 918, 1156, 920, 948)
t.test(x, alternative = "greater")
##
## One Sample t-test
##
## data: x
## t = 23.97, df = 9, p-value = 9.148e-10
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
## 920.8 Inf
## sample estimates:
## mean of x
## 997.1
# 單側區間估計的下限為920.8