Код Newtonraphson в R приводит к разным результатам

Question

Мне нужно приблизить параметры выборки из области Бирнбаум-Сондерс. Вот мой код:

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 

Проблема заключается в том, что хотя условие сходимости, кажется, встречались, мои приближения отличаются, на мой взгляд, значительно, в зависимости от тэта я выбираю в качестве отправной точки. Таким образом, я никогда не знаю, какие результаты я получу, выбрав еще немного другую отправную точку.

Любые идеи, почему метод настолько неустойчив?

Answer 1

Я бы не стал изобретать колесо для этой проблемы и использовать собственный алгоритм. Я бы использовал некоторую уже встроенную функцию в одном пакете R, который реализует алгоритм Ньютона-raphson.

Например, здесь с помощью 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).