【数据分析 R语言实战】学习笔记 第六章 参数估计与R实现（上）（1）

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

6.1点估计及R实现

6.1.1矩估计

R中的解方程函数:

函数及所在包：功能

uniroot()@stats：求解一元（非线性)方程

multiroot()@rootSolve：给定n个(非线性)方程，求解n个根

uniroot.all()@rootSolve：在一个区问内求解一个方程的多个根

BBsolve()@BB：使用Barzilai-Borwein步长求解非线性方程组

uniroot(f,interval, ...,lower = min(interval), upper = max(interval),f.lower = f(lower,...), f.upper = f(upper, ...),extendInt = c("no", "yes","downX", "upX"), check.conv = FALSE,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, positive = FALSE,

         jacfunc = NULL, jactype = "fullint",

         verbose = FALSE, bandup = 1, banddown = 1,

         parms = NULL, ...)

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

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

损失次数	0	1	2	3	4	5
保单数	1532	581	179	41	10	4

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

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

> k=0:5> ppois=dpois(k,lambda) > poisnum=ppois*length(num)#由poisson分布生成的损失次数> 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","poisson"),col=1:2,pch="0")

rootSolve包的函数multiroot()用于解方程组：

> 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){}> 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#均匀分布的两个参数值[0.75, 10.25]$f.root           f1            f2 -5.153211e-12  1.121688e-09 $iter[1] 4 $estim.precis[1] 5.634204e-10

验证一下：

> m1-sqrt(3*m2);m1+sqrt(3*m2)[1] -0.7523918[1] 10.25239

6.1.2极大似然估计

R中计算极值的函数（stats包）

optimize( )   计算单参数分布的极人似然估计值

optim()     计算多个参数分布的极大似然估计值

nlm()       计算非线性函数的最小值点

nlminb( )    非线性最小化函数

1.函数optimize()

当分布只包含一个参数时，我们可以使用R中计算极值的函数optimize()求极大似然估计值。

optimize(f = , interval = ,  ..., lower = min(interval),         upper = max(interval), maximum = FALSE,         tol = .Machine$double.eps^0.25)其中f是似然函数:interval指定参数的取值范围;lower/upper分别是参数的下界和上界:maximum默认为FALSE，表示求似然函数的极小值，若为TRUE则求极大值:tol表示计算的精度。

2.函数optim()和nlm()

当分布包含多个参数时，用函数optim()或nlm()计算似然函数的极大值点。

optim(par, fn, gr = NULL, ...,      method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN",                 "Brent"),      lower = -Inf, upper = Inf,      control = list(), hessian = FALSE)

par设置参数的初始值；fn为似然函数；method提供了5种计算极值的方法

nlm(f, p, ..., hessian = FALSE, typsize = rep(1, length(p)),    fscale = 1, print.level = 0, ndigit = 12, gradtol = 1e-6,    stepmax = max(1000 * sqrt(sum((p/typsize)^2)), 1000),steptol = 1e-6, iterlim = 100, check.analyticals = TRUE)nlm是非线性最小化函数，仅使用牛顿一拉夫逊算法，通过迭代计算函数的最小值点。一般只布要对前两个参数进行设置:f是需要最小化的函数:P设置参数初始值。

3.函数nlminb()

在实际应用中，上面这三个基本函数在遇到数据量较大或分布较复杂的计算时，就需要使用优化函数nlminb()

nlminb(start, objective, gradient = NULL, hessian = NULL, ...,       scale = 1, control = list(), lower = -Inf, upper = Inf)

参数start是数值向量，用于设置参数的初始值;objective指定要优化的函数:gradient和hess用于设置对数似然的梯度，通常采用默认状态;control是一个控制参数的列表:lower和upper设置参数的下限和上限，如果未指定，则假设所有参数都不受约束。

例：

> library(MASS)> head(geyser,5)  waiting duration1      80 4.0166672      71 2.1500003      57 4.0000004      80 4.0000005      75 4.000000> attach(geyser)> hist(waiting,freq=FALSE)#通过直方图了解数据分布的形态

猜测分布是两个正态分布的混合，需要估计出函数中的5个参数：p、μ1、μ2、σ1、σ2。

在R中编写对数似然函数时，5个参数都存放在向量para中，由于nlminb()是计算极小值的，因此函数function中最后返回的是对数似然函数的相反数。

> l1=function(para)+ {+ f1=dnorm(waiting,para[2],para[3])+ f2=dnorm(waiting,para[4],para[5])+ f=para[1]*f1+(1-para[1])*f2+ l1=sum(log(f))+ return(-11)+ }

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏5 回帖

关键词：R语言实战 数据分析 学习笔记 参数估计 R语言 方程组 容忍度

谢谢分享