Почему в mle2 возвращает subset() ошибку, а подмножество работает вручную

Question

Обновление: для минимального примера прокрутите вниз (не воспроизводит ту же ошибку, но разные результаты)

перед моим настоящим вопросом предостережение: я новичок в R, поэтому надеюсь, что в этом есть смысл. Я попытался выяснить проблему, используя Google и функции отладки R (хотя я все еще не слишком знаком с ними).

Я прочитал среди других:

Вопрос о том, когда использовать% в% и когда == во время поднабора

подмножество и адресация строк и столбцов через []

который привел меня в очень хороший блог Хэдли Уикхэма, где я буду дополнительно изучать возможности отладки R, его статьи на тему: Нестандартная оценка

Моя проблема

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

Я хочу использовать функцию mle2 из пакета bbmle, которая является оберткой для функции optim, чтобы оценить модель для структуры терминов, см. Код ниже. Тем не менее, является ли предположение о нормально распределенных остатках оправданным или нет, тем не менее остается открытым вопрос:

Подход 1: поднабор "вручную" работает просто отлично:

require("bbmle")

lifeBBB <- data$life[data$Rating == "BBB"]
yieldBBB <- data$yield[data$Rating == "BBB"]



LL <- function(b0, b1, b2, lambda, mu, sigma){ 

Score= yieldBBB - (b0+b1*((1-exp(-lambda*lifeBBB))/(lambda*lifeBBB))+ b2*((1-exp(-lambda*lifeBBB))/(lambda*lifeBBB)-exp(-lambda*lifeBBB)))

    Score = suppressWarnings(dnorm(Score, mu, sigma, log=TRUE))

    -sum(Score)

}

fit <- mle2(LL, start = list(b0 = 108, b1 = -85, b2=168, lambda= 0.12,  mu = -10, sigma=70),control=list(maxit = 5000))'

Call:
mle2(minuslogl = LL, start = list(b0 = 108, b1 = -85, b2 = 168, 
    lambda = 0.12, mu = -10, sigma = 70), control = list(maxit = 5000))

Coefficients:
         b0          b1          b2      lambda          mu 
110.2355968 -82.7010072 167.5960478   0.1410541  -7.7644032 
      sigma 
 69.4846302 

Log-likelihood: -1018.8 

Подход 2: Когда я пытаюсь вызвать подмножество из функции mle2, возникает ошибка:

LL2 <- function(b0, b1, b2, lmbd, mu, sigma){

    Score= yield - (b0+b1*((1-exp(-lmbd*life))/(lmbd*life))+ b2*((1-exp(-lmbd*life))/(lmbd*life)-exp(-lmbd*life)))

    Score = suppressWarnings(dnorm(Score, mu, sigma, log=TRUE)) 

    -sum(Score)    
}


 fit1 <- mle2(minuslogl=LL2, start = list(b0 = 108, b1 = -85, b2=168, lmbd= 0.12,  mu = -10, sigma=70),method="BFGS", data=data, subset= Rating=="BBB", control=list(maxit=5000))

Error in optim(par = c(108, -85, 168, 0.12, -10, 70), fn = function (p)  : 
  initial value in 'vmmin' is not finite

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

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

заранее спасибо

Б. Лёр

По запросу я попытался придумать минимальный пример, который использует примерные данные. Минимальный пример не выдает ошибку, а представляет собой два разных набора оценок. Это своеобразно, так как в моем понимании обе функции должны делать то же самое

#### minmial example

require("bbmle")
## approach 1:

set.seed(23456)

Rating <- c(rep("A",38),rep("BBB",39) )
yield <- rnorm(77,0.5,15)
life <- runif(77,1,20)

exdata <- data.frame(Rating,yield,life)


lifeBBB <- exdata$life[exdata$Rating == "BBB"]
yieldBBB <- exdata$yield[exdata$Rating == "BBB"]

LL <- function(b0, b1, b2, lambda, mu, sigma){

    Score= yieldBBB - (b0+b1*((1-exp(-lambda*lifeBBB))/(lambda*lifeBBB))+ b2*((1-exp(-lambda*lifeBBB))/(lambda*lifeBBB)-exp(-lambda*lifeBBB)))

    Score = suppressWarnings(dnorm(Score, mu, sigma, log=TRUE))

    -sum(Score)

}

fit <- mle2(LL, start = list(b0 = 100, b1 = -80, b2=160, lambda= 0.12,  mu = 1, sigma=14),method ="BFGS",control=list(maxit = 5000)) 


### approach 2:


LL2 <- function(b0, b1, b2, lmbd, mu, sigma){

    Score= exdata$yield - (b0+b1*((1-exp(-lmbd*exdata$life))/(lmbd*exdata$life))+ b2*((1-exp(-lmbd*exdata$life))/(lmbd*exdata$life)-
                                                                                          exp(-lmbd*exdata$life)))

    Score = suppressWarnings(dnorm(Score, mu, sigma, log=TRUE)) ## assumption is that residuals are normally distributed

    -sum(Score)    
}


fit1 <- mle2(minuslogl=LL2, start = list(b0 = 100, b1 = -80, b2=160, lmbd= 0.12,  mu = 1, sigma=14),method="BFGS", data=exdata, subset= Rating=="BBB", control=list(maxit=5000))

Call(fit):

Call:
mle2(minuslogl = LL, start = list(b0 = 100, b1 = -80, b2 = 160, 
    lambda = 0.12, mu = 1, sigma = 14), method = "BFGS", control = list(maxit = 5000))

Coefficients:
          b0           b1           b2       lambda           mu        sigma 
 94.34166416 -85.35582080 159.76952349  -0.00283408  -4.65833584  12.94283927 

Log-likelihood: -155.2 



Call (fit1):
all:
mle2(minuslogl = LL2, start = list(b0 = 100, b1 = -80, b2 = 160, 
    lmbd = 0.12, mu = 1, sigma = 14), method = "BFGS", data = exdata, 
    subset = Rating == "BBB", control = list(maxit = 5000))

Coefficients:
           b0            b1            b2          lmbd            mu         sigma 
 9.269496e+01 -8.600709e+01  1.586543e+02 -1.186371e-04 -6.305040e+00  1.458118e+01 

Log-likelihood: -315.6