Поставить mle с градиентами в R - PullRequest
Новая
       65

Поставить mle с градиентами в R

0 голосов
/ 06 марта 2020

Моя цель - передать градиенты моей целевой функции log (плотности) в mle (), так как это, по-видимому, ускоряет процесс сходимости и делает его более стабильным (в соответствии с https://stats.stackexchange.com/questions/282009/defining-gradient-function-argument-in-optim-function-r).

Моя упрощенная LL, немного данных и стандартное mle следующие: 

library(stats4)
LL <- function(n,s)
{
  V = (u(z1,n)-u(z2,n))*p + u(z2,n) 
  res = zce - u.inv(V,n)
  ll = dnorm(x=res, mean=0, sd=s,log=T)
  return(-sum(ll))
}

### Functions:
u <- function(x,n) 
{
  ifelse(n!=1, util <- x^(1-n)/(1-n), util <- log(x))
  return(util)
}
u.inv <- function(x,n)
{
  ifelse(n !=1, inv.util <- ((1-n)*(x))^(1/(1-n)), inv.util <- exp(x))
  return(inv.util)
}

### Data 
z1 <- c(0.1111111, 0.1037037, 0.1222222, 0.1111111, 0.1074074, 0.1666667, 0.1333333, 0.2000000, 0.1333333, 0.1074074,
        0.1037037, 0.1111111, 0.1333333, 0.2000000, 0.1222222, 0.1111111, 0.1666667, 0.1333333, 0.1111111, 0.1333333,
        0.1111111, 0.1666667, 0.1074074, 0.1333333, 0.1222222, 0.2000000, 0.1037037)

z2 <- c(0.08888889, 0.06666667, 0.07777778, 0.00000000, 0.03333333, 0.09259259, 0.09629630, 0.08888889, 0.06666667,
        0.03333333, 0.06666667, 0.08888889, 0.06666667, 0.08888889, 0.07777778, 0.00000000, 0.09259259, 0.09629630,
        0.00000000, 0.09629630, 0.08888889, 0.09259259, 0.03333333, 0.06666667, 0.07777778, 0.08888889, 0.06666667)

p <-  c(0.5, 0.9, 0.5, 0.9, 0.9, 0.1, 0.1, 0.1, 0.5, 0.9, 0.9, 0.5, 0.5, 0.1, 0.5, 0.9, 0.1, 0.1, 0.9, 0.1, 0.5, 0.1, 0.9, 0.5, 0.5, 0.1, 0.9)

zce <- c(0.11055556, 0.10277778, 0.11000000, 0.10833333, 0.10185185, 0.11666667, 0.13240741, 0.14166667, 0.13166667,
         0.07222222, 0.08796296, 0.09944444, 0.09500000,0.10833333, 0.09444444, 0.05277778, 0.10925926, 0.11759259,
         0.05833333, 0.10277778, 0.09277778, 0.10925926, 0.06111111, 0.08833333, 0.09222222, 0.12500000, 0.09166667)

### mle()
fit <- mle(LL,
           start = list(n = 0.1,s=0.1),
           method = "L-BFGS-B",
           lower = list(n=-Inf,s=0.0001),
           upper = list(n=0.9999,s=Inf),
           control = list(maxit = 500, ndeps = rep(0.000001,2),trace= 6),
           nobs=length(z1))

Из следующей ссылки я узнал, что мне следует 1) вычислить частичные производные мою функцию логарифмирования (плотности) относительно каждого параметра, который я оцениваю, и помещаю их в дополнительную функцию «градиента», которую я затем передаю mle () с аргументом «gr»: R optim () L-BFGS- B нуждается в конечных значениях 'fn' - Вейбулла . 

### Simplify log(dnorm)
# log(dnorm) = log((1/(s*sqrt(2*pi)))*exp(-((res-mu)^2/(2*s^2))))
# log(dnorm) = log((1/(s*sqrt(2*pi)))) - (res-mu)^2/(2*s^2)
# mu = 0
# log(dnorm) = log(1/s) + log(1/sqrt(2*pi)) - res^2/(2*s^2)

### Derivative of log(dnorm), that is, -res^2/(2*s^2) with respect to n
# = -1/(2*s^2)*2*res*dres/dn 

## Derivatives of u(x,n) and u.inv(x,n) with respect to n
d.u.n <- function(x,n){ifelse(n!=1, d.util <- x^(-n), d.util <- 1/x); return(d.util)}
d.u.inv.n <- function(x,n){ifelse(n!=1, d.inv.util <- x^(1/(1-n)-1), d.inv.util <- exp(x)); return(d.inv.util)}

## Derivative of res with respect to n: 
# dres/dn = d(zce - u.inv(V,n) )/dn
# dres/dn = -d.u.inv.n(V,n)*( (d.u.n(z1,n)-d.u.n(z2,n))*p + d.u.n(z2,n))

# Derivative of -res^2/(2*s^2) with respect to n
# dres/dn = -1/(2*s^2)*2*res*-d.u.inv.n(V,n)*( (d.u.n(z1,n)-d.u.n(z2,n))*p + d.u.n(z2,n))

# Derivative of log(dnorm) with respect to s:
# log(dnorm) = log(1/s) + log(1/sqrt(2*pi)) - res^2/(2*s^2)
# d(log(dnorm))/ds = s*-s^(-2) - res^2/2*s^(-3)*(-2)
# d(log(dnorm))/ds = -1/s + res^2*s^(-3)

### adopting the method from the link
pars <- c(1,1)
LL.gr <- function(pars, zz1=z1,zz2=z2,pp=p,zcee=zce)
{
  nn <- pars[1]
  sig <- pars[2]
  V = (u(z1,nn)-u(z2,nn))*pp + u(zz2,nn) 
  res = zcee - u.inv(V,nn)

  c(sum(-1/(2*sig^2)*2*res*-d.u.inv.n(V,nn)*( (d.u.n(zz1,nn)-d.u.n(zz2,nn))*pp + d.u.n(zz2,nn))),
    sum(-1/sig + res^2*sig^(-3))
  )
}

### mle() with gradient
fit <- mle(LL,
           start = list(n = 0.1,s=0.1),
           method = "L-BFGS-B",
           lower = list(n=-Inf,s=0.001),
           upper = list(n=0.9999,s=Inf),
           control = list(maxit = 500, ndeps = rep(0.000001,2),trace= 6),
           nobs=length(z1),
           gr = LL.gr)

Запуск этого приводит к ошибке: 

N = 2, M = 5 machine precision = 2.22045e-16
L = -inf 0.001 
X0 = 0.1 0.1 
U = 0.9999 inf 
At X0, 0 variables are exactly at the bounds
At iterate     0  f=      -36.749  |proj g|=       257.81
Iteration     0

---------------- CAUCHY entered-------------------

There are 0  breakpoints

GCP found in this segment
Piece      1 f1, f2 at start point         nan         nan
Distance to the stationary point =          nan
Cauchy X =  nan nan 

---------------- exit CAUCHY----------------------

2  variables are free at GCP on iteration 1
Error in optim(start, f, method = method, hessian = TRUE, ...) : 
  non-finite value supplied by optim

Спасибо за помощь!

ОБНОВЛЕНИЕ 1

Теперь я считаю, что способ, которым я передаю градиент, действительно является правильным способом предоставления градиента функции mle (), которая отвечает на главный вопрос моего поста.

Я все еще не совсем понимаю, как передать скрытые значения параметров, то есть значения параметров на любой итерации в ходе процедуры оптимизации, в функцию градиента. Я думал о том, чтобы поместить скрытые значения параметров в переменную "par", чтобы передать их функции градиента: 

LL <- function(n,s)
{
  V = (u(z1,n)-u(z2,n))*p + u(z2,n) 
  res = zce - u.inv(V,n)
  ll = dnorm(x=res, mean=0, sd=s,log=T)
  par <- c(n,s)
  return(-sum(ll))
}

LL.gr <- function(par, zz1=z1 ,zz2=z2,pp=p,zcee=zce)
{
  n <- par[1]
  s <- par[2]
  V = (u(zz1,n)-u(zz2,n))*pp + u(zz2,n) 
  res = zcee - u.inv(V,n)
  print(s)
  c(sum(-1/(2*s^2)*2*res*-d.u.inv.n(V,n)*( (d.u.n(zz1,n)-d.u.n(zz2,n))*pp + d.u.n(zz2,n))),
    sum(-1/s + res^2*s^(-3))
  )
}

Если я немного обману, чтобы решить проблему, которая вызвала "nan", которая результат d.u.n(zz2,n) равен бесконечности для значений z2 = zz2, которые равны 0, тогда процесс сходится, но (как ни странно?) он приближается к начальным значениям, которые я установил для n и s, таким образом, 0,1 в примере. 

z2 <- c(0.08888889, 0.06666667, 0.07777778, 0.0001, 0.03333333, 0.09259259, 0.09629630, 0.08888889, 0.06666667,
            0.03333333, 0.06666667, 0.08888889, 0.06666667, 0.08888889, 0.07777778, 0.0001, 0.09259259, 0.09629630,
            0.0001, 0.09629630, 0.08888889, 0.09259259, 0.03333333, 0.06666667, 0.07777778, 0.08888889, 0.06666667)

mle(LL,
           start = list(n = 0.1,s=0.1),
           method = "L-BFGS-B",
           lower = list(n=-Inf,s=0.001),
           upper = list(n=0.9999,s=Inf),
           control = list(maxit = 500, ndeps = rep(0.0001,2),trace= 6),
           nobs=length(z1),
           gr = LL.gr)

mle(minuslogl = LL, start = list(n = 0.1, s = 0.1), method = "L-BFGS-B", 
    nobs = length(z1), lower = list(n = -Inf, s = 0.001), upper = list(n = 0.9999, 
        s = Inf), control = list(maxit = 500, ndeps = rep(1e-04, 
        2), trace = 6), gr = LL.gr)

Coefficients:
  n   s 
0.1 0.1

Кто-нибудь знает, почему я получаю начальные значения? Является ли мой способ передачи скрытых переменных правильным? Еще раз спасибо.

...