# [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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

```