# Programming the error function (erf) in Python

1. Nov 13, 2007

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

PHP:
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?

2. Nov 13, 2007

### Astronuc

Staff Emeritus
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

3. Nov 13, 2007

### 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().

4. Nov 13, 2007

### Staff: Mentor

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.

5. Nov 13, 2007

### Tiiba

0

Mine does. Always did.

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

Last edited by a moderator: May 3, 2017
6. Nov 14, 2007