Integrating Gaussians
- To: mathgroup at yoda.ncsa.uiuc.edu
- Subject: Integrating Gaussians
- From: uunet!cello.hpl.hp.com!jacobson
- Date: Mon, 04 Jun 90 14:46:00 PDT
David Marchette reports problems with the integral
Integrate[Exp[-x^2/(2 sigma^2)],{x,-Infinity,Infinity}]
When we try it it comes out with
-Infinity
-(Sqrt[2] Sqrt[Pi] sigma Erf[-------------])
Sqrt[2] sigma
Out[1]= -------------------------------------------- +
2
Infinity
Sqrt[2] Sqrt[Pi] sigma Erf[-------------]
Sqrt[2] sigma
> -----------------------------------------
2
The problem is that Mathmatica does not know that Infinity/(Sqrt[2]
sigma) is still Infinity. However it does know that 3 Infinity or 1.5
Infinity is Infinity. So there must be some rules that let it know
this. Perhaps it is something like
Times[x_,Infinity] := Infinity /; x > 0
or maybe
Times[x_,Infinity] := Infinity /; Sign[x] == 1
(Of course these are equivalent, but Mathematica is really only a
pattern matching language and doesn't know that.)
Complaint:
The problem is that the Mathematica developers have chosen to not tell
us what the rules are, and consequently we don't know what facts to
give Mathematica about things like sigma so it can make the right
inferences.
It turns out that neither of the above is it, and I couldn't guess
any others, so I set about teaching Mathmatica about this stuff.
First I taught it the second form
DirectedInfinity/: Times[x_,Infinity] := Infinity /; Sign[x] == 1
and I taught it that square roots were positive :-)
Sign[Sqrt[x_]] := 1
and that sigma was positive
sigma/: Sign[sigma] := 1
Now that dealt with the sigma ok, even though it was in the denominator
(Sign[] is pretty smart!), but it couldn't deal with the Sqrt. The
reason is that the full form of 1/Sqrt[x] is Power[x, Rational[-1, 2]]
and I had only taught it about Sign[Power[x, Rational[1, 2]]].
Ok, so we have to teach it about Sqrt's in the denominator separately.
Sign[1/Sqrt[x_]] := 1 /; Sign[x] == 1
So now we try it out:
In[50]:= Integrate[Exp[-x^2/(2 sigma^2)],{x,-Infinity,Infinity}]
-Infinity
Sqrt[2] Sqrt[Pi] sigma Erf[-------------]
Sqrt[2] Sqrt[Pi] sigma Sqrt[2] sigma
Out[50]= ---------------------- - -----------------------------------------
2 2
Arrrrgggg!!!! It doesn't know about -Infinity.
So we change the teaching about Infinity:
DirectedInfinity/: Times[x_, DirectedInfinity[a_]] := DirectedInfinity[a] /;
Sign[x] == 1
And finally (drum roll):
In[50]:= Integrate[Exp[-x^2/(2 sigma^2)],{x,-Infinity,Infinity}]
Out[59]= Sqrt[2] Sqrt[Pi] sigma
Now really, shouldn't it be easier than all this?
-- David Jacobson