2

我需要近似来自 Birnbaum-Saunders 分布的样本的参数。这是我的代码：

    x =c(6.7508, 1.9345, 4.9612, 22.0232, 0.2665, 66.7933, 5.5582, 60.2324, 72.5214, 1.4188, 4.6318, 61.8093, 11.3845, 1.1587, 22.8475, 8.3223, 2.6085, 24.0875, 4.6762, 8.2369)
l.der1 = function(theta,x) {
gamma <- theta[1]
beta <- theta[2]
n <- length(x)
ausdruck1=sum((sqrt(x/beta)-sqrt(beta/x))^2)
ausdruck2=sqrt(x/beta)+sqrt(beta/x)
matrix(c(-n/gamma+ausdruck1/gamma^3, sum((1/(2*x*sqrt(beta/x))-x/(2*beta^2*sqrt(x/beta)))/ausdruck2)-1/(2*gamma^2)*sum(1/x-x/beta^2)),2, 1)
}

l.der2 = function(theta,x) {
gamma <- theta[1]
beta <- theta[2]
n <- length(x)
ausdruck1=sum((sqrt(x/beta)-sqrt(beta/x))^2)
ausdruck2=sqrt(x/beta)+sqrt(beta/x)
ausdruck3=(1/gamma^3)*sum(1/x-x/beta^2)
matrix(c(n/gamma^2-(3*ausdruck1)/gamma^4,ausdruck3,ausdruck3,sum((2-beta/x+x/beta)/(2*beta^2*ausdruck2^2))-(1/(2*gamma^2))*sum(2*x/beta^3)),2, 2, byrow=T)
 }

newtonraphson = function(theta,l.der1,l.der2,x,col=2,epsilon=10^(-6)) {
I <- l.der2(theta,x)
thetastar <- theta - solve(I) %*% l.der1(theta,x)
repeat {theta=thetastar 
  thetastar <- theta - solve(I) %*% l.der1(theta,x)
  if (((thetastar[1]-theta[1])^2)/thetastar[1]^2 < 10^(-6) && ((thetastar[2]-theta[2])^2)/thetastar[2]^2 < 10^(-6) )   #calculating relative convergence
  return(thetastar)
}
}

theta = c(1,4)  #starting point
theta= newtonraphson(theta,l.der1,l.der2,x=x)  
theta

问题是，尽管似乎满足了收敛条件，但在我看来，我的近似值差异很大，这取决于我选择哪个 theta 作为起点。因此，我永远不知道选择一个稍微不同的起点会得到什么结果。

关于为什么该方法如此不稳定的任何想法？

4

1 回答 1

4

我不会为这类问题重新发明轮子并使用自定义算法。我会在一个实现 Newton-raphson 算法的多个 R 包中使用一些已经内置的函数。

例如，这里使用rootSolve包：

library(rootSolve)
theta <- c(1,4) 
multiroot(l.der1,theta,jacfunc=l.der2,x=x)
$root
[1] 1.87116 6.83414

$f.root
             [,1]
[1,] 2.168992e-08
[2,] 6.425832e-09

$iter
[1] 8

$estim.precis
[1] 1.405788e-08

我得到了相同的结果theta2 <- c(1,3)。

于 2016-12-14T22:38:56.137 回答