✏️ Esercizio

import numpy as np
import pymc as pm
import scipy.stats as stats

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

Supponiamo di avere osservato il seguente campione:

# Observed data
y = np.array([
    16, 16.9, 14.5, 15.3, 16.8, 15.5, 9.9, 14.1, 16.1, 14.1, 15, 14.4, 17, 15.3,
    15, 16.9, 16.1, 13.6, 14.2, 14.3, 16.3, 15.3, 15.1, 14.9
])

Ipotizziamo che provenga da una popolazione che segue la legge Normale. Vogliamo fare inferenza sulla media \(\mu\) della popolazione. Supponiamo che la deviazione standard \(\sigma\) della popolazione sia nota e sia identica alla deviazione standard del campione (calcolata con \(n −1\) al denominatore). Imponiamo su \(\mu\) una distribuzione a priori Normale di media 12.7 e deviazione standard 0.9.

Si trovino la media e la deviazione standard della distribuzione a posteriori \(p(\mu \mid y)\) (a) usando la soluzione analitica, (b) usando PyMC.

Soluzione

mu_prior = 12.7
sigma_prior = 0.9
sigma = np.std(y, ddof=1)  # Standard deviation of the sample
n = len(y)
y_bar = np.mean(y)

# Calculating the posterior mean and standard deviation
tau_0 = sigma_prior
mu_posterior = ((1 / tau_0**2) * mu_prior + (n / sigma**2) * y_bar) / (1 / tau_0**2 + n / sigma**2)
sigma_posterior = np.sqrt(1 / ((1 / tau_0**2) + (n / sigma**2)))

print("Posterior mean:", mu_posterior)
print("Posterior standard deviation:", sigma_posterior)

Posterior mean: 14.86106928189399
Posterior standard deviation: 0.2883797105191614

with pm.Model() as model:
    # Priors
    mu = pm.Normal("mu", mu=mu_prior, sigma=sigma_prior)
    # Likelihood
    Y_obs = pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=y)

with model:
    idata = pm.sample()

Auto-assigning NUTS sampler...

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [mu]

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

100.00% [8000/8000 00:00<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 13 seconds.

# Summary of the MCMC results
pm.summary(idata)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
mu	14.859	0.295	14.288	15.405	0.007	0.005	1631.0	2656.0	1.0