[R] vectorizing a function
Bill.Venables@cmis.csiro.au
Bill.Venables at cmis.csiro.au
Wed Oct 23 05:55:05 CEST 2002
Robin Hankin asks:
> -----Original Message-----
> From: Robin Hankin [mailto:r.hankin at auckland.ac.nz]
> Sent: Wednesday, October 23, 2002 11:32 AM
> To: r-help at stat.math.ethz.ch
> Subject: [R] vectorizing a function
>
> Dear R-xperts
>
> I have just written a little hypergeometric function, included below
> [the hypergeometric function crops up when solving a common type of
> ODE].
>
> It works fine on single values of the primary argument z, but
> vectorizing it is getting confusing.
>
> The best I have come up with so far just tests for z being longer than
> 1 and if so, uses sapply() recursively. This is fine, except that it
> doesn't preserve the dimensions correctly if z is a matrix or an
> array. And the function doesn't work properly if z is a scalar but A
> is a vector.
>
> besselI() does The Right Thing (tm), but it is internal ; what is the
> best way to vectorize this type of function?
>
>
>
>
>
> hypergeo <- function(A,B,C,z,tol=1e-6){
> if(length(z) > 1) {
> return(sapply(z,hypergeo,A=A,B=B,C=C,tol=tol))
> } else {
> term <- tmp <- 1
>
> for(n in 1:100){
> term <- term*A*B/C
> term=term*z/n
>
> partial.sum <- tmp + term
> if ((abs(partial.sum-tmp)<tol) ||
> is.infinite(partial.sum)){return(partial.sum)}
> A <- A+1
> B <- B+1
> C <- C+1
> tmp <- partial.sum
> }
> return (NaN)
> }
> }
>
>
[WNV] Hmm. Not sure if this is the best way to calculate this
function, particularly if z is near 1, but here is one way to vectorize it.
I've added a few purist refinements, of course. Also, I'm really not sure
that the convergence check is all that secure. Use with extreme caution,
particularly if |z| is near 1.
(NOTE: this is also vectorized wrt A, B and C as a kind of bonus.)
hyper <- function(A, B, C, z, tol = sqrt(.Machine$double.eps), maxit
= 10000) {
term <- z
term[ ] <- 1
partial.sum <- term
n <- 0
while(max(abs(term)) > tol && (n <- n+1) < maxit) {
term <- term*A*B/C * z/n
partial.sum <- partial.sum + term
A <- A+1
B <- B+1
C <- C+1
}
if(n == maxit) warning("convergence may not have been achieved")
partial.sum
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list