On 8/3/07 22:12, "Gad Abraham" <g.abraham at ms.unimelb.edu.au>
wrote:
> Laura Hill wrote:
>>
>>
>> On 7/3/07 00:15, "Gad Abraham" <g.abraham at
ms.unimelb.edu.au> wrote:
>>
>>>> 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)
>>>>> }
>>
>>
>>
>> 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
Hi Gad,
Thanks for the tip about the empty vector, I ever knew you could do that. I
just have one problem,
Lets say Q is a 2x2 matrix
p is a 1x2 matrix
q is a 2x1 matrix
y is vector of times, say y = c(5, 10)
How do I multiply Q by each time y[i]?
I would like to get the answer to the equation
loglik[i] <- log(p %*% expm(Q * y[i]) %*% q)
Where first y=5 and then y=10 so that the answers to loglik for each i are
put into the empty vector.
I'm sure that I am missing something fairly obvious here but can't put
my
finger on it.
Thanks in advance
Laura
