How to implement an iterative Lucy algorithm in Mathematica?

Question

I am trying to implement the Lucy algorithm that the authors of this paper use to get a distribution of diffusion coefficients from the distribution of displacements. The equation looks like this.

$$P^{n+1}(D)=P^n(D)\int \frac{G_s(x)}{G_s^n(x)} g(x\vert D) dx$$

Where $G_s^n(x)=\int P^n(D) g(x \vert D) dx$.

I know $G_s(x)$, $g(x\vert D)$ and $P^1(D)$ so it seems pretty straightforward to iterate until reaching a solution inside the desired tolerance.

I have tried to define a function

loopfunction[x_, y_, z_, m_, d_] := x*NIntegrate[y*z/m, {t, -Infinity, Infinity}]

Taking as arguments $[P(D)$, $g(t|D)$,$G_s(t)$,$G_s^n(t)$,$d]$ and then using a While loop to try to find a value $P^*=P^{n+1}$ such that $\Vert P^* - P^n \Vert < \varepsilon$ for some tolerance.

I can't get the While loop to work mainly because of the integrations involved in the process (it seems to be some kind of problem when I try to consider as variables those variables over which I have not integrated). I think this is because NIntegrate is not capable of giving me a result if I don't specify a value for the free variables and I can't use Integrate because the expressions become quite difficult after a couple of iterations.

How can I implement this sort of iterative algorithm in Mathematica?

I can provide more code if it is necessary and clarification if the question is not clear.

Answer 1

I'm not sure that I get everything correctly, but I hope it would be helpful.

I did not found the form explicit form of $g(x|D)$ so I used some function (you can specify our own). This is the initialization step

dave = 1;(*Subscript[D, ave]*)
d0 = 1;(*Subscript[d, 0]*)

g[x_, d_] := d*Exp[-(x^2/2)];
p[1] := 1/dave Exp[-d/dave];
G[s_, x_][1] := Simplify[Integrate[p[1]*g[x, d], {d, 0, d0}]]

This is the recursive definition ($P(n)$ should be normalized as I understand from the paper):

 p[n_] := Module[{norm, pNotNorm},
   G[s, x][n - 1] = Simplify@Integrate[p[n - 1]*g[x, d], {d, 0, d0}];
   pNotNorm[n] = 
    p[n - 1]*
     Simplify[
      Integrate[
       G[s, x][1]/G[s, x][n - 1]*g[x, d], {x, -Infinity, Infinity}]];
   (p[n] = (1/Integrate[pNotNorm[n], {d, 0, d0}])*pNotNorm[n])];

We used the memorization technique that will save a lot of time

    p[10] // AbsoluteTiming
(*about 5 sec on my machine*)

Also, I prefer here using Integrate to NIntegrate, so that we will have analytical results. And then just to have some function that finds the number of iteration needed to have some tolerance:

ListLinePlot@Table[N[p[n] /. d -> 0.2], {n, 1, 11}]

findIterationNumber[eps_, D_, pstar_] := Module[{n},
  n = 1;
  While[Abs[(p[n] /. d -> D) - pstar] > eps,
   Print[Abs[(p[n] /. d -> D) - pstar]]; n++];
  n]

findIterationNumber[0.01, 0.2, 0]

Hope you will find it helpful.