Hesse-Matrix in optim() in R

Ich habe einige Probleme mit optim() in R, um für eine Wahrscheinlichkeit mit einem Integral zu lösen und die Hesse-Matrix für die 2 Parameter zu erhalten. Der Algorithmus konvergierte, aber ich erhalte einen Fehler, wenn ich die Option hessian = TRUE in optim() verwende. Der Fehler ist:Hesse-Matrix in optim() in R

Fehler bei integrieren (integrand1, lower = s1 [i] - 1 oben = s1 [i]): nicht-endlichen Funktionswert

auch eine Warnmeldung von NAs hatte

Hier ist mein Code:

s1=c(1384,1,1219,1597,2106,145,87,1535,290,1752,265,588,1188,160,745,237,479,39,99,56,1503,158,916,651,1064,166,635,19,553,51,79,155,85,1196,142,108,325 
,135,28,422,1032,1018,128,787,1704,307,854,6,896,902) 

LLL=function (par) { 

    integrand1 <- function(x){ (x-s1[i]+1)*dgamma(x, shape=par[1], rate=par[2]) } 
    integrand2 <- function(x){ (-x+s1[i]+1)*dgamma(x, shape=par[1],rate=par[2]) } 



    likelihood = vector() 

    for(i in 1:length(s1)) {likelihood[i] = 
    log(integrate(integrand1,lower=s1[i]-1,upper=s1[i])$value+ integrate(integrand2,lower=s1[i],upper=s1[i]+1)$value) 
    } 

    like= -sum(likelihood) 
    return(like) 

} 




optim(par=c(0.1,0.1),LLL,method="L-BFGS-B", lower=c(0,0)) 
optim(par=c(0.1,0.1),LLL,method="L-BFGS-B", lower=c(0,0), hessian=TRUE)

Vielen Dank für Ihre Hilfe!

Quelle

2017-01-09 Y. Ma

Ich glaube, dass Sie 'dgamma (0, shape = .1, scale = .1)', was 'Inf' ist, wenn Sie für' i == 2' von 's1 [2] -1' (0) zu 's1 [2]' (1). Was passiert, wenn Sie "1" von "s1" ausschließen? – Thales

@Thales Danke. –

Sie wurden gestern angewiesen, 'lower = c (.001, .001)' zu verwenden. Warum benutzen Sie das nicht für das Argument "lower"? Ändern Sie 'lower' in' c (0.01.0.01) 'und Sie sehen einen Hessian. – Bhas

optim minimiert die Funktion. Sie könnten die Likelihood-Funktion mit dem Argument rate plotten. Es braucht ein bisschen ein Fiedeln, um eine Handlung zu bekommen. Machen Sie es wie folgt aus:

z2 <- function(rate) { 
    par <- numeric(2) 
    par[1] <- .68 
    par[2] <- rate 
    y <- LLL(par) 
    y 
} 

z1 <- Vectorize(z2,vectorize.args="rate") 

curve(z1, from=.001,to=1)

Sie werden sehen, dass die Funktion für den niedrigsten Wert für rate minimal ist. Gleiches, wenn Sie from zu .1 ändern. Ich kann nicht beurteilen, ob die Schätzung gültig ist.

Quelle

2017-01-09 19:53:53 Bhas

Nochmals vielen Dank für Ihre großartige Hilfe. –

Hesse-Matrix in optim() in R

Antwort

Verwandte Themen