[R] Estimating parameters of 2 phase Coxian using optim

[Laura Hill]

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