[R] MLE, precision

Thu Jul 29 22:17:40 CEST 2004

      1.  Don't use "t" as a variable name.  It is the name of the 
matrix transpose function.  In most but not all contexts, R is smart 
enough to tell whether you want the system function or the local object. 

      2.  I can't tell from your question what you want.  "PLEASE do 
read the posting guide! http://www.R-project.org/posting-guide.html".  
You may be able to get answers to many of your questions by following 
that process.  If you follow that guide and still have a question for 
this listserve, the exercise will likely help you formulate a question 
that will be easier for people to understand, increasing thereby the 
likelihood that you will get an appropriate answer. 

      3.  I wonder if the following will help: 

 > (td <- outer(1:3, 4:5))
     [,1] [,2]
[1,]    4    5
[2,]    8   10
[3,]   12   15
 > rowSums(td)
[1]  9 18 27
 > colSums(td)
[1] 24 30

      hope this helps.  spencer graves

boshao zhang wrote:

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