Vectorising and loop (was Re: [R] optim "alog-likelihoodfunction")

[Zhen Pang]

Thank you for your kind help. my s does not depend on the input theta. Code 
without loop is really efficient.


>From: Gabor Grothendieck <ggrothendieck at myway.com>
>To: r-help at stat.math.ethz.ch
>Subject: Re: Vectorising and loop (was Re: [R] optim 
>"alog-likelihoodfunction")
>Date: Thu, 30 Sep 2004 17:59:21 +0000 (UTC)
>
>Zhen Pang <nusbj <at> hotmail.com> writes:
>
>:
>: Mr. Grothendieck does suggest me to paste the data here. I just show a 
>small
>: one here. I must metion that I made a mistake in my former email. The 
>first
>: column should be j and is greater than the second column. I have 
>corrected
>: ll.
>:
>: z is the matrix below
>:
>: 2 1 3
>: 3 1 1
>: 3 3 1
>: 5 2 4
>:
>: k<-max(z[,1])
>: ll <- function(theta)
>:   {t<-0
>:    for (ii in 1:k)
>:       {t<-t+exp(theta[ii])}
>:    lll<-0
>:    x00<-1/(1+t)
>:    x0<-x00*exp(theta)
>: for (m in 1:length(z[,1]))
>:    {j<-z[m,1]
>:     i<-z[m,2]
>:     a<-z[m,3]
>:     l<-i:(k-j+i)
>:     s<-rep(0,k)
>:     s[l]<-choose(j,i)*choose((k-j),(l-i))/choose(k,l)
>: # we only define some of s to be non-zero, since dim(l) might be smaller
>: than dim(s)
>:     ss<-sum(s*x0)  # ss is a weighted sum of x0
>:      lll<-lll+a*log(ss)
>:     }
>: -lll
>: # the negative sign is to find the maximum of the log-likelihood 
>function.
>: It can be omitted if we #use the finscale option in optim.
>: }
>:
>: Then I need to optim(b0,ll,hessian=T),
>: where b0<-c(0.8331934, 20.8009068, -7.0893623,  1.2109221, 18.7213273).
>:
>: optim(b0,ll,hessian=T)
>: $par
>: [1]  0.8331934 20.8009068 -7.0893623  1.2109221 18.7213273
>:
>: $value
>: [1] 5.182791
>:
>: $counts
>: function gradient
>:       52       NA
>:
>: $convergence
>: [1] 0
>:
>: $message
>: NULL
>:
>: $hessian
>:               [,1]          [,2]          [,3]          [,4]          
>[,5]
>: [1,]  1.065814e-08 -9.325873e-09  0.000000e+00 -3.330669e-10 
>-2.109424e-09
>: [2,] -9.325873e-09  8.887936e-01 -3.330669e-10 -1.620926e-08 
>-8.887936e-01
>: [3,]  0.000000e+00 -3.330669e-10 -6.661338e-10  0.000000e+00  
>0.000000e+00
>: [4,] -3.330669e-10 -1.620926e-08  0.000000e+00  7.549517e-09  
>7.105427e-09
>: [5,] -2.109424e-09 -8.887936e-01  0.000000e+00  7.105427e-09  
>8.887936e-01
>:
>:
>: I have tried to use eval() and modify my function, it seems to be able to
>: remove the m loop, however, optim() can not recognize it. So my main 
>concern
>: is to avoid the loop and optim() can works for my function. Thanks.
>
>
>
>I suspect your code may be wrong but taking it at face value
>s does not depend on the input theta so precalculate it
>as a matrix whose mth column is s[,m].  Also the only
>purpose of the loop indexed by m is to calculate lll and
>the final iteration of that loop calculates an lll which
>does not depend on the prior iterations so remove the loop
>and just run the final iteration.  Similarly we only need
>the final value of x0 that is calculated.  Note that the
>value of ll(b0), your loopy function, and ll2(b0) the one
>line non-loopy function using the precalculated s, give
>the same result:
>
>R> z <- matrix(c( 2,1,3, 3,1,1, 3,3,1, 5,2,4), 4, 3, byrow = TRUE)
>R> b0<-c(0.8331934, 20.8009068, -7.0893623, 1.2109221, 18.7213273)
>R>
>R> k<-max(z[,1])
>R> s <- apply(z, 1, function(z) {
>+ j<-z[1]; i<-z[2]
>+ l<-i:(k-j+i)
>+ s<-rep(0,k)
>+ s[l]<-choose(j,i)*choose((k-j),(l-i))/choose(k,l)
>+ s
>+ })
>
>R> ll2 <- function(theta) {
>+ - sum(z[,3]*log(crossprod(s, exp(theta)* 1/(1+sum(exp(theta))))))
>+ }
>R> ll(b0) # this is the ll function from your post
>[1] 5.182791
>R> ll2(b0)
>[1] 5.182791
>
>______________________________________________
>R-help at stat.math.ethz.ch mailing list
>https://stat.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide! 
>http://www.R-project.org/posting-guide.html