Функция R для правдоподобия

Question

Я пытаюсь проанализировать надежность отказоустойчивых систем с использованием моделей роста. Я уже установил модель Crow-Amsaa, но мне интересно, есть ли какой-либо пакет или какой-либо код для установки обобщенного процесса обновления (Kijima Model I) или типа II в R и найти его параметры Бета, Лямбда (или альфа) и кв. (или какая-то другая модель для средней кумулятивной функции MCF)

Уравнения номер 15 этой статьи дает выражение для логарифмического правдоподобия

http://arxiv.org/ftp/arxiv/papers/1006/1006.3718.pdf

Я пытался создать функцию, как это :

likelihood.G1=function(theta,x){ 
# x is a vector with the failure times, theta vector of parameters 
a=theta[1] #Alpha 
b=theta[2] #Beta 
q=theta[3] #q 

logl2=log(b/a) # First part of the equation 

for (i in 1:length(x)){ 
logl2=logl2 +(b-1)*log(x[i]/(a*(1+q)^(i-1))) -(x[i]/(a*(1+q)^(i-1)))^b 
} 
return(-logl2) #Negavite of the log-likelihood 
} 

И затем использовать некоторые rutine для минимизации -log (L)

theta=c(0.5,1.2,0.8) #Start parameters (lambda,beta,q) 

nlm(likelihood.G1,theta, x=Data) 

Или же

optim(theta,likelihood.G1,method="BFGS",x=Data) 

Однако это, кажется, какая-то ошибка, так как параметры, она возвращает не имеет смысла

Любые идеи, что я делаю не так?

Благодаря

Answer 1

Глядя на уравнение (16) бумаги вы ссылаетесь и сравнивая ее с кодом, похоже, вам не хватает одного термина в цикл. Кажется, что каждая точка данных способствует трем терминам логарифмического правдоподобия, но в вашем коде (внутри цикла) у вас есть только два термина (не считая срока обновления)

В частности, ваш код не включает 4-й термин в уравнении (16):

и ни это делает 7-й срок, и так далее. Это, по крайней мере, одна ошибка в коде. Дополнительным соображением будет то, что α и β должны быть больше нуля. Я не уверен, что решатель, который вы используете, рассматривает это ограничение.