$$ \mathbb{E}_p[f] = \int f(x) p(x)~dx $$
$$ \mathbb{E}_p[f] = \int f(x) p(x)~dx \approx \frac{1}{k}\sum_{j=1}^k f(x_j) $$ where $$ x_j \text{ is distributed like } p $$
$$ u \sim \operatorname{Uniform}(0, 1) $$ $$ \text{solve } \int_{-\infty}^x p(t)~dt = u $$
def metropolis_hastings(n_samples, pdf):
position = 0.
samples = []
for _ in range(n_samples):
proposal = np.random.normal(loc=position)
accept_prob = pdf(proposal) / pdf(position)
if np.random.rand() < accept_prob:
position = proposal
samples.append(position)
return samples
import arviz as az
import pymc3 as pm
with pm.Model():
x = pm.Normal('x')
trace = pm.sample(step=pm.Metropolis())
az.plot_trace(trace)