Programming the error function (erf) in Python

[Tiiba]

I'm having very unusual problems posting in the math forum.

I am trying to write a program that would calculate phi, the CDF of the standard normal curve, in Python. I tried calculating the erf using a series to do it, as described by Wikipedia. But mine diverges, when it should approach 1.

def erfterm(x, n):
   num = (-1.0)**n * x**(2*n+1)
   denum = fact(n) * (2*n+1)
   return num/denum
   
def erf(x):
   #By Maclaurin series
   def term(n):
      num = (-1.0)**n * x**(2*n+1)
      denum = fact(n) * (2*n+1)
      return num/denum
   return sum([erfterm(x, i) for i in range(25)]) *(2 * sqrt(pi))

Something seems to be wrong with erfterm, because it's giving very large results. What gives? Is the equation at Wikipedia wrong? Or did I copy it wrong? Or is my computer wrong?

 

[Astronuc]

Does n start at 0 or 1? If 1 then the series would start 2(1)+1 and one misses the first term 2(0)+1. One has to take enough terms.

By the time x gets out to 2 or -2, the erf(x) should be nearly 1.

In the last line (2 * sqrt(pi)) should be (2 / sqrt(pi)), but that's only a factor and would not contribute to divergence.

One could debug by writing out the successive terms.

The Maclaurin series on Wikipedia is correct. Confirm with
http://mathworld.wolfram.com/Erf.html

 

[mgb_phys]

Also be careful, python 2.x does not autmatically promote floats.
So 2 * sqrt() will give an integer answer, you need to write 2.0 * sqrt().

 

[jim mcnamara]

If it's erf-ing for learning programming, go for it. :)

But, erf() and erfc() have been available for python since 2003 at least. AFAIK. Other than as a fun exercise, why are you re-inventing the wheel? If this were going to be production code and I was involved in the project, I would insist that you use a widely tested implementation of something rather than a roll-your-own. It's called best practice.

 

[Tiiba]

Does n start at 0 or 1?

0

Also be careful, python 2.x does not autmatically promote floats.

Mine does. Always did.

But, erf() and erfc() have been available for python since 2003 at least.

http://docs.python.org/lib/module-math.html"

 

[jim mcnamara]

This:
http://pylab.sourceforge.net/packages/included_functions.html