Speeding up nested NIntegrate with convolution

Question

I'd like to calculate:

$ p(y_a, y_b | \phi_d) = \int_{-\infty}^{\infty} d\phi_d' \Big( \int_{-\pi/2}^{\pi/2} d\theta \exp(\frac{-(y_a-sin\theta)^2}{2v_a})\Big[\exp(\frac{-(y_b-sin(\theta+\phi_d')^2}{2v_b})+\exp(\frac{-(y_b-sin(\theta-\phi_d')^2}{2v_b}) \Big] \Big) \exp(\frac{-(\phi_d-\phi_d')^2}{2v_d})$

So there's an inner integral and then a convolution.

I'm currently doing this (changed the range of the convolution and results are fine). But I need to do this for many pairs of {ya,yb} to get a distribution as a function of $\phi_d$. Any suggestions for how to speed it up?

Thank you SO Much!!

ϕlist = Range[0.5, 0.35*Pi, 0.01];
Length[ϕlist]
pAllDP = {};

Timing[For[n = 1, n <= Length[yA], n++, pOneDP = {}; 
   For[i = 1, i <= Length[ϕlist], i++, 
    probIntegrand[θ_] := 
     Exp[(-(dpA - Sin[θ])^2/(2*varianceA))*
     (Exp[-(dpB - Sin[θ + ϕd])^2/(2*varianceB)] + 
      Exp[-(dpB - Sin[θ - ϕd])^2/(2*varianceB)])] /. dpA -> yA[[n]] /. 
       dpB -> yB[[n]] /. ϕd -> ϕlist[[i]];

pOffset[(z_)?NumericQ] := pOffest[z] = 
  NIntegrate[probIntegrand[θ], {θ, -Pi/2, Pi/2}, 
   MaxRecursion -> 20, AccuracyGoal -> 20, PrecisionGoal -> 7];

pDifferential[ϕ_] := Exp[-(ϕd - ϕ)^2/(2*varianceD)] /. ϕd -> ϕlist[[i]];
pTotal = NIntegrate[pOffset[z]*pDifferential[z], {z, 0, 2*Pi}];

AppendTo[pOneDP, pTotal];]; AppendTo[pAllDP, pOneDP];];]

Answer 1

I've benefited using ParallelDo or ParallelTable instead of For many times. You should check them out since they're easy to use.

Also, using compile can give you immense speed up. There are many examples and tips on using compile and parallelization in Mathematica across the internet.

I also recommend reading this post on 10 tips on writing faster Mathematica codes, a life-saver for me! Hope this helps.

Note: Adding a ParallelTable to @Bill 's comment should it get even faster:

\[Phi]list = Range[0.5, 0.35*Pi, 0.01]; 
Timing[
    pAllDP = ParallelTable[
        pOneDP = Table[
            probIntegrand[\[Theta]_] := 
                Exp[(-(dpA - Sin[\[Theta]])^2/(2*varianceA))* (Exp[-(dpB - Sin[\[Theta] + \[Phi]d])^2/(2*varianceB)] + Exp[-(dpB - Sin[\[Theta] - \[Phi]d])^2/(2*varianceB)])] /. dpA -> yA[[n]] /. dpB -> yB[[n]] /. \[Phi]d -> \[Phi]list[[i]]; 
           pOffset[(z_)?NumericQ] := pOffest[z] = NIntegrate[probIntegrand[\[Theta]], {\[Theta], -Pi/2, Pi/2}, MaxRecursion -> 20, AccuracyGoal -> 20, PrecisionGoal -> 7]; 
           pDifferential[\[Phi]_] := Exp[-(\[Phi]d - \[Phi])^2/(2*varianceD)] /. \[Phi]d -> \[Phi]list[[i]]; 
           pTotal = NIntegrate[pOffset[z]*pDifferential[z], {z, 0, 2*Pi}]
           ,{i,Length[\[Phi]list]}}]
        ,{n,Length[yA]]}}];
   ]