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

[Gabor Grothendieck]

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