Warum stoppt dieser Optimierungsalgorithmus in R nach einigen Funktionsauswertungen?

Question

Ich habe einen Code, der für etwas Papier verwendet wurde. Nach der Definition der zu optimierenden Funktion verwendete der Autor die Nelder-Mead-Methode, um die benötigten Parameter zu schätzen. Wenn ich den Code ausführe, friert er ein, nachdem 493 Funktionsauswertungen verwendet wurden, er zeigt keine Art von Fehlermeldung oder irgendetwas anderes. Ich habe versucht, ein paar Informationen zu finden, aber ich hatte kein Glück. Wie kann ich den Befehl optim ändern, um alle möglichen Kombinationen auszuwerten und/oder was verhindert, dass die Funktion optimiert wird?

Hier ist der Code. Es ist relativ lang, aber die vorletzte Zeile (system.time(stcopfit...)) ist die EINZIGE, die ich arbeiten/reparieren/ändern muss. Sie können also einfach kopieren & den Code einfügen (wie gesagt, vom Autor des erwähnten Papiers genommen) und lassen Sie es laufen, müssen Sie nicht durch den gesamten Code gehen, nur die letzten paar Zeilen. Dies ist die data, über die die Optimierung laufen soll, d.h. eine Matrix von [0,1] gleichförmigen Variablen der Dimension 2172x9.

Jede Hilfe wird geschätzt, danke!

Hier ist ein Screenshot in RStudio (es dauerte etwa 2 Minuten bei 493, zu kommen und dann es ist, wie dies für die letzten 30 Minuten stecken):

Code:

#download older version of "sn" package 
url <- "https://cran.r-project.org/src/contrib/Archive/sn/sn_1.0-0.tar.gz" 
install.packages(url, repos=NULL, type="source") 
install.packages(signal) 
library(sn) 
library(signal) 

#1. redefine qst function 
qst <- function (p, xi = 0, omega = 1, alpha = 0, nu = Inf, tol = 1e-08) 
{ 
    if (length(alpha) > 1) 
    stop("'alpha' must be a single value") 
    if (length(nu) > 1) 
    stop("'nu' must be a single value") 
    if (nu <= 0) 
    stop("nu must be non-negative") 
    if (nu == Inf) 
    return(qsn(p, xi, omega, alpha)) 
    if (nu == 1) 
    return(qsc(p, xi, omega, alpha)) 
    if (alpha == Inf) 
    return(xi + omega * sqrt(qf(p, 1, nu))) 
    if (alpha == -Inf) 
    return(xi - omega * sqrt(qf(1 - p, 1, nu))) 
    na <- is.na(p) | (p < 0) | (p > 1) 
    abs.alpha <- abs(alpha) 
    if (alpha < 0) 
    p <- (1 - p) 
    zero <- (p == 0) 
    one <- (p == 1) 
    x <- xa <- xb <- xc <- fa <- fb <- fc <- rep(NA, length(p)) 
    nc <- rep(TRUE, length(p)) 
    nc[(na | zero | one)] <- FALSE 
    fc[!nc] <- 0 
    xa[nc] <- qt(p[nc], nu) 
    xb[nc] <- sqrt(qf(p[nc], 1, nu)) 
    fa[nc] <- pst(xa[nc], 0, 1, abs.alpha, nu) - p[nc] 
    fb[nc] <- pst(xb[nc], 0, 1, abs.alpha, nu) - p[nc] 
    regula.falsi <- FALSE 
    while (sum(nc) > 0) { 
    xc[nc] <- if (regula.falsi) 
     xb[nc] - fb[nc] * (xb[nc] - xa[nc])/(fb[nc] - fa[nc]) 
    else (xb[nc] + xa[nc])/2 
    fc[nc] <- pst(xc[nc], 0, 1, abs.alpha, nu) - p[nc] 
    pos <- (fc[nc] > 0) 
    xa[nc][!pos] <- xc[nc][!pos] 
    fa[nc][!pos] <- fc[nc][!pos] 
    xb[nc][pos] <- xc[nc][pos] 
    fb[nc][pos] <- fc[nc][pos] 
    x[nc] <- xc[nc] 
    nc[(abs(fc) < tol)] <- FALSE 
    regula.falsi <- !regula.falsi 
    } 
    x <- replace(x, zero, -Inf) 
    x <- replace(x, one, Inf) 
    Sign <- function(x) sign(x)+ as.numeric(x==0) 
    q <- as.numeric(xi + omega * Sign(alpha)* x) 
    names(q) <- names(p) 
    return(q) 
} 

#2. initial parameter setting 
mkParam <- function(Omega, delta, nu){ 
    ndim <- length(delta)+1; 
    R <- diag(ndim); 
    for (i in 2:ndim){ 
    R[i,1] <- R[1,i] <- delta[i-1]; 
    if (i>=3){for (j in 2:(i-1)){R[i,j] <- R[j,i] <- Omega[i-1,j-1];}} 
    } 
    LTR <- t(chol(R)); 
    Mtheta <- matrix(0, nrow=ndim, ncol=ndim); 
    for (i in 2:ndim){ 
    Mtheta[i,1] <- acos(LTR[i,1]); 
    cumsin <- sin(Mtheta[i,1]); 
    if (i >=3){for (j in 2:(i-1)){ 
     Mtheta[i,j] <- acos(LTR[i,j]/cumsin); 
     cumsin <- cumsin*sin(Mtheta[i,j]);} 
    } 
    } 
    c(Mtheta[lower.tri(Mtheta)], log(nu-2)); 
} 

#3. from internal to original parameters 
paramToExtCorr <- function(param){ 
    ntheta <- dim*(dim+1)/2; 
    theta <- param[1:ntheta]; 
    ndim <- (1+sqrt(1+8*length(theta)))/2; 
    LTR <- diag(ndim); 
    for (i in 2:ndim){ 
    LTR[i,1] <- cos(theta[i-1]); 
    cumsin <- sin(theta[i-1]); 
    if (i >=3){for (j in 2:(i-1)){ 
     k <- i+ndim*(j-1)-j*(j+1)/2; 
     LTR[i,j] <- cumsin*cos(theta[k]); 
     cumsin <- cumsin*sin(theta[k]);} 
    } 
    LTR[i,i] <- cumsin; 
    } 
    R <- LTR %*% t(LTR); 
    R; 
} 

#4. show estimated parameters and log likelihood 
resultVec <- function(fit){ 
    R <- paramToExtCorr(fit$par); 
    logLik <- -fit$value; 
    Omega <- R[-1, -1]; 
    delta <- R[1, -1]; 
    ntheta <- dim*(dim+1)/2; 
    nu <- exp(fit$par[ntheta+1])+2; 
    c(Omega[lower.tri(Omega)], delta, nu, logLik); 
} 

#5. negative log likelihood for multivariate skew-t copula 
stcopn11 <- function(param){ 
    N <- nrow(udat); 
    mpoints <- 150; 
    npar <- length(param); 
    nu <- exp(param[npar])+2; 
    R <- paramToExtCorr(param); 
    Omega <- R[-1, -1]; 
    delta <- R[1, -1]; 
    zeta <- delta/sqrt(1-delta*delta); 
    iOmega <- solve(Omega); 
    alpha <- iOmega %*% delta/sqrt(1-(t(delta) %*% iOmega %*% delta)[1,1]); 
    ix <- matrix(0, nrow=N, ncol=dim); 
    lm <- matrix(0, nrow=N, ncol=dim); 
    for (j in 1:dim){ 
    minx <- qst(min(udat[,j]), alpha=zeta[j], nu=nu); 
    maxx <- qst(max(udat[,j]), alpha=zeta[j], nu=nu); 
    xx <- seq(minx, maxx, length=mpoints); 
    px <- sort(pst(xx, alpha=zeta[j], nu=nu)); 
    ix[,j] <- pchip(px, xx, udat[,j]); 
    lm[,j] <- dst(ix[,j], alpha=zeta[j], nu=nu, log=TRUE); 
    } 
    lc <- dmst(ix, Omega=Omega, alpha=alpha, nu=nu, log=TRUE); 
    -sum(lc)+sum(lm) 
} 

#6. sample setting 
dim <- 9; 
smdelta <- c(-0.36,-0.33,-0.48,-0.36,-0.33,-0.48,-0.36,-0.33,-0.48); 
smdf <- 5; 
smOmega <- cor(udat); 
smzeta <- smdelta/sqrt(1-smdelta*smdelta); 
iOmega <- solve(smOmega); 
smalpha <- iOmega %*% smdelta /sqrt(1-(t(smdelta) %*% iOmega %*% smdelta)[1,1]); 

#7. estimation 
iniPar <- mkParam(diag(dim),numeric(dim),6); 
system.time(stcopfit<-optim(iniPar,stcopn11,control=list(reltol=1e-8,trace=6))); 
resultVec(stcopfit); 

Answer 1

Die Parameter, die Sie bei Schritt 493 erreichen, führen zu einer Endlosschleife in Ihrer qst Funktion: ich habe keine Ahnung, was dieser sehr komplexe Code tatsächlich macht, ich fürchte, ich kann nicht weiter diagnostizieren. Hier ist, was ich tat, so weit kommen:

• Ich erklärte cur.params <- NULL im globalen Umfeld, legte dann cur.params <<- params innerhalb stcopn11; Dadurch wird der aktuelle Parametersatz in der globalen Umgebung gespeichert. Wenn Sie den optim()-Aufruf manuell (über Control-C oder ESC je nach Plattform) unterbrechen, können Sie den aktuellen Parametersatz überprüfen und problemlos von dort aus neu starten
• ich Altschule Debugging-Anweisungen setzen in (zB cat("entering stcopn11\n") und cat("leaving stcopn11\n") am Anfang und am nächsten zu letzten Zeile der Zielfunktion, um nur einige innerhalb stopc11 Fortschrittsindikatoren innerhalb anzuzeigen)
• einmal hatte ich das " schlechte "Parameter ich debug(stcopn11) und stcopn11(cur.param) verwendet, um durch die Funktion
• zu gehen Ich entdeckte, dass es auf Dimension 3 (j==3 in th hing e for Schleife innerhalb stcopn11) und vor allem auf der ersten qst() Anruf
• Ich fügte ein maxit=1e5 Argument zu qst; initialisiert it <- 1 vor der while Schleife; set it <- it+1 jedes Mal durch die Schleife; änderte das Stoppkriterium auf while (sum(nc) > 0 && it<maxit); und hinzugefügt if (it==maxit) stop("hit max number of iterations in qst") direkt nach der Schleife

1E5 Iterationen in qst 74 Sekunden dauerten; Ich habe keine Ahnung, ob es irgendwann aufhören könnte, aber ich wollte nicht warten, um es herauszufinden.

Dies war meine modifizierte Version von stcopn11:

cur.param <- NULL ## set parameter placeholder 

##5. negative log likelihood for multivariate skew-t copula 
stcopn11 <- function(param,debug=FALSE) { 
    cat("stcopn11\n") 
    cur.param <<- param ## record current params outside function 
    N <- nrow(udat) 
    mpoints <- 150 
    npar <- length(param) 
    nu <- exp(param[npar])+2 
    R <- paramToExtCorr(param) 
    Omega <- R[-1, -1] 
    delta <- R[1, -1] 
    zeta <- delta/sqrt(1-delta*delta) 
    cat("... solving iOmega") 
    iOmega <- solve(Omega) 
    alpha <- iOmega %*% delta/
      sqrt(1-(t(delta) %*% iOmega %*% delta)[1,1]) 
    ix <- matrix(0, nrow=N, ncol=dim) 
    lm <- matrix(0, nrow=N, ncol=dim) 
    cat("... entering dim loop\n") 
    for (j in 1:dim){ 
     if (debug) cat(j,"\n") 
     minx <- qst(min(udat[,j]), alpha=zeta[j], nu=nu) 
     maxx <- qst(max(udat[,j]), alpha=zeta[j], nu=nu) 
     xx <- seq(minx, maxx, length=mpoints) 
     px <- sort(pst(xx, alpha=zeta[j], nu=nu)) 
     ix[,j] <- pchip(px, xx, udat[,j]) 
     lm[,j] <- dst(ix[,j], alpha=zeta[j], nu=nu, log=TRUE) 
    } 
    lc <- dmst(ix, Omega=Omega, alpha=alpha, nu=nu, log=TRUE) 
    cat("leaving stcopn11\n") 
    -sum(lc)+sum(lm) 
}