From: Spencer Graves <spencer.graves_at_pdf.com>

Date: Fri 30 Jul 2004 - 06:17:40 EST

> rowSums(td)

[1] 9 18 27

> colSums(td)

[1] 24 30

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 06:22:31 2004

Date: Fri 30 Jul 2004 - 06:17:40 EST

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

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 06:22:31 2004

*
This archive was generated by hypermail 2.1.8
: Wed 03 Nov 2004 - 22:55:21 EST
*