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];
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];];]
  • 2
    $\begingroup$ Compare yours with ...Timing[ pAllDP=Table[ pOneDP=Table[ probIntegrand[... pTotal = NIntegrate[pOffset[z]*pDifferential[z], {z,0,2*Pi}], {i,Length[\[Phi]list]}], {n,Length[yA]}]; ] Replacing For and AppendTo with Table is supposed to be faster. $\endgroup$ – Bill Jun 21 '19 at 18:11

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]; 
    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}]
|improve this answer|||||
  • 2
    $\begingroup$ Thanks. This was actually the same time as doing the For loop but the notation is cleaner. I've gotten this down to 17s using: Method -> {Automatic, "SymbolicProcessing" -> 0} Any other suggestions for how I could manipulate this integral would be much appreciated. I've gotten the suggestion of maybe trying to pre-calculate parts, but not sure how to do this... $\endgroup$ – Megan Nantel Jun 21 '19 at 21:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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