# [R] MLE, precision

Spencer Graves spencer.graves at pdf.com
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
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
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.
>
>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
>>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
>>>
>>>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
>>
>>
>>>
>>>
>>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