Inverse stochastic problem solution

Question

I am trying to determine a probability density $p(\mu)$ such that, when $\mu$ is inputted into a forward simulation equation: $d \sim \mathcal{N}(\mu,0.1)$, I obtain a distribution on $d\sim \mathcal{N}(0, 1)$.

I can write down the joint density $p(d,\mu) = p(d|\mu) p(\mu)$, which explicitly is:

$p(d,\mu) = \frac{1}{\sqrt{2\pi\sigma^2}}\text{exp}(\frac{(d-\mu)^2}{2\times 0.1^2}) p(\mu).$

This, when integrated wrt $\mu$ should yield the marginal density of interest:

$\frac{1}{\sqrt{2\pi}}\text{exp}(\frac{d^2}{2}) = \int^{\infty}_{-\infty} \frac{1}{\sqrt{2\pi\times 0.1^2}}\text{exp}(\frac{(d-\mu)^2}{2\times 0.1^2}) p(\mu)\mathrm{d}\mu.$

I think this is a Fredholm equation and have tried solving it using,

eqn = PDF[NormalDistribution[0, 1], d] == Integrate[PDF[NormalDistribution[mu, 0.1], d] p[mu],
         {mu, -\[Infinity], \[Infinity]}]

sol = DSolveValue[eqn, p[x], x]

But this just returns unevaluated.

Does anyone have any idea of how to solve this inverse problem?

Answer 1

It's a bit disappointing that DSolve didn't get you the result you were looking for, but the answer is that p is distributed as NormalDistribution[0, Sqrt[99/100]:

Integrate[
  PDF[NormalDistribution[mu, 1/10], d] PDF[NormalDistribution[0, Sqrt[99/100]], mu],
  {mu, -\[Infinity], \[Infinity]}
] == PDF[NormalDistribution[], d]

True

The reason for this is that your integral is a convolution of a Gaussian with some unknown function which has to result in another Gaussian. Gaussian functions are closed under convolution, so p[mu] has to be another Gaussian. With that knowledge in hand, you can just do the convolution and figure out what the parameters of the distribution for p have to be.

Here is the general formula for the convolution of two Gaussians:

Integrate[
 PDF[NormalDistribution[mu, s1], d] PDF[NormalDistribution[mu1, s2], 
   mu],
 {mu, -\[Infinity], \[Infinity]},
 Assumptions -> s1 > 0 && s2 > 0
 ]

E^(-((d - mu1)^2/(2 (s1^2 + s2^2))))/(Sqrt[2 Pi] Sqrt[s1^2 + s2^2])

The reason your integral becomes a convolution is because mu in NormalDistribution[mu, 1/10] is a location parameter (i.e., it appears as mu -x in PDF[NormalDistribution[mu, 1/10], x]. If you're planning to do this kind of convolutions more frequently, I recommend you look up the characteristic function of the distribution you're working with (you can find it with CharacteristicFunction). The characteristic function is essentially the Fourier transform of the PDF of a distribution and the Fourier transforms turns convolutions into multiplication, so I think you'll find that useful.

You'll also find that the CharacteristicFunction of a normal distribution is a Gaussian function, which is a property unique to Gaussian functions. So when you convolve two Gaussians, you can multiply the characteristic functions (which are also Gaussians). Multiplying two Gaussians gives yet another Gaussian so now you have a characteristic function that again has to belong to a normal distribution. For other distributions you'll have to take the inverse Fourier transform to find out what distribution corresponds to the characteristic function you calculated.