[R] Estimating parameters of 2 phase Coxian using optim

[Gad Abraham]

Laura Hill wrote:
> 
> 
> On 7/3/07 00:15, "Gad Abraham" <g.abraham at ms.unimelb.edu.au> wrote:
> 
>> Andy Fugard wrote:
>>> Hi There,
>>>
>>> Perhaps the problem is the line
>>>
>>>          loglik<-log(p %*% expm(Q * y(i)) %*% q)
>>>
>>> You mention that y is a vector but here you're treating it as a
>>> function.  Maybe try
>>>
>>>          loglik<-log(p %*% expm(Q * y[i]) %*% q)
>>>
>>> ?
>>>
>>> Don't have a clue about the correctness of the contents of cox2.lik...
>>>
>>> Andy
>>>
>>>
>>> On 6 Mar 2007, at 08:54, Laura Hill wrote:
>>>
>>>> Hi,
>>>>
>>>> My name is Laura. I'm a PhD student at Queen's University Belfast
>>>> and have
>>>> just started learning R. I was wondering if somebody could help me
>>>> to see
>>>> where I am going wrong in my code for estimating the parameters
>>>> [mu1, mu2,
>>>> lambda1] of a 2-phase Coxian Distribution.
>>>>
>>>> cox2.lik<-function(theta, y){
>>>>     mu1<-theta[1]
>>>>
>>>>     mu2<-theta[2]
>>>>
>>>>     lambda1<-theta[3]
>>>>
>>>>     p<-Matrix(c(1, 0), nrow=1, ncol=2)
>>>>
>>>>     Q<-Matrix(c(-(lambda1 + mu1), 0, lambda1, -mu2), nrow=2, ncol=2)
>>>>
>>>>     q<-Matrix(c(mu1, mu2), nrow=2, ncol=1)
>>>>
>>>>     for (i in 1:length(y)){
>>>>         loglik<-log(p %*% expm(Q * y(i)) %*% q)
>>>>     return(loglik)}
>>>>
>>>>     sumloglik<-sum(loglik)
>>>>
>>>>     return(-sumloglik)
>>>>     }
>> Just to add my 2 AU cents regarding the for loop:
>>
>> You're trying to create a vector of log likelihoods to sum up later, but
>> that's not what's happening there. Instead, assign an empty vector of
>> same length as y, then assign the loglik from each iteration to a
>> different cell.
>>
>> Lastly, there's no need to return anything from a for loop, it's not a
>> function.
>>
>> HTH,
>> Gad
> 
> 
> 
> Hi Gad,
> 
> Yes that's exactly hat I am trying to do. If I gave you a simple example,
> could you perhaps tell me how I could create a vector of log likelihoods.
> 
> 
> 
> Lets say I have 1x1 matrices:
> 
> p=[1]
> Q=[0.05]      i.e. [mu1]
> q=[-0.05]     i.e. [-mu1]
> 
> Where mu1 is the parameter that I would like to estimate and I have chosen
> the initial value mu1=0.05
> 
> 
> Loglik<-p %*% expm(Q*y) %*% q
> 
> Where y=(5 10)
> 
> I want to sum the log likelihoods that I get for y=5 and y=10 using
> 
> Sumloglik<-sum(allloglik)
> 
> Where allloglik = vector of log likelihoods
> 
> 
> Any help would be greatly appreciated.
> 
> Thanks in advance
> Laura
> 

Hi Laura,

Make an empty vector of required length, then assign the loglik to each 
of its cells, and don't return() anything:

loglik <- rep(0, length(y))
for(i in 1:length(y)){
    loglik[i] <- log(p %*% expm(Q * y[i]) %*% q)
}

Then you can sum(loglik) like you did before.

Cheers,
Gad

-- 
Gad Abraham
Department of Mathematics and Statistics
The University of Melbourne
Parkville 3010, Victoria, Australia
email: g.abraham at ms.unimelb.edu.au
web: http://www.ms.unimelb.edu.au/~gabraham