Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Working with some iterative integral equations I have Gaussian density functions involved therein. Integrating the gaussian function I obrain 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

$$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$?

share|improve this question
    
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. –  Artes Oct 22 '13 at 23:34
    
hi @Artes thank you very much for your comment. I am looking at it right now. –  Seyhmus Güngören Oct 22 '13 at 23:36
    
I mean you can find only closed-form approximations, not an analytic closed-form formula. –  Artes Oct 22 '13 at 23:38
    
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? –  Seyhmus Güngören Oct 22 '13 at 23:42
    
@Artes by the way I checked the solutions over there and they are some functional approximations defined as a new function, a bit different. –  Seyhmus Güngören Oct 22 '13 at 23:43
show 4 more comments

2 Answers

up vote 4 down vote accepted

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.

share|improve this answer
    
It seems like this was what I was looking for. I will check and let you know. –  Seyhmus Güngören Oct 23 '13 at 10:35
add comment

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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