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.

I would like to solve the equations as given below in mathematica. I have written some codes before but I have no idea when there is sub indexing to indicate the series of functions. I will be very happy and grateful if you could help me. Here are the equations:

$$P[n]=\int_{(-\infty,B)\cup (A,\infty)} q_n (x) \mathrm{d}x$$

and

$$q_n(x)=\int_{B}^A q_{n-1} (\omega)f(x-\omega)\mathrm{d} \omega,\quad q_1=f,\quad n\geq 1.$$ and $f$ is some density function, say standard Gaussian density. Eventually, I want to obtain the probability mass $P[n]$.

Thank you very much!


I wrote this (from Hector's post):

f0[x_] := 1/Sqrt[2*Pi*1^2]*E^(-(x + 1)^2/(2*1^2))
q[n_ /; n > 1, B_, A_, f_][x_] := NIntegrate[q[n - 1, B, A, f][ω] f[x - ω], {ω, B, A}]
q[1, B_, A_, f_][x_] :=f[x]
p[n_, B_, A_, f_] := NIntegrate[q[n, B, A, f][x], {x, -∞, B}] +              
                     NIntegrate[q[n, B, A, f][x], {x, A, ∞}]

N[p[1, -2, 2, f0[x]]] --> 0.160005
N[p[2, -2, 2, f0[x]]] --> 0.35594
N[p[3, -2, 2, f0[x]]] --> 0.237045
N[p[4, -2, 2, f0[x]]] --> 0.229432
N[p[5, -2, 2, f0[x]]] --> 0.222063

...

these results seem incorrect because $p$ is the probability mass function and should add up to $1$.

Correct results for $f_0$ Gaussian with mean $-1$ and variance $1$ and $A=2$ and $B=-2$:

n=1-->0.1600
n=2-->0.3559
n=3-->0.2323
n=4-->0.1238
n=5-->0.0633
n=6-->0.0320
n=7-->0.0162
share|improve this question
    
You need the second line from my answer. That line ends the recursion defined in the first line. –  Hector Sep 20 '13 at 16:13
    
Also, the function must be a pure function (the type that uses # as the variable). I will update my answer to make it work numerically and with functions defined the way you do it. –  Hector Sep 20 '13 at 16:21

1 Answer 1

up vote 5 down vote accepted

Since OP seems to be interested in numerical calculations, this might be more suitable:

f0[x_]:= 1/Sqrt[2*Pi*1^2]*E^(-(x + 1)^2/(2*1^2))

q[1, A_, B_, f_] := f

q[n_ /; n > 1, A_, B_, f_][x_] :=
  Module[{ω},tempIntegrate[Evaluate[q[n - 1, A, B, f]][ω] f[x - ω], {ω, A, B}]]

p[n_, A_, B_, f_] :=
    tempIntegrate[q[n, A, B, f][x], {x, -∞, B}] + 
     tempIntegrate[q[n, A, B, f][x], {x, A, ∞}] //. {
      s_ tempIntegrate[a_, b__] :> tempIntegrate[s a, b], 
      tempIntegrate[tempIntegrate[a_, b__], c__] :> tempIntegrate[a, b, c] } /.
   tempIntegrate -> NIntegrate

Then, as a test case, p[3, -2, 2, f0] evaluates to 0.73581.


Original answer

The subindices are really parameters. One way to do it is to include the parameters along with the variables (see how Mathematica does it for the Legrendre polynomials for example). Subindices could also be part of the head of the function, for example f[a][x] can be thought of as the function $f_a(x)$.

The code you need is as below. Note that A and B are also parameters (and that I changed their order respect to your post).

q[n_ /; n > 1, A_, B_, f_][x_] := Integrate[q[n - 1, A, B, f][ω] f[x - ω], {ω, A, B}]
q[1, A_, B_, f_][x_] := HoldForm@f[x]
p[n_, A_, B_, f_] := Integrate[q[n, A, B, f][x], {x, -∞, A}] +
     Integrate[q[n, A, B, f][x], {x, B, ∞}]

I included a HoldForm to check whether the recursion is correct. The following input

p[2, 1, 3, Exp[-#^2] &]

produces

enter image description here

You might want to erase "HoldForm@" from the code above to get the actual answer. That is, change the second line to

q[1, A_, B_, f_][x_] := f[x]

Now, depending on the function f, Mathematica might take a long time to answer if you leave those integrals in analytic form. You might want to replace Integrate by NItegrate if you are interested in a numerical approximation.

share|improve this answer
    
what does the sign '#' mean? –  Seyhmus Güngören Sep 17 '13 at 20:22
    
It is the parameter of a pure function. Note that the final answer depends on what is the density function f. So, you need to supply it. In the example, f[x]=Exp[-x^2]. –  Hector Sep 17 '13 at 21:36
1  
You need to localize the variable in the integral Module[{ω}, Integrate[...]] should do it –  ssch Sep 17 '13 at 23:55
    
I tried the codes you posted. Although the codes seem correct, I was not able to get $P[n]$. I want to plot it an have some feeling what is going on with the computations but I cannot get it. I only need some numerical values for $P[n]$. What should I do? –  Seyhmus Güngören Sep 19 '13 at 19:05
    
Please see the code part which i appended to my question. Still I dont get any value. –  Seyhmus Güngören Sep 20 '13 at 11:42

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.