From: boshao zhang <zboshao_at_yahoo.com>

Date: Fri 30 Jul 2004 - 05:06:27 EST

Problem in f + g: needed atomic data, got an object of class "function"

Use traceback() to see the call stack

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

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

- Spencer Graves <spencer.graves@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@stat.math.ethz.ch mailing list
*

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

>

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
*