hello,大家好,上一篇分享了如何用R语言实现蒙特卡洛模拟,并用蒙特卡洛模拟计算了分布的均值和方差,今天给大家分享如何用R语言来进行矩估计和似然函数的求解。
因为在求解矩估计和似然函数时,可能会遇到非线性方程组,所以先给大家介绍一下如何用Newton法来求解非线性方程组。
本文所涉及的前两道例题来自于《R统计建模与R软件》——薛毅、陈立萍编著。
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.9129849 0.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},0<x<1,\theta>0. p(x;θ)=θxθ−1,0<x<1,θ>0.
按一般做法我们需要由 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)=∫01xθxθ−1dx推导出矩估计量 θ ^ \hat\theta θ^,现在我们不推导直接用uniroot求解;
因为这个密度函数是自定义的一个密度函数,因此我们需要先写一个服从该密度函数的随机数生成函数:
rdensity <- function(n, theta){
obj <- function(x){
sqrt(theta)*x^(sqrt(theta)-1)
}
u <- c()
while(length(u)<n){
x <- runif(1,0,1)
y <- runif(1,0,sqrt(theta)) #sqrt(theta)是当x=1是所对应的密度函数值
if(y<=obj(x)){
u <- c(u,x)
}
}
return(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π)211(1+pt2)2p+11,其中Γ(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∏nf(ti;p)=[Γ(2p)Γ(2p+1)(pπ)2p+11]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)−2nln(pπ)+2p+1i=1∑nln(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
dpdlnL(p;t)=2nΓ(2p+1)Γ′(2p+1)−2nΓ(2p)Γ′(2p)−2pn−21i=1∑nln(1+pti2)+2p+1i=1∑n1+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
注:
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,不然会报错。
获取代码
本文代码均已上传,关注公众号,回复“似然函数”,即可获得