Написание правильного нормального логарифмического правдоподобия в R

У меня есть проблема относительно следующей модели,

введите здесь описание изображения

где я хочу сделать вывод о и тау, u - известный вектор, а x - вектор данных. Логарифмическая вероятность

введите здесь описание изображения

У меня проблема с записью логарифмического правдоподобия в R.

x <- c(3.3569,1.9247,3.6156,1.8446,2.2196,6.8194,2.0820,4.1293,0.3609,2.6197)
mu <- seq(0,10,length=1000)

normal.lik1<-function(theta,x){ 
  u <- c(1,3,0.5,0.2,2,1.7,0.4,1.2,1.1,0.7)  
  mu<-theta[1] 
  tau<-theta[2] 
  n<-length(x) 

logl <-  sapply(c(mu,tau),function(mu,tau){logl<- -0.5*n*log(2*pi) -0.5*n*log(tau^2+u^2)- (1/(2*tau^2+u^2))*sum((x-mu)^2) } )

  return(logl) 
  } 

#test if it works for mu=1, tau=2
head(normal.lik1(c(1,2),x))
#Does not work..

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

theta.hat<-optim(c(1,1),loglike2,control=list(fnscale=-1),x=x,,method="BFGS")$par

Но это не работает. Любые предложения о том, как я могу написать вероятность?


r statistics log-likelihood normal-distribution mle
person Hans Christensen    schedule 09.01.2018    source источник
comment
sapply перебирает c(mu, tau). Вам не нужно sapply здесь. logl <- logl<- -0.5*n*log(2*pi) -0.5*n*log(tau^2+u^2)- (1/(2*tau^2+u^2))*sum((x-mu)^2)   -  person LyzandeR    schedule 09.01.2018
comment
Итак, как вы предлагаете мне записать эту вероятность? и построить его, например, по значениям mu? :)   -  person Hans Christensen    schedule 09.01.2018
comment
(я исправил вероятность)   -  person Hans Christensen    schedule 09.01.2018
comment
Какие оптимальные значения вы получаете для тау и мю?   -  person Hans Christensen    schedule 09.01.2018

Ответы (1)


arrow_upward
2
arrow_downward

Во-первых, как уже упоминалось в комментариях к вашему вопросу, нет необходимости использовать sapply(). Вы можете просто использовать sum() — так же, как в формуле logLikelihood.

Я изменил эту часть в normal.lik1() и умножил выражение, присвоенное logl, на минус 1, чтобы функция вычислила минус logLikelihood. Вы хотите найти минимум по тета, поскольку функция возвращает положительные значения.

x < c(3.3569,1.9247,3.6156,1.8446,2.2196,6.8194,2.0820,4.1293,0.3609,2.6197)
u <- c(1,3,0.5,0.2,2,1.7,0.4,1.2,1.1,0.7) 

normal.lik1 <- function(theta,x,u){ 
  mu <- theta[1] 
  tau <- theta[2] 
  n <- length(x) 
  logl <- - n/2 * log(2*pi) - 1/2 * sum(log(tau^2+u^2)) - 1/2 * sum((x-mu)^2/(tau^2+u^2))
  return(-logl) 
}

Это можно сделать с помощью nlm(), например

nlm(normal.lik1, c(0,1), hessian=TRUE, x=x,u=u)$estimate

где c(0,1) — начальные значения алгоритма.

Чтобы построить logLikelihood для диапазона значений mu и некоторого фиксированного tau, вы можете настроить функцию так, чтобы mu и tau были отдельными числовыми аргументами.

normal.lik2 <- function(mu,tau,x,u){ 
  n <- length(x) 
  logl <- - n/2 * log(2*pi) - 1/2 * sum(log(tau^2+u^2)) - 1/2 * sum((x-mu)^2/(tau^2+u^2))
  return(logl) 
}

Затем определите некоторый диапазон для mu, вычислите логарифмическую вероятность и используйте plot().

range.mu <- seq(-10,20,0.1)

loglik <- sapply(range.mu, function(m) normal.lik2(mu=m,tau=2,x=x,u=u))

plot(range.mu, loglik, type = "l")

введите здесь описание изображения

Я уверен, что есть более элегантные способы сделать это, но это помогает.

person Martin C. Arnold    schedule 09.01.2018
comment
имеет смысл! Если бы я изменил -logL на logL в normal.lik1, тогда бы: optim(c(1,1),normal.lik1,control=list(fnscale=-1),x=x, u=u,method =BFGS)$par делает то же самое, что и nlm? - person Hans Christensen; 10.01.2018
comment
По общему признанию, я не очень хорошо знаком с optim() и nlm(), не говоря уже об используемых алгоритмах, но для имеющихся данных оба дают практически одинаковые оценки (mu = 2,1863, tau = 0). - person Martin C. Arnold; 10.01.2018
comment
Я знаю, что меня беспокоит, так это то, что я не думаю, что тау должен быть равен 0, на самом деле я предполагаю, что он должен быть между 0,5 и 3. Исходная модель X_i=mu+lambda_i+epsilon_i. где лямбда — это N (0, тау ^ 2), а epsilon_i — это N (0, u ^ 2_i). Поэтому я предположил, что X - это N (mu, lambda ^ 2 + u ^ 2). Это означает, что тау не может быть отрицательным, и я бы предположил, что MLE для него должен принимать значение от 0 до 5... возможно, я ошибся в вероятности? - person Hans Christensen; 10.01.2018
comment
Я получаю те же мю и тау при использовании optim BTW. - person Hans Christensen; 10.01.2018
comment
Ваш вывод logLik правильный. Однако в normal.lik1() была другая ошибка (второе слагаемое, умноженное на n). Я отредактировал свой ответ. Теперь я получаю (mu=2,6525, tau=0,7692), используя nlm(). - person Martin C. Arnold; 10.01.2018