import os, sys
conda_lib = os.path.join(os.environ.get('CONDA_PREFIX',''), 'Library', 'bin')
if os.path.isdir(conda_lib) and conda_lib not in os.environ.get('PATH',''):
    os.environ['PATH'] = conda_lib + ';' + os.environ.get('PATH','')     # Windows DLL fix (PyMC needs this on PATH)
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import pymc as pm, arviz as az
from scipy.stats import t as t_dist, gaussian_kde
sys.path.insert(0, r'c:\Users\user\project')
from student_t_lm_gibbs import student_t_gibbs           # the from-scratch Gibbs, for side-by-side overlay
print('Import OK')

Import OK

# §1 synthetic data (identical to the from-scratch notebook) + run the discrete-ν Gibbs for comparison
rng1 = np.random.default_rng(0); n, k = 120, 3
beta_true = np.array([2.0, 0.8, -0.5]); sigma2_true, nu_true = 1.0, 4.0
X1 = np.column_stack([np.ones(n), rng1.standard_normal((n, k - 1))])
y1 = X1 @ beta_true + t_dist.rvs(df=nu_true, scale=np.sqrt(sigma2_true), size=n, random_state=rng1)
out_t = student_t_gibbs(y1, X1, R=15000, burn=3000, seed=1, nprint=0)
names1 = ['β[0]', 'β[1]', 'β[2]']
print('discrete-ν Gibbs: ν mean = %.1f' % out_t['nu'].mean())

Done in 0.0 min  |  kept draws: 12000
discrete-ν Gibbs: ν mean = 8.6

import pymc as pm
import arviz as az

# Use the same §1 data: y1, X1, beta_true, nu_true
with pm.Model() as model_t:
    # Priors
    beta_pm = pm.Flat('beta', shape=k)          # flat prior on β
    sigma   = pm.HalfNormal('sigma', sigma=2.0) # half-normal on σ (>0)
    nu_pm   = pm.Gamma('nu', alpha=2.0, beta=0.1)  # weakly informative; mean=20

    mu = pm.math.dot(X1, beta_pm)

    # Likelihood: Student-t with continuous ν
    obs = pm.StudentT('obs', nu=nu_pm, mu=mu, sigma=sigma, observed=y1)

    # Sample
    idata = pm.sample(
        draws=4000, tune=2000, chains=2,
        target_accept=0.9,
        random_seed=42, progressbar=False,
    )

print(az.summary(idata, var_names=['beta', 'sigma', 'nu']))

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\threadpoolctl.py:1226: RuntimeWarning: 
Found Intel OpenMP ('libiomp') and LLVM OpenMP ('libomp') loaded at
the same time. Both libraries are known to be incompatible and this
can cause random crashes or deadlocks on Linux when loaded in the
same Python program.
Using threadpoolctl may cause crashes or deadlocks. For more
information and possible workarounds, please see
    https://github.com/joblib/threadpoolctl/blob/master/multiple_openmp.md

  warnings.warn(msg, RuntimeWarning)
Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (2 chains in 2 jobs)

NUTS: [beta, sigma, nu]

Sampling 2 chains for 2_000 tune and 4_000 draw iterations (4_000 + 8_000 draws total) took 5 seconds.

We recommend running at least 4 chains for robust computation of convergence diagnostics

           mean     sd eti89_lb eti89_ub  ess_bulk  ess_tail r_hat mcse_mean  \
beta[0]   2.031  0.109      1.9      2.2      8991      6147  1.00    0.0012   
beta[1]   0.791  0.118      0.6     0.98      8638      5822  1.00    0.0013   
beta[2]  -0.664  0.107    -0.83    -0.49      7637      5495  1.00    0.0012   
sigma     1.062  0.126     0.87      1.3      5059      5468  1.00    0.0018   
nu          8.9      7      3.1       22      5132      5356  1.00     0.099   

         mcse_sd  
beta[0]  0.00081  
beta[1]  0.00094  
beta[2]  0.00089  
sigma     0.0012  
nu          0.19

# ── Side-by-side comparison: Gibbs (discrete) vs PyMC (continuous) ────────────
nu_gibbs = out_t['nu']
beta_gibbs = out_t['beta']

# PyMC posterior draws (flatten chains)
nu_pymc   = idata.posterior['nu'].values.ravel()
beta_pymc = idata.posterior['beta'].values.reshape(-1, k)
sigma_pymc = idata.posterior['sigma'].values.ravel()

# β comparison table
names1 = ['β[0]', 'β[1]', 'β[2]']
print('β posteriors  (true = [2.0, 0.8, -0.5])\n')
print(f'{"":6s} {"Gibbs mean":>12s} {"Gibbs sd":>10s} {"PyMC mean":>12s} {"PyMC sd":>10s}')
for j, nm in enumerate(names1):
    gm, gs = beta_gibbs[:, j].mean(), beta_gibbs[:, j].std()
    pm_, ps = beta_pymc[:, j].mean(), beta_pymc[:, j].std()
    print(f'{nm:6s} {gm:12.3f} {gs:10.3f} {pm_:12.3f} {ps:10.3f}')

print(f'\nσ² posterior:')
print(f'  Gibbs σ²  mean = {out_t["sigma2"].mean():.3f}  sd = {out_t["sigma2"].std():.3f}')
print(f'  PyMC  σ²  mean = {(sigma_pymc**2).mean():.3f}  sd = {(sigma_pymc**2).std():.3f}  (true σ²=1.0)')

print(f'\nν posterior:')
print(f'  Gibbs (discrete)  mean={nu_gibbs.mean():.1f}  median={np.median(nu_gibbs):.0f}'
      f'  mode={float(pd.Series(nu_gibbs).mode().iloc[0]):.0f}')
print(f'  PyMC  (continuous) mean={nu_pymc.mean():.1f}  median={np.median(nu_pymc):.1f}'
      f'  q025={np.quantile(nu_pymc, 0.025):.1f}  q975={np.quantile(nu_pymc, 0.975):.1f}')
print(f'  True ν = {nu_true:.0f}')

β posteriors  (true = [2.0, 0.8, -0.5])

         Gibbs mean   Gibbs sd    PyMC mean    PyMC sd
β[0]          2.028      0.111        2.031      0.109
β[1]          0.794      0.116        0.791      0.118
β[2]         -0.663      0.108       -0.664      0.107

σ² posterior:
  Gibbs σ²  mean = 1.102  sd = 0.288
  PyMC  σ²  mean = 1.143  sd = 0.270  (true σ²=1.0)

ν posterior:
  Gibbs (discrete)  mean=8.6  median=6  mode=4
  PyMC  (continuous) mean=8.9  median=6.6  q025=2.7  q975=28.4
  True ν = 4

fig, axes = plt.subplots(1, 4, figsize=(16, 4))

# β[0], β[1], β[2] — overlaid posteriors
for j, ax in enumerate(axes[:3]):
    lo = min(np.percentile(beta_gibbs[:, j], 0.5),
             np.percentile(beta_pymc[:, j], 0.5))
    hi = max(np.percentile(beta_gibbs[:, j], 99.5),
             np.percentile(beta_pymc[:, j], 99.5))
    bins = np.linspace(lo, hi, 60)
    ax.hist(beta_gibbs[:, j], bins=bins, density=True,
            alpha=0.55, color='orange', label='Gibbs (discrete ν)')
    ax.hist(beta_pymc[:, j],  bins=bins, density=True,
            alpha=0.45, color='steelblue', label='PyMC (continuous ν)', histtype='step', lw=2)
    ax.axvline(beta_true[j], color='k', ls='--', lw=1.5, label=f'True {beta_true[j]}')
    ax.set_title(names1[j]); ax.legend(fontsize=7)

# ν posteriors — discrete bar vs continuous KDE
ax = axes[3]
nu_vals, nu_cnt = np.unique(nu_gibbs, return_counts=True)
ax.bar(nu_vals, nu_cnt / len(nu_gibbs), width=0.7, color='orange',
       alpha=0.7, label='Gibbs (discrete)', zorder=2)

# KDE for PyMC continuous ν
from scipy.stats import gaussian_kde
kde = gaussian_kde(nu_pymc, bw_method=0.3)
nu_grid_plot = np.linspace(1, 40, 300)
ax.plot(nu_grid_plot, kde(nu_grid_plot), color='steelblue', lw=2,
        label='PyMC (continuous KDE)')
ax.axvline(nu_true, color='k', ls='--', lw=1.5, label=f'True ν={nu_true:.0f}')
ax.set_xlabel('ν'); ax.set_title('ν posterior: Gibbs vs PyMC')
ax.legend(fontsize=7); ax.set_xlim(1, 30)

fig.suptitle('§3 Comparison — Geweke Gibbs vs PyMC NUTS on §1 synthetic data\n'
             'β posteriors agree; ν — discrete Gibbs vs continuous KDE', fontsize=10)
plt.tight_layout(); plt.show()

# Nelson-Plosser load (focus series) + discrete-ν Gibbs modes for the comparison titles
np_df = pd.read_csv(r'c:\Users\user\project\nelplo.csv')
log_series = ['cpi','ip','gnp.nom','vel','emp','nom.wages','gnp.def','money.stock','gnp.real','stock.prices','gnp.capita','real.wages']
series_labels = {'cpi':'CPI','ip':'Industrial production','gnp.nom':'Nominal GNP','vel':'Velocity','emp':'Employment',
                 'int.rate':'Interest rate','nom.wages':'Nominal wages','gnp.def':'GNP deflator','money.stock':'Money stock',
                 'gnp.real':'Real GNP','stock.prices':'Stock prices','gnp.capita':'GNP per capita','real.wages':'Real wages','unemp':'Unemployment rate'}
np_df = np_df[(np_df['year']>=1909)&(np_df['year']<=1970)].copy()
for s in log_series: np_df[s] = np.log(np_df[s])
focus_series = ['gnp.real', 'int.rate', 'cpi']
results_np = {}
for s in focus_series:
    sub = np_df[['year', s]].dropna(); y_s = sub[s].values
    t_s = (sub['year'].values - sub['year'].mean()) / sub['year'].std(); X_s = np.column_stack([np.ones(len(y_s)), t_s])
    o = student_t_gibbs(y_s, X_s, R=12000, burn=3000, seed=hash(s) % (2**31), nprint=0)
    results_np[s] = {'nu_mode': float(pd.Series(o['nu']).mode().iloc[0])}
print('focus series Gibbs ν modes:', {s: results_np[s]['nu_mode'] for s in focus_series})

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000
focus series Gibbs ν modes: {'gnp.real': 2.0, 'int.rate': 6.0, 'cpi': 29.0}

focus_series = ['gnp.real', 'int.rate', 'cpi']

pymc_np_results = {}

for s in focus_series:
    sub = np_df[['year', s]].dropna()
    y_s = sub[s].values
    t_s = (sub['year'].values - sub['year'].mean()) / sub['year'].std()
    X_s = np.column_stack([np.ones(len(y_s)), t_s])
    n_s = len(y_s)

    with pm.Model() as m:
        beta_p  = pm.Flat('beta', shape=2)
        sigma_p = pm.HalfNormal('sigma', sigma=2.0)
        nu_p    = pm.Gamma('nu', alpha=2.0, beta=0.1)
        mu_p    = pm.math.dot(X_s, beta_p)
        pm.StudentT('obs', nu=nu_p, mu=mu_p, sigma=sigma_p, observed=y_s)
        idata_s = pm.sample(4000, tune=2000, chains=2, target_accept=0.9,
                            random_seed=hash(s) % (2**31),
                            progressbar=False)

    nu_draws = idata_s.posterior['nu'].values.ravel()
    pymc_np_results[s] = {
        'nu_mean':   nu_draws.mean(),
        'nu_median': np.median(nu_draws),
        'nu_q025':   np.quantile(nu_draws, 0.025),
        'nu_q975':   np.quantile(nu_draws, 0.975),
        'nu_draws':  nu_draws,
    }
    r = pymc_np_results[s]
    print(f'{series_labels[s]:20s}  '
          f'PyMC ν  mean={r["nu_mean"]:.1f}  median={r["nu_median"]:.1f}'
          f'  95% CrI [{r["nu_q025"]:.1f}, {r["nu_q975"]:.1f}]  |  '
          f'Gibbs mode={results_np[s]["nu_mode"]:.0f}')

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\threadpoolctl.py:1226: RuntimeWarning: 
Found Intel OpenMP ('libiomp') and LLVM OpenMP ('libomp') loaded at
the same time. Both libraries are known to be incompatible and this
can cause random crashes or deadlocks on Linux when loaded in the
same Python program.
Using threadpoolctl may cause crashes or deadlocks. For more
information and possible workarounds, please see
    https://github.com/joblib/threadpoolctl/blob/master/multiple_openmp.md

  warnings.warn(msg, RuntimeWarning)
Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (2 chains in 2 jobs)

NUTS: [beta, sigma, nu]

Sampling 2 chains for 2_000 tune and 4_000 draw iterations (4_000 + 8_000 draws total) took 5 seconds.

We recommend running at least 4 chains for robust computation of convergence diagnostics

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\threadpoolctl.py:1226: RuntimeWarning: 
Found Intel OpenMP ('libiomp') and LLVM OpenMP ('libomp') loaded at
the same time. Both libraries are known to be incompatible and this
can cause random crashes or deadlocks on Linux when loaded in the
same Python program.
Using threadpoolctl may cause crashes or deadlocks. For more
information and possible workarounds, please see
    https://github.com/joblib/threadpoolctl/blob/master/multiple_openmp.md

  warnings.warn(msg, RuntimeWarning)
Initializing NUTS using jitter+adapt_diag...

Real GNP              PyMC ν  mean=3.8  median=2.5  95% CrI [1.2, 16.5]  |  Gibbs mode=2

Multiprocess sampling (2 chains in 2 jobs)

NUTS: [beta, sigma, nu]

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

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

for ax, s in zip(axes, focus_series):
    # Gibbs discrete bars
    sub = np_df[['year', s]].dropna()
    y_s = sub[s].values
    t_s = (sub['year'].values - sub['year'].mean()) / sub['year'].std()
    X_s = np.column_stack([np.ones(len(y_s)), t_s])
    out_g = student_t_gibbs(y_s, X_s, R=12000, burn=3000,
                            seed=hash(s) % (2**31), nprint=0)
    nu_g = out_g['nu']
    nu_vals, nu_cnt = np.unique(nu_g, return_counts=True)
    ax.bar(nu_vals, nu_cnt / len(nu_g), width=0.7,
           color='orange', alpha=0.75, label='Gibbs (discrete)', zorder=2)

    # PyMC continuous KDE
    nu_p = pymc_np_results[s]['nu_draws']
    kde  = gaussian_kde(nu_p, bw_method=0.25)
    xg   = np.linspace(1, 42, 400)
    ax.plot(xg, kde(xg), color='steelblue', lw=2,
            label='PyMC (continuous KDE)')

    ax.set_xlim(1, 42); ax.set_xticks([5, 15, 25, 35])
    ax.set_xlabel('ν')
    ax.set_title(f'{series_labels[s]}\n'
                 f'Gibbs mode={results_np[s]["nu_mode"]:.0f}  '
                 f'PyMC median={pymc_np_results[s]["nu_median"]:.1f}')
    ax.legend(fontsize=8)

fig.suptitle('§3b Nelson-Plosser — ν posterior: Gibbs (orange bars) vs PyMC NUTS (blue KDE)\n'
             'Three representative series spanning heavy-tailed → near-Normal', fontsize=11)
plt.tight_layout(); plt.show()

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

	Discrete Gibbs (from-scratch)	PyMC / NUTS (here)
ν	discrete {2,…,40}	continuous $(0,\infty)$, Gamma(2, 0.1) prior
sampler	4-block Gibbs	NUTS (HMC)
latent λᵢ	explicit (outlier weights)	marginalised out

Bayesian Student-t Linear Model — PyMC (continuous ν)¶

1. Synthetic data: continuous-ν NUTS vs discrete-ν Gibbs¶

2. Nelson-Plosser: ν posterior, Gibbs vs PyMC (three representative series)¶

Results¶