MLE не сходится с «прописанным» dnorm в R - PullRequest
Новая
       20

MLE не сходится с «прописанным» dnorm в R

1 голос
/ 05 марта 2020

Я оцениваю параметры модели, используя функцию логарифмического правдоподобия. Для стандартной функции нормальной плотности я однажды использую встроенную функцию «dnorm» и однажды сам определяю эту функцию. Странно, но использование dnorm приводит к конвергенции, тогда как другой метод этого не делает: 

### 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)
}

v = function(x,n){return(1/(u(maxz,n)-u(minz,n))*(u(x,n)-u(minz,n)))}
v.inv = function(x,n){return(u.inv(x*(u(maxz,n)-u(minz,n))+u(minz,n),n))}

w <- function(p,a,b){return(exp(-b*(-log(p))^(1-a)))}

### 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)

maxz = 135
minz = 0

### Using dnorm:

LL <- function(n,a,b,s)
{
  V = (v(z1,n)-v(z2,n))*w(p,a,b) + v(z2,n) 
  res = zce - v.inv(V,n)
  ll = dnorm(res, 0, s,log=T)
  return(-sum(ll))
}

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

### Resulting coefficients saved in "fit"

Coefficients:
         n          a          b          s 
0.16533414 0.65254314 0.78727084 0.01475997

Теперь, используя записанный журнал стандартного нормаля вместо dnorm (..., log = T): 

ldens <- function(x,mu,sig){log((1/(sig*sqrt(2*pi)))*exp(-((x-mu)^2/(2*sig^2))))}

LL.ldens <- function(n,a,b,s)
{
  V = (v(z1,n)-v(z2,n))*w(p,a,b) + v(z2,n) 
  res = zce - v.inv(V,n)
  ll = ldens(x= res, mu=0, sig = s)
  return(-sum(ll))
}

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

, которое выдает сообщение «ошибка конечных значений»: 

Error in optim(start, f, method = method, hessian = TRUE, ...) : 
  L-BFGS-B needs finite values of 'fn'

Дело в том, что я не понимаю, почему. Если я возьму начальные значения для получения первого вектора «res», который будет использовать mle, я получу вектор плотностей журналов, используя собственную спецификацию. Более того, похоже, что это соответствует вектору, который я получаю при использовании dnorm (... log = T): 

n = a = b = s = 0.1
V = (v(z1,n)-v(z2,n))*w(p,a,b) + v(z2,n) 
res = zce - v.inv(V,n)

ldens(x= res, mu=0, sig = s)

[1] 1.383596 1.383637 1.379527 1.383579 1.382617 1.320485 1.381703 1.317098 1.383168 1.325277 1.372026 1.378537 1.327294 1.139934 1.353307 1.222810 1.291415 [18] 1.379966 1.252776 1.356281 1.369575 1.291415 1.281141 1.302690 1.347586 1.242405 1.376986

dnorm(res, 0, s, log=T)

[1] 1.383596 1.383637 1.379527 1.383579 1.382617 1.320485 1.381703 1.317098 1.383168 1.325277 1.372026 1.378537 1.327294 1.139934 1.353307 1.222810 1.291415 [18] 1.379966 1.252776 1.356281 1.369575 1.291415 1.281141 1.302690 1.347586 1.242405 1.376986

Интересно, что эти числа не совпадают при проверке равенства с "==" (кроме для одного): 

ldens(x= res,mu=0,sig = s) == dnorm(res, 0, s,log=T)

[1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE [26]  TRUE FALSE

Проверка чисел с более высокой точностью показывает, что они немного отличаются от других: 

sprintf("%.54f",ldens(res,0,s))[1]
"1.383596381246589235303190434933640062808990478515625000"
sprintf("%.54f",dnorm(res, 0, s,log=T))[1]
"1.383596381246589013258585509902331978082656860351562500"

Но это не может быть причиной того, что использование dnorm приводит к сходимости а другой нет?

1 Ответ

1 голос
/ 05 марта 2020

Изменяя вашу функцию для обнаружения ошибок с помощью: 

ldens <- function(x,mu,sig){v <- log((1/(sig*sqrt(2*pi)))*exp(-((x-mu)^2/(2*sig^2)))); if(is.infinite(sum(v))) browser(); v}

вы сможете увидеть значения параметров, которые вызывают проблему - экспоненциальная часть равна нулю, экспоненциальная оценка равна 0, и поэтому журнал возвращает -Inf. Внутренняя функция dnorm, вероятно, использует математически эквивалентную версию экспоненциального распределения, которая имеет лучшие свойства для арифметики с плавающей запятой c.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.

Похожие темы

...