Re: [R] MLE, precision

From: boshao zhang <zboshao_at_yahoo.com>
Date: Fri 30 Jul 2004 - 05:06:27 EST


Dear Spencer:

My problem get solved by using Matlab. It runs pretty quick(less than 5 seconds)and the result is stable with respect to the initial values. I was amaized. Here my t and are as long as 2390, sum the functions over t and d, the function becomes daunting. But I still like to try nlmb(I give up ms function in S or nlm in R).

But how can I sum the functions over t. To simplify the problem, I just need to know the following: t <- 1:2;
I would like to get f = sum(t*x + 1) over t. I tried the following:
 f<-0
 for (i in 1:2){

               g <- function(x){~t[i]*x+1}; f = f +g;
               }

Problem in f + g: needed atomic data, got an object of class "function"
Use traceback() to see the call stack

As you see, it refuse to work.

Please give me advice.
thank you.

Boshao

> Have you considered estimating ln.m1, ln.m2,
> and ln.b, which makes
> the negative log likelihood something like the
> following:
>
> l.ln<- function(ln.m1,ln.m2,ln.b){
> m1 <- exp(ln.m1); m2 <- exp(ln.m2); b <-
> exp(ln.b)
> lglk <- d*( ln.m1 + ln.m2
> + log1p(-exp(-(b+m2)*t)
> + (m1/b-d)*log(m2+b*exp(-b+m2)*t))
> + m1*t - m1/b*log(b+m2) )
>
> (-sum(lglk))
> }
> # NOT TESTED
>
> I don't know if I have this correct, but you
> should get the idea. Parameterizing in terms of
> logarithms automatically eliminates the constraints
> that m1, m2, and b must be positive.
>
> I also prefer to play with the function until I'm
> reasonably confident it will never produce NAs, and
> I use a few tricks to preserve numerical precision
> where I can. For example, log(b+m2) = log(b) +
> log1p(m2/b) = log(m2) + log1p(b/m2). If you use the
> first form when b is larger and the second when m1
> is larger, you should get more accurate answers.
> Using, e.g.:
>
> log.apb <- function(log.a, log.b){
> min.ab <- pmin(log.a, log.b)
> max.ab <- pmax(log.a, log.b)
> max.ab + log1p(exp(min.ab-max.ab))
> }
> # NOT TESTED
>
> If log.a and log.b are both really large, a and b
> could be Inf when log.a and log.b are finite.
> Computing log(a+b) like this eliminates that
> problem. The same problem occurs when log.a and
> log.b are so far negative that a and b are both
> numerically 0, even though log.a and log.b are very
> finite. This function eliminates that problem.
>
> Also, have you tried plotting your "l" vs. m1
> with m2 and b constant, and then vs. m2 with m2 and
> b constant and vs. b with m1 and m2 constant? Or
> (better) make contour plots of "l" vs. any 2 of
> these parameters with the other held constant. When
> I've done this in crudely similar situations, I've
> typically found that the log(likelihood) was more
> nearly parabolic in terms of ln.m1, ln.m2, and ln.b
> than in terms of the untransformed variables. This
> means that the traditional Wald confidence region
> procedures are more accurate, as they assume that
> the log(likelihood) is parabolic in the parameters
> estimated.
>
> hope this helps. spencer graves
> p.s. I suggest you avoid using "t" as a variable
> name: That's the name of the function for the
> transpose of a matrix. R and usually though not
> always tell from the context what you want.
> However, it's best to avoid that ambiguity. I often
> test at a command prompt variable names I want to
> use. If the response is "object not found", then I
> feel like I can use it.
>
> boshao zhang wrote:
>
> >Hi, everyone
> >
> >I am trying to estimate 3 parameters for my
> survival
> >function. It's very complicated. The negative
> >loglikelihood function is:
> >
> >l<- function(m1,m2,b) -sum( d*( log(m1) +
> log(m2)
> >+ log(1- exp(-(b + m2)*t)) ) + (m1/b - d)*log(m2 +
> >b*exp(-(b + m2)*t) ) + m1*t - m1/b*log(b+m2) )
> >
> >here d and t are given, "sum" means sum over these
> >two vairables.
> >the parameters are assumed small, m1, m2 in
> >thousandth, m2 in millionth.
> >
> >I used the function "nlm" to estimate m1,m2,b. But
> the
> >result is very bad. you can get more than 50
> warnings,
> >most of them are about "negative infinity"in log.
> And
> >the results are initial value dependent, or you
> will
> >get nothing when you choose some values.
> >
> >So I tried brutal force, i.e. evaluate the values
> of
> >grid point. It works well. Also, you can get the
> >correct answer of log(1e-12).
> >
> >My questions are:
> > What is the precision of a variable in R?
> > How to specify the constraint interval of
> parameters
> >in nlm? I tried lower, upper, it doesn't work.
> >any advice on MLE is appreciated.
> >
> >Thank you.
> >
> >Boshao
> >
> >______________________________________________
> >R-help@stat.math.ethz.ch mailing list

>
>https://www.stat.math.ethz.ch/mailman/listinfo/r-help

> >PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
> >
> >
>

>

R-help@stat.math.ethz.ch mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Fri Jul 30 05:13:43 2004

This archive was generated by hypermail 2.1.8 : Fri 18 Mar 2005 - 02:41:20 EST