Solving equation over TruncatedDistribution (FindRoot)

Question

I'm trying to solve the following equation :

Clear["Global`*"]
Val1 = 0.4582;
Val2 = Sqrt[0.0052];
FindRoot[{
  Mean[TruncatedDistribution[{0.35, 0.7},NormalDistribution[mu, sig]]] == Val1,
  StandardDeviation[TruncatedDistribution[{0.35, 0.7}, NormalDistribution[mu, sig]]] == Val2},
{{mu, Val1}, {sig, Val2}}]

which gives :

(*{mu -> 0.412198, sig -> 0.102817}*)

So far so good but I'm wondering how long it takes mathematica to solve the system... Ideas to speed up the calculation? I don't need a very precise result, 5 decimal digits are enough.

ps : I have already tried several variations of WorkingPrecision, PrecisionGoal... it doesn't really speed up anything

Answer 1

You could pre-calculate symbolic expressions for the mean and standard deviation of the truncated distribution:

ClearAll[mean, stdev]
mean[{a_, b_}, {mu_, sigma_}] = 
  Simplify[
    Mean@ TruncatedDistribution[{a, b}, NormalDistribution[mu, sigma]], 
    0 < a < b
  ];

stdev[{a_, b_}, {mu_, sigma_}] = 
  Simplify[
    StandardDeviation@ TruncatedDistribution[{a, b}, NormalDistribution[mu, sigma]], 
    0 < a < b
  ];

Note that, critically, these are defined with Set (=) and not SetDelayed (:=) so they are calculated once and then stored away for future use. These do take quite some time (a minute or two on my laptop), but then you get the nice bonus that you can change the truncation limits at will afterwards.

Once you have these, the FindRoot calculations becomes practically instantaneous:

Val1 = 0.4582;
Val2 = Sqrt[0.0052];

FindRoot[
  {mean[{0.35, 0.7}, {mu, sigma}] == Val1, 
   stdev[{0.35, 0.7}, {mu, sigma}] == Val2},
  {{mu, Val1}, {sigma, Val2}}
]

(* Out: {mu -> 0.412198, sigma -> 0.102817} *)

And you can easily change the truncation limits and the values at will, obtaining results very fast.

Answer 2

Saving the truncated distribution's properties would allow for this calculation to proceed much more quickly.

However, asking Mathematica to directly calculate the standard deviation of the truncated distribution in the general case takes a distressingly long time itself. Doing some of this by hand results in a much faster solution, so:

Let us grab the PDF of the interior distribution and manually truncate it ourselves:

fpdf[mu_, sig_, x_] = PDF[NormalDistribution[mu, sig], x];

And then calculate three related integrals where the limits are based on the truncation:

ex0 = Integrate[fpdf[mu, sig, x], {x, 35/100, 70/100}];
ex1 = Integrate[x fpdf[mu, sig, x], {x, 35/100, 70/100}]/ex0;
ex2 = Integrate[x^2 fpdf[mu, sig, x], {x, 35/100, 70/100}]/ex0;

ex0 captures the amount of normal distribution within the truncated region, ex1 captures the expected value of a variable sampled from the normal distribution, and ex2 capture $E(X^2)$. These three integrals take about 6 seconds total to evaluate on my machine.

From an answer over on math.SE, we can see that with this information we can calculate the mean and variance directly:

distMean[mu_, sig_] = ex1;
distVariance[mu_, sig_] = ex2 - ex1^2;

And relatedly, the standard deviation:

distStdDev[mu_, sig_] = Sqrt[ex2 - ex1^2];

If we plug in these distMean and distStdDev functions into the original FindRoot, we find that the solution is quite fast:

FindRoot[{distMean[mu, sig] == 0.4582, 
  distStdDev[mu, sig] == Sqrt[0.0052]}, {{mu, 0.5}, {sig, 0.1}}] // AbsoluteTiming

{0.002852, {mu -> 0.412198, sig -> 0.102817}}

This approach can potentially be generalized to other types of distributions, but it is somewhat dependent on being able to get a nice closed form for the 3 integrals involved for maximum speedup.

Answer 3

I'm assuming you want to estimate the parameters using the "method of moments" rather than using maximum likelihood. (If not, then you should consider maximum likelihood as it computes much faster and you can also easily obtain estimates of precision.)

(* Generate some data *)
d = TruncatedDistribution[{0.35, 0.7}, NormalDistribution[mu, sig]];
d2 = d /. mu -> 0.4 /. sig -> 0.1
SeedRandom[12345];
x = RandomVariate[d2, 1000];

(* Method of moments *)
AbsoluteTiming[
 FindDistributionParameters[x, d, {{mu, 0.4}, {sig, 0.1}},
  ParameterEstimator -> "MethodOfMoments"]]
(* {11.81 seconds, {mu -> 0.379214, sig -> 0.112361}} *)

(* Maximum likelihood *)
AbsoluteTiming[FindDistributionParameters[x, d, {{mu, 0.4}, {sig, 0.1}}]]
(* {0.37837 seconds, {mu -> 0.379214, sig -> 0.112361}} *)

For this particular sample the two sets of estimates match. However, that is not always the case.

Answer 4

A straightforward approach is as follows. Its advantage is that you find Mean and StandardDeviation only one time.

a = Mean[TruncatedDistribution[{7/20, 7/10},NormalDistribution[mu, sig]]]//AbsoluteTiming

{4.8678, -((E^(-(49/(200 sig^2)) - mu^2/( 2 sig^2)) (-Sqrt[2] E^(147/(800 sig^2) + (7 mu)/(20 sig^2)) sig Abs[7 - 20 mu] Abs[7 - 10 mu] + Sqrt[2] E^((7 mu)/(10 sig^2)) sig Abs[7 - 20 mu] Abs[7 - 10 mu] + 7 E^(49/(200 sig^2) + mu^2/(2 sig^2)) mu Sqrt[\[Pi]] Abs[7 - 10 mu] Erf[Abs[7 - 20 mu]/(20 Sqrt[2] sig)] - 20 E^(49/(200 sig^2) + mu^2/(2 sig^2)) mu^2 Sqrt[\[Pi]] Abs[7 - 10 mu] Erf[Abs[7 - 20 mu]/(20 Sqrt[2] sig)] - 7 E^(49/(200 sig^2) + mu^2/(2 sig^2)) mu Sqrt[\[Pi]] Abs[7 - 20 mu] Erf[Abs[7 - 10 mu]/(10 Sqrt[2] sig)] + 10 E^(49/(200 sig^2) + mu^2/(2 sig^2)) mu^2 Sqrt[\[Pi]] Abs[7 - 20 mu] Erf[Abs[7 - 10 mu]/( 10 Sqrt[2] sig)]))/(Sqrt[\[Pi]] Abs[7 - 20 mu] Abs[ 7 - 10 mu] (-Erf[(-(7/10) + mu)/(Sqrt[2] sig)] + Erf[(-(7/20) + mu)/(Sqrt[2] sig)])))}

b = StandardDeviation[TruncatedDistribution[{7/20, 7/10}, NormalDistribution[mu, sig]]]//AbsoluteTiming

{96.9828, (1/( 2^(3/4) Sqrt[ 5 \[Pi]]))(\[Sqrt](-((E^(-(49/(100 sig^2)) - mu^2/sig^2) sig (Erf[(7 - 20 mu)/(20 Sqrt[2] sig)] - Erf[(7 - 10 mu)/(10 Sqrt[2] sig)]) (-20 Sqrt[2] E^( 147/(400 sig^2) + (7 mu)/(10 sig^2)) sig + 40 Sqrt[2] E^(147/(800 sig^2) + (21 mu)/(20 sig^2)) sig - 20 Sqrt[2] E^((7 mu)/(5 sig^2)) sig - 7 E^(343/(800 sig^2) + (7 mu)/(20 sig^2) + mu^2/(2 sig^2)) Sqrt[\[Pi]] Erf[(7 - 20 mu)/(20 Sqrt[2] sig)] + 14 E^(49/(200 sig^2) + (7 mu)/(10 sig^2) + mu^2/(2 sig^2)) Sqrt[\[Pi]] Erf[(7 - 20 mu)/(20 Sqrt[2] sig)] + 20 E^(343/(800 sig^2) + (7 mu)/(20 sig^2) + mu^2/(2 sig^2)) mu Sqrt[\[Pi]] Erf[(7 - 20 mu)/(20 Sqrt[2] sig)] - 20 E^(49/(200 sig^2) + (7 mu)/(10 sig^2) + mu^2/(2 sig^2)) mu Sqrt[\[Pi]] Erf[(7 - 20 mu)/(20 Sqrt[2] sig)] + 10 Sqrt[2] E^( 49/(100 sig^2) + mu^2/sig^2) \[Pi] sig Erf[(7 - 20 mu)/( 20 Sqrt[2] sig)]^2 + 7 E^(343/(800 sig^2) + (7 mu)/(20 sig^2) + mu^2/(2 sig^2)) Sqrt[\[Pi]] Erf[(7 - 10 mu)/(10 Sqrt[2] sig)] - 14 E^(49/(200 sig^2) + (7 mu)/(10 sig^2) + mu^2/(2 sig^2)) Sqrt[\[Pi]] Erf[(7 - 10 mu)/(10 Sqrt[2] sig)] - 20 E^(343/(800 sig^2) + (7 mu)/(20 sig^2) + mu^2/(2 sig^2)) mu Sqrt[\[Pi]] Erf[(7 - 10 mu)/(10 Sqrt[2] sig)] + 20 E^(49/(200 sig^2) + (7 mu)/(10 sig^2) + mu^2/(2 sig^2)) mu Sqrt[\[Pi]] Erf[(7 - 10 mu)/(10 Sqrt[2] sig)] - 20 Sqrt[2] E^( 49/(100 sig^2) + mu^2/sig^2) \[Pi] sig Erf[(7 - 20 mu)/( 20 Sqrt[2] sig)] Erf[(7 - 10 mu)/(10 Sqrt[2] sig)] + 10 Sqrt[2] E^( 49/(100 sig^2) + mu^2/sig^2) \[Pi] sig Erf[(7 - 10 mu)/( 10 Sqrt[2] sig)]^2))/((Erf[(-(7/10) + mu)/( Sqrt[2] sig)] - Erf[(-(7/20) + mu)/( Sqrt[2] sig)])^2 (-Erf[(-(7/10) + mu)/(Sqrt[2] sig)] + Erf[(-(7/20) + mu)/(Sqrt[2] sig)])))))}

Now

Val1 = 0.4582;Val2 = Sqrt[0.0052];
FindRoot[{a[[2]] == Val1,  b[[2]] == Val2}, {{mu, Val1}, {sig, Val2}}] // AbsoluteTiming

{0.0036685, {mu -> 0.412198, sig -> 0.102817}}

For comparison, I could not execute your code during several minutes.