Почему этот код не оптимизирован для всех трех пунктов? - PullRequest
Новая
       11

Почему этот код не оптимизирован для всех трех пунктов?

3 голосов
/ 06 декабря 2010

Фон

Я пытаюсь приспособить распределение к 95% CI и режиму.Используемая мной функция стоимости решает три функции для 0: P (X = 2,5 | mu, сигма) = 0,025, P (X = 7,5 | mu, сигма) = 0,975 и режим log-N (mu, сигма).) = 3.3.примечание: режим логарифмического режима: = $ e ^ {\ mu- \ sigma ^ 2)} $:

Подход

Сначала я пишу функцию стоимости, prior

prior <- function(parms) {
  a <- abs(plnorm(2.5, parms[1], parms[2]) - 0.025)
  b <- abs(plnorm(7.5, parms[1], parms[2]) - 0.975)
  mode <- exp(parms[1] - parms[2]^2)
  c <- abs(mode-3.3)
  return(a + b + c)
}

И затем я ищу параметры, которые минимизируют функцию стоимости 

v = nlm(prior,c(log(3.3),0.14))

Очевидно, что функция максимизируется для режима LCL, но не для UCL. 

abs(plnorm(7.5, parms[1], parms[2]) - 0.975)
> [1] 0.02499989

Вот график с пунктирными линиями при желаемом 95% ДИ: 

x <- seq(0,10,0.1)
plot(x,dlnorm(x, v$estimate[1],v$estimate[2]),type='l')
abline(v=c(2.5,7.5), lty=2) #95%CI

Вопрос

Оптимизация по двум точкам и всеошибки в третьем.Тем не менее, я хотел бы, чтобы он соответствовал точкам равномерно.

Как я могу получить функцию, чтобы придать равный вес величине a , b , и c условия?Похоже, что функция подходит только a и c .

примечание: Этот вопрос основан на вопросе, заданном в [перекрестной проверке] [1], за исключением того, что эта версия конкретно относится к функции алгоритма оптимизации R nlm (), тогда как вопрос CVо поиске более подходящего дистрибутива.

Ответы [ 2 ]

1 голос
/ 14 декабря 2010

Причина, по которой ваша оптимизация "не работает", заключается в том, что шкала трех параметров, a, b и c, не совпадает.a и b измеряют разницу в вероятностях и всегда могут быть установлены не более 0,025 путем выбора действительно малого значения для стандартного отклонения (parms[2]), поскольку тогда plnorm(2.5, parms[1], parms[2]) будет равно 0 (то же самое)за 7,5).То же количество ошибок (0,025) было бы незаметным для c - это несоответствие масштабирования.

Вы можете переписать свою функцию оптимизации так, чтобы ошибки измерялись по шкале x для всех трех критериев.сравнивая квантили с 2,5 и 7,5: 

prior2 <- function(parms) { 
  a <- abs(qlnorm(0.025, parms[1], parms[2]) - 2.5) 
  b <- abs(qlnorm(0.975, parms[1], parms[2]) - 7.5) 
  mode <- exp(parms[1] - parms[2]^2) 
  c <- abs(mode-3.3) 
  return(a + b + c) 
}

Это похоже на то, что предложил Рамнатх, за исключением того, что не в логарифмическом масштабе.Этот подход не очень хорошо работает на левом хвосте, потому что распределение искажено вправо: небольшие изменения в расположении нижнего 2,5-го процентиля приводят к большим изменениям процентиля на 2,5, в то время как на 7,5 это не так.Предложение Рамната работать с логарифмической шкалой решает эту проблему, поскольку логарифмически нормальное распределение симметрично на логарифмической шкале.

Еще один способ улучшить ваше соответствие - изменить критерий оптимизации.Прямо сейчас вы минимизируете среднюю абсолютную ошибку.Это означает, что одна большая ошибка нормальна, если два других условия ошибки действительно малы.Вы можете наложить больший штраф на большие ошибки, уменьшив вместо этого среднеквадратическую ошибку (a^2+b^2_c^2).Эта последняя версия (в логарифмическом масштабе) дает наиболее привлекательную оценку с моей точки зрения. 

prior3 <- function(parms) { 
  a <- abs(parms[1] - 1.96*parms[2] - log(2.5))
  b <- abs(parms[1] + 1.96*parms[2] - log(7.5)) 
  c <- abs(parms[1] - parms[2]^2 - log(3.3)) 
  return(a^2 + b^2 + c^2) 
}
1 голос
/ 06 декабря 2010

попробуйте альтернативную формулировку вашей функции оптимизации.лог 95% доверительного интервала для логнормального распределения дается мю - 2 * сигма и мю + 2 * сигма.поэтому вы можете попытаться свести к минимуму abs (mu - 2 * sigma - log (2.5)) + abs (mu + 2 * sigma - log (7.5)) + abs (mu - sigma ^ 2 - log (3.3)).

Когда я минимизировал это, я обнаружил, что доверительные интервалы очень близки, в то время как режим немного отключен.в зависимости от характера вашей заявки, вы можете по-разному оценивать три термина

...