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

Question

Я пытаюсь оценить линейную модель с логарифмически нормальным распределенным сроком ошибки. У меня уже есть рабочий код для линейной модели с нормально распределенными ошибками:

library(Ecdat) 
library(assertthat) 
library(maxLik) 

# Load the data 
data(Wages1) 

# Check what R says 
summary(lm(wage ~ school + exper + sex, data = Wages1)) 


# Use maxLik from package maxLik 
# The likelihood function 
my_log_lik_pos <- function(theta, data){ 
    y <- data[, 1] 
    x <- data[, -1] 
    beta <- head(theta, -1) 
    sigma <- tail(theta, 1) 
    xb <- x%*%beta 
    are_equal(dim(xb), c(nrow(my_data), 1)) 
    return(sum(log(dnorm(y, mean = xb, sd = sigma)))) 
} 

# Bind the data 
my_data <- cbind(Wages1$wage, 1, Wages1$school, Wages1$exper, Wages1$sex) 

my_problem <- maxLik(my_log_lik_pos, data = my_data, 
       start = rep(1,5), method = "BFGS") 

summary(my_problem) 

Получаю примерно одинаковые результаты. Теперь я пытаюсь сделать то же самое, но используя log-normal вероятность. Для этого я должен первым имитировать некоторые данные:

true_beta <- c(0.1, 0.2, 0.3, 0.4, 0.5) 

ys <- my_data[, -1] %*% head(true_beta, -1) + 
     rlnorm(nrow(my_data), 0, tail(true_beta, 1)) 

my_data_2 <- cbind(ys, my_data[, -1]) 

И логарифмической функции правдоподобия:

my_log_lik_lognorm <- function(theta, data){ 
    y <- data[, 1] 
    x <- data[, -1] 
    beta <- head(theta, -1) 
    sigma <- tail(theta, 1) 
    xb <- x%*%beta 
    are_equal(dim(xb), c(nrow(data), 1)) 
    return(sum(log(dlnorm(y, mean = xb, sd = sigma)))) 
} 

my_problem2 <- maxLik(my_log_lik_lognorm, data = my_data_2, 
       start = rep(0.2,5), method = "BFGS") 

summary(my_problem2) 

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

Правильно ли я считать, что лог-правдоподобия логарифмически нормального распределения:

sum(log(dlnorm(y, mean = .., sd = ...)) 

Если я не ошибаюсь, то это определение логарифмической вероятности (сумма логарифмов плотности).

Answer 1

Я нашел проблему: кажется, проблема не в моей функции логарифмического правдоподобия. Когда я пытаюсь оценить модель с GLM:

summary(glm(ys ~ school + exper + sex, family=gaussian(link="log"), data=Wages1)) 

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

ys <- my_data[, -1] %*% head(true_beta, -1) + 
      rlnorm(nrow(my_data), 0, tail(true_beta, 1)) 

Правильный способ моделирования данных:

ys <- rlnorm(nrow(my_data), my_data[, -1] %*% head(true_beta, -1), tail(true_beta, 1)) 

Теперь все работает!