[R] MLE, precision

[boshao zhang]

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




--- Spencer Graves <spencer.graves at pdf.com> wrote:

>       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 at 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
> >  
> >
> 
>