Dear Arne,
Sorry for posting my mail twice and thanks a lot for your help.
Best regards,
Maram
Sent from my iPhone
> On Jul 7, 2015, at 9:55 PM, Arne Henningsen <arne.henning...@gmail.com> wrote:
>
> Dear Maram
>
> Please do NOT post your message twice!
>
> The warning messages occur each time, when maxLik() tries to calculate
> the logLik value for theta[1] <= 0, theta[1] + theta[2] <= 0, theta[3]
> <= 0 or something similar. According to the log-likelihood function,
> it seems that the parameters theta[1], theta[2], and theta[3] must be
> strictly positive. I suggest to re-parameterise your model so that the
> estimated parameters can take any values between minus infinity and
> infinity, e.g. by theta[1] <- exp( param[1] ); theta[2] <- exp(
> param[2] ); theta[3] <- exp( param[3] ) so that your estimated
> parameter vector 'param' consists of log( theta[1] ), log( theta[2] ),
> and log( theta[3] ). After the estimation, you can obtain the
> estimated values of the thetas by exp( param[1] ), exp( param[2] ),
> and exp( param[3] ) .
>
> Best regards,
> Arne
>
>
>
> 2015-07-06 2:29 GMT+02:00 Maram Salem <marammagdysa...@gmx.com>:
>> Dear All
>> I'm trying to find the maximum likelihood estimator of a certain
>> distribution based on the newton raphson method using maxLik package. I
>> wrote the log-likelihood , gradient, and hessian functionsusing the
>> following code.
>>
>> #Step 1: Creating the theta vector
>> theta <-vector(mode = "numeric", length = 3)
>> # Step 2: Setting the values of r and n
>> r<- 17
>> n <-30
>> # Step 3: Creating the T vector
>> T<-c(7.048,0.743,2.404,1.374,2.233,1.52,23.531,5.182,4.502,1.362,1.15,1.86,1.692,11.659,1.631,2.212,5.451)
>> # Step 4: Creating the C vector
>> C<-
>> c(0.562,5.69,12.603,3.999,6.156,4.004,5.248,4.878,7.122,17.069,23.996,1.538,7.792)
>> # The loglik. func.
>> loglik <- function(param) {
>> theta[1]<- param[1]
>> theta[2]<- param[2]
>> theta[3]<- param[3]
>> l<-(r*log(theta[3]))+(r*log(theta[1]+theta[2]))+(n*theta[3]*log(theta[1]))+(n*theta[3]*log(theta[2]))+
>> (-1*(theta[3]+1))*sum(log((T*(theta[1]+theta[2]))+(theta[1]*theta[2])))+
>> (-1*theta[3]*sum(log((C*(theta[1]+theta[2]))+(theta[1]*theta[2]))))
>> return(l)
>> }
>> # Step 5: Creating the gradient vector and calculating its inputs
>> U <- vector(mode="numeric",length=3)
>> gradlik<-function(param = theta,n, T,C)
>> {
>> U <- vector(mode="numeric",length=3)
>> theta[1] <- param[1]
>> theta[2] <- param[2]
>> theta[3] <- param[3]
>> r<- 17
>> n <-30
>> T<-c(7.048,0.743,2.404,1.374,2.233,1.52,23.531,5.182,4.502,1.362,1.15,1.86,1.692,11.659,1.631,2.212,5.451)
>> C<-
>> c(0.562,5.69,12.603,3.999,6.156,4.004,5.248,4.878,7.122,17.069,23.996,1.538,7.792)
>> U[1]<- (r/(theta[1]+theta[2]))+((n*theta[3])/theta[1])+(
>> -1*(theta[3]+1))*sum((T+theta[2])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+
>>
>> (-1*(theta[3]))*sum((C+theta[2])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
>> U[2]<-(r/(theta[1]+theta[2]))+((n*theta[3])/theta[2])+
>> (-1*(theta[3]+1))*sum((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+
>>
>> (-1*(theta[3]))*sum((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
>> U[3]<-(r/theta[3])+(n*log(theta[1]*theta[2]))+
>> (-1)*sum(log((T*(theta[1]+theta[2]))+(theta[1]*theta[2])))+(-1)*sum(log((C*(theta[1]+theta[2]))+(theta[1]*theta[2])))
>> return(U)
>> }
>> # Step 6: Creating the G (Hessian) matrix and Calculating its inputs
>> hesslik<-function(param=theta,n,T,C)
>> {
>> theta[1] <- param[1]
>> theta[2] <- param[2]
>> theta[3] <- param[3]
>> r<- 17
>> n <-30
>> T<-c(7.048,0.743,2.404,1.374,2.233,1.52,23.531,5.182,4.502,1.362,1.15,1.86,1.692,11.659,1.631,2.212,5.451)
>> C<-
>> c(0.562,5.69,12.603,3.999,6.156,4.004,5.248,4.878,7.122,17.069,23.996,1.538,7.792)
>> G<- matrix(nrow=3,ncol=3)
>> G[1,1]<-((-1*r)/((theta[1]+theta[2])^2))+((-1*n*theta[3])/(theta[1])^2)+
>> (theta[3]+1)*sum(((T+theta[2])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+(
>> theta[3])*sum(((C+theta[2])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
>> G[1,2]<-((-1*r)/((theta[1]+theta[2])^2))+
>> (theta[3]+1)*sum(((T)/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+
>> (theta[3])*sum(((C)/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
>> G[2,1]<-G[1,2]
>> G[1,3]<-(n/theta[1])+(-1)*sum(
>> (T+theta[2])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+(-1)*sum((C+theta[2])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
>> G[3,1]<-G[1,3]
>> G[2,2]<-((-1*r)/((theta[1]+theta[2])^2))+((-1*n*theta[3])/(theta[2])^2)+
>> (theta[3]+1)*sum(((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+(
>> theta[3])*sum(((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
>> G[2,3]<-(n/theta[2])+(-1)*sum((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+(-1)*sum((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
>> G[3,2]<-G[2,3]
>> G[3,3]<-((-1*r)/(theta[3])^2)
>> return(G)
>> }
>> mle<-maxLik(loglik, grad = gradlik, hess = hesslik, start=c(40,50,2))
>> There were 50 or more warnings (use warnings() to see the first 50)
>>
>> warnings ()
>> Warning messages:
>> 1: In log(theta[3]) : NaNs produced
>> 2: In log(theta[1] + theta[2]) : NaNs produced
>> 3: In log(theta[1]) : NaNs produced
>> 4: In log((T * (theta[1] + theta[2])) + (theta[1] * theta[2])) : NaNs
>> produced
>> and so on .......
>>
>> Although when I evaluate, for example, log(theta[3]) it gives me a number.
>> and the same applies for the other warnings.
>>
>> Then when I used summary (mle), I got
>>
>>
>> Maximum Likelihood estimation
>> Newton-Raphson maximisation, 7 iterations
>> Return code 1: gradient close to zero
>> Log-Likelihood: -55.89012
>> 3 free parameters
>> Estimates:
>> Estimate Std. error t value Pr(> t)
>> [1,] 11.132 Inf 0 1
>> [2,] 47.618 Inf 0 1
>> [3,] 1.293 Inf 0 1
>> --------------------------------------------
>>
>>
>> Where the estimates are far away from the starting values and they have
>> infinite standard errors. I think there is a problem with my gradlik or
>> hesslik functions, but I can't figure it out.
>> Any help?
>> Thank you in advance.
>>
>> Maram
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
>
>
> --
> Arne Henningsen
> http://www.arne-henningsen.name
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
