[S] Vectorization problem: Summary

Kent E. Holsinger (kent@darwin.eeb.uconn.edu)
07 Aug 1998 12:07:29 -0400


Last week I posed the following problem:

***begin***

> fst
function(s, mu, m)
{
p <- p.bar(s, mu, m)
const <- integrate(raw, 0, 1, sel = s, mut = mu, mig = m)$integral
p2 <- integrate(raw.x2, 0, 1, sel = s, mut = mu, mig = m)$integral/const
vp <- p2 - p^2
vp/(p * (1 - p))
}
> p.bar
function(s, mu, m)
{
const <- integrate(raw, 0, 1, sel = s, mut = mu, mig = m)$integral
p <- integrate(raw.x, 0, 1, sel = s, mut = mu, mig = m)$integral
p/const
}
> raw
function(x, sel, mut, mig)
{
s <- rep(sel, length(x))
mu <- rep(mut, length(x))
m <- rep(mig, length(x))
(exp( - s * (1 - x)) + exp( - s * x)) * (x * (1 - x))^(mu + (m/2) - 1)
}
> raw.x
function(x, sel, mut, mig)
{
s <- rep(sel, length(x))
mu <- rep(mut, length(x))
m <- rep(mig, length(x))
x * (exp( - s * (1 - x)) + exp( - s * x)) * (x * (1 - x))^(mu + (m/2) - 1)
}
> raw.x2
function(x, sel, mut, mig)
{
s <- rep(sel, length(x))
mu <- rep(mut, length(x))
m <- rep(mig, length(x))
(x^2) * (exp( - s * (1 - x)) + exp( - s * x)) * (x * (1 - x))^(mu + (m/2) - 1)
}

I also have vectors of data for s, mu, and m. What I'd like to do is
find a better way of doing this:

for (i in 1:length(s)) {
result[i] <- fst(s=s[i],mu=mu[i],m=m[i])
}

*** end ***

Bill Venables pointed out that I was evaluating the integrals more
frequently than necessary. His suggestion resulted in the following
improved set of functions.

> raw
function(x, sel, mut, mig)
{
s <- rep(sel, length(x))
mu <- rep(mut, length(x))
m <- rep(mig, length(x))
(exp( - s * (1 - x)) + exp( - s * x)) * (x * (1 - x))^(mu + (m/2) - 1)
}
> raw.x
function(x, sel, mut, mig)
{
s <- rep(sel, length(x))
mu <- rep(mut, length(x))
m <- rep(mig, length(x))
x * (exp( - s * (1 - x)) + exp( - s * x)) * (x * (1 - x))^(mu + (m/2) - 1)
}
> raw.x2
function(x, sel, mut, mig)
{
s <- rep(sel, length(x))
mu <- rep(mut, length(x))
m <- rep(mig, length(x))
(x^2) * (exp( - s * (1 - x)) + exp( - s * x)) * (x * (1 - x))^(mu + (m/2) - 1)
}
> fst
function(s, mu, m)
{
const <- integrate(raw, 0, 1, sel = s, mut = mu, mig = m)$integral
p <- integrate(raw.x, 0, 1, sel = s, mut = mu, mig = m)$integral
p2 <- integrate(raw.x2, 0, 1, sel = s, mut = mu, mig = m)$integral/const
vp <- p2 - p^2
vp/(p * (1 - p))
}

He also suggested the loop was not the problem, and he was right.

Pat Burns pointed out that I could use lapply() as follows:

> lapply(1:length(s) function(i, s, mu, m) fst(s[i], mu[i], m[i]),
s=s, mu=mu, m=m)

dos.time() for a for() loop with 100 elements on a P-166 and S-Plus
v4.5 was 44.18018s while the lapply() version finished about 0.9s
sooner (43.23999s). So the problem *is* in the integrations.

Thank you to Bill and Pat for their help.

Kent

P.S. If you haven't looked at Pat's Spoetry page, I recommend that you
do so (http://www.seanet.com/~pburns/Spoetry)

-- 
Kent E. Holsinger                Kent@Darwin.EEB.UConn.Edu
                                 http://darwin.eeb.uconn.edu
-- Department of Ecology & Evolutionary Biology          
-- University of Connecticut, U-43                                       
-- Storrs, CT   06269-3043                                               
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu.  To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message:  unsubscribe s-news