Marginal Posteriors

from pyinla import pyinla

# Gaussian approximation (fastest)
res_ga = pyinla(
    model=model, family='poisson', data=df,
    control={'inla': {'strategy': 'gaussian'}}
)

# Simplified Laplace (default when strategy='auto')
res_sla = pyinla(
    model=model, family='poisson', data=df,
    control={'inla': {'strategy': 'simplified.laplace'}}
)

# Full Laplace (most accurate)
res_la = pyinla(
    model=model, family='poisson', data=df,
    control={'inla': {'strategy': 'laplace'}}
)

# Adaptive (hybrid approach)
res_adaptive = pyinla(
    model=model, family='poisson', data=df,
    control={'inla': {'strategy': 'adaptive', 'adaptive_max': 15}}
)

import copy
import numpy as np
import pandas as pd
from pyinla import pyinla

# Synthetic data: y = 1 + 0.5 x + noise
rng = np.random.default_rng(0)
n = 200
x = rng.normal(size=n)
y = 1.0 + 0.5 * x + rng.normal(0, 0.3, size=n)
df = pd.DataFrame({"y": y, "x": x})

model = {"response": "y", "fixed": ["1", "x"]}

# Default: VB correction enabled with mean strategy
res = pyinla(
    model=copy.deepcopy(model), family='gaussian', data=df,
    control={
        'inla': {
            'control_vb': {
                'enable': True,           # Enable VB correction (default: auto)
                'strategy': 'mean',       # 'mean' or 'variance'
                'iter_max': 25,           # Maximum iterations
                'hessian_strategy': 'full'  # 'full', 'partial', 'diagonal'
            }
        }
    }
)

# Disable VB correction
res_no_vb = pyinla(
    model=copy.deepcopy(model), family='gaussian', data=df,
    control={'inla': {'control_vb': {'enable': False}}}
)

from pyinla import pyinla

# CCD integration (auto resolves to this for typical models)
res_ccd = pyinla(
    model=model, family='gaussian', data=df,
    control={'inla': {'int_strategy': 'ccd'}}
)

# Grid integration (thorough)
res_grid = pyinla(
    model=model, family='gaussian', data=df,
    control={'inla': {'int_strategy': 'grid'}}
)

# Empirical Bayes (fastest, mode only)
res_eb = pyinla(
    model=model, family='gaussian', data=df,
    control={'inla': {'int_strategy': 'eb'}}
)

# User-supplied design: each row is (theta_1, ..., theta_p, weight).
# Below is a minimal 3-point design for a single-hyperparameter model.
import numpy as np
design = np.array([
    [-1.0, 1.0],     # theta = -1 (log-precision), weight = 1
    [ 0.0, 2.0],     # theta =  0,                 weight = 2
    [ 1.0, 1.0],     # theta =  1,                 weight = 1
])
res_user = pyinla(
    model=model, family='gaussian', data=df,
    control={'inla': {'int_strategy': 'user', 'int_design': design}}
)

# Fast: Gaussian + Empirical Bayes
control_fast = {
    'inla': {
        'strategy': 'gaussian',
        'int_strategy': 'eb'
    }
}

# Accurate: Laplace + Grid
control_accurate = {
    'inla': {
        'strategy': 'laplace',
        'int_strategy': 'grid'
    }
}

# Balanced (default): Simplified Laplace + CCD
control_balanced = {
    'inla': {
        'strategy': 'simplified.laplace',
        'int_strategy': 'ccd'
    }
}

from pyinla import pyinla

res = pyinla(model=model, family='gaussian', data=df)

# Fixed effects (regression coefficients)
for name, marg in res.marginals_fixed.items():
    print(f"{name}: shape {marg.shape}")

# Hyperparameters (precision, variance, etc.)
for name, marg in res.marginals_hyperpar.items():
    print(f"{name}: shape {marg.shape}")

# Random effects (if model has them)
if res.marginals_random:
    for effect_name, marginals in res.marginals_random.items():
        print(f"{effect_name}: {len(marginals)} marginals")

# Linear predictor (fitted values)
if res.marginals_linear_predictor is not None:
    print(f"Linear predictor: {len(res.marginals_linear_predictor)} marginals")

import numpy as np
import pandas as pd
import pyinla

# --- Generate a small synthetic dataset --------------------------------
rng = np.random.default_rng(0)
n  = 200
x1 = rng.normal(size=n)
x2 = rng.normal(size=n)
y  = 1.0 + 2.0 * x1 - 0.5 * x2 + rng.normal(scale=0.3, size=n)
data = pd.DataFrame({"y": y, "x1": x1, "x2": x2})

# --- Fit a Gaussian linear model (dict-style spec) ---------------------
res = pyinla.pyinla(
    formula={'response': 'y', 'fixed': ['1', 'x1', 'x2']},
    family='gaussian',
    data=data,
)

# --- Access marginals --------------------------------------------------
marg_intercept = res.marginals_fixed["(Intercept)"]
marg_x1        = res.marginals_fixed["x1"]
# Hyperparameter keys use INLA's descriptive labels:
marg_precision = res.marginals_hyperpar["Precision for the Gaussian observations"]

# Each marginal is a pandas DataFrame with two columns: 'x' (parameter values) and 'y' (density)
print(marg_x1.shape)   # e.g., (43, 2)
print(marg_x1.head())  # First 5 rows of the DataFrame

import pyinla

# All functions available as pyinla.
pyinla.dmarginal(0.5, marg)   # Density at x=0.5
pyinla.pmarginal(0.5, marg)   # P(X <= 0.5)
pyinla.qmarginal(0.95, marg)  # 95th percentile

import pyinla
import numpy as np
import matplotlib.pyplot as plt

marg = res.marginals_fixed["x1"]   # DataFrame with columns ['x', 'y']

# Evaluate at a single point
d = pyinla.dmarginal(0.5, marg)
print(f"Density at 0.5: {d:.4f}")

# Log-density variant (useful for numerical stability and likelihood ratios)
log_d = pyinla.dmarginal(0.5, marg, log=True)

# Evaluate at multiple points for plotting
x_grid = np.linspace(marg["x"].min(), marg["x"].max(), 200)
y_grid = pyinla.dmarginal(x_grid, marg)

plt.figure(figsize=(8, 5))
plt.plot(x_grid, y_grid, 'b-', linewidth=2, label='Interpolated')
plt.scatter(marg["x"], marg["y"], s=20, c='red', alpha=0.5, label='Original points')
plt.xlabel(r'$\theta$')
plt.ylabel(r'$\pi(\theta \mid y)$')
plt.legend()
plt.title('Posterior Marginal Density')
plt.show()

import pyinla

marg = res.marginals_fixed["x1"]

# Probability that coefficient is negative
p_negative = pyinla.pmarginal(0, marg)
print(f"P(x1 < 0) = {p_negative:.4f}")

# Probability that coefficient is in interval [0.2, 0.8]
p_interval = pyinla.pmarginal(0.8, marg) - pyinla.pmarginal(0.2, marg)
print(f"P(0.2 < x1 < 0.8) = {p_interval:.4f}")

# One-sided credible interval
if p_negative < 0.05:
    print("Strong evidence that x1 > 0 (95% credible)")

import pyinla

marg = res.marginals_fixed["x1"]

# Median
median = pyinla.qmarginal(0.5, marg)
print(f"Median: {median:.4f}")

# 95% equal-tailed credible interval
q025 = pyinla.qmarginal(0.025, marg)
q975 = pyinla.qmarginal(0.975, marg)
print(f"95% CI: [{q025:.4f}, {q975:.4f}]")

# Multiple quantiles at once
quantiles = pyinla.qmarginal([0.025, 0.25, 0.5, 0.75, 0.975], marg)
print(f"Quantiles: {quantiles}")

import pyinla
import numpy as np

marg = res.marginals_fixed["x1"]

# Generate samples
samples = pyinla.rmarginal(10000, marg)

# For reproducibility
rng = np.random.default_rng(42)
samples_reproducible = pyinla.rmarginal(10000, marg, rng=rng)

# Monte Carlo estimates
print(f"Sample mean: {samples.mean():.4f}")
print(f"Sample std:  {samples.std():.4f}")
print(f"P(X > 0):    {(samples > 0).mean():.4f}")

import pyinla
import numpy as np

marg = res.marginals_fixed["x1"]

# Posterior mean (E[X])
mean = pyinla.emarginal(lambda x: x, marg)
print(f"Mean: {mean:.4f}")

# Posterior variance (E[X^2] - E[X]^2)
ex2 = pyinla.emarginal(lambda x: x**2, marg)
variance = ex2 - mean**2
print(f"Variance: {variance:.4f}")

# Expected value of exp(theta) - useful for log-scale parameters
e_exp = pyinla.emarginal(np.exp, marg)
print(f"E[exp(theta)]: {e_exp:.4f}")

import pyinla
import numpy as np

# pyinla returns marginals on the user (precision) scale, not log-scale,
# so to operate on log(precision) first transform back:
marg_prec     = res.marginals_hyperpar["Precision for the Gaussian observations"]
marg_log_prec = pyinla.tmarginal(np.log, marg_prec)

# Or transform precision to standard deviation directly: sd = 1 / sqrt(precision)
marg_sd = pyinla.tmarginal(lambda x: 1.0 / np.sqrt(x), marg_prec)

# Increase the internal sampling resolution for a smoother result
marg_sd_hi = pyinla.tmarginal(
    fun      = lambda x: 1.0 / np.sqrt(x),
    marginal = marg_prec,
    n        = 4096,
)

import pyinla

marg = res.marginals_fixed["x1"]

# Original might have ~50 points; smooth via spline interpolation
# factor controls density: ~factor x the number of input points.
marg_smooth = pyinla.smarginal(marg, factor=15)

print(f"Original points: {len(marg)}")
print(f"Smoothed points: {len(marg_smooth['x'])}")

import pyinla

marg = res.marginals_fixed["x1"]

# Mode (MAP estimate)
mode = pyinla.mmarginal(marg)

# Compare with mean and median
mean = pyinla.emarginal(lambda x: x, marg)
median = pyinla.qmarginal(0.5, marg)

print(f"Mode:   {mode:.4f}")
print(f"Mean:   {mean:.4f}")
print(f"Median: {median:.4f}")

import pyinla

marg = res.marginals_fixed["x1"]

# Get full summary
summary = pyinla.zmarginal(marg)

print("Posterior Summary for x1:")
print(f"  Mean:   {summary['mean']:.4f}")
print(f"  SD:     {summary['sd']:.4f}")
print(f"  95% CI: [{summary['quant0.025']:.4f}, {summary['quant0.975']:.4f}]")

# Summarize all fixed effects
for name, marg in res.marginals_fixed.items():
    s = pyinla.zmarginal(marg)
    print(f"{name:15s}: {s['mean']:8.4f} ({s['quant0.025']:8.4f}, {s['quant0.975']:8.4f})")

import numpy as np
import pandas as pd
import pyinla

# Synthetic data: y = 1 + 0.5 x1 + noise
rng = np.random.default_rng(0)
n = 200
x1 = rng.normal(size=n)
y = 1.0 + 0.5 * x1 + rng.normal(0, 0.3, size=n)
df = pd.DataFrame({"y": y, "x1": x1})

res = pyinla.pyinla(
    model={"response": "y", "fixed": ["1", "x1"]},
    family="gaussian", data=df,
)

marg = res.marginals_fixed["x1"]

# 95% HPD interval. hpdmarginal returns shape (len(p), 2); ravel for a flat [low, high].
hpd = pyinla.hpdmarginal(0.95, marg).ravel()
print(f"95% HPD:         [{hpd[0]:.4f}, {hpd[1]:.4f}]")

# 95% equal-tailed interval via the 2.5% / 97.5% quantiles.
q = pyinla.qmarginal([0.025, 0.975], marg)
print(f"95% Equal-tailed: [{q[0]:.4f}, {q[1]:.4f}]")

# HPD is ≤ equal-tailed in width; equal for symmetric marginals, shorter for skewed ones.
print(f"HPD width:         {hpd[1] - hpd[0]:.4f}")
print(f"Equal-tailed width: {q[1] - q[0]:.4f}")