# ParameterMixtureDistribution with custom distribution

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

ParameterMixtureDistribution[
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:

ParameterMixtureDistribution[
f,
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?

Clear["Global*"]

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

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


For comparison purposes,

AbsoluteTiming[
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).

AbsoluteTiming[
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])]}} *)

`