[R] vectorizing a function

[Bill.Venables@cmis.csiro.au]

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