[R] A strange behaviour in the graphical function "curve"

[Jeff Newmiller]

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:<jdnewmil at 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