5
$\begingroup$

Working with some iterative integral equations, I have Gaussian density functions involved therein. Integrating the Gaussian function I obtain the error function. When I take the second integration, Mathematica is not able to provide an analytical expression. Because apparently there is no closed-form expression for the integral of the error function.

Question: Is it possible the tell Mathematica to consider the error function as some given exact expression which is analytically integrable?

For example

$$\mathrm{erf}(x)=\frac{1}{(1+a_1x+a_2x^2+a_3x^3+a_4x^4)^4}$$

for some appropriate constants $a_1$, $a_2$, $a_3$, $a_4$?

$\endgroup$
9
  • $\begingroup$ Yes, it's possible. Take a look at a closely related question Numerical underflow for a scaled error function. I guess therein you'll find answer to your question. $\endgroup$
    – Artes
    Oct 22, 2013 at 23:34
  • $\begingroup$ hi @Artes thank you very much for your comment. I am looking at it right now. $\endgroup$ Oct 22, 2013 at 23:36
  • $\begingroup$ I mean you can find only closed-form approximations, not an analytic closed-form formula. $\endgroup$
    – Artes
    Oct 22, 2013 at 23:38
  • $\begingroup$ yes but my question is that, whenever a result comes as $erf$ can I tell mathematica to use a predefined function instead of it automatically. Can I interrupt the general understanding of mathematica from the error function? $\endgroup$ Oct 22, 2013 at 23:42
  • 1
    $\begingroup$ And you can't adapt Integrate[Erf[t],{t,0,x}] -> x*Erf[x] + (Exp[-x^2] - 1)/Sqrt[Pi]?? $\endgroup$ Oct 24, 2013 at 6:03

2 Answers 2

6
$\begingroup$

You could unprotect Erf and redefine it as you like it.

Unprotect[Erf];
 (* Here you define the constants and your erf approximation *)
   a1=...; etc.
   Erf[x_] := 1/(1 + a1 x + a2 x^2 + a3 x^3 + a4 x^4)^4;
 (* Here you perform your calculations with the approximate value of erf *)
   myComplexCalculation = Erf[x]/2
 (* When you are done, clear your definition and reprotect Erf *)
  Clear[Erf];
 Protect[Erf];

This can be dangerous but, as long as you know what you are doing, it will work.

A more local workaround is to enclose your 'approximate' calculation within a Block that blocks the definition of Erf:

Block[ {Erf},
  (* Here you define the constants and your erf approximation *)
  a1=...; etc.
  Erf[x_] := 1/(1 + a1 x + a2 x^2 + a3 x^3 + a4 x^4)^4;

  (* Here you perform your calculations with the approximate value of erf *)
  myComplexCalculation = Erf[x]/2
]

Is this what you were looking for or was it the actual expression for the approximation?

As a more general approach to temporarily modifying built-in functions within the walls of a scoping construct, you might also want to have a look at the undocumented function InheritedBlock, described in this post on Stack Overflow. It appears that it will copy all definitions associated with the symbol you want to block, so that you can change, within the block only the part you need to. Basically, while

Block[ {builtinFun},
    code redefining builtinFun from scratch
    and using the new builtinFun
]

would give you a completely clean symbol for builtinFun to code from scratch, using

Internal`InheritedBlock[ {builtinFun},
    code redefining parts of the inherited builtinFun
    and using the modified builtinFun
]

would give you a completely modifiable copy of builtinFun, so that you can modify only the parts you need to modify.

$\endgroup$
1
  • $\begingroup$ It seems like this was what I was looking for. I will check and let you know. $\endgroup$ Oct 23, 2013 at 10:35
4
$\begingroup$

Altrnate to redefining Erf you might define your own approximate probability distribution function and work with that , eg..

cvals = {a1 -> .278393, a2 -> .230389, a3 -> .000972, a4 -> .078108}

(* erf approximation from wikipedia , take for what its worth *)
approxerf[x_] := 
       Sign[x] ( 
            1 - 1/(1 + a1 Abs[x] + a2 x^2 + a3 Abs[x]^3 + a4 x^4)^4) /. cvals

approxpdf[mean_, sd_] = (512/sd) Piecewise[{
     {(  
        4 g a2 + 4 a4 g^3  + Sqrt[2] 2 a1 + 
          3 Sqrt[2] a3 (g^2)  )/((4 + 2 Sqrt[2] a1 g + 2 a2 g^2 + 
          Sqrt[2] a3 g^3 + a4 g^4)^5)  /. g -> (mean - x)/sd,
        ( mean - x)/sd >= 0},
     {(  
        4 g a2 + 4 a4 g^3  + Sqrt[2] 2 a1 + 
         3 Sqrt[2] a3 (g^2)  )/((4 + 2 Sqrt[2] a1 g + 2 a2 g^2 + 
         Sqrt[2] a3 g^3 + a4 g^4)^5)  /. g -> -((mean - x)/sd), 
        True}}] /. cvals;

mynormal[mean_, sd_] = 
     ProbabilityDistribution[
        approxpdf[mean, sd], {x, -Infinity, Infinity}];

Plot[ { CDF[mynormal[2, 5], x] , 
        CDF[NormalDistribution[2, 5], x]}  , {x, -20, 20}]

This should be equivalent - more a matter of semantics. Do you want an approximate solution to an exact problem or an exact solution to an approximate problem.

Note this performs pretty slowly..but I assume the other approach has the same issue as that analytic expression isn't particularly nice to integrate.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.