Estimating probabilities from MCMC samples

[hdelecki]

Hello all,

I'm relatively new to probabilistic programming, so this question might be silly. As a simple example, say I want to estimate the probability that a sample from the 1D normal distribution is above a certain threshold using MCMC samples. A real application would involve higher-dimensional, non-Gaussian distributions.

Maybe a naïve way to generate some samples from the 1D distribution would be

@model function example(threshold)
	x ~ Normal()
	y ~ Dirac(x > threshold)
end

m = example(2.0) | (;y=true)
example_chain = sample(m, MH(), 100000)

Is there a way that I could use these samples with another Turing model to infer p(y=true)?

Thanks!