参数估计(Parameter Estimation)是指用样本指标(称为统计量)估计总体指标(称为参数)。参数估计有点估计(point estimation)和区间估计(interval estimation)两种。

点估计

设总体\(X\)的分布函数\(F(x;\theta)\)形式已知，其中\(\theta\)是待估计的参数，点估计就是利用样本(\(x_{1},x_{2},...,x_{n}\))，构造一个统计量\(\hat{\theta }=\hat{\theta}(x_{1},x_{2},...,x_{n})\)来估计\(\theta\)，称\(\hat{\theta}(x_{1},x_{2},...,x_{n})\)为\(\theta\)的点估计量，它是一个随机变量。将样本观测值(\(x_{1},x_{2},...,x_{n}\))代入估计量\(\hat{\theta}(x_{1},x_{2},...,x_{n})\)，就得到它的一个具体数值\(\hat{\theta}(x_{1},x_{2},...,x_{n})\)，这个数值成为\(\theta\)的点估计值。 点估计是依据样本估计总体分布中所含的未知参数或未知参数的函数。通常它是总体的某个特征值，如数学期望、方差和相关系数等。点估计问题就是要构造一个只依赖于样本的量，作为未知参数或未知参数的函数的估计值。 构造点估计常用的方法是：

矩估计法

设\((x_{1},x_{2},...,x_{n})\)是来自总体\(X\)的一个样本，根据大数定律，对任意\(\varepsilon >0\),有 \[\lim_{n\rightarrow \infty}P\left \{\right |\bar{X}-E(X)|\geq \varepsilon\}=0\] 并且对于任何\(k\)，只要\(E(X^{k})\)存在，同样有 \[\lim_{n\rightarrow \infty}P\left \{ \right|\frac{1}{n}\sum_{i=1}^{n}X_{i}^{k}-E(X^{k})|\geq \varepsilon \}=0,k=1,2,...\] 因此用样本矩估计总体矩，从而得到总体分布中参数的一种估计。如用样本均值估计总体均值。矩法的优点是简单易行，并不需要事先知道总体是什么分布，缺点是当总体类型已知时，没有充分利用分布提供的信息，且矩估计量不具有唯一性。

例1 设某药厂一天中发生着火现象的次数X服从参数为\(\lambda\)的Poisson分布，\(\lambda\)未知，有以下样本值，试用矩法估计参数\(\lambda\)。

着火的次数 0 1 2 3 4 5 6
发生k次着火的天数\(n_{k}\) 75 90 54 22 6 2 1

解 \(EX=\lambda,A_{1}=\frac{1}{n}\sum_{i=1}^{n}X_{i}=\bar{X}\),

令\(\bar{X}=\lambda\)则\(\bar{\lambda }=\bar{x}=\frac{1}{250}(0\times 57+1\times 90+...+6\times 1)=1.22\)，

所以\(\bar{X}=\lambda\),估计值\(\hat{\lambda }=1.22\)

例2 正态分布N（0,1）的矩估计

x<-rnorm(100) #产生N（0,1）的100个随机数
mu<-mean(x)   #对N(mu,sigma)中的mu做矩估计
sigma<-var(x) #这里的var并不是样本方差的计算函数，而是修正的样本方差，其实也就是x的总体方差
mu
## [1] -0.08923805
sigma
## [1] 1.286883

极大似然估计法（MLE）

它是建立在极大似然原理的基础上的一个统计方法，极大似然原理的直观想法是：一个随机试验如有若干个可能的结果A，B，C，…。若在一次试验中，结果A出现，则一般认为试验条件对A出现有利，也即A出现的概率很大。当从模型总体随机抽取n组样本观测值后，最合理的参数估计量应该是使得从模型中抽取该n组样本观测值的概率最大。在任一次随机抽取中，样本观测值都以一定的概率出现。如果已经知道总体的参数，当然由变量的频率函数可以计算其概率。如果只知道总体服从某种分布，但不知道其分布参数，通过随机样本可以求出总体的参数估计。

例3 对MASS包中的geyser数据，该数据采集自美国黄石公园内的一个名叫Old Faithful 的喷泉。“waiting”就是喷泉两次喷发的间隔时间，“duration”当然就是指每次喷发的持续时间。在这里，我们只用到“waiting”数据

panderOptions('table.split.table', Inf)
pander(head(geyser))
waiting duration
80 4.017
71 2.15
57 4
80 4
75 4
77 2 
hist(geyser$waiting,freq = F) #从图中可以看出，其分布是两个正态分布的混合。

用如下的分布函数来描述该数据 \[f(x)=pN(x_i;\mu_1,\sigma_1)+(1-p)N(x_i;\mu_2,\sigma_2)\] 该函数中有5个参数\(p、\mu_1、\sigma_1、\mu_2、\sigma_2\)需要确定。上述分布函数的对数极大似然函数为： \[l=\sum_{i=1}^n\log \{pN(x_i;\mu_1,\sigma_1)+(1-p)N(x_i;\mu_2,\sigma_2)\}\] 在R中定义对数似然函数

LL<-function(params,data) #定义log-likelihood函数,参数"params"是一个向量，依次包含了五个参数：p,mu1,sigma1,#mu2,sigma2.#参数"data"，是观测数据。
{
t1<-dnorm(data,params[2],params[3])  #这里的dnorm()函数是用来生成正态密度函数的。
t2<-dnorm(data,params[4],params[5])
f<-params[1]*t1+(1-params[1])*t2
ll<-sum(log(f))  #混合密度函数,log-likelihood函数
return(-ll) #nlminb()函数是最小化一个函数的值，但我们是要最大化log-likeilhood函数，所以需要在“ll”前加个“-”号。
}
#参数估计
hist(geyser$waiting,freq = F)
lines(density(geyser$waiting)) #初始值为p=0.5,mu1=50,sigma1=10,mu2=80,sigma2=10

geyser.res<-nlminb(c(0.5,50,10,80,10),LL,data=geyser$waiting,lower=c(0.0001,-Inf,0.0001,-Inf,-Inf,0.0001),upper=c(0.9999,Inf,Inf,Inf,Inf))

估计结果

geyser.res$par #查看拟合的参数
## [1]  0.3075937 54.2026518  4.9520026 80.3603085  7.5076330
X<-seq(40,120,length=100)
p<-geyser.res$par[1]
mu1<-geyser.res$par[2]
sig1<-geyser.res$par[3]
mu2<-geyser.res$par[4]
sig2<-geyser.res$par[5]
f<-p*dnorm(X,mu1,sig1)+(1-p)*dnorm(X,mu2,sig2) #将估计的参数函数代入原密度函数。
hist(geyser$waiting,probability=T,col=0,ylab="Density", #作出数据的直方图
     ylim=c(0,0.04),xlab="Eruption waiting times")
lines(X,f) #画出拟合的曲线

最小二乘法

当从模型总体随机抽取n组样本观测值后，最合理的参数估计量应该使得模型能最好地拟合样本数据，即实际值与估计值的距离最小，主要用于线性统计模型中的参数估计问题。

例4 用最小二乘法估计线性回归模型

x <- c(5.05, 6.75, 3.21, 2.66)
y <- c(1.65, 26.5, -5.93, 7.96)
lsfit(x, y)$coefficients #或者lm(y ~ x)$coefficients
##  Intercept          X 
## -16.281128   5.393577
plot(x, y)
abline(lsfit(x, y)$coefficients, col="red")

EM算法

EM算法是一种在观测到数据后，用迭代法估计未知参数的方法。可以证明EM算法得到的序列是稳定单调递增的。这种算法对于截尾数据或参数中有一些不感兴趣的参数时特别有效。EM算法的步骤为：E-step（求期望）：利用对隐藏变量的现有估计值，计算其最大似然估计值。M-step（求极值）：最大化在 E 步上求得的最大似然值来计算参数的值，重复以上两步，直至收敛即可得到theta的MLE。可以看到对于一个参数的情况，EM仅仅只是求解MLE的一个迭代算法。

sim.x <- c()  
sim.y <- c()  

# 用循环产生2000个点
for (i in 1:2000) {
    # first draw to determine which normal distribution is used
  first.draw = rmultinom(1, 1, c(0.1, 0.2, 0.7))[, 1]
  y = which(first.draw == 1)
  sim.y[i] = y
  
  # second draw to generate X from corresponding distribution
  if (y == 1) {
    x = rnorm(1, mean = 0, sd = 1)
    sim.x[i] = x
  }
  if (y == 2) {
    x = rnorm(1, mean = 10, sd = 5)
    sim.x[i] = x
  }
  if (y == 3) {
    x = rnorm(1, mean = -10, sd = 1)
    sim.x[i] = x
  }
}
plot(density(sim.x), main = "Density plot of sim.x")

mix.model <- normalmixEM(sim.x, lambda = c(0.3, 0.3, 0.4), mu = c(-20, 0, 20), sigma = c(1, 1, 1), k = 3)
## number of iterations= 35
summary(mix.model)
## summary of normalmixEM object:
##           comp 1     comp 2   comp 3
## lambda   0.71040  0.0987118 0.190888
## mu     -10.02994 -0.0328829 9.778046
## sigma    0.97781  1.1686164 4.769408
## loglik at estimate:  -4934.399
plot(mix.model, which = 2, density = TRUE)

Bootstrap法

以原始数据为基础的模拟抽样统计推断法,可用于研究一组数据的某统计量的分布特征,特别适用于那些难以用常规方法导出对参数的区间估计、假设检验等问题。“Bootstrap”的基本思想是:在原始数据的围内作有放回的再抽样,样本容量仍为n,原始数据中每个观察单位每次被抽到的概率相等,为1，…,n,所得样本称为bootstrap样本。于是可得到参数Η的一个估计值Η(b),这样重复若干次,记为B。设B=1000,就得到该参数的1000个估计值,则参数Η的标准误的bootstrap估计。简而言之就是：就是从样本中重复抽样。

gauss<-rnorm(1000,4,10)
boot<-0
for(i in 1:1000){boot[i]=mean(sample(gauss,replace=T))}
summary(boot)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.477   3.493   3.726   3.728   3.951   4.702
summary(gauss)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -32.610  -3.252   4.181   3.744  10.680  35.870
sd(boot)
## [1] 0.3392106

区间估计

由于点估计不能说明估计值与真实值的偏差到底有多大，也不能说明这个估计有多大的可行度，这些问题需要区间估计予以解决。区间估计是依据抽取的样本，根据一定的正确度与精确度的要求，构造出适当的区间，作为总体分布的未知参数或参数的函数的真值所在范围的估计。求置信区间常用的三种方法：1.利用已知的抽样分布.利用区间估计与假设检验的联系。3.利用大样本理论。

设总体X的分布中含有未知参数\(\theta,\alpha\)是任意给定的正数\((0< \alpha < 1)\)，如果能从样本除服确定出两个统计量\(\hat{\theta }_{1}(x_{1},x_{2},...,x_{n}),\hat{\theta }_{2}(x_{1},x_{2},...,x_{n})\)，使得 \[P\left \{ \hat{\theta }_{1}<\theta <\hat{\theta }_{2} \right \}=1-\alpha \] 成立，我们称\(1-\alpha\)为置信度或置信概率，区间\(\hat{\theta }_{1},\hat{\theta }_{2}\)为参数\(\theta\)的置信度为\(1-\alpha\)的置信区间。分别称\(\hat{\theta }_{1},\hat{\theta }_{2}\)为置信上线和置信下线。 置信度为0.95是指100组样本值所得置信区间的实现中，约有95个能覆盖\(\theta\)，而不是说一个实现以0.95的概率覆盖了\(\theta\)。区间的宽度反应了估计的精度，区间越小，精度越高。区间估计中精确性和可靠性是相互矛盾的。当样本容量一定时，提供估计的可靠性，将降低估计的精度，相反，提高估计的精确性，将降低估计的可靠性。实际使用中，总是在保证一定的可靠度的情况下尽可能地提高其精度。 区间估计的基本步骤 1.选取一个合适的随机变量T,这个随机变量一方面包括了待估参数\(\theta\)，另一方面，它的分布是已知的； 2.根据实际需要，选取合适的置信度\(1-\alpha\)； 3.根据相应分布的分位数的概念，写出如下形式的概率表达式 \(P\left \{ T_{1}<T <T_{2} \right \}=1-\alpha\) 4.将上式表达形式变为 \(P\left \{ \hat{\theta }_{1}<\theta <\hat{\theta }_{2} \right \}=1-\alpha\) 5.写出参数\(\theta\)的置信区间\(\hat{\theta }_{1},\hat{\theta }_{2}\)

单正态总体参数的区间估计

方差已知时的均值的区间估计

总体方差已知，均值的置信度为\(1-\alpha\)的单侧置信上限\(\bar{X}+\frac{\sigma }{\sqrt{n}}z_{1-\alpha }\),单侧置信下线\(\bar{X}-\frac{\sigma }{\sqrt{n}}z_{1-\alpha }\)。

例5 某单位随机抽样的15位员工的身高分别为: 159 158 164 169 161 161 160 157 158 163 161 154 166 168 159, 假设身高服从方差为4的正态分布, 要求估计该单位员工身高均值的置信区间，置信水平为95%。

z.test<-function(x,sigma,conf.level=0.95,u0=0,alt="two.sided"){
  result<-list()
  mean<-mean(x)
  a=1-conf.level
  n <- length(x)
  z<-(mean-u0)/(sigma/sqrt(n))
  p<-pnorm(z,lower.tail=F)
  result$z<-z
  result$p.value<-p
  if(alt=="two.sided"){
    result$p.value<-2*p
  }
  else if (alt == "greater"|alt =="less" ){
    result$p.value<p
  }
  result$interval<-c(mean-sigma*qnorm(1-a/2,0,1)/sqrt(n),mean+sigma*qnorm(1-a/2,0,1)/sqrt(n))
  result
}

x<-c(159,158,164,169,161,161,160,157,158,163,161,154,166,168,159)
result<-z.test(x,4) #默认95%的置信区间
result
## $z
## [1] 156.0812
## 
## $p.value
## [1] 0
## 
## $interval
## [1] 159.1758 163.2242

方差未知时的均值的区间估计

总体方差未知，均值的置信度为\(1-\alpha\)的置信区间为\(\bar{X}+\frac{\sigma }{\sqrt{n}}t_{1-\alpha },\bar{X}-\frac{\sigma }{\sqrt{n}}t_{1-\alpha }\)。 方差未知时我们直接利用R语言的t.test( )来求置信区间。 例6 假设不知道例5中总体的方差，要求估计该单位员工身高均值的置信区间，置信水平为95%。

t.test(x)
## 
##  One Sample t-test
## 
## data:  x
## t = 150.04, df = 14, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  158.8957 163.5043
## sample estimates:
## mean of x 
##     161.2

方差的区间估计

方差置信水平\(1-\alpha\)的置信区间为\(\left ( \frac{(n-1)S^{2}}{\chi_{1-\alpha /2}^{2}(n-1)},\frac{(n-1)S^{2}}{\chi_{\alpha/2}^{2}(n-1)} \right )\) 例7 假设不知道例5中总体的方差，要求估计该单位员工身高方差的置信区间，置信水平为95%。

chisq.var.test<-function(x,conf.level=0.95,alt="two.sided",sigma0=1)   #默认95%的置信区间 双侧检验
{
  result<-list()
  n <- length(x)
  a=1-conf.level
  v<-var(x)
  result$interval<-c((n-1)*v/qchisq(1-a/2,n-1,lower.tail=T),(n-1)*v/qchisq(a/2,n-1,lower.tail=T))
  chi2<-(n-1)*v/sigma0
  result$chi2<-chi2
  p<-pchisq(chi2,n-1)
  if(alt=="two.sided")
    result$p.value<-2*min(pchisq(chi2,n-1),pchisq(chi2,n-1,lower.tail=F))
  else
    result$p.value<-pchisq(chi2,n-1,lower.tail=F)
  result
}

chisq.var.test(x)
## $interval
## [1]  9.280619 43.064806
## 
## $chi2
## [1] 242.4
## 
## $p.value
## [1] 2.138353e-43

两正态总体参数的区间估计

均值差\(\mu _{1}-\mu_{2}\)的置信区间

两方差都已知时两均值差的置信区间

两方差已知\(\mu _{1}-\mu_{2}\)的置信水平\(1-\alpha\)的置信区间为 \((\bar{X}-\bar{Y}-z_{1-\alpha /2}\sqrt{\frac{\sigma_{1}^{2} }{n_{1}}+\frac{\sigma_{2}^{2} }{n_{2}}},\bar{X}-\bar{Y}+z_{1-\alpha /2}\sqrt{\frac{\sigma_{1}^{2} }{n_{1}}+\frac{\sigma_{2}^{2} }{n_{2}}})\)

例8 为比较两种药品的降糖效果, 选择20名条件相似的受试者, 分别服用甲、乙两种药品后，测得甲组的血糖为6.48 4.00 5.54 6.89 5.14 4.60 3.67 4.32 4.80 7.50,乙组的血糖为4.16 10.77 9.08 5.95 6.36 3.77 5.18 6.76 3.86 3.63。假定甲乙两组血糖均服从正态分布，甲组的方差为5.0，乙组的方差为5.2，试求这两组平均血糖差的置信区间(取\(\alpha=0.05\))

two.sample.ci<-function(x,y,sigma1,sigma2,conf.level=0.95){
  m= length(x)
  n = length(y)
  xbar=mean(x)-mean(y) 
  alpha = 1 - conf.level
  zstar= qnorm(1-alpha/2)* (sigma1/m+sigma2/n)^(1/2)
  xbar +c(-zstar, +zstar)
}

x <- c(6.48,4.00,5.54,6.89,5.14,4.60,3.67,4.32,4.80,7.50)
y <- c(4.16,10.77,9.08,5.95,6.36,3.77,5.18,6.76,3.86,3.63)
sigma1<-5.0
sigma2<-5.2

two.sample.ci(x,y,sigma1,sigma2)
## [1] -2.637467  1.321467
两方差都未知但相等时两均值差的置信区间

两方差未知\(\mu _{1}-\mu_{2}\)的置信水平\(1-\alpha\)的置信区间为 \((\bar{X}-\bar{Y}-t_{1-\alpha /2}\sqrt{\frac{\sigma_{1}^{2} }{n_{1}}+\frac{\sigma_{2}^{2} }{n_{2}}},\bar{X}-\bar{Y}+t_{1-\alpha /2}\sqrt{\frac{\sigma_{1}^{2} }{n_{1}}+\frac{\sigma_{2}^{2} }{n_{2}}})\)

例9 假设例8中两组的方差未知，试求这两组平均血糖差的置信区间(取\(\alpha=0.05\))

t.test(x,y,var.equal=TRUE) #两方差都未知但相等时两均值差的置信区间
## 
##  Two Sample t-test
## 
## data:  x and y
## t = -0.76237, df = 18, p-value = 0.4557
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.471296  1.155296
## sample estimates:
## mean of x mean of y 
##     5.294     5.952

两方差比的置信区间

\(\sigma _{1}^{2}/\sigma _{1}^{2}\)的置信水平\(1-\alpha\)的置信区间为\(\begin{pmatrix} \frac{S_{1}^{2}}{S_{2}^{2}}\frac{1}{F_{1-\alpha/2}(n_{1}-1,n_{2}-1)},\frac{S_{1}^{2}}{S_{2}^{2}}\frac{1}{F_{\alpha/2}(n_{1}-1,n_{2}-1)} \end{pmatrix}\)

例10 假设例9中两组的方差未知，试求这两组平均血糖方差比的置信区间(取\(\alpha=0.05\))

var.test(x,y)
## 
##  F test to compare two variances
## 
## data:  x and y
## F = 0.28552, num df = 9, denom df = 9, p-value = 0.07585
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.07091824 1.14948739
## sample estimates:
## ratio of variances 
##          0.2855164

单总体比率\(p\)的区间估计

\(x\)为容量为n的样本中具有某种特征的个体数量，则样本比例为\(x/n\)。当总体中的样品数足够多时，\(x\)近似服从二项分布\(b(n,p)\) (实际上它是超几何分布),这时总体比例可用样本比例来估计,总体比例为\(p\)的置信水平\(1-\alpha\)的置信区间为\(\left ( \hat{p}-z_{1-\alpha/2}\sqrt{\hat{p}(1-\hat{p})/n},\hat{p}+z_{1-\alpha/2}\sqrt{\hat{p}(1-\hat{p})/n} \right )\)

例10 在某小学随机抽取了120人，发现其中34人有不同程度的视力下降，假定样本的数量服从正态分布，以95%的置信度, 估计这个小学视力下降比例。

prop.test(34,120,correct=TRUE) #correct选项为是否做连续性校正
## 
##  1-sample proportions test with continuity correction
## 
## data:  34 out of 120, null probability 0.5
## X-squared = 21.675, df = 1, p-value = 3.23e-06
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.2067136 0.3740761
## sample estimates:
##         p 
## 0.2833333
binom.test(34,120)
## 
##  Exact binomial test
## 
## data:  34 and 120
## number of successes = 34, number of trials = 120, p-value =
## 2.275e-06
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.2048557 0.3728205
## sample estimates:
## probability of success 
##              0.2833333

两总体比率差\(p_{1}-p{2}\)的区间估计

在近似正态性下\(p_{1}-p{2}\)置信水平为\(1-\alpha\)的区间估计\((\hat{p_{1}}-\hat{p_{2}})\pm z_{1-\alpha/2}\sqrt{\frac{\hat{p_{1}}(1-\hat{p_{1}})}{n_{1}}+\frac{\hat{p_{2}}(1-\hat{p_{2}})}{n_{2}}}\)

例11 对某疾病进行调查。在甲地区调查了160人，有98人符合诊断标准，在乙地区调查了206人，有132人符合诊断标准。试以95%的可靠性对该病在两地差别作出区间估计。

s <- c(98,132)
t <- c(160,206)
prop.test(s,t)
## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  s out of t
## X-squared = 0.19915, df = 1, p-value = 0.6554
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.13378294  0.07722954
## sample estimates:
##    prop 1    prop 2 
## 0.6125000 0.6407767

基于Bootstrap的区间估计

对于大多数的统计量而言，计算置信区间的数学公式太过复制或者根本没有数学公式，因而没有计算置信区间的已知的解析形式的公式。采用Bootstrap方法，解析形式的数据公式是不必知道的。当样本不符合理论分布假设时，求样本统计量的置信区间就成为一个难题。而自助法(Bootstrap)的思路是对原始样本重复抽样产生多个新样本，针对每个样本求取统计量，然后得到它的经验分布，再通过求经验分布的分位数来得到统计量的置信区间，这种方法不需要对统计量有任何理论分布的假设。一般认为，只要样本具有代表性，采用自助法需要的原始样本只要20-30个,重复抽样1000次就能达到满意的结果。在R中进行自助法是利用boot扩展包，其流程如下：1.编写一个求取统计量的自定义函数 2.将上面的函数放入boot（）函数中进行运算，得到自助法的结果 3.用boot.ci()函数求取置信区间

例12 对mtcars数据集中的mpg变量估算其中位数的置信区间，将wt和disp作为自变量，mpg 作为因变量，进行回归后对判定系数R-square，用自助法求它的95%置信区间。

myfun <- function(data,indices){
  d <- data[indices,]
  rs <- median(d$mpg)
  return (rs)
}

results=boot(data=mtcars,statistic=myfun ,R=1000) #results这个数据结构中包括了原始样本的统计(results$t0)和再抽样样本的统计量(results$t0)，
boot.ci(results,conf=0.95,type=c('perc','bca')) 
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = results, conf = 0.95, type = c("perc", "bca"))
## 
## Intervals : 
## Level     Percentile            BCa          
## 95%   (16.4, 21.4 )   (15.8, 21.2 )  
## Calculations and Intervals on Original Scale
## Some BCa intervals may be unstable
#其中conf表示置信水平，type表示了用何种算法来求区间，perc即使用百分位方法，bca表示adjusted bootstrap percentile，即对偏差进行了调整

rsq=function(data,indices){
d=data[indices,]
fit=lm(formula=mpg~wt+disp,data=d)
return(summary(fit)$r.square)
}
results=boot(data=mtcars,statistic=rsq,R=1000)
print(results)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = mtcars, statistic = rsq, R = 1000)
## 
## 
## Bootstrap Statistics :
##      original    bias    std. error
## t1* 0.7809306 0.0110107  0.04997466
plot(results) 

#左侧的直方图表示了再抽样样本的统计量的经验分布，其中的虚线表示了原始样本的统计量，从中可以观察到偏差。右侧QQ图有助于判断经验分布是否正态。
boot.ci(results,conf=0.95,type=c('perc','bca'))
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = results, conf = 0.95, type = c("perc", "bca"))
## 
## Intervals : 
## Level     Percentile            BCa          
## 95%   ( 0.6826,  0.8757 )   ( 0.6187,  0.8516 )  
## Calculations and Intervals on Original Scale
## Some BCa intervals may be unstable

例12 假定例10的样本数量的分布状况未知，以95%的置信度, 估计这个小学视力下降比例。

p.test <- function(x,n)
{
  n <- rep(0,each=n)
  n[seq(1,x)] <- 1
  data <- as.data.frame(n)
  
  p <- function(data,indices)
  {
    d <- data[indices,]
    ratio<- sum(d)/length(d)
    return (ratio)
  }
  results=boot(data=data,statistic=p,R=1000)
  boot.ci(results,conf=0.95,type=c('perc','bca'))
}
  
p.test(34,120)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = results, conf = 0.95, type = c("perc", "bca"))
## 
## Intervals : 
## Level     Percentile            BCa          
## 95%   ( 0.2083,  0.3667 )   ( 0.2000,  0.3583 )  
## Calculations and Intervals on Original Scale