This page presents Python code for implementing Bayesian conjugate models using standard prior–likelihood formulations and their associated posterior distributions.
- Binomial Beta model
- Poisson Gamma model
- Normal Normal model
- Normal Gamma model
- Normal Inverse gamma model
For preliminary theoretical part, refer
The choice of parameters should always follow the native argument convention of the corresponding R function. Several probability distributions, such as the Gamma distribution, permit alternative but equivalent parameterizations (e.g., shape–scale or shape–rate). Consequently, we should exercise appropriate care while specifying parameters to ensure consistency with the underlying R implementation.
The implementations use closed-form analytical expressions for posterior summaries wherever available, thereby avoiding numerical approximation except where necessary (e.g., posterior quantiles).
Numerical output includes posterior summaries such as the mean, median, and mode (where defined and computable), with user-selectable posterior quantiles and credible intervals of any desired probability level. Appropriate safeguards have been incorporated for cases in which the posterior mode does not exist or cannot be determined analytically.
Graphical output includes the corresponding prior and posterior density curves to facilitate visual comparison. All prior distribution parameters are fully user-configurable, allowing the models to be adapted to different prior beliefs and applications.
Python Code is
Bayesian Conjugate Model – Binomial-Beta
Show code
# =============================================================
# Bayesian Conjugate Model – Binomial-Beta
# Closed-form analytical solution (standalone script)
#
# Likelihood : X | theta ~ Binomial(n, theta)
# Prior : theta ~ Beta(alpha, beta)
# Posterior : theta | x ~ Beta(alpha + x, beta + n – x)
#
# Input : scalar / list / vector OR DataFrame + var + group_var
# Output : clean console tables, Plotly HTML saved to disk
# Plots : Plotly only — saved as HTML, never auto-rendered
# =============================================================
import os
import math
import numpy as np
import pandas as pd
from scipy.stats import beta as sp_beta
import plotly.graph_objects as go
from plotly.subplots import make_subplots
# =============================================================
# INTERNAL UTILITIES
# =============================================================
def _print_table(title: str, df: pd.DataFrame) -> None:
“””Plain text table to console — no color, no clutter.”””
col_names = list(df.columns)
col_widths = [
max(len(nm), df[nm].astype(str).map(len).max())
for nm in col_names
]
col_widths = [max(w, 10) for w in col_widths]
total_width = sum(col_widths) + 3 * (len(col_widths) – 1) + 2
divider = “=” * total_width
thin_div = “-” * total_width
header = ” “.join(nm.ljust(w) for nm, w in zip(col_names, col_widths))
lines = []
lines.append(divider)
lines.append(f” {title}”)
lines.append(divider)
lines.append(f” {header}”)
lines.append(f” {thin_div}”)
for _, row in df.iterrows():
row_str = ” “.join(str(val).ljust(w) for val, w in zip(row, col_widths))
lines.append(f” {row_str}”)
lines.append(divider)
lines.append(“”)
print(“\n”.join(lines))
def _fmt_mode(v) -> str:
return “undefined” if v is None else f”{v:.6f}”
def _expand_cri(cri_levels, quantiles):
“””Expand cri_levels to symmetric tail quantile pairs; merge with raw quantiles.”””
q_all = sorted(set(quantiles)) if quantiles else []
if cri_levels:
for lev in sorted(set(cri_levels)):
lo = (1 – lev / 100) / 2
hi = 1 – lo
q_all = sorted(set(q_all + [lo, hi]))
return q_all
def _build_cri_table(cri_levels, all_q, q_vals, cri_type):
“””Build CrI table from cri_levels only — symmetric pairs.”””
if not cri_levels:
return None
rows = []
for lev in sorted(set(cri_levels)):
lo = (1 – lev / 100) / 2
hi = 1 – lo
idx_lo = next((i for i, q in enumerate(all_q) if abs(q – lo) < 1e-10), None)
idx_hi = next((i for i, q in enumerate(all_q) if abs(q - hi) < 1e-10), None)
if idx_lo is not None and idx_hi is not None:
rows.append({
"Level": f"{lev}% {cri_type}",
"Lower": f"{q_vals[idx_lo]:.6f}",
"Upper": f"{q_vals[idx_hi]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _build_quantile_table(quantiles, all_q, q_vals):
"""Build individual quantile table from quantiles arg only."""
if not quantiles:
return None
rows = []
for q in sorted(set(quantiles)):
idx = next((i for i, v in enumerate(all_q) if abs(v - q) < 1e-10), None)
if idx is not None:
rows.append({
"Quantile": f"Q({q:.4f})",
"Value": f"{q_vals[idx]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _validate_df_var(data: pd.DataFrame, var: str) -> None:
“””Validate DataFrame and column for the Binomial-Beta model.”””
if not isinstance(data, pd.DataFrame):
raise TypeError(“`data` must be a pandas DataFrame”)
if not isinstance(var, str):
raise TypeError(“`var` must be a single column name (string)”)
if var not in data.columns:
raise ValueError(f”Column ‘{var}’ not found in data”)
col = data[var].dropna()
if col.empty:
raise ValueError(f”Column ‘{var}’ is entirely NA”)
uvals = col.unique()
if len(uvals) > 2:
raise ValueError(
f”Column ‘{var}’ has {len(uvals)} unique values. ”
f”Binomial model requires a binary variable (exactly 2 levels).”
)
if len(uvals) < 2:
raise ValueError(
f"Column '{var}' has only one unique non-NA value — cannot fit Binomial model"
)
def _validate_group(data: pd.DataFrame, group_var) -> None:
if group_var is not None and group_var not in data.columns:
raise ValueError(f”group_var column ‘{group_var}’ not found in data”)
def _split_entries(data: pd.DataFrame, var: str, group_var):
“””Split data into labelled entries by group.”””
if group_var is not None:
entries = []
for g in data[group_var].unique():
sub = data[data[group_var] == g]
entries.append({
“label”: f”{group_var} = {g}”,
“col”: sub[var].dropna().values
})
return entries
else:
return [{“label”: “All”, “col”: data[var].dropna().values}]
def _resolve_plot_path(plot_path, plot_filename, default_name: str) -> str:
“””Resolve final HTML save path.”””
directory = plot_path if plot_path is not None else os.getcwd()
filename = (plot_filename if plot_filename is not None else default_name)
if not filename.endswith(“.html”):
filename += “.html”
return os.path.join(directory, filename)
def _save_plot(fig, path: str) -> None:
fig.write_html(path, auto_open=False)
print(f” [Plot saved] {path}\n”)
# =============================================================
# BINOMIAL – BETA
# =============================================================
def conjugate_binomial_beta(
data = None,
var = None,
group_var = None,
success_level = None,
x = None,
n = None,
alpha = None,
beta = None,
cri_levels = (90, 95),
quantiles = None,
cri_type = “CrI”,
plot_path = None,
plot_filename = None,
):
cri_levels = list(cri_levels) if cri_levels else []
quantiles = list(quantiles) if quantiles else []
all_q = _expand_cri(cri_levels, quantiles)
# — Build entries —
if data is not None:
_validate_df_var(data, var)
_validate_group(data, group_var)
entries_raw = _split_entries(data, var, group_var)
entries = []
for e in entries_raw:
col = e[“col”]
uvals = sorted(set(col.astype(str)))
sl = str(success_level) if success_level is not None else uvals[0]
if sl not in [str(v) for v in col]:
raise ValueError(f”success_level ‘{sl}’ not found in column ‘{var}'”)
xv = int(np.sum(col.astype(str) == sl))
nv = len(col)
entries.append({“label”: e[“label”], “x”: xv, “n”: nv, “sl”: sl})
else:
if x is None or n is None:
raise ValueError(“Provide `data` + `var` OR both `x` and `n`”)
x_arr = np.atleast_1d(np.asarray(x, dtype=float))
n_arr = np.atleast_1d(np.asarray(n, dtype=float))
if len(n_arr) == 1:
n_arr = np.repeat(n_arr, len(x_arr))
if len(x_arr) != len(n_arr):
raise ValueError(“`x` and `n` must be the same length”)
entries = [
{“label”: str(i + 1), “x”: x_arr[i], “n”: n_arr[i], “sl”: “1”}
for i in range(len(x_arr))
]
results = []
for e in entries:
xv, nv = e[“x”], e[“n”]
a_post = alpha + xv
b_post = beta + (nv – xv)
pr_mean = alpha / (alpha + beta)
pr_var = (alpha * beta) / ((alpha + beta)**2 * (alpha + beta + 1))
pr_mode = (alpha – 1) / (alpha + beta – 2) if (alpha > 1 and beta > 1) else None
pr_median = sp_beta.ppf(0.5, alpha, beta)
po_mean = a_post / (a_post + b_post)
po_var = (a_post * b_post) / ((a_post + b_post)**2 * (a_post + b_post + 1))
po_mode = (a_post – 1) / (a_post + b_post – 2) if (a_post > 1 and b_post > 1) else None
po_median = sp_beta.ppf(0.5, a_post, b_post)
all_q_int = sorted(set(all_q + [0.5]))
q_vals = [sp_beta.ppf(q, a_post, b_post) for q in all_q_int]
_print_table(
f”Binomial-Beta | {e[‘label’]} | Data Summary”,
pd.DataFrame([{
“Successes”: int(xv),
“Trials”: int(nv),
“Observed Prop”: f”{xv/nv:.6f}”,
“Success Level”: str(e[“sl”]),
}])
)
_print_table(
f”Binomial-Beta | {e[‘label’]} | Parameters & Moments”,
pd.DataFrame({
“Quantity”: [“alpha”, “beta”, “Mean”, “Median”, “Mode”, “Variance”, “SD”],
“Prior”: [f”{alpha:.6f}”, f”{beta:.6f}”, f”{pr_mean:.6f}”,
f”{pr_median:.6f}”, _fmt_mode(pr_mode),
f”{pr_var:.6f}”, f”{math.sqrt(pr_var):.6f}”],
“Posterior”: [f”{a_post:.6f}”, f”{b_post:.6f}”, f”{po_mean:.6f}”,
f”{po_median:.6f}”, _fmt_mode(po_mode),
f”{po_var:.6f}”, f”{math.sqrt(po_var):.6f}”],
})
)
cri_df = _build_cri_table(cri_levels, all_q_int, q_vals, cri_type)
if cri_df is not None:
_print_table(f”Binomial-Beta | {e[‘label’]} | Credible Intervals ({cri_type})”, cri_df)
q_df = _build_quantile_table(quantiles, all_q_int, q_vals)
if q_df is not None:
_print_table(f”Binomial-Beta | {e[‘label’]} | Posterior Quantiles”, q_df)
# — Plot —
th_lo = 1e-4 if min(alpha, beta, a_post, b_post) < 1 else 0.0
th_hi = 1 - th_lo
th = np.linspace(th_lo, th_hi, 1000)
fig = make_subplots(rows=1, cols=2, subplot_titles=["Prior", "Posterior"])
fig.add_trace(go.Scatter(
x=th, y=sp_beta.pdf(th, alpha, beta),
mode="lines", line=dict(color="#2C7BB6", width=1.5), showlegend=False
), row=1, col=1)
fig.add_trace(go.Scatter(
x=th, y=sp_beta.pdf(th, a_post, b_post),
mode="lines", line=dict(color="#2C7BB6", width=1.5), showlegend=False
), row=1, col=2)
fig.update_xaxes(title_text="θ", row=1, col=1)
fig.update_xaxes(title_text="θ", row=1, col=2)
fig.update_yaxes(title_text="Density", row=1, col=1)
fig.update_layout(
title_text=(
f"Binomial-Beta [{e['label']}]
”
f”Prior: Beta({alpha:.2f}, {beta:.2f}) | ”
f”Posterior: Beta({a_post:.2f}, {b_post:.2f})”
),
template=”plotly_white”,
)
default_name = f”binomial_beta_{e[‘label’].replace(‘ ‘, ‘_’)}”
fpath = _resolve_plot_path(plot_path, plot_filename, default_name)
_save_plot(fig, fpath)
results.append({
“label”: e[“label”],
“data”: {“x”: xv, “n”: nv, “success_level”: e[“sl”]},
“prior”: {“alpha”: alpha, “beta”: beta,
“mean”: pr_mean, “median”: pr_median, “mode”: pr_mode,
“var”: pr_var, “sd”: math.sqrt(pr_var)},
“posterior”: {“alpha”: a_post, “beta”: b_post,
“mean”: po_mean, “median”: po_median, “mode”: po_mode,
“var”: po_var, “sd”: math.sqrt(po_var),
“quantiles”: dict(zip(all_q_int, q_vals))},
“plot”: fig,
})
return results[0] if len(results) == 1 else results
# =============================================================
# USAGE EXAMPLES
# =============================================================
if __name__ == “__main__”:
# scalar
# res = conjugate_binomial_beta(x=14, n=20, alpha=2, beta=0.2,
# cri_levels=(80, 90, 95, 99), quantiles=(0.75,))
# df + var (binary 0/1)
# df = pd.DataFrame({“response”: [1,0,1,1,0,1,0,0,1,1]})
# res = conjugate_binomial_beta(data=df, var=”response”,
# alpha=2, beta=2, cri_levels=(90, 95))
# df + var (character binary) + group + success_level
# df = pd.DataFrame({
# “response”: [“Yes”,”No”,”Yes”,”Yes”,”No”,”Yes”,”No”,”No”,”Yes”,”Yes”],
# “group”: [“A”]*5 + [“B”]*5
# })
# res = conjugate_binomial_beta(data=df, var=”response”, group_var=”group”,
# success_level=”Yes”, alpha=2, beta=2,
# cri_levels=(90, 95))
# —- access plot object and enhance —-
# res = conjugate_binomial_beta(x=14, n=20, alpha=2, beta=2)
# fig = res[“plot”]
# fig.update_layout(title_text=”My Custom Title”)
# fig.write_html(“my_custom_plot.html”, auto_open=False)
# —- custom plot path and filename —-
res = conjugate_binomial_beta(
x=14, n=20, alpha=2, beta=2,
cri_levels=[80, 90, 95, 99], quantiles=[0.2, 0.75],
plot_path=”.”, plot_filename=”binom_result”
)
# =============================================================
# Bayesian Conjugate Model - Binomial-Beta
# Closed-form analytical solution (standalone script)
#
# Likelihood : X | theta ~ Binomial(n, theta)
# Prior : theta ~ Beta(alpha, beta)
# Posterior : theta | x ~ Beta(alpha + x, beta + n - x)
#
# Input : scalar / list / vector OR DataFrame + var + group_var
# Output : clean console tables, Plotly HTML saved to disk
# Plots : Plotly only — saved as HTML, never auto-rendered
# =============================================================
import os
import math
import numpy as np
import pandas as pd
from scipy.stats import beta as sp_beta
import plotly.graph_objects as go
from plotly.subplots import make_subplots
# =============================================================
# INTERNAL UTILITIES
# =============================================================
def _print_table(title: str, df: pd.DataFrame) -> None:
"""Plain text table to console — no color, no clutter."""
col_names = list(df.columns)
col_widths = [
max(len(nm), df[nm].astype(str).map(len).max())
for nm in col_names
]
col_widths = [max(w, 10) for w in col_widths]
total_width = sum(col_widths) + 3 * (len(col_widths) - 1) + 2
divider = "=" * total_width
thin_div = "-" * total_width
header = " ".join(nm.ljust(w) for nm, w in zip(col_names, col_widths))
lines = []
lines.append(divider)
lines.append(f" {title}")
lines.append(divider)
lines.append(f" {header}")
lines.append(f" {thin_div}")
for _, row in df.iterrows():
row_str = " ".join(str(val).ljust(w) for val, w in zip(row, col_widths))
lines.append(f" {row_str}")
lines.append(divider)
lines.append("")
print("\n".join(lines))
def _fmt_mode(v) -> str:
return "undefined" if v is None else f"{v:.6f}"
def _expand_cri(cri_levels, quantiles):
"""Expand cri_levels to symmetric tail quantile pairs; merge with raw quantiles."""
q_all = sorted(set(quantiles)) if quantiles else []
if cri_levels:
for lev in sorted(set(cri_levels)):
lo = (1 - lev / 100) / 2
hi = 1 - lo
q_all = sorted(set(q_all + [lo, hi]))
return q_all
def _build_cri_table(cri_levels, all_q, q_vals, cri_type):
"""Build CrI table from cri_levels only — symmetric pairs."""
if not cri_levels:
return None
rows = []
for lev in sorted(set(cri_levels)):
lo = (1 - lev / 100) / 2
hi = 1 - lo
idx_lo = next((i for i, q in enumerate(all_q) if abs(q - lo) < 1e-10), None)
idx_hi = next((i for i, q in enumerate(all_q) if abs(q - hi) < 1e-10), None)
if idx_lo is not None and idx_hi is not None:
rows.append({
"Level": f"{lev}% {cri_type}",
"Lower": f"{q_vals[idx_lo]:.6f}",
"Upper": f"{q_vals[idx_hi]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _build_quantile_table(quantiles, all_q, q_vals):
"""Build individual quantile table from quantiles arg only."""
if not quantiles:
return None
rows = []
for q in sorted(set(quantiles)):
idx = next((i for i, v in enumerate(all_q) if abs(v - q) < 1e-10), None)
if idx is not None:
rows.append({
"Quantile": f"Q({q:.4f})",
"Value": f"{q_vals[idx]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _validate_df_var(data: pd.DataFrame, var: str) -> None:
"""Validate DataFrame and column for the Binomial-Beta model."""
if not isinstance(data, pd.DataFrame):
raise TypeError("`data` must be a pandas DataFrame")
if not isinstance(var, str):
raise TypeError("`var` must be a single column name (string)")
if var not in data.columns:
raise ValueError(f"Column '{var}' not found in data")
col = data[var].dropna()
if col.empty:
raise ValueError(f"Column '{var}' is entirely NA")
uvals = col.unique()
if len(uvals) > 2:
raise ValueError(
f"Column '{var}' has {len(uvals)} unique values. "
f"Binomial model requires a binary variable (exactly 2 levels)."
)
if len(uvals) < 2:
raise ValueError(
f"Column '{var}' has only one unique non-NA value — cannot fit Binomial model"
)
def _validate_group(data: pd.DataFrame, group_var) -> None:
if group_var is not None and group_var not in data.columns:
raise ValueError(f"group_var column '{group_var}' not found in data")
def _split_entries(data: pd.DataFrame, var: str, group_var):
"""Split data into labelled entries by group."""
if group_var is not None:
entries = []
for g in data[group_var].unique():
sub = data[data[group_var] == g]
entries.append({
"label": f"{group_var} = {g}",
"col": sub[var].dropna().values
})
return entries
else:
return [{"label": "All", "col": data[var].dropna().values}]
def _resolve_plot_path(plot_path, plot_filename, default_name: str) -> str:
"""Resolve final HTML save path."""
directory = plot_path if plot_path is not None else os.getcwd()
filename = (plot_filename if plot_filename is not None else default_name)
if not filename.endswith(".html"):
filename += ".html"
return os.path.join(directory, filename)
def _save_plot(fig, path: str) -> None:
fig.write_html(path, auto_open=False)
print(f" [Plot saved] {path}\n")
# =============================================================
# BINOMIAL - BETA
# =============================================================
def conjugate_binomial_beta(
data = None,
var = None,
group_var = None,
success_level = None,
x = None,
n = None,
alpha = None,
beta = None,
cri_levels = (90, 95),
quantiles = None,
cri_type = "CrI",
plot_path = None,
plot_filename = None,
):
cri_levels = list(cri_levels) if cri_levels else []
quantiles = list(quantiles) if quantiles else []
all_q = _expand_cri(cri_levels, quantiles)
# --- Build entries ---
if data is not None:
_validate_df_var(data, var)
_validate_group(data, group_var)
entries_raw = _split_entries(data, var, group_var)
entries = []
for e in entries_raw:
col = e["col"]
uvals = sorted(set(col.astype(str)))
sl = str(success_level) if success_level is not None else uvals[0]
if sl not in [str(v) for v in col]:
raise ValueError(f"success_level '{sl}' not found in column '{var}'")
xv = int(np.sum(col.astype(str) == sl))
nv = len(col)
entries.append({"label": e["label"], "x": xv, "n": nv, "sl": sl})
else:
if x is None or n is None:
raise ValueError("Provide `data` + `var` OR both `x` and `n`")
x_arr = np.atleast_1d(np.asarray(x, dtype=float))
n_arr = np.atleast_1d(np.asarray(n, dtype=float))
if len(n_arr) == 1:
n_arr = np.repeat(n_arr, len(x_arr))
if len(x_arr) != len(n_arr):
raise ValueError("`x` and `n` must be the same length")
entries = [
{"label": str(i + 1), "x": x_arr[i], "n": n_arr[i], "sl": "1"}
for i in range(len(x_arr))
]
results = []
for e in entries:
xv, nv = e["x"], e["n"]
a_post = alpha + xv
b_post = beta + (nv - xv)
pr_mean = alpha / (alpha + beta)
pr_var = (alpha * beta) / ((alpha + beta)**2 * (alpha + beta + 1))
pr_mode = (alpha - 1) / (alpha + beta - 2) if (alpha > 1 and beta > 1) else None
pr_median = sp_beta.ppf(0.5, alpha, beta)
po_mean = a_post / (a_post + b_post)
po_var = (a_post * b_post) / ((a_post + b_post)**2 * (a_post + b_post + 1))
po_mode = (a_post - 1) / (a_post + b_post - 2) if (a_post > 1 and b_post > 1) else None
po_median = sp_beta.ppf(0.5, a_post, b_post)
all_q_int = sorted(set(all_q + [0.5]))
q_vals = [sp_beta.ppf(q, a_post, b_post) for q in all_q_int]
_print_table(
f"Binomial-Beta | {e['label']} | Data Summary",
pd.DataFrame([{
"Successes": int(xv),
"Trials": int(nv),
"Observed Prop": f"{xv/nv:.6f}",
"Success Level": str(e["sl"]),
}])
)
_print_table(
f"Binomial-Beta | {e['label']} | Parameters & Moments",
pd.DataFrame({
"Quantity": ["alpha", "beta", "Mean", "Median", "Mode", "Variance", "SD"],
"Prior": [f"{alpha:.6f}", f"{beta:.6f}", f"{pr_mean:.6f}",
f"{pr_median:.6f}", _fmt_mode(pr_mode),
f"{pr_var:.6f}", f"{math.sqrt(pr_var):.6f}"],
"Posterior": [f"{a_post:.6f}", f"{b_post:.6f}", f"{po_mean:.6f}",
f"{po_median:.6f}", _fmt_mode(po_mode),
f"{po_var:.6f}", f"{math.sqrt(po_var):.6f}"],
})
)
cri_df = _build_cri_table(cri_levels, all_q_int, q_vals, cri_type)
if cri_df is not None:
_print_table(f"Binomial-Beta | {e['label']} | Credible Intervals ({cri_type})", cri_df)
q_df = _build_quantile_table(quantiles, all_q_int, q_vals)
if q_df is not None:
_print_table(f"Binomial-Beta | {e['label']} | Posterior Quantiles", q_df)
# --- Plot ---
th_lo = 1e-4 if min(alpha, beta, a_post, b_post) < 1 else 0.0
th_hi = 1 - th_lo
th = np.linspace(th_lo, th_hi, 1000)
fig = make_subplots(rows=1, cols=2, subplot_titles=["Prior", "Posterior"])
fig.add_trace(go.Scatter(
x=th, y=sp_beta.pdf(th, alpha, beta),
mode="lines", line=dict(color="#2C7BB6", width=1.5), showlegend=False
), row=1, col=1)
fig.add_trace(go.Scatter(
x=th, y=sp_beta.pdf(th, a_post, b_post),
mode="lines", line=dict(color="#2C7BB6", width=1.5), showlegend=False
), row=1, col=2)
fig.update_xaxes(title_text="θ", row=1, col=1)
fig.update_xaxes(title_text="θ", row=1, col=2)
fig.update_yaxes(title_text="Density", row=1, col=1)
fig.update_layout(
title_text=(
f"Binomial-Beta [{e['label']}]<br>"
f"<sup>Prior: Beta({alpha:.2f}, {beta:.2f}) | "
f"Posterior: Beta({a_post:.2f}, {b_post:.2f})</sup>"
),
template="plotly_white",
)
default_name = f"binomial_beta_{e['label'].replace(' ', '_')}"
fpath = _resolve_plot_path(plot_path, plot_filename, default_name)
_save_plot(fig, fpath)
results.append({
"label": e["label"],
"data": {"x": xv, "n": nv, "success_level": e["sl"]},
"prior": {"alpha": alpha, "beta": beta,
"mean": pr_mean, "median": pr_median, "mode": pr_mode,
"var": pr_var, "sd": math.sqrt(pr_var)},
"posterior": {"alpha": a_post, "beta": b_post,
"mean": po_mean, "median": po_median, "mode": po_mode,
"var": po_var, "sd": math.sqrt(po_var),
"quantiles": dict(zip(all_q_int, q_vals))},
"plot": fig,
})
return results[0] if len(results) == 1 else results
# =============================================================
# USAGE EXAMPLES
# =============================================================
if __name__ == "__main__":
# scalar
# res = conjugate_binomial_beta(x=14, n=20, alpha=2, beta=0.2,
# cri_levels=(80, 90, 95, 99), quantiles=(0.75,))
# df + var (binary 0/1)
# df = pd.DataFrame({"response": [1,0,1,1,0,1,0,0,1,1]})
# res = conjugate_binomial_beta(data=df, var="response",
# alpha=2, beta=2, cri_levels=(90, 95))
# df + var (character binary) + group + success_level
# df = pd.DataFrame({
# "response": ["Yes","No","Yes","Yes","No","Yes","No","No","Yes","Yes"],
# "group": ["A"]*5 + ["B"]*5
# })
# res = conjugate_binomial_beta(data=df, var="response", group_var="group",
# success_level="Yes", alpha=2, beta=2,
# cri_levels=(90, 95))
# ---- access plot object and enhance ----
# res = conjugate_binomial_beta(x=14, n=20, alpha=2, beta=2)
# fig = res["plot"]
# fig.update_layout(title_text="My Custom Title")
# fig.write_html("my_custom_plot.html", auto_open=False)
# ---- custom plot path and filename ----
res = conjugate_binomial_beta(
x=14, n=20, alpha=2, beta=2,
cri_levels=[80, 90, 95, 99], quantiles=[0.2, 0.75],
plot_path=".", plot_filename="binom_result"
)Bayesian Conjugate Model – Poisson-Gamma
Show code
# =============================================================
# Bayesian Conjugate Model – Poisson-Gamma
# Closed-form analytical solution (standalone script)
#
# Likelihood : X_i ~ Poisson(lambda)
# Prior : lambda ~ Gamma(alpha, rate = beta)
# Posterior : lambda | x ~ Gamma(alpha + sum(x), rate = beta + n)
#
# Input : scalar / list / vector OR DataFrame + var + group_var
# Output : clean console tables, Plotly HTML saved to disk
# Plots : Plotly only — saved as HTML, never auto-rendered
# =============================================================
import os
import math
import numpy as np
import pandas as pd
from scipy.stats import gamma as sp_gamma
import plotly.graph_objects as go
from plotly.subplots import make_subplots
# =============================================================
# INTERNAL UTILITIES
# =============================================================
def _print_table(title: str, df: pd.DataFrame) -> None:
“””Plain text table to console — no color, no clutter.”””
col_names = list(df.columns)
col_widths = [
max(len(nm), df[nm].astype(str).map(len).max())
for nm in col_names
]
col_widths = [max(w, 10) for w in col_widths]
total_width = sum(col_widths) + 3 * (len(col_widths) – 1) + 2
divider = “=” * total_width
thin_div = “-” * total_width
header = ” “.join(nm.ljust(w) for nm, w in zip(col_names, col_widths))
lines = []
lines.append(divider)
lines.append(f” {title}”)
lines.append(divider)
lines.append(f” {header}”)
lines.append(f” {thin_div}”)
for _, row in df.iterrows():
row_str = ” “.join(str(val).ljust(w) for val, w in zip(row, col_widths))
lines.append(f” {row_str}”)
lines.append(divider)
lines.append(“”)
print(“\n”.join(lines))
def _fmt_mode(v) -> str:
return “undefined” if v is None else f”{v:.6f}”
def _expand_cri(cri_levels, quantiles):
“””Expand cri_levels to symmetric tail quantile pairs; merge with raw quantiles.”””
q_all = sorted(set(quantiles)) if quantiles else []
if cri_levels:
for lev in sorted(set(cri_levels)):
lo = (1 – lev / 100) / 2
hi = 1 – lo
q_all = sorted(set(q_all + [lo, hi]))
return q_all
def _build_cri_table(cri_levels, all_q, q_vals, cri_type):
“””Build CrI table from cri_levels only — symmetric pairs.”””
if not cri_levels:
return None
rows = []
for lev in sorted(set(cri_levels)):
lo = (1 – lev / 100) / 2
hi = 1 – lo
idx_lo = next((i for i, q in enumerate(all_q) if abs(q – lo) < 1e-10), None)
idx_hi = next((i for i, q in enumerate(all_q) if abs(q - hi) < 1e-10), None)
if idx_lo is not None and idx_hi is not None:
rows.append({
"Level": f"{lev}% {cri_type}",
"Lower": f"{q_vals[idx_lo]:.6f}",
"Upper": f"{q_vals[idx_hi]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _build_quantile_table(quantiles, all_q, q_vals):
"""Build individual quantile table from quantiles arg only."""
if not quantiles:
return None
rows = []
for q in sorted(set(quantiles)):
idx = next((i for i, v in enumerate(all_q) if abs(v - q) < 1e-10), None)
if idx is not None:
rows.append({
"Quantile": f"Q({q:.4f})",
"Value": f"{q_vals[idx]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _validate_df_var(data: pd.DataFrame, var: str) -> None:
“””Validate DataFrame and column for the Poisson-Gamma model.”””
if not isinstance(data, pd.DataFrame):
raise TypeError(“`data` must be a pandas DataFrame”)
if not isinstance(var, str):
raise TypeError(“`var` must be a single column name (string)”)
if var not in data.columns:
raise ValueError(f”Column ‘{var}’ not found in data”)
col = data[var].dropna()
if col.empty:
raise ValueError(f”Column ‘{var}’ is entirely NA”)
if not pd.api.types.is_numeric_dtype(col):
raise TypeError(f”Column ‘{var}’ must be numeric for Poisson-Gamma model”)
if (col < 0).any():
raise ValueError(
f"Column '{var}' contains negative values — Poisson model requires non-negative counts"
)
if not (col == col.apply(math.floor)).all():
raise ValueError(
f"Column '{var}' contains non-integer values — Poisson model requires integer counts"
)
def _validate_group(data: pd.DataFrame, group_var) -> None:
if group_var is not None and group_var not in data.columns:
raise ValueError(f”group_var column ‘{group_var}’ not found in data”)
def _split_entries(data: pd.DataFrame, var: str, group_var):
“””Split data into labelled entries by group.”””
if group_var is not None:
entries = []
for g in data[group_var].unique():
sub = data[data[group_var] == g]
entries.append({
“label”: f”{group_var} = {g}”,
“col”: sub[var].dropna().values
})
return entries
else:
return [{“label”: “All”, “col”: data[var].dropna().values}]
def _resolve_plot_path(plot_path, plot_filename, default_name: str) -> str:
“””Resolve final HTML save path.”””
directory = plot_path if plot_path is not None else os.getcwd()
filename = (plot_filename if plot_filename is not None else default_name)
if not filename.endswith(“.html”):
filename += “.html”
return os.path.join(directory, filename)
def _save_plot(fig, path: str) -> None:
fig.write_html(path, auto_open=False)
print(f” [Plot saved] {path}\n”)
# =============================================================
# POISSON – GAMMA
# =============================================================
def conjugate_poisson_gamma(
data = None,
var = None,
group_var = None,
x = None,
alpha = None,
beta = None,
cri_levels = (90, 95),
quantiles = None,
cri_type = “CrI”,
plot_path = None,
plot_filename = None,
):
cri_levels = list(cri_levels) if cri_levels else []
quantiles = list(quantiles) if quantiles else []
all_q = _expand_cri(cri_levels, quantiles)
if data is not None:
_validate_df_var(data, var)
_validate_group(data, group_var)
entries_raw = _split_entries(data, var, group_var)
entries = [{“label”: e[“label”], “x”: e[“col”]} for e in entries_raw]
else:
if x is None:
raise ValueError(“Provide `data` + `var` OR `x`”)
entries = [{“label”: “All”, “x”: np.atleast_1d(np.asarray(x, dtype=float))}]
results = []
for e in entries:
xv = np.asarray(e[“x”], dtype=float)
nv = len(xv)
sum_x = float(np.sum(xv))
a_post = alpha + sum_x
b_post = beta + nv
pr_mean = alpha / beta
pr_var = alpha / beta**2
pr_mode = (alpha – 1) / beta if alpha >= 1 else None
pr_median = sp_gamma.ppf(0.5, a=alpha, scale=1/beta)
po_mean = a_post / b_post
po_var = a_post / b_post**2
po_mode = (a_post – 1) / b_post if a_post >= 1 else None
po_median = sp_gamma.ppf(0.5, a=a_post, scale=1/b_post)
all_q_int = sorted(set(all_q + [0.5]))
q_vals = [sp_gamma.ppf(q, a=a_post, scale=1/b_post) for q in all_q_int]
_print_table(
f”Poisson-Gamma | {e[‘label’]} | Data Summary”,
pd.DataFrame([{
“n”: nv,
“sum(x)”: int(sum_x),
“mean(x)”: f”{float(np.mean(xv)):.6f}”,
}])
)
_print_table(
f”Poisson-Gamma | {e[‘label’]} | Parameters & Moments”,
pd.DataFrame({
“Quantity”: [“alpha (shape)”, “beta (rate)”, “Mean”, “Median”, “Mode”, “Variance”, “SD”],
“Prior”: [f”{alpha:.6f}”, f”{beta:.6f}”, f”{pr_mean:.6f}”,
f”{pr_median:.6f}”, _fmt_mode(pr_mode),
f”{pr_var:.6f}”, f”{math.sqrt(pr_var):.6f}”],
“Posterior”: [f”{a_post:.6f}”, f”{b_post:.6f}”, f”{po_mean:.6f}”,
f”{po_median:.6f}”, _fmt_mode(po_mode),
f”{po_var:.6f}”, f”{math.sqrt(po_var):.6f}”],
})
)
cri_df = _build_cri_table(cri_levels, all_q_int, q_vals, cri_type)
if cri_df is not None:
_print_table(f”Poisson-Gamma | {e[‘label’]} | Credible Intervals ({cri_type})”, cri_df)
q_df = _build_quantile_table(quantiles, all_q_int, q_vals)
if q_df is not None:
_print_table(f”Poisson-Gamma | {e[‘label’]} | Posterior Quantiles”, q_df)
# — Plot —
pr_lo = max(1e-6, sp_gamma.ppf(0.001, a=alpha, scale=1/beta))
pr_hi = sp_gamma.ppf(0.999, a=alpha, scale=1/beta)
po_lo = max(1e-6, sp_gamma.ppf(0.001, a=a_post, scale=1/b_post))
po_hi = sp_gamma.ppf(0.999, a=a_post, scale=1/b_post)
lam_pr = np.linspace(pr_lo, pr_hi, 1000)
lam_po = np.linspace(po_lo, po_hi, 1000)
fig = make_subplots(rows=1, cols=2, subplot_titles=[“Prior”, “Posterior”])
fig.add_trace(go.Scatter(
x=lam_pr, y=sp_gamma.pdf(lam_pr, a=alpha, scale=1/beta),
mode=”lines”, line=dict(color=”#2C7BB6″, width=1.5), showlegend=False
), row=1, col=1)
fig.add_trace(go.Scatter(
x=lam_po, y=sp_gamma.pdf(lam_po, a=a_post, scale=1/b_post),
mode=”lines”, line=dict(color=”#2C7BB6″, width=1.5), showlegend=False
), row=1, col=2)
fig.update_xaxes(title_text=”λ”, row=1, col=1)
fig.update_xaxes(title_text=”λ”, row=1, col=2)
fig.update_yaxes(title_text=”Density”, row=1, col=1)
fig.update_layout(
title_text=(
f”Poisson-Gamma [{e[‘label’]}]
”
f”Prior: Gamma({alpha:.2f}, rate={beta:.2f}) | ”
f”Posterior: Gamma({a_post:.2f}, rate={b_post:.2f})”
),
template=”plotly_white”,
)
default_name = f”poisson_gamma_{e[‘label’].replace(‘ ‘, ‘_’)}”
fpath = _resolve_plot_path(plot_path, plot_filename, default_name)
_save_plot(fig, fpath)
results.append({
“label”: e[“label”],
“data”: {“n”: nv, “sum_x”: sum_x},
“prior”: {“alpha”: alpha, “beta”: beta,
“mean”: pr_mean, “median”: pr_median, “mode”: pr_mode,
“var”: pr_var, “sd”: math.sqrt(pr_var)},
“posterior”: {“alpha”: a_post, “beta”: b_post,
“mean”: po_mean, “median”: po_median, “mode”: po_mode,
“var”: po_var, “sd”: math.sqrt(po_var),
“quantiles”: dict(zip(all_q_int, q_vals))},
“plot”: fig,
})
return results[0] if len(results) == 1 else results
# =============================================================
# USAGE EXAMPLES
# =============================================================
if __name__ == “__main__”:
# scalar / list
# res = conjugate_poisson_gamma(x=[3,5,2,4,6], alpha=2, beta=1,
# cri_levels=(80, 90, 95, 99))
# df + var + group
# df = pd.DataFrame({“counts”: [3,5,2,4,6,1,3,2,5,4],
# “site”: [“S1”]*5 + [“S2″]*5})
# res = conjugate_poisson_gamma(data=df, var=”counts”, group_var=”site”,
# alpha=2, beta=1, cri_levels=(90, 95))
res = conjugate_poisson_gamma(
x=[3, 5, 2, 4, 6], alpha=2, beta=1,
cri_levels=(80, 90, 95, 99),
plot_path=”.”, plot_filename=”poisson_result”
)
# =============================================================
# Bayesian Conjugate Model - Poisson-Gamma
# Closed-form analytical solution (standalone script)
#
# Likelihood : X_i ~ Poisson(lambda)
# Prior : lambda ~ Gamma(alpha, rate = beta)
# Posterior : lambda | x ~ Gamma(alpha + sum(x), rate = beta + n)
#
# Input : scalar / list / vector OR DataFrame + var + group_var
# Output : clean console tables, Plotly HTML saved to disk
# Plots : Plotly only — saved as HTML, never auto-rendered
# =============================================================
import os
import math
import numpy as np
import pandas as pd
from scipy.stats import gamma as sp_gamma
import plotly.graph_objects as go
from plotly.subplots import make_subplots
# =============================================================
# INTERNAL UTILITIES
# =============================================================
def _print_table(title: str, df: pd.DataFrame) -> None:
"""Plain text table to console — no color, no clutter."""
col_names = list(df.columns)
col_widths = [
max(len(nm), df[nm].astype(str).map(len).max())
for nm in col_names
]
col_widths = [max(w, 10) for w in col_widths]
total_width = sum(col_widths) + 3 * (len(col_widths) - 1) + 2
divider = "=" * total_width
thin_div = "-" * total_width
header = " ".join(nm.ljust(w) for nm, w in zip(col_names, col_widths))
lines = []
lines.append(divider)
lines.append(f" {title}")
lines.append(divider)
lines.append(f" {header}")
lines.append(f" {thin_div}")
for _, row in df.iterrows():
row_str = " ".join(str(val).ljust(w) for val, w in zip(row, col_widths))
lines.append(f" {row_str}")
lines.append(divider)
lines.append("")
print("\n".join(lines))
def _fmt_mode(v) -> str:
return "undefined" if v is None else f"{v:.6f}"
def _expand_cri(cri_levels, quantiles):
"""Expand cri_levels to symmetric tail quantile pairs; merge with raw quantiles."""
q_all = sorted(set(quantiles)) if quantiles else []
if cri_levels:
for lev in sorted(set(cri_levels)):
lo = (1 - lev / 100) / 2
hi = 1 - lo
q_all = sorted(set(q_all + [lo, hi]))
return q_all
def _build_cri_table(cri_levels, all_q, q_vals, cri_type):
"""Build CrI table from cri_levels only — symmetric pairs."""
if not cri_levels:
return None
rows = []
for lev in sorted(set(cri_levels)):
lo = (1 - lev / 100) / 2
hi = 1 - lo
idx_lo = next((i for i, q in enumerate(all_q) if abs(q - lo) < 1e-10), None)
idx_hi = next((i for i, q in enumerate(all_q) if abs(q - hi) < 1e-10), None)
if idx_lo is not None and idx_hi is not None:
rows.append({
"Level": f"{lev}% {cri_type}",
"Lower": f"{q_vals[idx_lo]:.6f}",
"Upper": f"{q_vals[idx_hi]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _build_quantile_table(quantiles, all_q, q_vals):
"""Build individual quantile table from quantiles arg only."""
if not quantiles:
return None
rows = []
for q in sorted(set(quantiles)):
idx = next((i for i, v in enumerate(all_q) if abs(v - q) < 1e-10), None)
if idx is not None:
rows.append({
"Quantile": f"Q({q:.4f})",
"Value": f"{q_vals[idx]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _validate_df_var(data: pd.DataFrame, var: str) -> None:
"""Validate DataFrame and column for the Poisson-Gamma model."""
if not isinstance(data, pd.DataFrame):
raise TypeError("`data` must be a pandas DataFrame")
if not isinstance(var, str):
raise TypeError("`var` must be a single column name (string)")
if var not in data.columns:
raise ValueError(f"Column '{var}' not found in data")
col = data[var].dropna()
if col.empty:
raise ValueError(f"Column '{var}' is entirely NA")
if not pd.api.types.is_numeric_dtype(col):
raise TypeError(f"Column '{var}' must be numeric for Poisson-Gamma model")
if (col < 0).any():
raise ValueError(
f"Column '{var}' contains negative values — Poisson model requires non-negative counts"
)
if not (col == col.apply(math.floor)).all():
raise ValueError(
f"Column '{var}' contains non-integer values — Poisson model requires integer counts"
)
def _validate_group(data: pd.DataFrame, group_var) -> None:
if group_var is not None and group_var not in data.columns:
raise ValueError(f"group_var column '{group_var}' not found in data")
def _split_entries(data: pd.DataFrame, var: str, group_var):
"""Split data into labelled entries by group."""
if group_var is not None:
entries = []
for g in data[group_var].unique():
sub = data[data[group_var] == g]
entries.append({
"label": f"{group_var} = {g}",
"col": sub[var].dropna().values
})
return entries
else:
return [{"label": "All", "col": data[var].dropna().values}]
def _resolve_plot_path(plot_path, plot_filename, default_name: str) -> str:
"""Resolve final HTML save path."""
directory = plot_path if plot_path is not None else os.getcwd()
filename = (plot_filename if plot_filename is not None else default_name)
if not filename.endswith(".html"):
filename += ".html"
return os.path.join(directory, filename)
def _save_plot(fig, path: str) -> None:
fig.write_html(path, auto_open=False)
print(f" [Plot saved] {path}\n")
# =============================================================
# POISSON - GAMMA
# =============================================================
def conjugate_poisson_gamma(
data = None,
var = None,
group_var = None,
x = None,
alpha = None,
beta = None,
cri_levels = (90, 95),
quantiles = None,
cri_type = "CrI",
plot_path = None,
plot_filename = None,
):
cri_levels = list(cri_levels) if cri_levels else []
quantiles = list(quantiles) if quantiles else []
all_q = _expand_cri(cri_levels, quantiles)
if data is not None:
_validate_df_var(data, var)
_validate_group(data, group_var)
entries_raw = _split_entries(data, var, group_var)
entries = [{"label": e["label"], "x": e["col"]} for e in entries_raw]
else:
if x is None:
raise ValueError("Provide `data` + `var` OR `x`")
entries = [{"label": "All", "x": np.atleast_1d(np.asarray(x, dtype=float))}]
results = []
for e in entries:
xv = np.asarray(e["x"], dtype=float)
nv = len(xv)
sum_x = float(np.sum(xv))
a_post = alpha + sum_x
b_post = beta + nv
pr_mean = alpha / beta
pr_var = alpha / beta**2
pr_mode = (alpha - 1) / beta if alpha >= 1 else None
pr_median = sp_gamma.ppf(0.5, a=alpha, scale=1/beta)
po_mean = a_post / b_post
po_var = a_post / b_post**2
po_mode = (a_post - 1) / b_post if a_post >= 1 else None
po_median = sp_gamma.ppf(0.5, a=a_post, scale=1/b_post)
all_q_int = sorted(set(all_q + [0.5]))
q_vals = [sp_gamma.ppf(q, a=a_post, scale=1/b_post) for q in all_q_int]
_print_table(
f"Poisson-Gamma | {e['label']} | Data Summary",
pd.DataFrame([{
"n": nv,
"sum(x)": int(sum_x),
"mean(x)": f"{float(np.mean(xv)):.6f}",
}])
)
_print_table(
f"Poisson-Gamma | {e['label']} | Parameters & Moments",
pd.DataFrame({
"Quantity": ["alpha (shape)", "beta (rate)", "Mean", "Median", "Mode", "Variance", "SD"],
"Prior": [f"{alpha:.6f}", f"{beta:.6f}", f"{pr_mean:.6f}",
f"{pr_median:.6f}", _fmt_mode(pr_mode),
f"{pr_var:.6f}", f"{math.sqrt(pr_var):.6f}"],
"Posterior": [f"{a_post:.6f}", f"{b_post:.6f}", f"{po_mean:.6f}",
f"{po_median:.6f}", _fmt_mode(po_mode),
f"{po_var:.6f}", f"{math.sqrt(po_var):.6f}"],
})
)
cri_df = _build_cri_table(cri_levels, all_q_int, q_vals, cri_type)
if cri_df is not None:
_print_table(f"Poisson-Gamma | {e['label']} | Credible Intervals ({cri_type})", cri_df)
q_df = _build_quantile_table(quantiles, all_q_int, q_vals)
if q_df is not None:
_print_table(f"Poisson-Gamma | {e['label']} | Posterior Quantiles", q_df)
# --- Plot ---
pr_lo = max(1e-6, sp_gamma.ppf(0.001, a=alpha, scale=1/beta))
pr_hi = sp_gamma.ppf(0.999, a=alpha, scale=1/beta)
po_lo = max(1e-6, sp_gamma.ppf(0.001, a=a_post, scale=1/b_post))
po_hi = sp_gamma.ppf(0.999, a=a_post, scale=1/b_post)
lam_pr = np.linspace(pr_lo, pr_hi, 1000)
lam_po = np.linspace(po_lo, po_hi, 1000)
fig = make_subplots(rows=1, cols=2, subplot_titles=["Prior", "Posterior"])
fig.add_trace(go.Scatter(
x=lam_pr, y=sp_gamma.pdf(lam_pr, a=alpha, scale=1/beta),
mode="lines", line=dict(color="#2C7BB6", width=1.5), showlegend=False
), row=1, col=1)
fig.add_trace(go.Scatter(
x=lam_po, y=sp_gamma.pdf(lam_po, a=a_post, scale=1/b_post),
mode="lines", line=dict(color="#2C7BB6", width=1.5), showlegend=False
), row=1, col=2)
fig.update_xaxes(title_text="λ", row=1, col=1)
fig.update_xaxes(title_text="λ", row=1, col=2)
fig.update_yaxes(title_text="Density", row=1, col=1)
fig.update_layout(
title_text=(
f"Poisson-Gamma [{e['label']}]<br>"
f"<sup>Prior: Gamma({alpha:.2f}, rate={beta:.2f}) | "
f"Posterior: Gamma({a_post:.2f}, rate={b_post:.2f})</sup>"
),
template="plotly_white",
)
default_name = f"poisson_gamma_{e['label'].replace(' ', '_')}"
fpath = _resolve_plot_path(plot_path, plot_filename, default_name)
_save_plot(fig, fpath)
results.append({
"label": e["label"],
"data": {"n": nv, "sum_x": sum_x},
"prior": {"alpha": alpha, "beta": beta,
"mean": pr_mean, "median": pr_median, "mode": pr_mode,
"var": pr_var, "sd": math.sqrt(pr_var)},
"posterior": {"alpha": a_post, "beta": b_post,
"mean": po_mean, "median": po_median, "mode": po_mode,
"var": po_var, "sd": math.sqrt(po_var),
"quantiles": dict(zip(all_q_int, q_vals))},
"plot": fig,
})
return results[0] if len(results) == 1 else results
# =============================================================
# USAGE EXAMPLES
# =============================================================
if __name__ == "__main__":
# scalar / list
# res = conjugate_poisson_gamma(x=[3,5,2,4,6], alpha=2, beta=1,
# cri_levels=(80, 90, 95, 99))
# df + var + group
# df = pd.DataFrame({"counts": [3,5,2,4,6,1,3,2,5,4],
# "site": ["S1"]*5 + ["S2"]*5})
# res = conjugate_poisson_gamma(data=df, var="counts", group_var="site",
# alpha=2, beta=1, cri_levels=(90, 95))
res = conjugate_poisson_gamma(
x=[3, 5, 2, 4, 6], alpha=2, beta=1,
cri_levels=(80, 90, 95, 99),
plot_path=".", plot_filename="poisson_result"
)Bayesian Conjugate Model – Normal-Normal (Mean, variance known)
Show code
# =============================================================
# Bayesian Conjugate Model – Normal-Normal (Mean, variance known)
# Closed-form analytical solution (standalone script)
#
# Likelihood : X_i ~ N(mu, sigma2) sigma2 known
# Prior : mu ~ N(mu0, tau2)
# Posterior : mu | x ~ N(mu_post, tau2_post)
# tau2_post = 1 / (1/tau2 + n/sigma2)
# mu_post = tau2_post * (mu0/tau2 + n*xbar/sigma2)
#
# Input : scalar / list / vector OR DataFrame + var + group_var
# Output : clean console tables, Plotly HTML saved to disk
# Plots : Plotly only — saved as HTML, never auto-rendered
# =============================================================
import os
import math
import numpy as np
import pandas as pd
from scipy.stats import norm as sp_norm
import plotly.graph_objects as go
from plotly.subplots import make_subplots
# =============================================================
# INTERNAL UTILITIES
# =============================================================
def _print_table(title: str, df: pd.DataFrame) -> None:
“””Plain text table to console — no color, no clutter.”””
col_names = list(df.columns)
col_widths = [
max(len(nm), df[nm].astype(str).map(len).max())
for nm in col_names
]
col_widths = [max(w, 10) for w in col_widths]
total_width = sum(col_widths) + 3 * (len(col_widths) – 1) + 2
divider = “=” * total_width
thin_div = “-” * total_width
header = ” “.join(nm.ljust(w) for nm, w in zip(col_names, col_widths))
lines = []
lines.append(divider)
lines.append(f” {title}”)
lines.append(divider)
lines.append(f” {header}”)
lines.append(f” {thin_div}”)
for _, row in df.iterrows():
row_str = ” “.join(str(val).ljust(w) for val, w in zip(row, col_widths))
lines.append(f” {row_str}”)
lines.append(divider)
lines.append(“”)
print(“\n”.join(lines))
def _expand_cri(cri_levels, quantiles):
“””Expand cri_levels to symmetric tail quantile pairs; merge with raw quantiles.”””
q_all = sorted(set(quantiles)) if quantiles else []
if cri_levels:
for lev in sorted(set(cri_levels)):
lo = (1 – lev / 100) / 2
hi = 1 – lo
q_all = sorted(set(q_all + [lo, hi]))
return q_all
def _build_cri_table(cri_levels, all_q, q_vals, cri_type):
“””Build CrI table from cri_levels only — symmetric pairs.”””
if not cri_levels:
return None
rows = []
for lev in sorted(set(cri_levels)):
lo = (1 – lev / 100) / 2
hi = 1 – lo
idx_lo = next((i for i, q in enumerate(all_q) if abs(q – lo) < 1e-10), None)
idx_hi = next((i for i, q in enumerate(all_q) if abs(q - hi) < 1e-10), None)
if idx_lo is not None and idx_hi is not None:
rows.append({
"Level": f"{lev}% {cri_type}",
"Lower": f"{q_vals[idx_lo]:.6f}",
"Upper": f"{q_vals[idx_hi]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _build_quantile_table(quantiles, all_q, q_vals):
"""Build individual quantile table from quantiles arg only."""
if not quantiles:
return None
rows = []
for q in sorted(set(quantiles)):
idx = next((i for i, v in enumerate(all_q) if abs(v - q) < 1e-10), None)
if idx is not None:
rows.append({
"Quantile": f"Q({q:.4f})",
"Value": f"{q_vals[idx]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _validate_df_var(data: pd.DataFrame, var: str) -> None:
“””Validate DataFrame and column for the Normal-Normal model.”””
if not isinstance(data, pd.DataFrame):
raise TypeError(“`data` must be a pandas DataFrame”)
if not isinstance(var, str):
raise TypeError(“`var` must be a single column name (string)”)
if var not in data.columns:
raise ValueError(f”Column ‘{var}’ not found in data”)
col = data[var].dropna()
if col.empty:
raise ValueError(f”Column ‘{var}’ is entirely NA”)
if not pd.api.types.is_numeric_dtype(col):
raise TypeError(f”Column ‘{var}’ must be numeric for Normal conjugate models”)
def _validate_group(data: pd.DataFrame, group_var) -> None:
if group_var is not None and group_var not in data.columns:
raise ValueError(f”group_var column ‘{group_var}’ not found in data”)
def _split_entries(data: pd.DataFrame, var: str, group_var):
“””Split data into labelled entries by group.”””
if group_var is not None:
entries = []
for g in data[group_var].unique():
sub = data[data[group_var] == g]
entries.append({
“label”: f”{group_var} = {g}”,
“col”: sub[var].dropna().values
})
return entries
else:
return [{“label”: “All”, “col”: data[var].dropna().values}]
def _resolve_plot_path(plot_path, plot_filename, default_name: str) -> str:
“””Resolve final HTML save path.”””
directory = plot_path if plot_path is not None else os.getcwd()
filename = (plot_filename if plot_filename is not None else default_name)
if not filename.endswith(“.html”):
filename += “.html”
return os.path.join(directory, filename)
def _save_plot(fig, path: str) -> None:
fig.write_html(path, auto_open=False)
print(f” [Plot saved] {path}\n”)
# =============================================================
# NORMAL – NORMAL (estimating MEAN | variance known)
# =============================================================
def conjugate_normal_normal(
data = None,
var = None,
group_var = None,
x = None,
mu0 = None,
tau2 = None,
sigma2 = None,
cri_levels = (90, 95),
quantiles = None,
cri_type = “CrI”,
plot_path = None,
plot_filename = None,
):
cri_levels = list(cri_levels) if cri_levels else []
quantiles = list(quantiles) if quantiles else []
all_q = _expand_cri(cri_levels, quantiles)
if data is not None:
_validate_df_var(data, var)
_validate_group(data, group_var)
entries_raw = _split_entries(data, var, group_var)
entries = [{“label”: e[“label”], “x”: e[“col”]} for e in entries_raw]
else:
if x is None:
raise ValueError(“Provide `data` + `var` OR `x`”)
entries = [{“label”: “All”, “x”: np.atleast_1d(np.asarray(x, dtype=float))}]
results = []
for e in entries:
xv = np.asarray(e[“x”], dtype=float)
nv = len(xv)
xbar = float(np.mean(xv))
tau2_post = 1 / (1/tau2 + nv/sigma2)
mu_post = tau2_post * (mu0/tau2 + nv*xbar/sigma2)
# Normal is symmetric: mean = median = mode
pr_median = mu0
pr_mode = mu0
po_median = mu_post
po_mode = mu_post
all_q_int = sorted(set(all_q + [0.5]))
q_vals = [sp_norm.ppf(q, loc=mu_post, scale=math.sqrt(tau2_post)) for q in all_q_int]
_print_table(
f”Normal-Normal | {e[‘label’]} | Data Summary”,
pd.DataFrame([{
“n”: nv,
“xbar”: f”{xbar:.6f}”,
“sigma2 (known)”: f”{sigma2:.6f}”,
}])
)
_print_table(
f”Normal-Normal | {e[‘label’]} | Parameters & Moments”,
pd.DataFrame({
“Quantity”: [“mu (mean=median=mode)”, “tau2 (variance)”, “SD”],
“Prior”: [f”{mu0:.6f}”, f”{tau2:.6f}”, f”{math.sqrt(tau2):.6f}”],
“Posterior”: [f”{mu_post:.6f}”, f”{tau2_post:.6f}”, f”{math.sqrt(tau2_post):.6f}”],
})
)
cri_df = _build_cri_table(cri_levels, all_q_int, q_vals, cri_type)
if cri_df is not None:
_print_table(f”Normal-Normal | {e[‘label’]} | Credible Intervals ({cri_type})”, cri_df)
q_df = _build_quantile_table(quantiles, all_q_int, q_vals)
if q_df is not None:
_print_table(f”Normal-Normal | {e[‘label’]} | Posterior Quantiles”, q_df)
# — Plot —
pr_lo = sp_norm.ppf(0.001, loc=mu0, scale=math.sqrt(tau2))
pr_hi = sp_norm.ppf(0.999, loc=mu0, scale=math.sqrt(tau2))
po_lo = sp_norm.ppf(0.001, loc=mu_post, scale=math.sqrt(tau2_post))
po_hi = sp_norm.ppf(0.999, loc=mu_post, scale=math.sqrt(tau2_post))
mu_pr = np.linspace(pr_lo, pr_hi, 1000)
mu_po = np.linspace(po_lo, po_hi, 1000)
fig = make_subplots(rows=1, cols=2, subplot_titles=[“Prior”, “Posterior”])
fig.add_trace(go.Scatter(
x=mu_pr, y=sp_norm.pdf(mu_pr, loc=mu0, scale=math.sqrt(tau2)),
mode=”lines”, line=dict(color=”#2C7BB6″, width=1.5), showlegend=False
), row=1, col=1)
fig.add_trace(go.Scatter(
x=mu_po, y=sp_norm.pdf(mu_po, loc=mu_post, scale=math.sqrt(tau2_post)),
mode=”lines”, line=dict(color=”#2C7BB6″, width=1.5), showlegend=False
), row=1, col=2)
fig.update_xaxes(title_text=”μ”, row=1, col=1)
fig.update_xaxes(title_text=”μ”, row=1, col=2)
fig.update_yaxes(title_text=”Density”, row=1, col=1)
fig.update_layout(
title_text=(
f”Normal-Normal (Mean) [{e[‘label’]}]
”
f”Prior: N({mu0:.4f}, {tau2:.4f}) | ”
f”Posterior: N({mu_post:.4f}, {tau2_post:.4f})”
),
template=”plotly_white”,
)
default_name = f”normal_normal_{e[‘label’].replace(‘ ‘, ‘_’)}”
fpath = _resolve_plot_path(plot_path, plot_filename, default_name)
_save_plot(fig, fpath)
results.append({
“label”: e[“label”],
“data”: {“n”: nv, “xbar”: xbar, “sigma2”: sigma2},
“prior”: {“mu0”: mu0, “tau2”: tau2,
“mean”: mu0, “median”: pr_median, “mode”: pr_mode,
“sd”: math.sqrt(tau2)},
“posterior”: {“mu”: mu_post, “tau2”: tau2_post,
“mean”: mu_post, “median”: po_median, “mode”: po_mode,
“sd”: math.sqrt(tau2_post),
“quantiles”: dict(zip(all_q_int, q_vals))},
“plot”: fig,
})
return results[0] if len(results) == 1 else results
# =============================================================
# USAGE EXAMPLES
# =============================================================
if __name__ == “__main__”:
# scalar / list
# res = conjugate_normal_normal(x=[2.1,1.9,2.3,2.0], mu0=2, tau2=1,
# sigma2=0.5, cri_levels=(90, 95))
# df + var + group
# df = pd.DataFrame({“measure”: [2.1,1.9,2.3,2.0,2.5,2.2,1.8,2.4],
# “batch”: [“B1”]*4 + [“B2″]*4})
# res = conjugate_normal_normal(data=df, var=”measure”, group_var=”batch”,
# mu0=2, tau2=1, sigma2=0.5,
# cri_levels=(90, 95), quantiles=(0.10, 0.75))
res = conjugate_normal_normal(
x=[2.1, 1.9, 2.3, 2.0], mu0=2, tau2=1, sigma2=0.5,
cri_levels=(90, 95),
plot_path=”.”, plot_filename=”normal_normal_result”
)
# =============================================================
# Bayesian Conjugate Model - Normal-Normal (Mean, variance known)
# Closed-form analytical solution (standalone script)
#
# Likelihood : X_i ~ N(mu, sigma2) sigma2 known
# Prior : mu ~ N(mu0, tau2)
# Posterior : mu | x ~ N(mu_post, tau2_post)
# tau2_post = 1 / (1/tau2 + n/sigma2)
# mu_post = tau2_post * (mu0/tau2 + n*xbar/sigma2)
#
# Input : scalar / list / vector OR DataFrame + var + group_var
# Output : clean console tables, Plotly HTML saved to disk
# Plots : Plotly only — saved as HTML, never auto-rendered
# =============================================================
import os
import math
import numpy as np
import pandas as pd
from scipy.stats import norm as sp_norm
import plotly.graph_objects as go
from plotly.subplots import make_subplots
# =============================================================
# INTERNAL UTILITIES
# =============================================================
def _print_table(title: str, df: pd.DataFrame) -> None:
"""Plain text table to console — no color, no clutter."""
col_names = list(df.columns)
col_widths = [
max(len(nm), df[nm].astype(str).map(len).max())
for nm in col_names
]
col_widths = [max(w, 10) for w in col_widths]
total_width = sum(col_widths) + 3 * (len(col_widths) - 1) + 2
divider = "=" * total_width
thin_div = "-" * total_width
header = " ".join(nm.ljust(w) for nm, w in zip(col_names, col_widths))
lines = []
lines.append(divider)
lines.append(f" {title}")
lines.append(divider)
lines.append(f" {header}")
lines.append(f" {thin_div}")
for _, row in df.iterrows():
row_str = " ".join(str(val).ljust(w) for val, w in zip(row, col_widths))
lines.append(f" {row_str}")
lines.append(divider)
lines.append("")
print("\n".join(lines))
def _expand_cri(cri_levels, quantiles):
"""Expand cri_levels to symmetric tail quantile pairs; merge with raw quantiles."""
q_all = sorted(set(quantiles)) if quantiles else []
if cri_levels:
for lev in sorted(set(cri_levels)):
lo = (1 - lev / 100) / 2
hi = 1 - lo
q_all = sorted(set(q_all + [lo, hi]))
return q_all
def _build_cri_table(cri_levels, all_q, q_vals, cri_type):
"""Build CrI table from cri_levels only — symmetric pairs."""
if not cri_levels:
return None
rows = []
for lev in sorted(set(cri_levels)):
lo = (1 - lev / 100) / 2
hi = 1 - lo
idx_lo = next((i for i, q in enumerate(all_q) if abs(q - lo) < 1e-10), None)
idx_hi = next((i for i, q in enumerate(all_q) if abs(q - hi) < 1e-10), None)
if idx_lo is not None and idx_hi is not None:
rows.append({
"Level": f"{lev}% {cri_type}",
"Lower": f"{q_vals[idx_lo]:.6f}",
"Upper": f"{q_vals[idx_hi]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _build_quantile_table(quantiles, all_q, q_vals):
"""Build individual quantile table from quantiles arg only."""
if not quantiles:
return None
rows = []
for q in sorted(set(quantiles)):
idx = next((i for i, v in enumerate(all_q) if abs(v - q) < 1e-10), None)
if idx is not None:
rows.append({
"Quantile": f"Q({q:.4f})",
"Value": f"{q_vals[idx]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _validate_df_var(data: pd.DataFrame, var: str) -> None:
"""Validate DataFrame and column for the Normal-Normal model."""
if not isinstance(data, pd.DataFrame):
raise TypeError("`data` must be a pandas DataFrame")
if not isinstance(var, str):
raise TypeError("`var` must be a single column name (string)")
if var not in data.columns:
raise ValueError(f"Column '{var}' not found in data")
col = data[var].dropna()
if col.empty:
raise ValueError(f"Column '{var}' is entirely NA")
if not pd.api.types.is_numeric_dtype(col):
raise TypeError(f"Column '{var}' must be numeric for Normal conjugate models")
def _validate_group(data: pd.DataFrame, group_var) -> None:
if group_var is not None and group_var not in data.columns:
raise ValueError(f"group_var column '{group_var}' not found in data")
def _split_entries(data: pd.DataFrame, var: str, group_var):
"""Split data into labelled entries by group."""
if group_var is not None:
entries = []
for g in data[group_var].unique():
sub = data[data[group_var] == g]
entries.append({
"label": f"{group_var} = {g}",
"col": sub[var].dropna().values
})
return entries
else:
return [{"label": "All", "col": data[var].dropna().values}]
def _resolve_plot_path(plot_path, plot_filename, default_name: str) -> str:
"""Resolve final HTML save path."""
directory = plot_path if plot_path is not None else os.getcwd()
filename = (plot_filename if plot_filename is not None else default_name)
if not filename.endswith(".html"):
filename += ".html"
return os.path.join(directory, filename)
def _save_plot(fig, path: str) -> None:
fig.write_html(path, auto_open=False)
print(f" [Plot saved] {path}\n")
# =============================================================
# NORMAL - NORMAL (estimating MEAN | variance known)
# =============================================================
def conjugate_normal_normal(
data = None,
var = None,
group_var = None,
x = None,
mu0 = None,
tau2 = None,
sigma2 = None,
cri_levels = (90, 95),
quantiles = None,
cri_type = "CrI",
plot_path = None,
plot_filename = None,
):
cri_levels = list(cri_levels) if cri_levels else []
quantiles = list(quantiles) if quantiles else []
all_q = _expand_cri(cri_levels, quantiles)
if data is not None:
_validate_df_var(data, var)
_validate_group(data, group_var)
entries_raw = _split_entries(data, var, group_var)
entries = [{"label": e["label"], "x": e["col"]} for e in entries_raw]
else:
if x is None:
raise ValueError("Provide `data` + `var` OR `x`")
entries = [{"label": "All", "x": np.atleast_1d(np.asarray(x, dtype=float))}]
results = []
for e in entries:
xv = np.asarray(e["x"], dtype=float)
nv = len(xv)
xbar = float(np.mean(xv))
tau2_post = 1 / (1/tau2 + nv/sigma2)
mu_post = tau2_post * (mu0/tau2 + nv*xbar/sigma2)
# Normal is symmetric: mean = median = mode
pr_median = mu0
pr_mode = mu0
po_median = mu_post
po_mode = mu_post
all_q_int = sorted(set(all_q + [0.5]))
q_vals = [sp_norm.ppf(q, loc=mu_post, scale=math.sqrt(tau2_post)) for q in all_q_int]
_print_table(
f"Normal-Normal | {e['label']} | Data Summary",
pd.DataFrame([{
"n": nv,
"xbar": f"{xbar:.6f}",
"sigma2 (known)": f"{sigma2:.6f}",
}])
)
_print_table(
f"Normal-Normal | {e['label']} | Parameters & Moments",
pd.DataFrame({
"Quantity": ["mu (mean=median=mode)", "tau2 (variance)", "SD"],
"Prior": [f"{mu0:.6f}", f"{tau2:.6f}", f"{math.sqrt(tau2):.6f}"],
"Posterior": [f"{mu_post:.6f}", f"{tau2_post:.6f}", f"{math.sqrt(tau2_post):.6f}"],
})
)
cri_df = _build_cri_table(cri_levels, all_q_int, q_vals, cri_type)
if cri_df is not None:
_print_table(f"Normal-Normal | {e['label']} | Credible Intervals ({cri_type})", cri_df)
q_df = _build_quantile_table(quantiles, all_q_int, q_vals)
if q_df is not None:
_print_table(f"Normal-Normal | {e['label']} | Posterior Quantiles", q_df)
# --- Plot ---
pr_lo = sp_norm.ppf(0.001, loc=mu0, scale=math.sqrt(tau2))
pr_hi = sp_norm.ppf(0.999, loc=mu0, scale=math.sqrt(tau2))
po_lo = sp_norm.ppf(0.001, loc=mu_post, scale=math.sqrt(tau2_post))
po_hi = sp_norm.ppf(0.999, loc=mu_post, scale=math.sqrt(tau2_post))
mu_pr = np.linspace(pr_lo, pr_hi, 1000)
mu_po = np.linspace(po_lo, po_hi, 1000)
fig = make_subplots(rows=1, cols=2, subplot_titles=["Prior", "Posterior"])
fig.add_trace(go.Scatter(
x=mu_pr, y=sp_norm.pdf(mu_pr, loc=mu0, scale=math.sqrt(tau2)),
mode="lines", line=dict(color="#2C7BB6", width=1.5), showlegend=False
), row=1, col=1)
fig.add_trace(go.Scatter(
x=mu_po, y=sp_norm.pdf(mu_po, loc=mu_post, scale=math.sqrt(tau2_post)),
mode="lines", line=dict(color="#2C7BB6", width=1.5), showlegend=False
), row=1, col=2)
fig.update_xaxes(title_text="μ", row=1, col=1)
fig.update_xaxes(title_text="μ", row=1, col=2)
fig.update_yaxes(title_text="Density", row=1, col=1)
fig.update_layout(
title_text=(
f"Normal-Normal (Mean) [{e['label']}]<br>"
f"<sup>Prior: N({mu0:.4f}, {tau2:.4f}) | "
f"Posterior: N({mu_post:.4f}, {tau2_post:.4f})</sup>"
),
template="plotly_white",
)
default_name = f"normal_normal_{e['label'].replace(' ', '_')}"
fpath = _resolve_plot_path(plot_path, plot_filename, default_name)
_save_plot(fig, fpath)
results.append({
"label": e["label"],
"data": {"n": nv, "xbar": xbar, "sigma2": sigma2},
"prior": {"mu0": mu0, "tau2": tau2,
"mean": mu0, "median": pr_median, "mode": pr_mode,
"sd": math.sqrt(tau2)},
"posterior": {"mu": mu_post, "tau2": tau2_post,
"mean": mu_post, "median": po_median, "mode": po_mode,
"sd": math.sqrt(tau2_post),
"quantiles": dict(zip(all_q_int, q_vals))},
"plot": fig,
})
return results[0] if len(results) == 1 else results
# =============================================================
# USAGE EXAMPLES
# =============================================================
if __name__ == "__main__":
# scalar / list
# res = conjugate_normal_normal(x=[2.1,1.9,2.3,2.0], mu0=2, tau2=1,
# sigma2=0.5, cri_levels=(90, 95))
# df + var + group
# df = pd.DataFrame({"measure": [2.1,1.9,2.3,2.0,2.5,2.2,1.8,2.4],
# "batch": ["B1"]*4 + ["B2"]*4})
# res = conjugate_normal_normal(data=df, var="measure", group_var="batch",
# mu0=2, tau2=1, sigma2=0.5,
# cri_levels=(90, 95), quantiles=(0.10, 0.75))
res = conjugate_normal_normal(
x=[2.1, 1.9, 2.3, 2.0], mu0=2, tau2=1, sigma2=0.5,
cri_levels=(90, 95),
plot_path=".", plot_filename="normal_normal_result"
)Bayesian Conjugate Model – Normal-Gamma (Precision, mean known)
Show code
# =============================================================
# Bayesian Conjugate Model – Normal-Gamma (Precision, mean known)
# Closed-form analytical solution (standalone script)
#
# Likelihood : X_i ~ N(mu_known, 1/phi) phi = precision
# Prior : phi ~ Gamma(alpha, rate = beta)
# Posterior : phi | x ~ Gamma(alpha + n/2, rate = beta + SS/2)
# SS = sum((x – mu_known)^2)
#
# Input : scalar / list / vector OR DataFrame + var + group_var
# Output : clean console tables, Plotly HTML saved to disk
# Plots : Plotly only — saved as HTML, never auto-rendered
# =============================================================
import os
import math
import numpy as np
import pandas as pd
from scipy.stats import gamma as sp_gamma
import plotly.graph_objects as go
from plotly.subplots import make_subplots
# =============================================================
# INTERNAL UTILITIES
# =============================================================
def _print_table(title: str, df: pd.DataFrame) -> None:
“””Plain text table to console — no color, no clutter.”””
col_names = list(df.columns)
col_widths = [
max(len(nm), df[nm].astype(str).map(len).max())
for nm in col_names
]
col_widths = [max(w, 10) for w in col_widths]
total_width = sum(col_widths) + 3 * (len(col_widths) – 1) + 2
divider = “=” * total_width
thin_div = “-” * total_width
header = ” “.join(nm.ljust(w) for nm, w in zip(col_names, col_widths))
lines = []
lines.append(divider)
lines.append(f” {title}”)
lines.append(divider)
lines.append(f” {header}”)
lines.append(f” {thin_div}”)
for _, row in df.iterrows():
row_str = ” “.join(str(val).ljust(w) for val, w in zip(row, col_widths))
lines.append(f” {row_str}”)
lines.append(divider)
lines.append(“”)
print(“\n”.join(lines))
def _fmt_mode(v) -> str:
return “undefined” if v is None else f”{v:.6f}”
def _expand_cri(cri_levels, quantiles):
“””Expand cri_levels to symmetric tail quantile pairs; merge with raw quantiles.”””
q_all = sorted(set(quantiles)) if quantiles else []
if cri_levels:
for lev in sorted(set(cri_levels)):
lo = (1 – lev / 100) / 2
hi = 1 – lo
q_all = sorted(set(q_all + [lo, hi]))
return q_all
def _build_cri_table(cri_levels, all_q, q_vals, cri_type):
“””Build CrI table from cri_levels only — symmetric pairs.”””
if not cri_levels:
return None
rows = []
for lev in sorted(set(cri_levels)):
lo = (1 – lev / 100) / 2
hi = 1 – lo
idx_lo = next((i for i, q in enumerate(all_q) if abs(q – lo) < 1e-10), None)
idx_hi = next((i for i, q in enumerate(all_q) if abs(q - hi) < 1e-10), None)
if idx_lo is not None and idx_hi is not None:
rows.append({
"Level": f"{lev}% {cri_type}",
"Lower": f"{q_vals[idx_lo]:.6f}",
"Upper": f"{q_vals[idx_hi]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _build_quantile_table(quantiles, all_q, q_vals):
"""Build individual quantile table from quantiles arg only."""
if not quantiles:
return None
rows = []
for q in sorted(set(quantiles)):
idx = next((i for i, v in enumerate(all_q) if abs(v - q) < 1e-10), None)
if idx is not None:
rows.append({
"Quantile": f"Q({q:.4f})",
"Value": f"{q_vals[idx]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _validate_df_var(data: pd.DataFrame, var: str) -> None:
“””Validate DataFrame and column for the Normal-Gamma model.”””
if not isinstance(data, pd.DataFrame):
raise TypeError(“`data` must be a pandas DataFrame”)
if not isinstance(var, str):
raise TypeError(“`var` must be a single column name (string)”)
if var not in data.columns:
raise ValueError(f”Column ‘{var}’ not found in data”)
col = data[var].dropna()
if col.empty:
raise ValueError(f”Column ‘{var}’ is entirely NA”)
if not pd.api.types.is_numeric_dtype(col):
raise TypeError(f”Column ‘{var}’ must be numeric for Normal conjugate models”)
def _validate_group(data: pd.DataFrame, group_var) -> None:
if group_var is not None and group_var not in data.columns:
raise ValueError(f”group_var column ‘{group_var}’ not found in data”)
def _split_entries(data: pd.DataFrame, var: str, group_var):
“””Split data into labelled entries by group.”””
if group_var is not None:
entries = []
for g in data[group_var].unique():
sub = data[data[group_var] == g]
entries.append({
“label”: f”{group_var} = {g}”,
“col”: sub[var].dropna().values
})
return entries
else:
return [{“label”: “All”, “col”: data[var].dropna().values}]
def _resolve_plot_path(plot_path, plot_filename, default_name: str) -> str:
“””Resolve final HTML save path.”””
directory = plot_path if plot_path is not None else os.getcwd()
filename = (plot_filename if plot_filename is not None else default_name)
if not filename.endswith(“.html”):
filename += “.html”
return os.path.join(directory, filename)
def _save_plot(fig, path: str) -> None:
fig.write_html(path, auto_open=False)
print(f” [Plot saved] {path}\n”)
# =============================================================
# NORMAL – GAMMA (estimating PRECISION | mean known)
# =============================================================
def conjugate_normal_gamma(
data = None,
var = None,
group_var = None,
x = None,
mu_known = None,
alpha = None,
beta = None,
cri_levels = (90, 95),
quantiles = None,
cri_type = “CrI”,
plot_path = None,
plot_filename = None,
):
cri_levels = list(cri_levels) if cri_levels else []
quantiles = list(quantiles) if quantiles else []
all_q = _expand_cri(cri_levels, quantiles)
if data is not None:
_validate_df_var(data, var)
_validate_group(data, group_var)
entries_raw = _split_entries(data, var, group_var)
entries = [{“label”: e[“label”], “x”: e[“col”]} for e in entries_raw]
else:
if x is None:
raise ValueError(“Provide `data` + `var` OR `x`”)
entries = [{“label”: “All”, “x”: np.atleast_1d(np.asarray(x, dtype=float))}]
results = []
for e in entries:
xv = np.asarray(e[“x”], dtype=float)
nv = len(xv)
ss = float(np.sum((xv – mu_known)**2))
a_post = alpha + nv / 2
b_post = beta + ss / 2
pr_mean = alpha / beta
pr_var = alpha / beta**2
pr_mode = (alpha – 1) / beta if alpha >= 1 else None
pr_median = sp_gamma.ppf(0.5, a=alpha, scale=1/beta)
po_mean = a_post / b_post
po_var = a_post / b_post**2
po_mode = (a_post – 1) / b_post if a_post >= 1 else None
po_median = sp_gamma.ppf(0.5, a=a_post, scale=1/b_post)
all_q_int = sorted(set(all_q + [0.5]))
q_vals = [sp_gamma.ppf(q, a=a_post, scale=1/b_post) for q in all_q_int]
_print_table(
f”Normal-Gamma | {e[‘label’]} | Data Summary”,
pd.DataFrame([{
“n”: nv,
“mu_known”: f”{mu_known:.6f}”,
“SS = sum((x-mu)^2)”: f”{ss:.6f}”,
}])
)
_print_table(
f”Normal-Gamma | {e[‘label’]} | Parameters & Moments (Precision)”,
pd.DataFrame({
“Quantity”: [“alpha (shape)”, “beta (rate)”, “Mean”, “Median”, “Mode”, “Variance”, “SD”],
“Prior”: [f”{alpha:.6f}”, f”{beta:.6f}”, f”{pr_mean:.6f}”,
f”{pr_median:.6f}”, _fmt_mode(pr_mode),
f”{pr_var:.6f}”, f”{math.sqrt(pr_var):.6f}”],
“Posterior”: [f”{a_post:.6f}”, f”{b_post:.6f}”, f”{po_mean:.6f}”,
f”{po_median:.6f}”, _fmt_mode(po_mode),
f”{po_var:.6f}”, f”{math.sqrt(po_var):.6f}”],
})
)
cri_df = _build_cri_table(cri_levels, all_q_int, q_vals, cri_type)
if cri_df is not None:
_print_table(
f”Normal-Gamma | {e[‘label’]} | Credible Intervals ({cri_type}) [precision scale]”,
cri_df
)
q_df = _build_quantile_table(quantiles, all_q_int, q_vals)
if q_df is not None:
_print_table(
f”Normal-Gamma | {e[‘label’]} | Posterior Quantiles [precision scale]”,
q_df
)
# — Plot —
pr_lo = max(1e-6, sp_gamma.ppf(0.001, a=alpha, scale=1/beta))
pr_hi = sp_gamma.ppf(0.999, a=alpha, scale=1/beta)
po_lo = max(1e-6, sp_gamma.ppf(0.001, a=a_post, scale=1/b_post))
po_hi = sp_gamma.ppf(0.999, a=a_post, scale=1/b_post)
phi_pr = np.linspace(pr_lo, pr_hi, 1000)
phi_po = np.linspace(po_lo, po_hi, 1000)
fig = make_subplots(rows=1, cols=2, subplot_titles=[“Prior”, “Posterior”])
fig.add_trace(go.Scatter(
x=phi_pr, y=sp_gamma.pdf(phi_pr, a=alpha, scale=1/beta),
mode=”lines”, line=dict(color=”#2C7BB6″, width=1.5), showlegend=False
), row=1, col=1)
fig.add_trace(go.Scatter(
x=phi_po, y=sp_gamma.pdf(phi_po, a=a_post, scale=1/b_post),
mode=”lines”, line=dict(color=”#2C7BB6″, width=1.5), showlegend=False
), row=1, col=2)
fig.update_xaxes(title_text=”φ (precision)”, row=1, col=1)
fig.update_xaxes(title_text=”φ (precision)”, row=1, col=2)
fig.update_yaxes(title_text=”Density”, row=1, col=1)
fig.update_layout(
title_text=(
f”Normal-Gamma (Precision) [{e[‘label’]}]
”
f”Prior: Gamma({alpha:.2f}, rate={beta:.2f}) | ”
f”Posterior: Gamma({a_post:.2f}, rate={b_post:.2f})”
),
template=”plotly_white”,
)
default_name = f”normal_gamma_{e[‘label’].replace(‘ ‘, ‘_’)}”
fpath = _resolve_plot_path(plot_path, plot_filename, default_name)
_save_plot(fig, fpath)
results.append({
“label”: e[“label”],
“data”: {“n”: nv, “mu_known”: mu_known, “ss”: ss},
“prior”: {“alpha”: alpha, “beta”: beta,
“mean”: pr_mean, “median”: pr_median, “mode”: pr_mode,
“var”: pr_var, “sd”: math.sqrt(pr_var)},
“posterior”: {“alpha”: a_post, “beta”: b_post,
“mean”: po_mean, “median”: po_median, “mode”: po_mode,
“var”: po_var, “sd”: math.sqrt(po_var),
“quantiles”: dict(zip(all_q_int, q_vals))},
“plot”: fig,
})
return results[0] if len(results) == 1 else results
# =============================================================
# USAGE EXAMPLES
# =============================================================
if __name__ == “__main__”:
# scalar / list
# res = conjugate_normal_gamma(x=[2.1,1.9,2.3,2.0], mu_known=2,
# alpha=2, beta=1, cri_levels=(90, 95))
# df + var
# df = pd.DataFrame({“measure”: [2.1,1.9,2.3,2.0,2.5,2.2,1.8,2.4]})
# res = conjugate_normal_gamma(data=df, var=”measure”, mu_known=2,
# alpha=2, beta=1, cri_levels=(90, 95))
res = conjugate_normal_gamma(
x=[2.1, 1.9, 2.3, 2.0], mu_known=2, alpha=2, beta=1,
cri_levels=(90, 95),
plot_path=”.”, plot_filename=”normal_gamma_result”
)
# =============================================================
# Bayesian Conjugate Model - Normal-Gamma (Precision, mean known)
# Closed-form analytical solution (standalone script)
#
# Likelihood : X_i ~ N(mu_known, 1/phi) phi = precision
# Prior : phi ~ Gamma(alpha, rate = beta)
# Posterior : phi | x ~ Gamma(alpha + n/2, rate = beta + SS/2)
# SS = sum((x - mu_known)^2)
#
# Input : scalar / list / vector OR DataFrame + var + group_var
# Output : clean console tables, Plotly HTML saved to disk
# Plots : Plotly only — saved as HTML, never auto-rendered
# =============================================================
import os
import math
import numpy as np
import pandas as pd
from scipy.stats import gamma as sp_gamma
import plotly.graph_objects as go
from plotly.subplots import make_subplots
# =============================================================
# INTERNAL UTILITIES
# =============================================================
def _print_table(title: str, df: pd.DataFrame) -> None:
"""Plain text table to console — no color, no clutter."""
col_names = list(df.columns)
col_widths = [
max(len(nm), df[nm].astype(str).map(len).max())
for nm in col_names
]
col_widths = [max(w, 10) for w in col_widths]
total_width = sum(col_widths) + 3 * (len(col_widths) - 1) + 2
divider = "=" * total_width
thin_div = "-" * total_width
header = " ".join(nm.ljust(w) for nm, w in zip(col_names, col_widths))
lines = []
lines.append(divider)
lines.append(f" {title}")
lines.append(divider)
lines.append(f" {header}")
lines.append(f" {thin_div}")
for _, row in df.iterrows():
row_str = " ".join(str(val).ljust(w) for val, w in zip(row, col_widths))
lines.append(f" {row_str}")
lines.append(divider)
lines.append("")
print("\n".join(lines))
def _fmt_mode(v) -> str:
return "undefined" if v is None else f"{v:.6f}"
def _expand_cri(cri_levels, quantiles):
"""Expand cri_levels to symmetric tail quantile pairs; merge with raw quantiles."""
q_all = sorted(set(quantiles)) if quantiles else []
if cri_levels:
for lev in sorted(set(cri_levels)):
lo = (1 - lev / 100) / 2
hi = 1 - lo
q_all = sorted(set(q_all + [lo, hi]))
return q_all
def _build_cri_table(cri_levels, all_q, q_vals, cri_type):
"""Build CrI table from cri_levels only — symmetric pairs."""
if not cri_levels:
return None
rows = []
for lev in sorted(set(cri_levels)):
lo = (1 - lev / 100) / 2
hi = 1 - lo
idx_lo = next((i for i, q in enumerate(all_q) if abs(q - lo) < 1e-10), None)
idx_hi = next((i for i, q in enumerate(all_q) if abs(q - hi) < 1e-10), None)
if idx_lo is not None and idx_hi is not None:
rows.append({
"Level": f"{lev}% {cri_type}",
"Lower": f"{q_vals[idx_lo]:.6f}",
"Upper": f"{q_vals[idx_hi]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _build_quantile_table(quantiles, all_q, q_vals):
"""Build individual quantile table from quantiles arg only."""
if not quantiles:
return None
rows = []
for q in sorted(set(quantiles)):
idx = next((i for i, v in enumerate(all_q) if abs(v - q) < 1e-10), None)
if idx is not None:
rows.append({
"Quantile": f"Q({q:.4f})",
"Value": f"{q_vals[idx]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _validate_df_var(data: pd.DataFrame, var: str) -> None:
"""Validate DataFrame and column for the Normal-Gamma model."""
if not isinstance(data, pd.DataFrame):
raise TypeError("`data` must be a pandas DataFrame")
if not isinstance(var, str):
raise TypeError("`var` must be a single column name (string)")
if var not in data.columns:
raise ValueError(f"Column '{var}' not found in data")
col = data[var].dropna()
if col.empty:
raise ValueError(f"Column '{var}' is entirely NA")
if not pd.api.types.is_numeric_dtype(col):
raise TypeError(f"Column '{var}' must be numeric for Normal conjugate models")
def _validate_group(data: pd.DataFrame, group_var) -> None:
if group_var is not None and group_var not in data.columns:
raise ValueError(f"group_var column '{group_var}' not found in data")
def _split_entries(data: pd.DataFrame, var: str, group_var):
"""Split data into labelled entries by group."""
if group_var is not None:
entries = []
for g in data[group_var].unique():
sub = data[data[group_var] == g]
entries.append({
"label": f"{group_var} = {g}",
"col": sub[var].dropna().values
})
return entries
else:
return [{"label": "All", "col": data[var].dropna().values}]
def _resolve_plot_path(plot_path, plot_filename, default_name: str) -> str:
"""Resolve final HTML save path."""
directory = plot_path if plot_path is not None else os.getcwd()
filename = (plot_filename if plot_filename is not None else default_name)
if not filename.endswith(".html"):
filename += ".html"
return os.path.join(directory, filename)
def _save_plot(fig, path: str) -> None:
fig.write_html(path, auto_open=False)
print(f" [Plot saved] {path}\n")
# =============================================================
# NORMAL - GAMMA (estimating PRECISION | mean known)
# =============================================================
def conjugate_normal_gamma(
data = None,
var = None,
group_var = None,
x = None,
mu_known = None,
alpha = None,
beta = None,
cri_levels = (90, 95),
quantiles = None,
cri_type = "CrI",
plot_path = None,
plot_filename = None,
):
cri_levels = list(cri_levels) if cri_levels else []
quantiles = list(quantiles) if quantiles else []
all_q = _expand_cri(cri_levels, quantiles)
if data is not None:
_validate_df_var(data, var)
_validate_group(data, group_var)
entries_raw = _split_entries(data, var, group_var)
entries = [{"label": e["label"], "x": e["col"]} for e in entries_raw]
else:
if x is None:
raise ValueError("Provide `data` + `var` OR `x`")
entries = [{"label": "All", "x": np.atleast_1d(np.asarray(x, dtype=float))}]
results = []
for e in entries:
xv = np.asarray(e["x"], dtype=float)
nv = len(xv)
ss = float(np.sum((xv - mu_known)**2))
a_post = alpha + nv / 2
b_post = beta + ss / 2
pr_mean = alpha / beta
pr_var = alpha / beta**2
pr_mode = (alpha - 1) / beta if alpha >= 1 else None
pr_median = sp_gamma.ppf(0.5, a=alpha, scale=1/beta)
po_mean = a_post / b_post
po_var = a_post / b_post**2
po_mode = (a_post - 1) / b_post if a_post >= 1 else None
po_median = sp_gamma.ppf(0.5, a=a_post, scale=1/b_post)
all_q_int = sorted(set(all_q + [0.5]))
q_vals = [sp_gamma.ppf(q, a=a_post, scale=1/b_post) for q in all_q_int]
_print_table(
f"Normal-Gamma | {e['label']} | Data Summary",
pd.DataFrame([{
"n": nv,
"mu_known": f"{mu_known:.6f}",
"SS = sum((x-mu)^2)": f"{ss:.6f}",
}])
)
_print_table(
f"Normal-Gamma | {e['label']} | Parameters & Moments (Precision)",
pd.DataFrame({
"Quantity": ["alpha (shape)", "beta (rate)", "Mean", "Median", "Mode", "Variance", "SD"],
"Prior": [f"{alpha:.6f}", f"{beta:.6f}", f"{pr_mean:.6f}",
f"{pr_median:.6f}", _fmt_mode(pr_mode),
f"{pr_var:.6f}", f"{math.sqrt(pr_var):.6f}"],
"Posterior": [f"{a_post:.6f}", f"{b_post:.6f}", f"{po_mean:.6f}",
f"{po_median:.6f}", _fmt_mode(po_mode),
f"{po_var:.6f}", f"{math.sqrt(po_var):.6f}"],
})
)
cri_df = _build_cri_table(cri_levels, all_q_int, q_vals, cri_type)
if cri_df is not None:
_print_table(
f"Normal-Gamma | {e['label']} | Credible Intervals ({cri_type}) [precision scale]",
cri_df
)
q_df = _build_quantile_table(quantiles, all_q_int, q_vals)
if q_df is not None:
_print_table(
f"Normal-Gamma | {e['label']} | Posterior Quantiles [precision scale]",
q_df
)
# --- Plot ---
pr_lo = max(1e-6, sp_gamma.ppf(0.001, a=alpha, scale=1/beta))
pr_hi = sp_gamma.ppf(0.999, a=alpha, scale=1/beta)
po_lo = max(1e-6, sp_gamma.ppf(0.001, a=a_post, scale=1/b_post))
po_hi = sp_gamma.ppf(0.999, a=a_post, scale=1/b_post)
phi_pr = np.linspace(pr_lo, pr_hi, 1000)
phi_po = np.linspace(po_lo, po_hi, 1000)
fig = make_subplots(rows=1, cols=2, subplot_titles=["Prior", "Posterior"])
fig.add_trace(go.Scatter(
x=phi_pr, y=sp_gamma.pdf(phi_pr, a=alpha, scale=1/beta),
mode="lines", line=dict(color="#2C7BB6", width=1.5), showlegend=False
), row=1, col=1)
fig.add_trace(go.Scatter(
x=phi_po, y=sp_gamma.pdf(phi_po, a=a_post, scale=1/b_post),
mode="lines", line=dict(color="#2C7BB6", width=1.5), showlegend=False
), row=1, col=2)
fig.update_xaxes(title_text="φ (precision)", row=1, col=1)
fig.update_xaxes(title_text="φ (precision)", row=1, col=2)
fig.update_yaxes(title_text="Density", row=1, col=1)
fig.update_layout(
title_text=(
f"Normal-Gamma (Precision) [{e['label']}]<br>"
f"<sup>Prior: Gamma({alpha:.2f}, rate={beta:.2f}) | "
f"Posterior: Gamma({a_post:.2f}, rate={b_post:.2f})</sup>"
),
template="plotly_white",
)
default_name = f"normal_gamma_{e['label'].replace(' ', '_')}"
fpath = _resolve_plot_path(plot_path, plot_filename, default_name)
_save_plot(fig, fpath)
results.append({
"label": e["label"],
"data": {"n": nv, "mu_known": mu_known, "ss": ss},
"prior": {"alpha": alpha, "beta": beta,
"mean": pr_mean, "median": pr_median, "mode": pr_mode,
"var": pr_var, "sd": math.sqrt(pr_var)},
"posterior": {"alpha": a_post, "beta": b_post,
"mean": po_mean, "median": po_median, "mode": po_mode,
"var": po_var, "sd": math.sqrt(po_var),
"quantiles": dict(zip(all_q_int, q_vals))},
"plot": fig,
})
return results[0] if len(results) == 1 else results
# =============================================================
# USAGE EXAMPLES
# =============================================================
if __name__ == "__main__":
# scalar / list
# res = conjugate_normal_gamma(x=[2.1,1.9,2.3,2.0], mu_known=2,
# alpha=2, beta=1, cri_levels=(90, 95))
# df + var
# df = pd.DataFrame({"measure": [2.1,1.9,2.3,2.0,2.5,2.2,1.8,2.4]})
# res = conjugate_normal_gamma(data=df, var="measure", mu_known=2,
# alpha=2, beta=1, cri_levels=(90, 95))
res = conjugate_normal_gamma(
x=[2.1, 1.9, 2.3, 2.0], mu_known=2, alpha=2, beta=1,
cri_levels=(90, 95),
plot_path=".", plot_filename="normal_gamma_result"
)Bayesian Conjugate Model – Normal-Inverse-Gamma (Variance, mean known)
Show code
# =============================================================
# Bayesian Conjugate Model – Normal-Inverse-Gamma (Variance, mean known)
# Closed-form analytical solution (standalone script)
#
# Likelihood : X_i ~ N(mu_known, sigma2)
# Prior : sigma2 ~ InvGamma(alpha, scale = beta)
# Posterior : sigma2 | x ~ InvGamma(alpha + n/2, scale = beta + SS/2)
#
# InvGamma(a, b) moments:
# Mean = b / (a – 1) for a > 1
# Var = b^2 / ((a-1)^2 * (a-2)) for a > 2
#
# Quantile: 1 / qgamma(1 – p, shape = a, rate = b)
# = 1 / Gamma(a, scale=1/b).ppf(1 – p)
#
# Input : scalar / list / vector OR DataFrame + var + group_var
# Output : clean console tables, Plotly HTML saved to disk
# Plots : Plotly only — saved as HTML, never auto-rendered
# =============================================================
import os
import math
import numpy as np
import pandas as pd
from scipy.stats import gamma as sp_gamma
import plotly.graph_objects as go
from plotly.subplots import make_subplots
# =============================================================
# INTERNAL UTILITIES
# =============================================================
def _print_table(title: str, df: pd.DataFrame) -> None:
“””Plain text table to console — no color, no clutter.”””
col_names = list(df.columns)
col_widths = [
max(len(nm), df[nm].astype(str).map(len).max())
for nm in col_names
]
col_widths = [max(w, 10) for w in col_widths]
total_width = sum(col_widths) + 3 * (len(col_widths) – 1) + 2
divider = “=” * total_width
thin_div = “-” * total_width
header = ” “.join(nm.ljust(w) for nm, w in zip(col_names, col_widths))
lines = []
lines.append(divider)
lines.append(f” {title}”)
lines.append(divider)
lines.append(f” {header}”)
lines.append(f” {thin_div}”)
for _, row in df.iterrows():
row_str = ” “.join(str(val).ljust(w) for val, w in zip(row, col_widths))
lines.append(f” {row_str}”)
lines.append(divider)
lines.append(“”)
print(“\n”.join(lines))
def _fmt(v, digits: int = 6) -> str:
“””Format float; return ‘undefined’ for None/NaN.”””
if v is None or (isinstance(v, float) and math.isnan(v)):
return “undefined”
return f”{v:.{digits}f}”
def _expand_cri(cri_levels, quantiles):
“””Expand cri_levels to symmetric tail quantile pairs; merge with raw quantiles.”””
q_all = sorted(set(quantiles)) if quantiles else []
if cri_levels:
for lev in sorted(set(cri_levels)):
lo = (1 – lev / 100) / 2
hi = 1 – lo
q_all = sorted(set(q_all + [lo, hi]))
return q_all
def _build_cri_table(cri_levels, all_q, q_vals, cri_type):
“””Build CrI table from cri_levels only — symmetric pairs.”””
if not cri_levels:
return None
rows = []
for lev in sorted(set(cri_levels)):
lo = (1 – lev / 100) / 2
hi = 1 – lo
idx_lo = next((i for i, q in enumerate(all_q) if abs(q – lo) < 1e-10), None)
idx_hi = next((i for i, q in enumerate(all_q) if abs(q - hi) < 1e-10), None)
if idx_lo is not None and idx_hi is not None:
rows.append({
"Level": f"{lev}% {cri_type}",
"Lower": f"{q_vals[idx_lo]:.6f}",
"Upper": f"{q_vals[idx_hi]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _build_quantile_table(quantiles, all_q, q_vals):
"""Build individual quantile table from quantiles arg only."""
if not quantiles:
return None
rows = []
for q in sorted(set(quantiles)):
idx = next((i for i, v in enumerate(all_q) if abs(v - q) < 1e-10), None)
if idx is not None:
rows.append({
"Quantile": f"Q({q:.4f})",
"Value": f"{q_vals[idx]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _validate_df_var(data: pd.DataFrame, var: str) -> None:
“””Validate DataFrame and column for the Normal-InvGamma model.”””
if not isinstance(data, pd.DataFrame):
raise TypeError(“`data` must be a pandas DataFrame”)
if not isinstance(var, str):
raise TypeError(“`var` must be a single column name (string)”)
if var not in data.columns:
raise ValueError(f”Column ‘{var}’ not found in data”)
col = data[var].dropna()
if col.empty:
raise ValueError(f”Column ‘{var}’ is entirely NA”)
if not pd.api.types.is_numeric_dtype(col):
raise TypeError(f”Column ‘{var}’ must be numeric for Normal conjugate models”)
def _validate_group(data: pd.DataFrame, group_var) -> None:
if group_var is not None and group_var not in data.columns:
raise ValueError(f”group_var column ‘{group_var}’ not found in data”)
def _split_entries(data: pd.DataFrame, var: str, group_var):
“””Split data into labelled entries by group.”””
if group_var is not None:
entries = []
for g in data[group_var].unique():
sub = data[data[group_var] == g]
entries.append({
“label”: f”{group_var} = {g}”,
“col”: sub[var].dropna().values
})
return entries
else:
return [{“label”: “All”, “col”: data[var].dropna().values}]
def _resolve_plot_path(plot_path, plot_filename, default_name: str) -> str:
“””Resolve final HTML save path.”””
directory = plot_path if plot_path is not None else os.getcwd()
filename = (plot_filename if plot_filename is not None else default_name)
if not filename.endswith(“.html”):
filename += “.html”
return os.path.join(directory, filename)
def _save_plot(fig, path: str) -> None:
fig.write_html(path, auto_open=False)
print(f” [Plot saved] {path}\n”)
# =============================================================
# NORMAL – INVERSE GAMMA (estimating VARIANCE | mean known)
# =============================================================
def conjugate_normal_igamma(
data = None,
var = None,
group_var = None,
x = None,
mu_known = None,
alpha = None,
beta = None,
cri_levels = (90, 95),
quantiles = None,
cri_type = “CrI”,
plot_path = None,
plot_filename = None,
):
cri_levels = list(cri_levels) if cri_levels else []
quantiles = list(quantiles) if quantiles else []
all_q = _expand_cri(cri_levels, quantiles)
def _ig_mean(a, b): return b / (a – 1) if a > 1 else None
def _ig_var (a, b): return b**2 / ((a-1)**2 * (a-2)) if a > 2 else None
def _ig_q (p, a, b): return 1 / sp_gamma.ppf(1 – p, a=a, scale=1/b)
def _ig_pdf (s, a, b): return sp_gamma.pdf(1/s, a=a, scale=1/b) / s**2
if data is not None:
_validate_df_var(data, var)
_validate_group(data, group_var)
entries_raw = _split_entries(data, var, group_var)
entries = [{“label”: e[“label”], “x”: e[“col”]} for e in entries_raw]
else:
if x is None:
raise ValueError(“Provide `data` + `var` OR `x`”)
entries = [{“label”: “All”, “x”: np.atleast_1d(np.asarray(x, dtype=float))}]
results = []
for e in entries:
xv = np.asarray(e[“x”], dtype=float)
nv = len(xv)
ss = float(np.sum((xv – mu_known)**2))
a_post = alpha + nv / 2
b_post = beta + ss / 2
pr_mean = _ig_mean(alpha, beta)
pr_var = _ig_var (alpha, beta)
pr_mode = beta / (alpha + 1)
pr_median = _ig_q(0.5, alpha, beta)
po_mean = _ig_mean(a_post, b_post)
po_var = _ig_var (a_post, b_post)
po_mode = b_post / (a_post + 1)
po_median = _ig_q(0.5, a_post, b_post)
all_q_int = sorted(set(all_q + [0.5]))
q_vals = [_ig_q(q, a_post, b_post) for q in all_q_int]
_print_table(
f”Normal-IGamma | {e[‘label’]} | Data Summary”,
pd.DataFrame([{
“n”: nv,
“mu_known”: f”{mu_known:.6f}”,
“SS = sum((x-mu)^2)”: f”{ss:.6f}”,
}])
)
_print_table(
f”Normal-IGamma | {e[‘label’]} | Parameters & Moments (Variance)”,
pd.DataFrame({
“Quantity”: [“alpha (shape)”, “beta (scale)”, “Mean”, “Median”, “Mode”, “Variance”, “SD”],
“Prior”: [f”{alpha:.6f}”, f”{beta:.6f}”, _fmt(pr_mean),
f”{pr_median:.6f}”, f”{pr_mode:.6f}”,
_fmt(pr_var),
_fmt(math.sqrt(pr_var) if pr_var is not None else None)],
“Posterior”: [f”{a_post:.6f}”, f”{b_post:.6f}”, _fmt(po_mean),
f”{po_median:.6f}”, f”{po_mode:.6f}”,
_fmt(po_var),
_fmt(math.sqrt(po_var) if po_var is not None else None)],
})
)
cri_df = _build_cri_table(cri_levels, all_q_int, q_vals, cri_type)
if cri_df is not None:
_print_table(
f”Normal-IGamma | {e[‘label’]} | Credible Intervals ({cri_type}) [variance scale]”,
cri_df
)
q_df = _build_quantile_table(quantiles, all_q_int, q_vals)
if q_df is not None:
_print_table(
f”Normal-IGamma | {e[‘label’]} | Posterior Quantiles [variance scale]”,
q_df
)
# — Plot —
pr_lo = max(_ig_q(0.001, alpha, beta) * 0.5, 1e-6)
pr_hi = _ig_q(0.999, alpha, beta)
po_lo = max(_ig_q(0.001, a_post, b_post) * 0.5, 1e-6)
po_hi = _ig_q(0.999, a_post, b_post)
s2_pr = np.linspace(pr_lo, pr_hi, 1000)
s2_po = np.linspace(po_lo, po_hi, 1000)
fig = make_subplots(rows=1, cols=2, subplot_titles=[“Prior”, “Posterior”])
fig.add_trace(go.Scatter(
x=s2_pr, y=_ig_pdf(s2_pr, alpha, beta),
mode=”lines”, line=dict(color=”#2C7BB6″, width=1.5), showlegend=False
), row=1, col=1)
fig.add_trace(go.Scatter(
x=s2_po, y=_ig_pdf(s2_po, a_post, b_post),
mode=”lines”, line=dict(color=”#2C7BB6″, width=1.5), showlegend=False
), row=1, col=2)
fig.update_xaxes(title_text=”σ² (variance)”, row=1, col=1)
fig.update_xaxes(title_text=”σ² (variance)”, row=1, col=2)
fig.update_yaxes(title_text=”Density”, row=1, col=1)
fig.update_layout(
title_text=(
f”Normal-InvGamma (Variance) [{e[‘label’]}]
”
f”Prior: IG({alpha:.2f}, {beta:.2f}) | ”
f”Posterior: IG({a_post:.2f}, {b_post:.2f})”
),
template=”plotly_white”,
)
default_name = f”normal_igamma_{e[‘label’].replace(‘ ‘, ‘_’)}”
fpath = _resolve_plot_path(plot_path, plot_filename, default_name)
_save_plot(fig, fpath)
results.append({
“label”: e[“label”],
“data”: {“n”: nv, “mu_known”: mu_known, “ss”: ss},
“prior”: {“alpha”: alpha, “beta”: beta,
“mean”: pr_mean, “median”: pr_median, “mode”: pr_mode,
“var”: pr_var,
“sd”: math.sqrt(pr_var) if pr_var is not None else None},
“posterior”: {“alpha”: a_post, “beta”: b_post,
“mean”: po_mean, “median”: po_median, “mode”: po_mode,
“var”: po_var,
“sd”: math.sqrt(po_var) if po_var is not None else None,
“quantiles”: dict(zip(all_q_int, q_vals))},
“plot”: fig,
})
return results[0] if len(results) == 1 else results
# =============================================================
# USAGE EXAMPLES
# =============================================================
if __name__ == “__main__”:
# scalar / list
# res = conjugate_normal_igamma(x=[2.1,1.9,2.3,2.0], mu_known=2,
# alpha=3, beta=1, cri_levels=(90, 95))
# df + var
# df = pd.DataFrame({“measure”: [2.1,1.9,2.3,2.0,2.5,2.2,1.8,2.4]})
# res = conjugate_normal_igamma(data=df, var=”measure”, mu_known=2,
# alpha=3, beta=1, cri_levels=(90, 95),
# quantiles=(0.01, 0.10))
res = conjugate_normal_igamma(
x=[2.1, 1.9, 2.3, 2.0], mu_known=2, alpha=3, beta=1,
cri_levels=(90, 95),
plot_path=”.”, plot_filename=”normal_igamma_result”
)
# =============================================================
# Bayesian Conjugate Model - Normal-Inverse-Gamma (Variance, mean known)
# Closed-form analytical solution (standalone script)
#
# Likelihood : X_i ~ N(mu_known, sigma2)
# Prior : sigma2 ~ InvGamma(alpha, scale = beta)
# Posterior : sigma2 | x ~ InvGamma(alpha + n/2, scale = beta + SS/2)
#
# InvGamma(a, b) moments:
# Mean = b / (a - 1) for a > 1
# Var = b^2 / ((a-1)^2 * (a-2)) for a > 2
#
# Quantile: 1 / qgamma(1 - p, shape = a, rate = b)
# = 1 / Gamma(a, scale=1/b).ppf(1 - p)
#
# Input : scalar / list / vector OR DataFrame + var + group_var
# Output : clean console tables, Plotly HTML saved to disk
# Plots : Plotly only — saved as HTML, never auto-rendered
# =============================================================
import os
import math
import numpy as np
import pandas as pd
from scipy.stats import gamma as sp_gamma
import plotly.graph_objects as go
from plotly.subplots import make_subplots
# =============================================================
# INTERNAL UTILITIES
# =============================================================
def _print_table(title: str, df: pd.DataFrame) -> None:
"""Plain text table to console — no color, no clutter."""
col_names = list(df.columns)
col_widths = [
max(len(nm), df[nm].astype(str).map(len).max())
for nm in col_names
]
col_widths = [max(w, 10) for w in col_widths]
total_width = sum(col_widths) + 3 * (len(col_widths) - 1) + 2
divider = "=" * total_width
thin_div = "-" * total_width
header = " ".join(nm.ljust(w) for nm, w in zip(col_names, col_widths))
lines = []
lines.append(divider)
lines.append(f" {title}")
lines.append(divider)
lines.append(f" {header}")
lines.append(f" {thin_div}")
for _, row in df.iterrows():
row_str = " ".join(str(val).ljust(w) for val, w in zip(row, col_widths))
lines.append(f" {row_str}")
lines.append(divider)
lines.append("")
print("\n".join(lines))
def _fmt(v, digits: int = 6) -> str:
"""Format float; return 'undefined' for None/NaN."""
if v is None or (isinstance(v, float) and math.isnan(v)):
return "undefined"
return f"{v:.{digits}f}"
def _expand_cri(cri_levels, quantiles):
"""Expand cri_levels to symmetric tail quantile pairs; merge with raw quantiles."""
q_all = sorted(set(quantiles)) if quantiles else []
if cri_levels:
for lev in sorted(set(cri_levels)):
lo = (1 - lev / 100) / 2
hi = 1 - lo
q_all = sorted(set(q_all + [lo, hi]))
return q_all
def _build_cri_table(cri_levels, all_q, q_vals, cri_type):
"""Build CrI table from cri_levels only — symmetric pairs."""
if not cri_levels:
return None
rows = []
for lev in sorted(set(cri_levels)):
lo = (1 - lev / 100) / 2
hi = 1 - lo
idx_lo = next((i for i, q in enumerate(all_q) if abs(q - lo) < 1e-10), None)
idx_hi = next((i for i, q in enumerate(all_q) if abs(q - hi) < 1e-10), None)
if idx_lo is not None and idx_hi is not None:
rows.append({
"Level": f"{lev}% {cri_type}",
"Lower": f"{q_vals[idx_lo]:.6f}",
"Upper": f"{q_vals[idx_hi]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _build_quantile_table(quantiles, all_q, q_vals):
"""Build individual quantile table from quantiles arg only."""
if not quantiles:
return None
rows = []
for q in sorted(set(quantiles)):
idx = next((i for i, v in enumerate(all_q) if abs(v - q) < 1e-10), None)
if idx is not None:
rows.append({
"Quantile": f"Q({q:.4f})",
"Value": f"{q_vals[idx]:.6f}",
})
return pd.DataFrame(rows) if rows else None
def _validate_df_var(data: pd.DataFrame, var: str) -> None:
"""Validate DataFrame and column for the Normal-InvGamma model."""
if not isinstance(data, pd.DataFrame):
raise TypeError("`data` must be a pandas DataFrame")
if not isinstance(var, str):
raise TypeError("`var` must be a single column name (string)")
if var not in data.columns:
raise ValueError(f"Column '{var}' not found in data")
col = data[var].dropna()
if col.empty:
raise ValueError(f"Column '{var}' is entirely NA")
if not pd.api.types.is_numeric_dtype(col):
raise TypeError(f"Column '{var}' must be numeric for Normal conjugate models")
def _validate_group(data: pd.DataFrame, group_var) -> None:
if group_var is not None and group_var not in data.columns:
raise ValueError(f"group_var column '{group_var}' not found in data")
def _split_entries(data: pd.DataFrame, var: str, group_var):
"""Split data into labelled entries by group."""
if group_var is not None:
entries = []
for g in data[group_var].unique():
sub = data[data[group_var] == g]
entries.append({
"label": f"{group_var} = {g}",
"col": sub[var].dropna().values
})
return entries
else:
return [{"label": "All", "col": data[var].dropna().values}]
def _resolve_plot_path(plot_path, plot_filename, default_name: str) -> str:
"""Resolve final HTML save path."""
directory = plot_path if plot_path is not None else os.getcwd()
filename = (plot_filename if plot_filename is not None else default_name)
if not filename.endswith(".html"):
filename += ".html"
return os.path.join(directory, filename)
def _save_plot(fig, path: str) -> None:
fig.write_html(path, auto_open=False)
print(f" [Plot saved] {path}\n")
# =============================================================
# NORMAL - INVERSE GAMMA (estimating VARIANCE | mean known)
# =============================================================
def conjugate_normal_igamma(
data = None,
var = None,
group_var = None,
x = None,
mu_known = None,
alpha = None,
beta = None,
cri_levels = (90, 95),
quantiles = None,
cri_type = "CrI",
plot_path = None,
plot_filename = None,
):
cri_levels = list(cri_levels) if cri_levels else []
quantiles = list(quantiles) if quantiles else []
all_q = _expand_cri(cri_levels, quantiles)
def _ig_mean(a, b): return b / (a - 1) if a > 1 else None
def _ig_var (a, b): return b**2 / ((a-1)**2 * (a-2)) if a > 2 else None
def _ig_q (p, a, b): return 1 / sp_gamma.ppf(1 - p, a=a, scale=1/b)
def _ig_pdf (s, a, b): return sp_gamma.pdf(1/s, a=a, scale=1/b) / s**2
if data is not None:
_validate_df_var(data, var)
_validate_group(data, group_var)
entries_raw = _split_entries(data, var, group_var)
entries = [{"label": e["label"], "x": e["col"]} for e in entries_raw]
else:
if x is None:
raise ValueError("Provide `data` + `var` OR `x`")
entries = [{"label": "All", "x": np.atleast_1d(np.asarray(x, dtype=float))}]
results = []
for e in entries:
xv = np.asarray(e["x"], dtype=float)
nv = len(xv)
ss = float(np.sum((xv - mu_known)**2))
a_post = alpha + nv / 2
b_post = beta + ss / 2
pr_mean = _ig_mean(alpha, beta)
pr_var = _ig_var (alpha, beta)
pr_mode = beta / (alpha + 1)
pr_median = _ig_q(0.5, alpha, beta)
po_mean = _ig_mean(a_post, b_post)
po_var = _ig_var (a_post, b_post)
po_mode = b_post / (a_post + 1)
po_median = _ig_q(0.5, a_post, b_post)
all_q_int = sorted(set(all_q + [0.5]))
q_vals = [_ig_q(q, a_post, b_post) for q in all_q_int]
_print_table(
f"Normal-IGamma | {e['label']} | Data Summary",
pd.DataFrame([{
"n": nv,
"mu_known": f"{mu_known:.6f}",
"SS = sum((x-mu)^2)": f"{ss:.6f}",
}])
)
_print_table(
f"Normal-IGamma | {e['label']} | Parameters & Moments (Variance)",
pd.DataFrame({
"Quantity": ["alpha (shape)", "beta (scale)", "Mean", "Median", "Mode", "Variance", "SD"],
"Prior": [f"{alpha:.6f}", f"{beta:.6f}", _fmt(pr_mean),
f"{pr_median:.6f}", f"{pr_mode:.6f}",
_fmt(pr_var),
_fmt(math.sqrt(pr_var) if pr_var is not None else None)],
"Posterior": [f"{a_post:.6f}", f"{b_post:.6f}", _fmt(po_mean),
f"{po_median:.6f}", f"{po_mode:.6f}",
_fmt(po_var),
_fmt(math.sqrt(po_var) if po_var is not None else None)],
})
)
cri_df = _build_cri_table(cri_levels, all_q_int, q_vals, cri_type)
if cri_df is not None:
_print_table(
f"Normal-IGamma | {e['label']} | Credible Intervals ({cri_type}) [variance scale]",
cri_df
)
q_df = _build_quantile_table(quantiles, all_q_int, q_vals)
if q_df is not None:
_print_table(
f"Normal-IGamma | {e['label']} | Posterior Quantiles [variance scale]",
q_df
)
# --- Plot ---
pr_lo = max(_ig_q(0.001, alpha, beta) * 0.5, 1e-6)
pr_hi = _ig_q(0.999, alpha, beta)
po_lo = max(_ig_q(0.001, a_post, b_post) * 0.5, 1e-6)
po_hi = _ig_q(0.999, a_post, b_post)
s2_pr = np.linspace(pr_lo, pr_hi, 1000)
s2_po = np.linspace(po_lo, po_hi, 1000)
fig = make_subplots(rows=1, cols=2, subplot_titles=["Prior", "Posterior"])
fig.add_trace(go.Scatter(
x=s2_pr, y=_ig_pdf(s2_pr, alpha, beta),
mode="lines", line=dict(color="#2C7BB6", width=1.5), showlegend=False
), row=1, col=1)
fig.add_trace(go.Scatter(
x=s2_po, y=_ig_pdf(s2_po, a_post, b_post),
mode="lines", line=dict(color="#2C7BB6", width=1.5), showlegend=False
), row=1, col=2)
fig.update_xaxes(title_text="σ² (variance)", row=1, col=1)
fig.update_xaxes(title_text="σ² (variance)", row=1, col=2)
fig.update_yaxes(title_text="Density", row=1, col=1)
fig.update_layout(
title_text=(
f"Normal-InvGamma (Variance) [{e['label']}]<br>"
f"<sup>Prior: IG({alpha:.2f}, {beta:.2f}) | "
f"Posterior: IG({a_post:.2f}, {b_post:.2f})</sup>"
),
template="plotly_white",
)
default_name = f"normal_igamma_{e['label'].replace(' ', '_')}"
fpath = _resolve_plot_path(plot_path, plot_filename, default_name)
_save_plot(fig, fpath)
results.append({
"label": e["label"],
"data": {"n": nv, "mu_known": mu_known, "ss": ss},
"prior": {"alpha": alpha, "beta": beta,
"mean": pr_mean, "median": pr_median, "mode": pr_mode,
"var": pr_var,
"sd": math.sqrt(pr_var) if pr_var is not None else None},
"posterior": {"alpha": a_post, "beta": b_post,
"mean": po_mean, "median": po_median, "mode": po_mode,
"var": po_var,
"sd": math.sqrt(po_var) if po_var is not None else None,
"quantiles": dict(zip(all_q_int, q_vals))},
"plot": fig,
})
return results[0] if len(results) == 1 else results
# =============================================================
# USAGE EXAMPLES
# =============================================================
if __name__ == "__main__":
# scalar / list
# res = conjugate_normal_igamma(x=[2.1,1.9,2.3,2.0], mu_known=2,
# alpha=3, beta=1, cri_levels=(90, 95))
# df + var
# df = pd.DataFrame({"measure": [2.1,1.9,2.3,2.0,2.5,2.2,1.8,2.4]})
# res = conjugate_normal_igamma(data=df, var="measure", mu_known=2,
# alpha=3, beta=1, cri_levels=(90, 95),
# quantiles=(0.01, 0.10))
res = conjugate_normal_igamma(
x=[2.1, 1.9, 2.3, 2.0], mu_known=2, alpha=3, beta=1,
cri_levels=(90, 95),
plot_path=".", plot_filename="normal_igamma_result"
)if name == “main“:
# scalar / list
# res = conjugate_normal_igamma(x=[2.1,1.9,2.3,2.0], mu_known=2,
# alpha=3, beta=1, cri_levels=(90, 95))
The output can be stored as an R object, allowing users to further enhance and customize the generated plots using the extensive capabilities of plotly, including options for graphical modifications.
The provided script generates and displays a basic plotly object on screen, which can be subsequently modified according to the desired presentation style, formatting preferences, or reporting requirements.