Using MLE on a somewhat unusual likelihood function

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Using MLE on a somewhat unusual likelihood function

Hugo Andre
So I am trying to use the mle command (from stats4 package) to estimate a
number of parameters using data but it keeps throwing up this error message:

Error in solve.default(oout$hessian) :
  Lapack routine dgesv: system is exactly singular: U[1,1] = 0

This error sometimes indicates that the list of starting values is too far
from optimum but this is unlikely since I picked values close to where the
parameters usually end up. I have also tried switching these around a bit.

Here is the code:

  xhat = c(statemw-(1-alpha)*rval)
  survivalf <- function(x) {(1-plnorm(statemw,mean=mu,sd=logalpha))}

wagefn <- function(lam, eta, alpha, xhat, mu, logalpha)  {
  n=nrow(cpsdata2)
  wagevec = matrix(nrow=n,ncol=1)
    for (i in 1:n) {

    if (cpsdata2[i,2] > 0){
     wagevec[i,] <-
c(eta*lam*survivalf(statemw)*exp(-lam*survivalf(statemw)*cpsdata2[i,2,]))
    } else if (cpsdata2[i,1,]==statemw) {
      wagevec[i,] <-
c(lam*(survivalf(statemw)-survivalf((statemw-(1-alpha)*xhat)/alpha))/(eta+lam*survivalf(statemw)))
    } else if (cpsdata2[i,1,]>statemw) {
      wagevec[i,] <-
c(lam*plnorm((cpsdata2[i,1,]-(1-alpha)*xhat)/alpha,mean=mu,sd=logalpha)/(alpha*(eta+lam*survivalf(statemw))))
    }
      else {
       wagevec[i,] <- NA
      }
    }
  lnwagevec <- log(wagevec)
  -sum(lnwagevec>-200 & ln2wagevec<200, na.rm=TRUE)
}

fit <- mle(wagefn, start=listmat, method= "L-BFGS-B",lower=
c(-Inf,0),upper=c(Inf,Inf)


I know the likelihood function is a handful but it does return a reasonable
looking vector of values. The "lnwagevec>-200" etc is an inelegant way of
preventing values of Inf and -Inf from entering the sum, the actual values
rarely go as high as 8 or low as -5.

Thank you in advance to anyone responding!
/Hugo

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Using MLE on a somewhat unusual likelihood function

Peter Dalgaard-2
(inline)

> On 7 Nov 2017, at 06:02 , Hugo André <[hidden email]> wrote:
>
> So I am trying to use the mle command (from stats4 package) to estimate a
> number of parameters using data but it keeps throwing up this error message:
>
> Error in solve.default(oout$hessian) :
>  Lapack routine dgesv: system is exactly singular: U[1,1] = 0
>
> This error sometimes indicates that the list of starting values is too far
> from optimum but this is unlikely since I picked values close to where the
> parameters usually end up. I have also tried switching these around a bit.
>
> Here is the code:
>
>  xhat = c(statemw-(1-alpha)*rval)
>  survivalf <- function(x) {(1-plnorm(statemw,mean=mu,sd=logalpha))}
>
> wagefn <- function(lam, eta, alpha, xhat, mu, logalpha)  {
>  n=nrow(cpsdata2)
>  wagevec = matrix(nrow=n,ncol=1)
>    for (i in 1:n) {
>
>    if (cpsdata2[i,2] > 0){
>     wagevec[i,] <-
> c(eta*lam*survivalf(statemw)*exp(-lam*survivalf(statemw)*cpsdata2[i,2,]))
>    } else if (cpsdata2[i,1,]==statemw) {
>      wagevec[i,] <-
> c(lam*(survivalf(statemw)-survivalf((statemw-(1-alpha)*xhat)/alpha))/(eta+lam*survivalf(statemw)))
>    } else if (cpsdata2[i,1,]>statemw) {
>      wagevec[i,] <-
> c(lam*plnorm((cpsdata2[i,1,]-(1-alpha)*xhat)/alpha,mean=mu,sd=logalpha)/(alpha*(eta+lam*survivalf(statemw))))
>    }
>      else {
>       wagevec[i,] <- NA
>      }
>    }
>  lnwagevec <- log(wagevec)
>  -sum(lnwagevec>-200 & ln2wagevec<200, na.rm=TRUE)
> }
>
> fit <- mle(wagefn, start=listmat, method= "L-BFGS-B",lower=
> c(-Inf,0),upper=c(Inf,Inf)
>
>
> I know the likelihood function is a handful but it does return a reasonable
> looking vector of values. The "lnwagevec>-200" etc is an inelegant way of
> preventing values of Inf and -Inf from entering the sum, the actual values
> rarely go as high as 8 or low as -5.
>

If you're getting infinite likelihoods, something may be needing a rethink, but first this:

I don't think

 -sum(lnwagevec>-200 & ln2wagevec<200, na.rm=TRUE)

does what I think you think it does....

Presumably you wanted -sum(lnagevec[lnwagevec>-200 & ln2wagevec<200], na.rm=TRUE)

-pd


> Thank you in advance to anyone responding!
> /Hugo
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] 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.

--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: [hidden email]  Priv: [hidden email]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Using MLE on a somewhat unusual likelihood function

Peter Dalgaard-2

> On 7 Nov 2017, at 11:38 , peter dalgaard <[hidden email]> wrote:
>
> Presumably you wanted -sum(lnagevec[lnwagevec>-200 & ln2wagevec<200], na.rm=TRUE)

Well, after correcting the obvious spelling error. Sorry about that...

-pd

--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: [hidden email]  Priv: [hidden email]

______________________________________________
[hidden email] 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.