See below
On Fri, 12 Apr 2013, Julio Sergio wrote:
Berend Hasselman <bhh <at> xs4all.nl> writes:
Your function miBeta returns a scalar when the argument mu is a vector.
Use Vectorize to vectorize it. Like this
VmiBeta <- Vectorize(miBeta,vectorize.args=c("mu"))
VmiBeta(c(420,440))
and draw the curve with this
curve(VmiBeta,xlim=c(370,430), xlab="mu", ylab="L(mu)")
Berend
Taking into account what you have pointed out, I reprogrammed my function
as follows, as an alternative solution to yours:
zetas <- function(alpha) {z <- qnorm(alpha/2); c(z,-z)}
# First transformation function
Tzx <- function(z, sigma_p, mu_p) sigma_p*z + mu_p
# Second transformation function
Txz <- function(x, sigma_p, mu_p) (x - mu_p)/sigma_p
BetaG <- function(mu, alpha, n, sigma, mu_0) {
lasZ <- zetas(alpha) # Zs corresponding to alpha
sigma_M <- sigma/sqrt(n) # sd of my distribution
lasX <- Tzx(lasZ, sigma_M, mu_0) # Transformed Zs into Xs
# Now I consider mu to be a vector composed of m's
NewZ <- lapply(mu, function(m) Txz(lasX, sigma_M, m))
# NewZ is a list, the same length as mu, with 2D vectors
# The result will be a vector, the same length as mu (and NewZ)
sapply(NewZ, function(zz) pnorm(zz[2]) - pnorm(zz[1]))
}
miBeta <- function(mu) BetaG(mu, 0.05, 36, 48, 400)
plot(miBeta,xlim=c(370,430), xlab="mu", ylab="L(mu)")
I hope this is useful to people following this discussion,
-Sergio.
Adding to what you have defined above, consider the benefit that true
vectorization provides:
BetaGv <- function(mu, alpha, n, sigma, mu_0) {
lasZ <- zetas( alpha ) # Zs corresponding to alpha
sigma_M <- sigma/sqrt( n ) # sd of my distribution
lasX <- Tzx( lasZ, sigma_M, mu_0 ) # Transformed Zs into Xs
# Now fold lasX and mu into matrices where columns are defined by lasX
# and rows are defined by mu
lasX_M <- matrix( lasX, nrow=length(mu), ncol=2, byrow=TRUE )
mu_M <- matrix( mu, nrow=length(mu), ncol=2 ) )
# Compute newZ
NewZ <- Txz( lasX_M, sigma_M, mu_M )
# NewZ is a matrix, where the first column corresponds to lower Z, and
# the second column corresponds to upper Z
# The result will be a vector, the same length as mu
pnorm(NewZ[,2]) - pnorm(NewZ[,1])
}
miBetav <- function(mu) BetaGv(mu, 0.05, 36, 48, 400)
system.time( miBeta( seq( 370, 430, length.out=1e5 ) ) )
system.time( miBetav( seq( 370, 430, length.out=1e5 ) ) )
On my machine, I get:
system.time( miBeta( seq( 370, 430, length.out=1e5 ) ) )
user system elapsed
1.300 0.024 1.476
system.time( miBetav( seq( 370, 430, length.out=1e5 ) ) )
user system elapsed
0.076 0.000 0.073
So a bit of matrix manipulation speeds things up considerably.
(That being said, I have a feeling that some algorithmic optimization
could speed things up even more, using theoretical considerations.)
---------------------------------------------------------------------------
Jeff Newmiller The ..... ..... Go Live...
DCN:<jdnew...@dcn.davis.ca.us> Basics: ##.#. ##.#. Live Go...
Live: OO#.. Dead: OO#.. Playing
Research Engineer (Solar/Batteries O.O#. #.O#. with
/Software/Embedded Controllers) .OO#. .OO#. rocks...1k
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
