Newtonraphson-Code in R führt zu unterschiedlichen Ergebnissen

Question

Ich brauche die Parameter einer Probe aus Birnbaum-Saunders Distr. hier ist mein Code:

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 

Das Problem ist, dass, obwohl der Zustand der Konvergenz erfüllt zu sein scheint, meine Annäherungen unterscheiden, meiner Meinung nach deutlich, je nachdem, welche Theta ich als Ausgangspunkt wählen. Somit weiß ich nie, welche Ergebnisse ich erzielen werde, wenn ich einen etwas anderen Startpunkt wähle.

Irgendwelche Ideen, warum die Methode so instabil ist?

Answer 1

Ich würde das Rad für diese Art von Problem nicht neu erfinden und einen benutzerdefinierten Algorithmus verwenden. Ich würde eine bereits eingebaute Funktion in einem der mehreren R-Pakete verwenden, die den Newton-Raphson-Algorithmus implementieren.

Zum Beispiel, hier mit rootSolve Paket:

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 

bekam ich das gleiche Ergebnis mit theta2 <- c(1,3).