Animate a FindMinimum with NIntegrate

Question

The objective is to investigate how a fitted normal distribution approximates the sum of two normal distributions with a difference in $\mu$, $\sigma$, and amplitude. The purpose is finding out how the presence of undifferentiated peaks (that we know a priori are there) affect the center of the fitted distribution in order to compensate for it.

So we have normal distributions $N_i$, which we sample unevenly from with probabilities $p_i$ such that $\sum_i p_i = 1$. So the pdf of the total distribution is $\sum_i p_i * \text{pdf}(N_i)$.

f[μ_, σ_] := 
 1/(σ Sqrt[2 π]) Exp[-(1/(2 σ^2)) (x - μ)^2]
weightedsum[specs_] := Total[Map[#1[[1]] f @@ #1[[2 ;; 3]] &, specs]]/
 Total[Transpose[specs][[1]]]

We can then test this and find the for the case of two distributions $N_1 = N(0, 1)$, $N_2 = N(0.3, 1)$, $p_1 = 1/3$, $p_2 = 2/3$, we get the following pdf:

spec = {{1, 0, 1}, {2, 0.3, 1}};
weightedsum[spec]

1/3 (E^(-(1/2) (-0.3 + x)^2) Sqrt[2/π] + E^(-(x^2/2))/Sqrt[
   2 π])

Which is what we would expect.

Next, we need to find the fit that minimizes the L2 distance between the functions. We must defer evaluation of the internal of the FindMinimum, because Mathematica will attempt to evaluate the NIntegrate first, and we don't want that. The method I'm using to deter evaluation is to place a constraint on the function inputs.

fn[μ_?NumericQ, σ_?NumericQ] := 
 NIntegrate[(f[μ, σ] - weightedsum[spec])^2, {x, -6, 6}]
FindMinimum[fn[μ, σ], {{μ, 0}, {σ, 1}}]


{2.39408*10^-8, {μ -> 0.200499, σ -> 1.01001}}

(Note, the bounds are -6 and 6 because I was getting precision warnings with -Infinity and Infinity, eventually I will handle them better, but this was just a quick and dirty solution.)

We can plot this to make sure everything is functioning correctly.

Plot[{
  weightedsum[spec],
  f[μ, σ] /. {FindMinimum[
      fn[μ, σ], {{μ, 0}, {σ, 1}}][[2]]}
  }, {x, -3, 3}]

As mentioned previously, the point is to get a feel for this data, so that I can determine the importance of 2nd order effects and try to get a feel on how to model this, is second order sufficient? Is first order sufficient? Etc.

So I'm trying to set up Animate to rapidly scrub through many values in the interesting range. However, now I've got two places where I have to make sure values are being replaced before being passed into the NIntegrate function, and I don't think I'm doing a good enough job. I've tried several things, but here's the attempt which feels the closest to working.

fn[μ_?NumericQ, σ_?NumericQ, 
  spec_?(MatrixQ[#, NumericQ]) &] := 
 NIntegrate[(f[μ, σ] - weightedsum[spec])^2, {x, -6, 6}]
Animate[With[{
   spec = {
     {1, 0, 1},
     {pow, mu, 1}
     }
   },
  Plot[{
    weightedsum[spec],
    f[μ, σ] /. { 
      ReleaseHold[
        Hold[FindMinimum[
           fn[μ, σ, foo],
           {{μ, 0}, {σ, 1}}]] /. {foo -> spec}][[2]]
      }
    }, 
   {x, -7 , 7}]],
 {pow, 1, 5}, {mu, 0, 5}]

Which gives me the error:

ReplaceAll::reps: {μ,0} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.
ReplaceAll::reps: {σ,1} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

(Please note that mu and sigma are different variables from μ and σ, I'm using mu and sigma to replace into spec and μ and σ are my objective variables for the FindMinimum function.)

The error message stack trace indicates that this error is occurring within the FindMinimum function, and I don't quite understand why.

Here's what I think should be happening but this is my first time really trying to make Hold and ReleaseHold work so while I know this is wrong, this could also be funny:

FindMinimum will return {error, {μ -> best_mu, σ -> best_sigma}}

However, I want to pass spec to FindMinimum, and spec is specified by Animate so FindMinimum needs to be delayed until after that happens.

So I wrap it in Hold, and I perform the replacement outside of Hold. spec is replaced while the evaluation is being held. Now ReleaseHold occurs. FindMinimum is evaluated entirely with numerical values, and returns the result, which is passed through Hold and ReleaseHold, then I take the second index to get the replacement list.

Obviously, something in this sequence I described is not what happens, but I have no idea what it could be.

Answer 1

As a higher-level comment, most programming languages use "functions" as building blocks. You explicitly specify inputs and outputs.

Mathematica allows you to use "expressions" which reduces boilerplate as inputs are implicit. This helps for simple code, but for more complicated stuff you have to battle the evaluation order which increases the boilerplate.

As someone who used to know evaluation tricks and is still in top 10 Mathematica answerers on SO, I'll tell you from experience -- you can figure out ReleaseHold tricks to make the expression approach work in your case, but you'll forget them unless you are a full-time Mathematica dev.

More maintainable way to to deal with complexity is to structure your code in terms of functions.

Looking at your code, you want to animate things by changing the spec, hence wrap FindMinimum into a helper function getFit that takes specs as an argument

ClearAll["Global`*"];
f[\[Mu]_, \[Sigma]_] = 
  1/(\[Sigma]  Sqrt[
       2  \[Pi]])  Exp[-(1/(2  \[Sigma]^2))  (x - \[Mu])^2];
weightedSum[specs_] := 
  Total[Map[#1[[1]]  f @@ #1[[2 ;; 3]] &, specs]]/
   Total[Transpose[specs][[1]]];
fn[\[Mu]_?NumericQ, \[Sigma]_?NumericQ, specs_] := 
 NIntegrate[(f[\[Mu], \[Sigma]] - weightedSum[specs])^2, {x, -6, 6}]
fittedSum[m_, s_, specs_] := 
  f[\[Mu], \[Sigma]] /. 
   Last@FindMinimum[
     fn[\[Mu], \[Sigma], specs], {{\[Mu], m}, {\[Sigma], s}}];
specs = {{1, 0, 1}, {2, 0.3, 1}};
fn1 = weightedSum[specs];
fn2 = fittedSum[0, 1, specs];
Plot[{fn1, fn2}, {x, -7, 7}]

Since you only vary two argument of the spec, refactor that part further

compoundPlot[pow_, mu_] := (
   specs = {{1, 0, 1}, {pow, mu, 1}};
   fn1 = weightedSum[specs];
   fn2 = fittedSum[0, 1, specs];
   Plot[{fn1, fn2}, {x, -7, 7}]
   );
compoundPlot[1, 2]

Now you can animate compoundPlot by varying its two arguments