optimering - Problem med 64 bit R er optimalt under Windows 7

Indlæg af Hanne Mølgaard Plasc

Problem



Jeg er i færd med at udvikle min første R-pakke, som formodes at passe til multinomial processing tree (MPT) modeller (se hjemmeside for den aktuelle version). Modelleringen opnås ved R 's optim funktion.

I dag spillede jeg med det for første gang på en Windows 7-maskine og bemærkede noget virkelig mærkeligt: ​​optim konvergerer ikke succesfuldt, når du bruger 64 bitversionen af ​​R. Dette ligner en fejl for mig (især som nlminb konvergerer for begge R versioner). Som optim er kernen i min pakke, er enhver hjælp på dette spørgsmål meget værdsat. [18]


Her kommer et minimalt reproducerbart eksempel (normalt er modellen angivet via udtryk og ikke angivet i objektivfunktionen, men for enkelhed sætter jeg alt i objektivfunktionen her):


# The objective function:
llk.tree <- function (Q, data) 
{
    p <- Q[1]
    q <- Q[2]
    r <- Q[3]
    e <- rep(NA,4)
    e[1] <- p * q * r
    e[2] <- p * q * (1-r)
    e[3] <- p * (1-q) * r
    e[4] <- p * (1-q) * (1-r) + (1-p)

    llk <- sum(data * log(e))
    if (is.na(llk)) 
        llk <- -1e+19
    if (llk == -Inf) 
        llk <- -1e+19
    return(-llk)
}

# The call to optim:
optim(runif(3, 0.1, 0.9), llk.tree, data = c(24, 65, 30, 61), method = "L-BFGS-B", lower = rep(0, 3), upper = rep(1, 3))


Dette eksempel gengiver et eksempel fra et seminal papir på MPT'er af Riefer & Batchelder, nemlig række 1 i tabel 1 s. 327 (forventede parameterværdier ville være p=1, q=.49 og r=.30). [19]


Kører det på en 32 bit R altid giver det rigtige resultat (prøvet med versioner 2.12.2 og 2.13.0):


$par
[1] 1.0000000 0.4944449 0.3000001

$value
[1] 234.7110

$counts
function gradient 
      11       11 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL\_REDUCTION\_OF\_F <= FACTR*EPSMCH"


(Bemærk at tællingen kan afvige på grund af tilfældige startværdier.)


At køre det på en 64 bit R på den anden side kan frembringe et sådant (forkert) resultat:


$par
[1] 0.8668081 0.6326655 0.1433857

$value
[1] 257.7328

$counts
function gradient 
       3        3 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL\_REDUCTION\_OF\_F <= FACTR*EPSMCH"


Den returnerede værdi af objektivfunktionen og de returnerede parameterværdier afviger hver gang, men tællingen er altid 3!


Bemærk at kørsel nlminb producerer korrekte resultater på 32 bit og 64 bit R:


nlminb(runif(3, 0.1, 0.9), llk.tree, data = c(24, 65, 30, 61), lower = 0, upper = 1)

$par
[1] 1.0000000 0.4944445 0.3000000

$objective
[1] 234.711

$convergence
[1] 0

$iterations
[1] 14

$evaluations
function gradient 
      19       55 

$message
[1] "relative convergence (4)"


En sidste note: Vi har eksempler (dette er vores enkleste eksempelmodel), der fungerede på 64 bit R og optim, men flere eksempler (som den her viste) fungerede ikke.


Og tælle er altid 3 ...


EDIT:


Når man fastsætter startværdierne (takket være Joshua Ulrich) optim, bevæger man sig ikke væk fra de faste værdier under 64 bit R:


optim(c(0.5, 0.5, 0.5), llk.tree, data = c(24, 65, 30, 61), method = "L-BFGS-B", lower = rep(0, 3), upper = rep(1, 3))

$par
[1] 0.5 0.5 0.5

$value
[1] 276.1238

$counts
function gradient 
       3        3 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL\_REDUCTION\_OF\_F <= FACTR*EPSMCH"

Bedste reference


Vi gjorde nogle flere test i dag og fandt det samme problem som beskrevet i spørgsmålet under Linux ved hjælp af 64 bit R.


Men takket være Joachim Vandekerckhove for denne geniale idé forsøgte vi en simpel ændring, der løste problemet (selv om spørgsmålet er mistænkt). Ved afslutningen af ​​objektivfunktionen, hvis llk er Inf, satte vi det til en ekstremt høj værdi (var 1e19.


Ved at bruge en mindre værdi (f.eks. 1e10) fjerner problemet på 64 bit-maskiner (hidtil testet på Linux): [20]


llk.tree <- function (Q, data) 
{
    p <- Q[1]
    q <- Q[2]
    r <- Q[3]
    e <- rep(NA,4)
    e[1] <- p * q * r
    e[2] <- p * q * (1-r)
    e[3] <- p * (1-q) * r
    e[4] <- p * (1-q) * (1-r) + (1-p)

    llk <- sum(data * log(e))
    if (is.na(llk)) 
        llk <- -1e+10
    if (llk == -Inf) 
        llk <- -1e+10
    return(-llk)
}

# The call to optim:
optim(runif(3, 0.1, 0.9), llk.tree, data = c(24, 65, 30, 61), method = "L-BFGS-B", lower = rep(0, 3), upper = rep(1, 3))


Dette giver det rigtige resultat!