锐客网

  1. 首页 > it技术 > >

R语言|R语言（Newton法、似然函数）

数据分析r语言似然函数Newton法

hello，大家好，上一篇分享了如何用R语言实现蒙特卡洛模拟，并用蒙特卡洛模拟计算了分布的均值和方差，今天给大家分享如何用R语言来进行矩估计和似然函数的求解。
因为在求解矩估计和似然函数时，可能会遇到非线性方程组，所以先给大家介绍一下如何用Newton法来求解非线性方程组。

本文所涉及的前两道例题来自于《R统计建模与R软件》——薛毅、陈立萍编著。

目录
  • Newton法
    • 例1：求解方程组
  • 矩估计
    • 例2：设总体 X X X服从二项分布 B ( k , p ) B(k,p) B(k,p)，其中 k , p k,p k,p为未知参数， X 1 , X 2 , … , X n X_1,X_2,\dots,X_n X1?,X2?,…,Xn?是总体 X X X的一个样本，求参数 k , p k,p k,p的矩估计 k ^ , p ^ \hat{k},\hat{p} k^,p^?.
    • 例3：设总体密度函数如下， x 1 , … , x n x_1,\dots,x_n x1?,…,xn?是样本，试求未知参数的矩估计
  • 似然函数
    • 例3：设总体 X X X服从自由度为p的t分布，其概率密度函数为
  • 获取代码

Newton法 牛顿迭代法（Newton’s method）是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法。
设 r r r是 f ( x ) = 0 f(x)=0 f(x)=0的根，选取 x 0 x_0 x0?作为 r r r的初始近似值，过点 ( x 0 , f ( x 0 ) ) (x_0,f(x_0)) (x0?,f(x0?))做曲线 y = f ( x ) y=f(x) y=f(x)的切线 L L L， L : y = f ( x 0 ) + f ′ ( x 0 ) ( x ? x 0 ) L: y=f(x_0)+f'(x_0)(x-x_0) L:y=f(x0?)+f′(x0?)(x?x0?)，则 L L L与 x x x轴的交点的横坐标 x 1 = x 0 ? f ( x 0 ) f ′ ( x 0 ) x_1=x_0-\frac{f(x_0)}{f'(x_0)} x1?=x0??f′(x0?)f(x0?)?，称 x 1 x_1 x1?为 r r r的一次近似值。过点 ( x 1 , f ( x 1 ) ) (x_1,f(x_1)) (x1?,f(x1?))做曲线 y = f ( X ) y=f(X) y=f(X)的切线，并求该切线与 x x x轴交点的横坐标 x 2 = x 1 ? f ( x 1 ) f ′ ( x 1 ) x_2=x_1-\frac{f(x_1)}{f'(x_1)} x2?=x1??f′(x1?)f(x1?)?称 x 2 x_2 x2?为 r r r的二次近似值。重复以上过程，得 r r r的近似值序列，其中， x n + 1 = x n ? f ( x n ) f ′ ( x n ) x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)} xn+1?=xn??f′(xn?)f(xn?)?称为 r r r的 n + 1 n+1 n+1次近似值，上式称为牛顿迭代公式。
利用牛顿迭代算法的基本思路：确定迭代变量->建立迭代关系式->对迭代过程进行控制，接下来我们用一个例子来讲解。
例1：求解方程组 { x 1 2 + x 2 2 ? 5 = 0 ( x 1 + 1 ) x 2 ? ( 3 x 1 + 1 ) = 0 \begin{cases} x_1^2+x_2^2-5=0\\ (x_1+1)x_2-(3x_1+1)=0 \end{cases} {x12?+x22??5=0(x1?+1)x2??(3x1?+1)=0?
1、确定迭代变量 x = ( x 1 , x 2 ) T x=(x_1,x_2)^T x=(x1?,x2?)T，设定初始值 x ( 0 ) = ( 0 , 1 ) T x^{(0)}=(0,1)^T x(0)=(0,1)T；
2、建立迭代关系式：
x ( k + 1 ) = x ( k ) ? [ J ( x ( k ) ) ? 1 ] f ( x k ) x^{(k+1)}=x^{(k)}-[J(x^{(k)})^{-1}]f(x^{k}) x(k+1)=x(k)?[J(x(k))?1]f(xk)
其中 J ( x ) J(x) J(x)为函数 f ( x ) f(x) f(x)的Jacobi矩阵，即
J = ( ? f 1 ? x 1 ? f 1 ? x 2 … ? f 1 ? x n ? f 2 ? x 1 ? f 2 ? x 2 … ? f 2 ? x n ? ? ? ? f n ? x 1 ? f n ? x 2 … ? f n ? x n ) J=\begin{pmatrix} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \dots & \frac{\partial f_1}{\partial x_n}\\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \dots & \frac{\partial f_2}{\partial x_n}\\ \vdots & \vdots & & \vdots\\ \frac{\partial f_n}{\partial x_1} & \frac{\partial f_n}{\partial x_2} & \dots & \frac{\partial f_n}{\partial x_n} \end{pmatrix} J=????????x1??f1???x1??f2????x1??fn????x2??f1???x2??f2????x2??fn???………??xn??f1???xn??f2????xn??fn??????????
3、对迭代过程进行控制，即精度要求 ε = 1 0 ? 5 \varepsilon = 10^{-5} ε=10?5。
Newtons <- function(fun, x, ep = 1e-5, it_max = 100){ index <- 0; k <- 1 while (k <= it_max){ x1 <- x; obj <- fun(x) x <- x - solve(obj$J, obj$f) norm <- sqrt((x - x1) %*% (x - x1)) if (norm < ep){ index <- 1; break } k <- k + 1 } obj <- fun(x) list(root = x, it = k, index = index, FunVal = obj$f) }

funs <- function(x){ f <- c(x[1]^2 + x[2]^2 - 5, (x[1] + 1)*x[2] - (3*x[1] + 1)) J <- matrix(c(2*x[1], 2*x[2], x[2]-3, x[1]+1), nrow = 2, byrow = T) list(f = f, J = J) }

Newtons(funs,c(0,1)) ## $root ## [1] 1 2 ## ## $it ## [1] 6 ## ## $index ## [1] 1 ## ## $FunVal ## [1] 1.598721e-14 6.217249e-15

所以方程的解为 x ? = ( 1 , 2 ) T x^*=(1,2)^T x?=(1,2)T，总共迭代了6次。
矩估计 例2：设总体 X X X服从二项分布 B ( k , p ) B(k,p) B(k,p)，其中 k , p k,p k,p为未知参数， X 1 , X 2 , … , X n X_1,X_2,\dots,X_n X1?,X2?,…,Xn?是总体 X X X的一个样本，求参数 k , p k,p k,p的矩估计 k ^ , p ^ \hat{k},\hat{p} k^,p^?. 由二项分布的均值（一阶原点矩）和方差（二阶中心矩）可得方程组
{ k p ? X ˉ = 0 k p ( 1 ? p ) ? M 2 = 0 \begin{cases} kp-\bar{X}=0\\ kp(1-p)-M_2=0 \end{cases} {kp?Xˉ=0kp(1?p)?M2?=0?
moment_fun <- function(p){ f <- c(p[1]*p[2]-A1,p[1]*p[2]-p[1]*p[2]^2-M2) J <- matrix(c(p[2],p[1],p[2]-p[2]^2,p[1]-2*p[1]*p[2]),nrow=2,byrow=T) list(f=f,J=J) }

x <- rbinom(100,20,0.7) n <- length(x) A1 <- mean(x) M2 <- (n-1)/n*var(x) p <- c(10,0.5) Newtons(moment_fun,p) ## $root ## [1] 19.91298490.7221419 ## ## $it ## [1] 6 ## ## $index ## [1] 1 ## ## $FunVal ## [1] 1.776357e-15 1.776357e-15

从结果可以看出，误差非常的小，但是也发现了一个弊端，在能用这个算法的情况下，计算往往也比较简单，所以它的效率相对较低。
接下来再给大家分享一个新的工具——uniroot函数，在遇到一些较为复杂的一元方程时可以用uniroot函数进行求解。
例3：设总体密度函数如下， x 1 , … , x n x_1,\dots,x_n x1?,…,xn?是样本，试求未知参数的矩估计 p ( x ; θ ) = θ x θ ? 1 , 0 < x < 1 , θ > 0. p(x; \theta)=\sqrt\theta x^{\sqrt\theta-1},00. p(x; θ)=θ ?xθ ??1,00.
按一般做法我们需要由 E ( X ) = ∫ 0 1 x θ x θ ? 1 d x E(X)=\int_0^1x\sqrt\theta x^{\sqrt\theta-1}\mathrm{d}x E(X)=∫01?xθ ?xθ ??1dx推导出矩估计量 θ ^ \hat\theta θ^，现在我们不推导直接用uniroot求解；
因为这个密度函数是自定义的一个密度函数，因此我们需要先写一个服从该密度函数的随机数生成函数：
rdensity <- function(n, theta){ obj <- function(x){ sqrt(theta)*x^(sqrt(theta)-1) } u <- c() while(length(u)

注：
这里的随机数生成采用的是随机投点的方式，取落在密度函数内的值。
x <- rdensity(100, 5) fun <- function(theta){ obj <- function(x) x*sqrt(theta)*x^(sqrt(theta)-1) integrate(obj, 0, 1)$value-mean(x) } uniroot(fun,c(1,10)) ## $root ## [1] 5.214032 ## ## $f.root ## [1] 2.338409e-07 ## ## $iter ## [1] 6 ## ## $init.it ## [1] NA ## ## $estim.prec ## [1] 6.103516e-05

注：
integrate()为定积分函数。
从结果可以看出，准确率还是非常高的。
似然函数 最后，我们再在似然函数上进行下实验，以t分布为例；
例3：设总体 X X X服从自由度为p的t分布，其概率密度函数为 f ( t ; p ) = Γ ( p + 1 2 ) Γ ( p 2 ) 1 ( p π ) 1 2 1 ( 1 + t 2 p ) p + 1 2 , 其 中 Γ ( x ) = ∫ 0 ∞ t x ? 1 e ? t d t f(t; p)=\frac{\varGamma(\frac{p+1}{2})}{\varGamma(\frac{p}{2})}\frac{1}{(p\pi)^\frac12}\frac{1}{(1+\frac{t^2}{p})^{\frac{p+1}2}},其中\varGamma(x)=\int_0^\infty t^{x-1}e^{-t}\mathrm{d}t f(t; p)=Γ(2p?)Γ(2p+1?)?(pπ)21?1?(1+pt2?)2p+1?1?,其中Γ(x)=∫0∞?tx?1e?tdt
其中 p p p为未知参数. X 1 , X 2 , … , X n X_1,X_2,\dots,X_n X1?,X2?,…,Xn?是来自总体 X X X的样本，求 p p p的极大似然估计.
解：t分布的似然函数为
L ( p ; t ) = ∏ i = 1 n f ( t i ; p ) = [ Γ ( p + 1 2 ) Γ ( p 2 ) 1 ( p π ) p + 1 2 ] n 1 ∏ i = 1 n [ ( 1 + t i 2 p ) p + 1 2 ] L(p; t)=\prod_{i=1}^nf(t_i; p)=[\frac{\varGamma(\frac{p+1}{2})}{\varGamma(\frac{p}{2})}\frac{1}{(p\pi)^\frac{p+1}{2}}]^n\frac1{\prod_{i=1}^n[(1+\frac{t_i^2}p)^\frac{p+1}2]} L(p; t)=i=1∏n?f(ti?; p)=[Γ(2p?)Γ(2p+1?)?(pπ)2p+1?1?]n∏i=1n?[(1+pti2??)2p+1?]1?
相应的对数似然函数为
ln ? L ( p ; t ) = n ln ? Γ ( p + 1 2 ) ? n ln ? Γ ( p 2 ) ? n 2 ln ? ( p π ) + p + 1 2 ∑ i = 1 n ln ? ( 1 + t i 2 p ) \ln L(p; t)=n\ln\varGamma(\frac{p+1}{2})-n\ln{\varGamma(\frac{p}{2})}-\frac n2\ln(p\pi)+\frac{p+1}2\sum_{i=1}^n\ln(1+\frac{t_i^2}p) lnL(p; t)=nlnΓ(2p+1?)?nlnΓ(2p?)?2n?ln(pπ)+2p+1?i=1∑n?ln(1+pti2??)
得到对数似然方程
d d p ln ? L ( p ; t ) = n 2 Γ ′ ( p + 1 2 ) Γ ( p + 1 2 ) ? n 2 Γ ′ ( p 2 ) Γ ( p 2 ) ? n 2 p ? 1 2 ∑ i = 1 n ln ? ( 1 + t i 2 p ) + p + 1 2 ∑ i = 1 n ( t i p ) 2 1 + t i 2 p = 0 \frac{\mathrm{d}}{\mathrm{d}p}\ln L(p; t)=\frac n2\frac{\varGamma'(\frac{p+1}2)}{\varGamma(\frac{p+1}2)}-\frac n2\frac{\varGamma'(\frac p2)}{\varGamma(\frac p2)}-\frac{n}{2p}-\frac12\sum_{i=1}^n\ln(1+\frac{t_i^2}p)+\frac{p+1}2\sum_{i=1}^n\frac{(\frac{t_i}{p})^2}{1+\frac{t_i^2}p}=0 dpd?lnL(p; t)=2n?Γ(2p+1?)Γ′(2p+1?)??2n?Γ(2p?)Γ′(2p?)??2pn??21?i=1∑n?ln(1+pti2??)+2p+1?i=1∑n?1+pti2??(pti??)2?=0
其中 Γ ′ ( x ) = ∫ 0 ∞ t x ? 1 ln ? t e ? t d t \varGamma'(x)=\int_0^\infty t^{x-1}\ln te^{-t}\mathrm{d}t Γ′(x)=∫0∞?tx?1lnte?tdt
n <- 100000 t <- rt(n, 6) lnL <- function(p){ (n/2)*digamma((p+1)/2)-(n/2)*digamma(p/2)-n/(2*p)-0.5*sum(log(1+t^2/p))+(p+1)/2*sum((t/p)^2/(1+t^2/p)) } uniroot(lnL,c(1,10)) ## $root ## [1] 6.092751 ## ## $f.root ## [1] -0.004239805 ## ## $iter ## [1] 8 ## ## $init.it ## [1] NA ## ## $estim.prec ## [1] 6.103516e-05

注：
【R语言|R语言（Newton法、似然函数）】1、 d i g a m m a ( x ) = Γ ′ ( x ) Γ ( x ) digamma(x)=\frac{\varGamma'(x)}{\varGamma(x)} digamma(x)=Γ(x)Γ′(x)?；
2、因为方程中p有充当分母，因此给定的范围不能包含0，不然会报错。
获取代码 本文代码均已上传，关注公众号，回复“似然函数”，即可获得
R语言|R语言（Newton法、似然函数）
文章图片

    推荐阅读


    上一篇：反射望远镜的发展历程(1)

    下一篇：使用Google|使用Google Cloud虚拟机做学习测试