I am trying to use ParameterMixtureDistribution with a custom distribution as the first argument. It does not evaluate. In attempting to solve my problem, I reverted to a much simpler scenario (below), which I am puzzled I cannot get to work. My assumption is that I am missing something terribly obvious.

If I do

  NormalDistribution[m, s], 
  m \[Distributed] NormalDistribution[m2, s2]

I get

NormalDistribution[m2, Sqrt[s^2 + s2^2]]

However, if I try and use a self-defined function as the first parameter, here just a Gaussian, it fails to evaluate.

First define a Gaussian:

f = 1/(s Sqrt[2 Pi]) Exp[-((x - m)^2/(2 s^2))]

Now try the same mixture:

  m \[Distributed] NormalDistribution[m2, s2]

and it does not evaluate. Which is puzzling as I thought the two would be equivalent.

I have been playing around with a few things, none of which seem to work.

  • Make f a ProbabilityDistribution
  • Make f a PDF
  • Add Assumptions to ParameterMixtureDistribution

I am starting to think I am just missing something very obvious. What am I missing?


1 Answer 1


dist = ParameterMixtureDistribution[NormalDistribution[m, s], 
  m \[Distributed] NormalDistribution[m2, s2]]

(* NormalDistribution[m2, Sqrt[s^2 + s2^2]] *)

For comparison purposes,

 stats = #[dist] & /@ {Mean, StandardDeviation, PDF[#, x] &, CDF[#, x] &}]

(* {0.000525, {m2, Sqrt[s^2 + s2^2], E^(-((-m2 + x)^2/(2 (s^2 + s2^2))))/(
  Sqrt[2 π] Sqrt[s^2 + s2^2]), 
  1/2 Erfc[(m2 - x)/(Sqrt[2] Sqrt[s^2 + s2^2])]}} *)

For a user-defined function,

f = 1/(s Sqrt[2 π]) Exp[-((x - m)^2/(2 s^2))];

For Mathematica to consider f as a distribution you must use ProbabilityDistribution

fdist = ProbabilityDistribution[f, {x, -Infinity, Infinity},
   Assumptions -> m ∈ Reals && s > 0];

dist2 = ParameterMixtureDistribution[fdist, 
  m \[Distributed] NormalDistribution[m2, s2]]

(* ParameterMixtureDistribution[
 ProbabilityDistribution[E^(-((\[FormalX] - m)^2/(2 s^2)))/(
  Sqrt[2 π] s), {\[FormalX], -∞, ∞}, 
  Assumptions -> m ∈ Reals && s > 0], m \[Distributed] NormalDistribution[m2, s2]] *)

While Mathematica does not recognize nor provide a concise notation, dist2 is a valid distribution and produces the same results (albeit extremely slowly).

 stats2 = #[dist2] & /@ {Mean, StandardDeviation, PDF[#, x] &, CDF[#, x] &}]

(* {58.3517, {m2, Sqrt[s^2 + s2^2], E^(-((m2 - x)^2/(2 (s^2 + s2^2))))/(
  Sqrt[2 π] Sqrt[s^2 + s2^2]), 
  1/2 Erfc[(m2 - x)/(Sqrt[2] Sqrt[s^2 + s2^2])]}} *)

Thread[stats == stats2] // Simplify

(* {True, True, True, True} *)
  • $\begingroup$ Thanks Bob, that's great. $\endgroup$
    – flyingmind
    Commented May 16, 2022 at 8:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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