矩法是K.Pearson于1894年基于某些生物学方面的数据并非来自正态分布这一发现而提出的，这种估计方法在实际中有着广泛的应用。矩估计法的主要思想是用样本矩去估计总体矩，我们首先得知道样本矩和总体矩是什么。

  最常用的矩有两种，定义：对于样本\(X_1,...X_n\)及任意一正整数\(k\)，我们称\[a_k=\frac{1}{n}\sum_{i=1}^nX_i^k,m_k=\frac{1}{n}\sum_{i=1}^n(X_i-\bar X)^k\]为样本的\(k\)阶原点矩和\(k\)阶中心矩。

  而一个随机变量\(X\)的\(k\)阶原点矩和中心矩就叫作总体的\(k\)阶原点矩和中心矩，定义为：\[\mu_k=EX^k,v_k=E(X-\mu_k)^k\]

  对比定义可以发现，样本矩不依赖于总体中的参数，它们是统计量；而总体矩则与分布中的未知参数有关，而在极限情况下，样本矩又是总体矩的一个很好的估计，于是，我们可以利用上述想法来估计总体中的未知参数。

  设总体X的分布中的未知参数为\(\theta=(\theta_1,\theta_2,...,\theta_m)'\),并且总体的k阶原点矩存在，我们令总体的k阶原点矩等于样本的k阶原点矩，即\[a_k(\theta_1,\theta_2,...\theta_m)=E(X^k)=\frac{1}{k}\sum\limits_{i=1}^nX_i^k,k=1,2,...m\]  由上式可以得到关于\(\theta\)的解，此解则可以作为\(\theta\)的矩估计。

  举一个具体的例子，设总体\(X\)的均值为\(\mu\)，方差为\(\sigma\)，\(X_1,...X_n\)是来自总体\(X\)的一个样本，试用矩法估计均值\(\mu\)和方差\(\sigma\)。

  计算总体\(X\)的一阶、二阶原点矩，有\[\alpha_1=E(X)=\mu,\\\alpha_2=E(X^2)=E(X-E(X)+E(X))^2\\=E(X-E(X))^2+(E(X))^2+2E(X-E(X))E(X)\\=E(X-E(X))^2+(E(X))^2=\sigma^2+\mu^2\]

  计算总体的一阶、二阶原点矩，有\[A_1=\frac{1}{n}\sum_{i=1}^nX_i=\bar X,A_2=\frac{1}{n}\sum_{i=1}^nX_i^2\]

  令样本矩和总体矩对应相等，可得到方程组：\[\begin{cases}\mu=\bar X,\\\sigma^2+\mu^2=\frac{1}{n}\sum_{i=1}^nX_i^2\end{cases}\]

  解上述方程组得到均值\(\mu\)和方差\(\sigma^2\)的矩估计如下：\[\hat{\mu}=\bar X,\\\hat{\sigma}^2=\frac{1}{n}\sum_{i=1}^nX_i^2-{\bar X}^2=\frac{1}{n}\sum_{i=1}^n(X_i-\bar X)^2\]

  可能会有人想问，为什么方差的矩估计并不等于样本方差\(S_n^2\)呢？参数的估计方法有很多，之所以用\(S_n^2\)来定义样本方差，是因为它是好的，是经得起推敲的，那为什么它是好的呢？在后边估计量的优良性准则中会具体讲到。

  由于不同分布造成函数f形态各异，因此在R中不可能有固定的R程序来直接得到矩估计结果，只能利用R的计算功能具体问题具体分析。首先介绍几个R中常用的解方程函数，如下表所示

  其中BBsolve()与给定的参数初值有关，在实际中并不常用，因此这里终点介绍前两个函数的使用方法。首先是解一次方程的函数uniroot()：

uniroot(f, interval, ...,
        lower = min(interval), upper = max(interval),
        f.lower = f(lower, ...), f.upper = f(upper, ...),...
        tol = .Machine$double.eps^0.25, maxiter = 1000, trace = 0)

  其中f是指定所要求解方程的函数；interval是一个数值向量，指定要求解的根的区间范围；或者用lower和upper分别指定区间的两个端点；tol表示所需的精度（收敛容忍度）；maxiter为最大迭代次数。

  如果遇到多元方程的求解，就需要利用rootSolve包的函数multiroot()来解方程组。multiroot()用于对n个非线性方程组求解n个根，其要求完整的雅可比矩阵，采用Newton-Raphson方法。调用格式为

multiroot(f,start,maxiter=100,rtol=1e-6,atol=1e-8,ctol=1e-8,useFortran=TRUE,...)

  f指定所要求解的函数；由于使用的是牛顿迭代法，因而必须通过start给定根的初始值，其中的name属性还可以标记输出变量的名称；maxiter是允许的最大迭代次数；rtol和atol分别为绝对误差和相对误差，一般保持默认值即可；ctol也是一个用于控制迭代次数的标量，如果两次迭代的最大变化值小于col，那么迭代停止，得到方程组的根。

  接下来用实例说明在R中求矩估计的方法，由于样本分布的多样性，矩估计必须“具体问题具体分析”，并没有一个固定的求解套路。

  例如，已知某种保险产品在一个保单年度内的损失情况如下表所示，其中给出了不同损失次数下的保单数，我们对损失次数的分布进行估计。已知分布类型是Poisson分布，其样本均值即为参数\(\lambda\)的矩估计。

num=c(rep(0:5,c(1532,581,179,41,10,4)))#用rep函数生成样本，样本值由0-5个数字构成，函数中第二个向量对应表示每个数字的重复次数
lambda=mean(num)
lambda

## [1] 0.4780571

  根据估计的参数值，画图比较损失次数的估计值和样本值之间的差别

k=0:5
ppois=dpois(k,lambda)
poisnum=ppois*length(num)#由possion分布生成的损失次数
plot(k,poisnum,ylim=c(0,1600))#画图比较，为图形效果更好，用参数ylim设置纵轴的范围，最小值为0，最大值要大于样本的最值，选取1600
samplenum=as.vector(table(num))#样本的损失次数
points(k,samplenum,type="p",col=2)
legend(4,1000,legend=c("num","possion"),col=1:2,pch="o")#设置图例的位置颜色形状和内容

#install.packages("rootSolve")#安装软件包
x=c(4,5,4,3,9,9,5,7,9,8,0,3,8,0,8,7,2,1,1,2)
m1=mean(x)#样本均值
m2=var(x)#样本方差
model=function(x,m1,m2){
  c(f1=x[1]+x[2]-2*m1,
    f2=(x[2]-x[1])^2/12-m2)
}
library(rootSolve)#加载软件包
multiroot(f=model,start=c(0,10),m1=m1,m2=m2)

## $root
## [1] -0.7523918 10.2523918
## 
## $f.root
##            f1            f2 
## -5.153211e-12  1.121688e-09 
## 
## $iter
## [1] 4
## 
## $estim.precis
## [1] 5.634204e-10