Re: Vectorising and loop (was Re: [R] optim "alog-likelihoodfunction")
Thank you for your kind help. my s does not depend on the input theta. Code
without loop is really efficient.
From: Gabor Grothendieck <[EMAIL PROTECTED]>
To: [EMAIL PROTECTED]
Subject: Re: Vectorising and loop (was Re: [R] optim
"alog-likelihoodfunction")
Date: Thu, 30 Sep 2004 17:59:21 + (UTC)
Zhen Pang 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.00e+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.00e+00 -3.330669e-10 -6.661338e-10 0.00e+00
0.00e+00
: [4,] -3.330669e-10 -1.620926e-08 0.00e+00 7.549517e-09
7.105427e-09
: [5,] -2.109424e-09 -8.887936e-01 0.00e+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
__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: Vectorising and loop (was Re: [R] optim "alog-likelihoodfunction")
Zhen Pang 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.00e+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.00e+00 -3.330669e-10 -6.661338e-10 0.00e+00 0.00e+00
: [4,] -3.330669e-10 -1.620926e-08 0.00e+00 7.549517e-09 7.105427e-09
: [5,] -2.109424e-09 -8.887936e-01 0.00e+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
__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: Vectorising and loop (was Re: [R] optim "alog-likelihoodfunction")
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.00e+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.00e+00 -3.330669e-10 -6.661338e-10 0.00e+00 0.00e+00
[4,] -3.330669e-10 -1.620926e-08 0.00e+00 7.549517e-09 7.105427e-09
[5,] -2.109424e-09 -8.887936e-01 0.00e+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.
__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
