Skip to content

Existence of non-vectorised functions

6 messages · Uwe Ligges, Ravi Varadhan, Martin Maechler +2 more

Original
Ravi Varadhan
# 
Dear R-Group:

Recently, I ran into a problem.  I was using a function called "I.1", 
which evaluates the first-order modified Bessel function of the first 
kind, in the package "CircStats". This function is not vectorized, 
since it uses a couple of "if" conditions.  However, when I called this 
function with a vector argument, I got no error/warning messages in 
R1.5.1 (under Windows 2000). Furthermore, the results were generally 
incorect (please see the attached results).  However, if I vectorized 
it using "ifelse" everything is fine (I created a function 
called "myI.1" to do this).  I was wondering if there are functions in 
other packages, where such problems exist due to non-vectorization?

thanks,
Ravi.
 
#######################################################
[1]    0.5651591    1.5906369    3.9533702    9.7594643   24.3351955
 [6]   61.3232698  155.7162813  396.4825232 1005.2747089 2517.6347571
[1]    0.5651591    1.5906369    3.9533702    9.7594652   24.3356418
 [6]   61.3419369  156.0390965  399.8731348 1030.9147087 2670.9883206



-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Uwe Ligges
#
Ravi Varadhan wrote:
That's extremely probable. It depends on the package author and the
complexity of the problem, thus you have to check it yourself whether it
works vectorized or not, e.g. in the help files.

BTW: You might want to send your improvements of I.1() to the package
author as well (or even at first - hint, hint) ... 

Uwe Ligges
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Ravi Varadhan
# 
Thanks to Uwe Ligges and Ben Bolker for their responses.  My main concern
was with the fact that R doesn't issue a warning when the user attempts an
"if" test with a vector argument.  Splus (at least Splus 6) does do this,
saying that only the first element is used in testing the condition. I feel
that R should also give this warning.

thanks again,
Ravi.

----- Original Message -----
From: "Uwe Ligges" <ligges at statistik.uni-dortmund.de>
To: "Ravi Varadhan" <rvaradha at jhsph.edu>
Cc: <r-help at stat.math.ethz.ch>
Sent: Tuesday, November 26, 2002 12:43 PM
Subject: Re: [R] Existence of non-vectorised functions
-.-.-
http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
_._
-.-.-
http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
_._

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Martin Maechler
# 
Ravi> Dear R-Group: Recently, I ran into a problem.  I was
    Ravi> using a function called "I.1", which evaluates the
    Ravi> first-order modified Bessel function of the first
    Ravi> kind, in the package "CircStats". This function is not
    Ravi> vectorized, since it uses a couple of "if" conditions....

    Ravi> <...........>

If you use  library(help = "CircStats")
(something you should do for any package when you find problems!)
you see that it has been ported from S (probably S-plus).

In R, the Bessel functions have been available for long time,
in fast and vectorized form  --- although for "numeric" (i.e
non-complex) arguments only, and help(bessel) {e.g.} gives you

     besselI(x, nu, expon.scaled = FALSE)
     besselK(x, nu, expon.scaled = FALSE)
     besselJ(x, nu)
     besselY(x, nu)

i.e. I.0(x) = besselI(x,0) \
     I.1(x) = besselI(x,1)  > when x is numeric
     I.p(x) = besselI(x,p) /

Only if you are interested in the non-numeric complex case,
you (and even more the maintainer of CircStats) should
investigate improving I.1 (and I.p) to work vectorized and you
are right,  ifelse() might a good way to do so.

Martin Maechler <maechler at stat.math.ethz.ch>	http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum  LEO C16	Leonhardstr. 27
ETH (Federal Inst. Technology)	8092 Zurich	SWITZERLAND
phone: x-41-1-632-3408		fax: ...-1228			<><
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Ben Bolker
# 
I have hacked together some pre-existing code (using C code from TOMS
written by D. E. Amos) for computing complex Bessel functions.  I 
considered sending it in, but (1) it's part of a larger package I haven't 
cleaned up yet, (2) I wasn't sure about the licensing/redistribution 
implications of the ACM license (permission required for commercial use: 
http://www.acm.org/pubs/copyright_policy/softwareCRnotice.html).
  If anyone would find it useful, drop me a line.

  Ben Bolker
On Wed, 27 Nov 2002, Martin Maechler wrote:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Brian Ripley
    

    

      
      #
    

  


  

      
On Tue, 26 Nov 2002, Ravi Varadhan wrote:

        

            

      
      
      
    

        
It now does in R-devel (and also for while()).