how can I solve integrals including recursion of some sequences of functions

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

Answer

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

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.