R-код ручной прокатки для Poisson MLE

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

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

n <- 10000 # sample size
b0 <- 1.0 # intercept
b1 <- 0.2 # coefficient
x <- runif(n=n, min=0, max=1.5) # generate covariate values
lp <- b0+b1*x # linear predictor
lambda <- exp(lp) # compute lamda
y <- rpois(n=n, lambda=lambda) # generate y-values
dta <- data.frame(y=y, x=x) # generate dataset
negloglike <- function(lambda) {n*lambda-sum(x)*log(lambda) + sum(log(factorial(y)))} # build negative log-likelihood
starting.vals <- c(0,0) # one starting value for each parameter
pars <- c(b0, b1)
maxLike <- optim(par=pars,fn=negloglike, data = dta) # optimize

Мой вывод R, когда я ввожу maxLike, выглядит следующим образом:

Error in fn(par, ...) : unused argument (data = list(y = c(2, 4....

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


r poisson glm mle
person Austin    schedule 04.10.2014    source источник
comment
ты не показываешь свое определение negloglike   -  person Dason    schedule 04.10.2014
comment
@Dason, извиняюсь, это уже там.   -  person Austin    schedule 04.10.2014
comment
Вам нужно negloglike <- function(lambda) {n*lambda-sum(data$x)*log(lambda) + sum(log(factorial(data$y)))}. Рассмотрим также dpois(...,log=TRUE) и bbmle::mle2.   -  person Ben Bolker    schedule 04.10.2014
comment
@BenBolker - вы, вероятно, захотите добавить туда параметр данных   -  person Dason    schedule 04.10.2014
comment
вас также может заинтересовать ms.mcmaster.ca/~bolker/emdbook/chap6A. pdf   -  person Ben Bolker    schedule 04.10.2014
comment
@BenBolker Вы уверены, что доверяете этому автору? ;)   -  person Dason    schedule 04.10.2014
comment
@BenBolker - спасибо за ссылку, а также за то, что ваши исследования доступны для социологов!   -  person Austin    schedule 04.10.2014

Ответы (2)


arrow_upward
4
arrow_downward

optim может использовать вашу функцию только определенным образом. Предполагается, что первый параметр в вашей функции принимает параметры в виде вектора. Если вам нужно передать в эту функцию другую информацию (в вашем случае данные), вам необходимо указать ее в качестве параметра вашей функции. Ваша функция negloglike не имеет параметра data, и именно на это она жалуется. То, как вы это закодировали, вам не нужно, поэтому вы, вероятно, могли бы решить свою проблему, просто удалив часть данных = dat вашего вызова для оптимизации, но я не проверял это. Вот небольшой пример выполнения простого MLE только для пуассона (не glm)

negloglike_pois <- function(par, data){
  x <- data$x
  lambda <- par[1]

  -sum(dpois(x, lambda, log = TRUE))
}

dat <- data.frame(x = rpois(30, 5))
optim(par = 4, fn = negloglike_pois, data = dat)
mean(dat$x)

> optim(par = 4, fn = negloglike_pois, data = dat)
$par
[1] 4.833594

$value
[1] 65.7394

$counts
function gradient 
      22       NA 

$convergence
[1] 0

$message
NULL

Warning message:
In optim(par = 4, fn = negloglike_pois, data = dat) :
  one-dimensional optimization by Nelder-Mead is unreliable:
use "Brent" or optimize() directly
> # The "true" MLE. We didn't hit it exactly but came really close
> mean(dat$x)
[1] 4.833333
person Dason    schedule 04.10.2014
comment
Спасибо, @Dason. Как бы это выглядело, если бы я захотел параметризовать лямбду, чтобы понять, как это работает в среде GLM? - Остин - person Austin; 05.10.2014
comment
Сделайте так, чтобы ваша неглогподобная функция имела первый параметр, принимающий вектор, содержащий b0 и b1, а затем вытягивайте эти части в верхней части функции (например, как я вытаскиваю лямбда из номинала в верхней части примера) и используйте их, как вы хотите после этого. - person Dason; 05.10.2014

arrow_upward
0
arrow_downward

Реализация комментариев из ответа Дейсона довольно проста, но на всякий случай:

library("data.table")

d <- data.table(id = as.character(1:100), 
                x1 = runif(100, 0, 1),
                x2 = runif(100, 0, 1))

#' the assumption is that lambda can be written as
#' log(lambda) = b1*x1 + b2*x2 
#' (In addition, could add a random component)
d[, mean := exp( 1.57*x1 + 5.86*x2 )]
#' draw a y for each of the observations
#' (rpois is not vectorized, need to use sapply)
d[, y := sapply(mean, function(x)rpois(1,x)) ]

negloglike_pois <- function(par, data){
  data <- copy(d)
  # update estimate of the mean
  data[, mean_tmp := exp( par[1]*x1 + par[2]*x2 )]
  # calculate the contribution of each observation to the likelihood
  data[, log_p := dpois(y, mean_tmp, log = T)]
  #' Now we can sum up the probabilities
  data[, -sum(log_p)]
}

optim(par = c(1,1), fn = negloglike_pois, data = d)
$par
[1] 1.554759 5.872219

$value
[1] 317.8094

$counts
function gradient 
      95       NA 

$convergence
[1] 0

$message
NULL
person desval    schedule 08.04.2020