Functions of Distributions vs. Distributions of Function Values

[RationalityEnhancement]

I am trying to use Squiggle to compute E[f(P(X))] where P(x) is a probability distribution that depends on the value x of the random variable X. What I want Squiggle to do is compute f(P(x)) for each sample x and then average the resulting scalars. Instead, Squiggle seems to compute f(mx(P(x1), P(x2), ...)). This leads to an entirely different result from what I expected.

To illustrate the problem, I have created the following example. The code is supposed to compute the distribution of the width of the 95% credible interval for the posterior distribution that you get by combining a standard normal prior with an extremely precise likelihood function. The correct answer is that it should always be extremely small. However, the answer I get from Squiggle is 2.3. I think this happens because Squiggle computes the quantile range of the mixture distribution instead of computing the quantile range for each distribution separately and then returning those value as a distribution.

I would very much appreciate if someone could either help me understand how I can get Squiggle to do what I want it to do or point out that this is not possible and that I should use a different tool instead.

posterior(prior,observation,precision)={
likelihood = normal(observation,pow(precision,-1/2))
normalize(prior.*likelihood)
}

posterior_quantile_range(prior,precision,observation)={
  quantile(posterior(prior,observation,precision),0.95)-quantile(posterior(prior,observation,precision),0.05)
}

prior = normal(0,1)
observation=normal(0,1)
precision = 100000
posterior_quantile_range(prior,precision,observation)

[Hazelfire]

Hey Falk!

I might be confused, but I don't know why your observation is a distribution, shouldn't it be a scalar? Turning this into a scalar fixes this issue.

If I'm wrong, and it makes sense for your observation to be a distribution. Then I'm not sure what you mean by that. Do you mean you want to run it as if you had x many observations that were sampled from the observation dist? If so, then how many? That would change the answer dramatically. If you want I can show you the code for this.

[RationalityEnhancement]

I want to evaluate f(P(X)) separately for many samples (x_1, ..., x_10000 ~ N(0,1)). In a procedural programming language, the code would be something like

for i=1:nr_samples
x ~ N(0,1) // a scalar value is sampled from a standard normal distribution
values(i) <- posterior_quantile_range(prior, precision, x)
end

mean(values)