If you need to loop over some data set, you essentially have two options. Control flow iteration, which uses your common for loops and the like, and recursion. At first blush, recursion is very mysterious. The canonical example is the factorial() function where you multiply an integer by all the integers less than it and stop at 1. You remember from high school: 5! = 5x4x3x2x1.
I've written the following R implementation of it below:
factorial <- function(x){
if(x==1)
return( 1)
else
return(x*factorial(x-1))
}If you need to stare at it for a while, you're not alone. How does a function call itself? Shouldn't this simply blow up with a laundry list of errors at minimum and at worse a complete crash of your computer? Well, actually it works. Try it. This is how it flows:
factorial(5) 5*factorial(4) 5*(4*factorial(3)) 5*(4*(3*factorial(2))) 5*(4*(3*(2*factorial(1)))) 5*(4*(3*(2*1))) 120
I thought to see how R prefers to implement the factorial() function and found my answer in the gregmisc package. This is the code for its factorial() function:
function (x) gamma(x + 1)
Okay, so what is the gamma() function? From the R help files, I got this definition:
Γ(x) = integral_0^Inf t^(x-1) exp(-t) dt
At some point I need to have the ability to blog functions in common mathematical notation and add some clarity to a dense concept, but for now that's all I got. To be honest, if I had to explain the gamma function for 5 minutes, I would fail. At minimum though, I think mine is simpler.
UPDATE: Best wishes to Matz and his native Japan. May damage be minimized and recovery swift.
Hello!
ReplyDeleteInteresting post. For fun, I've been trying to write my own recursive factorial function using a 'while' loop, but have failed so far. Your method is much simpler, but not as flexible as I would like.
For example:
> factorial(1)
[1] 1
> factorial(2)
[1] 2
> factorial(1:10)
[1] 1
Warning message:
In if (x == 1) return(1) else return(x * factorial(x - 1)) :
the condition has length > 1 and only the first element will be used
>
Any ideas how to fix this? Please keep in mind that I'm not criticizing your method--on the contrary, it's very interesting and I'd like to know how to improve it...
R's native, gamma-based factorial() function works fine for all normal vector values, as in:
> factorial(1:10)
[1] 1 2 6 24 120 720 5040 40320 362880 3628800
>
Hello again,
ReplyDeleteAlso, LaTeX script is a great way to render equations (see http://en.wikibooks.org/wiki/LaTeX/Mathematics for notation info). Most of the equations presented in Wikipedia articles are rendered in TeX. There are many ways to publish equations on the web, but for Mac I prefer LaTeXiT (bundled with LaTeX, but available separately), which allows you to export in several formats (include .GIF, .PDF, and MathML). MathML is great if you don't want to insert an image, and it works in all modern browsers (even explorer, with an add-on).
TeX script would read:
\Gamma(x)=\int_{0}^{\infty}t^{x-1}e^{-t}dx