∑MathIntermediate

Exponential Family Distributions

Key Points

•
Exponential family distributions express many common probability models in a single template p(x| $η)$ = h(x) exp( $η^{T}$ T(x) − A( $η)) .$
•
They are powerful because T(x) are sufficient statistics, which means data can be compressed without losing information for inference.
•
The log-partition function A( $η)$ normalizes the density and its gradients give moments: ∇A( $η)$ = E[T(X)] and ∇^2A( $η)$ = Var[T(X)].
•
Many familiar distributions (Bernoulli, Poisson, Gaussian, Categorical, Multinomial, Exponential, Gamma) are in the exponential family.
•
MLE for exponential families reduces to matching observed sufficient statistics to their expectations, often solvable with convex optimization.
•
Conjugate priors exist and make Bayesian updating simple: posterior hyperparameters are prior hyperparameters plus data sufficient statistics.
•
Computations are numerically stable when using tricks like log-sum-exp and log1p/expit to avoid overflow/underflow.
•
In C++, fitting and updating exponential-family models typically runs in O(nk) time where n is data size and k is the number of sufficient statistics.

Prerequisites

→Basic probability and random variables — Understanding PMFs/PDFs, expectations, and independence is required to interpret distributions.
→Calculus (derivatives, gradients, Hessians) — MLE and properties of A(η) rely on differentiation and convexity.
→Linear algebra and vectors — Natural parameters and sufficient statistics are often vectors; gradients and Hessians are matrix-valued.
→Convex optimization basics — MLE in exponential families leverages concavity/convexity and Newton-like methods.
→Bayesian inference basics — Conjugate priors and posterior updates are central to exponential families.
→Numerical stability techniques — Stable evaluation of log-sum-exp, log1p, and logistic functions prevents overflow/underflow.
→C++ programming fundamentals — Required to implement estimators, iterate, and handle numeric precision.

Detailed Explanation

Tap terms for definitions

01Overview

Exponential family distributions form a unifying framework that captures a wide range of commonly used probability distributions through a single mathematical template. A distribution belongs to this family if its probability density (or mass) function can be written as p(x\mid\eta) = h(x) \exp(\eta^{\top} T(x) - A(\eta)), where T(x) are the sufficient statistics, \eta are the natural (canonical) parameters, h(x) is the base measure, and A(\eta) is the log-partition function that ensures normalization. This representation provides beautiful and useful properties: the gradient of A(\eta) equals the mean of T(x), the Hessian equals their covariance, and maximum likelihood estimation (MLE) becomes a convex optimization problem in the natural parameters.

Many standard models fit this template. For example, Bernoulli has T(x) = x and A(\eta) = \log(1+e^{\eta}); Poisson has T(x) = x and A(\eta) = e^{\eta}; Gaussian with known variance has T(x) = x and x^{2} terms depending on parameterization. In Bayesian analysis, exponential families admit conjugate priors, enabling closed-form posterior updates by simply adding data’s sufficient statistics to prior hyperparameters. These properties make exponential families central to statistics, machine learning (e.g., GLMs), and probabilistic modeling, particularly for scalable and stable algorithms.

02Intuition & Analogies

Imagine you run a factory that produces data-cookies. Each batch of cookies (your dataset) can be summarized by a few key ingredients: total sugar used, number of chocolate chips, and the average size. Rather than storing every cookie, you track just these summaries. In exponential families, those summaries are called sufficient statistics T(x). They’re like perfect recipe summaries: for the purposes of fitting or updating your model, no additional details from the batch are needed.

Now think of the control panel in your factory. You have a set of knobs that determine how sweet or crunchy your cookies turn out. Those knobs are the natural parameters \eta. Turn the sweetness knob up (increase \eta) and the expected sugar in the cookies (the mean of T(x)) also goes up. This relationship is governed by A(\eta), the log-partition function, which acts like an automatic calibration system ensuring everything remains a valid probability distribution. The first derivative of A tells you how the average summary (mean statistics) change with the knobs; the second derivative tells you how variable the summaries will be.

Finally, consider ordering supplies. If you know how many chips and how much sugar you used in previous batches, updating your plan for the next batch is simple—just add the latest batch summary to your current inventory. That’s conjugacy: the prior acts like your current inventory, the data contribute new sufficient statistics, and the posterior is the updated inventory. This is why exponential families are computation-friendly: learning from data reduces to accumulating and transforming a small set of summary numbers, leading to fast, stable algorithms.

03Formal Definition

A family of probability distributions on sample space

X

is an exponential family if there exist functions T:

X

\to

R^{k}

, h:

X

\to

R_{+}

, a log-partition A:

Ω

\subseteq

R^{k}

\to

R

, and a natural parameter

η

\in

Ω

such that the densities (or PMFs) have the form p(x

∣ η

) = h(x)

\exp\big

(

η^{⊤}

T(x) - A(

η

)

\big

\in

X

. The natural parameter space

Ω

is the set of

η

for which A(

η

) <

\infty

, where A(

η

) =

lo g

\int_{X}

h(x)

\exp\big

(

η^{⊤}

T(x)

\big

) \, d

μ

(x) with

μ

the base measure (counting or Lebesgue). If T(x) are linearly independent a.s., the family is minimal. The mean parameter is

μ

E_{η}

[T(X)] =

\nabla

η

), and the covariance of T(X) is

\nabla^{2}

η

), which is positive semidefinite, so A is convex. For IID samples

x_{1 : n}

, the log-likelihood is

ℓ

(

η

) =

η^{⊤}

\sum_{i = 1}^{n}

x_{i}

) - n A(

η

) +

\sum_{i = 1}^{n}

lo g

x_{i}

), showing that data enter only through the aggregate sufficient statistics

\sum_{i = 1}^{n}

x_{i}

). A prior is conjugate if it has the form p(

η ∣ χ

ν

)

\propto

exp

(

η^{⊤}

χ

ν

η

) - B(

χ

ν

)), which yields a posterior with updated hyperparameters (

χ

\sum

x_{i}

ν

+ n).

04When to Use

Use exponential family models when your data’s distribution is one of the many members (Bernoulli, Binomial, Categorical/Multinomial, Poisson, Exponential, Gamma, Gaussian variants, etc.) or when you need tractable learning and inference.

Generalized Linear Models (GLMs): Logistic regression (Bernoulli), Poisson regression (Poisson), and softmax regression (Multinomial) rely on exponential-family likelihoods and canonical links.
Bayesian updating at scale: Conjugate priors let you stream data and update posteriors by adding sufficient statistics, ideal for online learning and A/B testing.
Maximum entropy modeling: If you want the least-committal distribution subject to fixed expectations of features, the solution lies in the exponential family with those features as T(x).
Variational inference and message passing: Many algorithms exploit \nabla A and \nabla^{2}A identities to compute expectations and uncertainties efficiently.
Numerical stability and convexity: When you want MLE that is convex in natural parameters, exponential families provide guarantees and enable fast Newton or quasi-Newton methods with reliable convergence.

⚠️Common Mistakes

• Confusing natural and mean parameters: The natural parameter \eta is not usually the same as intuitive parameters like p or \lambda. For Bernoulli, \eta = \log\frac{p}{1-p}. Always use link functions to convert between spaces. • Ignoring the base measure h(x): Dropping h(x) can lead to wrong likelihoods, especially when comparing non-identical distributions or parameterizations. • Domain violations: Not all \eta are valid. The natural parameter space \Omega is defined by A(\eta) < \infty. For example, Gaussian with unknown variance requires constraints; ensure parameters lie in \Omega. • Numerical overflow/underflow: Naively computing \log(1+e^{\eta}) or \log \sum e^{\eta_i} can overflow. Use log1p/expit and log-sum-exp tricks for stability. • Non-minimal representations: Redundant components in T(x) can make inference ill-conditioned. Prefer minimal sufficient statistics when possible. • Misusing conjugacy: Conjugate priors depend on parameterization. Ensure your prior corresponds to the natural parameters (e.g., Beta prior is conjugate to Bernoulli in p but corresponds to linear natural parameters in \eta). • Confusing IID with exponential-family structure: Sufficient statistics add across IID samples, but dependence (time series, graphs) breaks simple additivity unless the model factors appropriately.

Key Formulas

Exponential Family Form

p (x ∣ η) = h (x) exp (η^{⊤} T (x) - A (η))

Explanation: Defines the family using sufficient statistics T(x), natural parameters $η,$ base measure h(x), and log-partition A( $η) .$ Use this to recognize whether a distribution is in the exponential family.

Log-Partition

A (η) = lo g \int h (x) exp (η^{⊤} T (x)) d μ (x)

Explanation: Computes the normalizing constant ensuring probabilities sum/integrate to 1. It is convex and central to moment computations.

Mean Parameters

\nabla A (η) = E_{η} [T (X)]

Explanation: The gradient of A equals the expected sufficient statistics. This lets you map natural parameters to mean parameters.

Covariance / Fisher Information

\nabla^{2} A (η) = Var_{η} [T (X)]

Explanation: The Hessian of A equals the covariance matrix of sufficient statistics, giving curvature and information content.

Log-Likelihood

ℓ (η) = i = 1 \sum n lo g p (x_{i} ∣ η) = η^{⊤} i = 1 \sum n T (x_{i}) - n A (η) + i = 1 \sum n lo g h (x_{i})

Explanation: The log-likelihood depends on data only through aggregated sufficient statistics. This is the basis for efficient MLE.

Score Function

\nabla ℓ (η) = i = 1 \sum n T (x_{i}) - n \nabla A (η)

Explanation: The gradient of the log-likelihood is observed statistics minus their expectations. Setting it to zero matches moments.

Concavity of Log-Likelihood

\nabla^{2} ℓ (η) = - n \nabla^{2} A (η) ⪯ 0

Explanation: For IID data, the log-likelihood is concave in η (or linear minus convex), enabling fast convex optimization methods.

Conjugate Prior

p (η ∣ χ, ν) \propto exp (η^{⊤} χ - ν A (η) - B (χ, ν))

Explanation: Conjugate priors maintain the exponential-family form in $η,$ yielding closed-form posterior updates.

Posterior Hyperparameters

χ^{'} = χ + i = 1 \sum n T (x_{i}), ν^{'} = ν + n

Explanation: Conjugate posterior updates are simple additions of sufficient statistics and counts to prior hyperparameters.

Bernoulli Canonical Form

p = σ (η) = \frac{1}{1 + e ^{- η}}, A (η) = lo g (1 + e^{η})

Explanation: Bernoulli’s mean is the logistic of $η,$ and A( $η)$ is log-sum-exp for two outcomes. Useful for logistic regression and MLE.

Poisson Canonical Form

λ = e^{η}, A (η) = e^{η}

Explanation: Poisson’s rate is the exponential of $η .$ This mapping is used in Poisson regression and count modeling.

Sufficient-Statistic Complexity

O (nk) time, O (k) space

Explanation: Computing and storing sufficient statistics for n samples with k-dimensional T(x) takes linear time in n and memory in k.

Complexity Analysis

Core computations with exponential families scale with the number of data points n and the dimension k of the sufficient statistics. Aggregating sufficient statistics S =

\sum_{i = 1}^{n}

x_{i}

) requires O(nk) time and O(k) space, since we need only a running sum. Evaluating the log-likelihood ℓ(

η)

or its gradient ∇ℓ(

η)

then costs O(k) once S is known, plus any cost to evaluate A(

η)

and ∇A(

η),

which is typically O(k) for common families (Bernoulli, Poisson, Categorical/Multinomial with

K = k + 1

). For MLE in natural parameters, Newton–Raphson or quasi-Newton methods are effective because ℓ is concave in

η .

Each Newton step requires computing the gradient and Hessian: ∇ℓ(

η)

= S − n∇A(

η)

and ∇^2ℓ(

η)

= −n∇^2A(

η) .

For scalar η (e.g., Bernoulli, Poisson), each iteration is O(n) to form S if not precomputed, but with S precomputed it becomes O(1) per iteration. For vector η ∈

R^{k}

, a Newton step also involves solving a k×k linear system, costing O(

k^{3}

) in general, but often k is small (e.g., category counts) and the Hessian has structure (diagonal for independent components or softmax structure for Multinomial), reducing cost to O(k)–O(

k^{2}

). Bayesian conjugate updates are especially efficient: updating hyperparameters (χ, ν) to (χ + S, ν + n) is O(k) time and O(k) space. In online/streaming contexts, maintaining S is O(k) memory regardless of n. Numerical stability may require log-domain computations (e.g., log-sum-exp), which add constant-factor overhead but do not change asymptotic complexity. Overall, exponential-family inference is near-optimal in n, typically linear, and modest in k, enabling scalability.

Code Examples

Bernoulli MLE in Natural Parameter via Newton–Raphson (Stable)

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 // Stable logistic function sigma(eta) = 1 / (1 + exp(-eta))
5 static inline double logistic(double eta) {
6     if (eta >= 0) {
7         double z = exp(-eta);
8         return 1.0 / (1.0 + z);
9     } else {
10         double z = exp(eta);
11         return z / (1.0 + z);
12     }
13 }
14 
15 // Stable log(1 + exp(eta))
16 static inline double log1pexp(double eta) {
17     if (eta > 0) return eta + log1p(exp(-eta));
18     return log1p(exp(eta));
19 }
20 
21 struct BernoulliMLE {
22     // Returns log-likelihood for given natural parameter eta
23     static double log_likelihood(const vector<int>& x, double eta) {
24         long long s = 0; // sum of x_i
25         for (int xi : x) s += xi;
26         long long n = (long long)x.size();
27         // l(eta) = eta * sum x - n * log(1 + exp(eta))
28         return eta * (double)s - (double)n * log1pexp(eta);
29     }
30 
31     // One Newton step: eta <- eta - g / H
32     static double newton_update(const vector<int>& x, double eta) {
33         long long s = 0; // sufficient statistic
34         for (int xi : x) s += xi;
35         double n = (double)x.size();
36         double p = logistic(eta);
37         double g = (double)s - n * p;                 // gradient
38         double H = - n * p * (1.0 - p);               // Hessian (negative)
39         // If H is near zero (data all 0 or all 1), damp to avoid huge steps
40         if (fabs(H) < 1e-12) H = (H < 0 ? -1e-12 : 1e-12);
41         return eta - g / H;
42     }
43 
44     // Fit eta using Newton's method with damping and bounds
45     static double fit_eta(const vector<int>& x, int max_iters = 50, double tol = 1e-10) {
46         // Initialize with logit of sample mean (with small smoothing)
47         long long s = 0; for (int xi : x) s += xi; double n = (double)x.size();
48         double p0 = (s + 0.5) / (n + 1.0); // Laplace smoothing to avoid 0 or 1
49         double eta = log(p0 / (1.0 - p0));
50         for (int it = 0; it < max_iters; ++it) {
51             double eta_new = newton_update(x, eta);
52             // Dampen if needed to ensure increase in log-likelihood
53             double ll_old = log_likelihood(x, eta);
54             double ll_new = log_likelihood(x, eta_new);
55             double step = 1.0;
56             while (ll_new < ll_old && step > 1e-6) {
57                 step *= 0.5;
58                 eta_new = eta + step * (eta_new - eta);
59                 ll_new = log_likelihood(x, eta_new);
60             }
61             if (fabs(eta_new - eta) < tol) return eta_new;
62             eta = eta_new;
63         }
64         return eta;
65     }
66 };
67 
68 int main() {
69     ios::sync_with_stdio(false);
70     cin.tie(nullptr);
71 
72     // Example data: 0/1 outcomes (e.g., coin flips)
73     vector<int> x = {1,0,1,1,0,1,0,1,1,0,1,1};
74 
75     double eta_hat = BernoulliMLE::fit_eta(x);
76     double p_hat = logistic(eta_hat);
77 
78     cout.setf(std::ios::fixed); cout<<setprecision(6);
79     cout << "Estimated natural parameter eta: " << eta_hat << "\n";
80     cout << "Estimated mean parameter p: " << p_hat << "\n";
81     cout << "Log-likelihood at eta_hat: " << BernoulliMLE::log_likelihood(x, eta_hat) << "\n";
82 
83     return 0;
84 }
85

This program fits a Bernoulli model using the exponential-family natural parameter η. It uses Newton–Raphson with line search to maximize the concave log-likelihood ℓ(η) = η Σ x_i − n log(1+e^{η}). The code includes numerically stable logistic and log1pexp helpers, Laplace-smoothed initialization, and damping to ensure monotonic improvement.

Time: O(n) per Newton iteration (O(1) if the sufficient statistic Σx_i is precomputed).Space: O(1) extra space beyond the input data.

Poisson MLE in Natural Parameter and Closed-Form Rate; Stable Computations

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 // Clamp to avoid overflow in exp
5 static inline double safe_exp(double x) {
6     const double MAX_ARG = 709.0; // ~ log(DBL_MAX)
7     if (x > MAX_ARG) return exp(MAX_ARG);
8     if (x < -1000.0) return 0.0; // underflow to 0
9     return exp(x);
10 }
11 
12 struct PoissonMLE {
13     static double log_likelihood(const vector<int>& x, double eta) {
14         long long s = 0; for (int xi : x) s += xi; double n = (double)x.size();
15         // l(eta) = eta * sum x - n * exp(eta) - sum log(x!)  (last term is const in eta)
16         return eta * (double)s - n * safe_exp(eta);
17     }
18 
19     static double newton_update(const vector<int>& x, double eta) {
20         long long s = 0; for (int xi : x) s += xi; double n = (double)x.size();
21         double e = safe_exp(eta);
22         double g = (double)s - n * e;     // gradient
23         double H = - n * e;               // Hessian (negative)
24         if (fabs(H) < 1e-18) H = (H < 0 ? -1e-18 : 1e-18);
25         return eta - g / H;
26     }
27 
28     static double fit_eta(const vector<int>& x, int max_iters = 30, double tol = 1e-12) {
29         // Initialize with log of sample mean (with small epsilon)
30         double mean_x = 0.0; for (int xi : x) mean_x += xi; mean_x /= max<size_t>(1, x.size());
31         double eta = log(max(1e-12, mean_x));
32         for (int it = 0; it < max_iters; ++it) {
33             double eta_new = newton_update(x, eta);
34             double ll_old = log_likelihood(x, eta);
35             double ll_new = log_likelihood(x, eta_new);
36             double step = 1.0;
37             while (ll_new < ll_old && step > 1e-7) {
38                 step *= 0.5;
39                 eta_new = eta + step * (eta_new - eta);
40                 ll_new = log_likelihood(x, eta_new);
41             }
42             if (fabs(eta_new - eta) < tol) return eta_new;
43             eta = eta_new;
44         }
45         return eta;
46     }
47 };
48 
49 int main() {
50     ios::sync_with_stdio(false);
51     cin.tie(nullptr);
52 
53     // Example count data
54     vector<int> x = {3, 2, 0, 4, 1, 2, 3, 5, 2, 1};
55 
56     // Closed-form MLE for rate: lambda_hat = mean(x)
57     double mean_x = 0.0; for (int xi : x) mean_x += xi; mean_x /= (double)x.size();
58     double lambda_hat = mean_x;
59     double eta_closed = log(max(1e-12, lambda_hat));
60 
61     // Fit eta via Newton as a cross-check
62     double eta_hat = PoissonMLE::fit_eta(x);
63 
64     cout.setf(std::ios::fixed); cout<<setprecision(6);
65     cout << "Closed-form lambda_hat: " << lambda_hat << "\n";
66     cout << "Closed-form eta (log lambda): " << eta_closed << "\n";
67     cout << "Newton-estimated eta: " << eta_hat << "\n";
68     cout << "Log-likelihood at eta_hat: " << PoissonMLE::log_likelihood(x, eta_hat) << "\n";
69 
70     return 0;
71 }
72

This program fits a Poisson model using both the closed-form MLE for the rate λ (the sample mean) and Newton’s method for the natural parameter η = log λ. The log-likelihood ℓ(η) = η Σx − n e^{η} + const is concave, and Newton converges in a few steps. safe_exp prevents overflows/underflows.

Time: O(n) to compute sufficient statistics; O(1) per Newton iteration with precomputed sums.Space: O(1) extra space beyond the input data.

Conjugate Posterior Updates: Beta–Bernoulli and Gamma–Poisson

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 struct BetaBernoulli {
5     // Prior: Beta(alpha, beta) on p
6     // Data: x_i in {0,1}
7     // Posterior: Beta(alpha + s, beta + f), where s = sum x_i, f = n - s
8     static pair<double,double> posterior(double alpha, double beta, const vector<int>& x) {
9         long long s = 0; for (int xi : x) s += xi; long long n = (long long)x.size();
10         return {alpha + (double)s, beta + (double)(n - s)};
11     }
12     static double posterior_mean(double alpha, double beta, const vector<int>& x) {
13         auto [a, b] = posterior(alpha, beta, x);
14         return a / (a + b);
15     }
16 };
17 
18 struct GammaPoisson {
19     // Prior: Gamma(a, b) with shape a > 0 and rate b > 0 on lambda
20     // Likelihood: Poisson(x | lambda)
21     // Posterior: Gamma(a + sum x_i, b + n)
22     static pair<double,double> posterior(double a, double b, const vector<int>& x) {
23         long long s = 0; for (int xi : x) s += xi; long long n = (long long)x.size();
24         return {a + (double)s, b + (double)n};
25     }
26     static double posterior_mean(double a, double b, const vector<int>& x) {
27         auto [ap, bp] = posterior(a, b, x);
28         return ap / bp; // mean of Gamma(shape, rate)
29     }
30 };
31 
32 int main() {
33     ios::sync_with_stdio(false);
34     cin.tie(nullptr);
35 
36     // Bernoulli example: prior Beta(2,2), data with 7 successes, 3 failures
37     vector<int> bern_x = {1,1,0,1,1,0,1,0,1,1};
38     auto [alpha_post, beta_post] = BetaBernoulli::posterior(2.0, 2.0, bern_x);
39     cout.setf(std::ios::fixed); cout<<setprecision(6);
40     cout << "Beta-Bernoulli posterior: alpha'=" << alpha_post << ", beta'=" << beta_post << "\n";
41     cout << "Posterior mean p: " << BetaBernoulli::posterior_mean(2.0, 2.0, bern_x) << "\n";
42 
43     // Poisson example: prior Gamma(3, 1), count data
44     vector<int> pois_x = {3,2,0,4,1,2,3,5,2,1};
45     auto [a_post, b_post] = GammaPoisson::posterior(3.0, 1.0, pois_x);
46     cout << "Gamma-Poisson posterior: a'=" << a_post << ", b'=" << b_post << "\n";
47     cout << "Posterior mean lambda: " << GammaPoisson::posterior_mean(3.0, 1.0, pois_x) << "\n";
48 
49     return 0;
50 }
51

This program demonstrates fast Bayesian updating using conjugate priors in exponential families. For Bernoulli, the Beta prior yields posterior Beta(α+s, β+f). For Poisson, the Gamma prior yields posterior Gamma(a+Σx, b+n). Posterior means are simple functions of updated hyperparameters.

Time: O(n) to aggregate sufficient statistics for each dataset.Space: O(1) extra space beyond the input.

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	// Stable logistic function sigma(eta) = 1 / (1 + exp(-eta))
5	static inline double logistic(double eta) {
6	if (eta >= 0) {
7	double z = exp(-eta);
8	return 1.0 / (1.0 + z);
9	} else {
10	double z = exp(eta);
11	return z / (1.0 + z);
12	}
13	}
14
15	// Stable log(1 + exp(eta))
16	static inline double log1pexp(double eta) {
17	if (eta > 0) return eta + log1p(exp(-eta));
18	return log1p(exp(eta));
19	}
20
21	struct BernoulliMLE {
22	// Returns log-likelihood for given natural parameter eta
23	static double log_likelihood(const vector<int>& x, double eta) {
24	long long s = 0; // sum of x_i
25	for (int xi : x) s += xi;
26	long long n = (long long)x.size();
27	// l(eta) = eta * sum x - n * log(1 + exp(eta))
28	return eta * (double)s - (double)n * log1pexp(eta);
29	}
30
31	// One Newton step: eta <- eta - g / H
32	static double newton_update(const vector<int>& x, double eta) {
33	long long s = 0; // sufficient statistic
34	for (int xi : x) s += xi;
35	double n = (double)x.size();
36	double p = logistic(eta);
37	double g = (double)s - n * p; // gradient
38	double H = - n * p * (1.0 - p); // Hessian (negative)
39	// If H is near zero (data all 0 or all 1), damp to avoid huge steps
40	if (fabs(H) < 1e-12) H = (H < 0 ? -1e-12 : 1e-12);
41	return eta - g / H;
42	}
43
44	// Fit eta using Newton's method with damping and bounds
45	static double fit_eta(const vector<int>& x, int max_iters = 50, double tol = 1e-10) {
46	// Initialize with logit of sample mean (with small smoothing)
47	long long s = 0; for (int xi : x) s += xi; double n = (double)x.size();
48	double p0 = (s + 0.5) / (n + 1.0); // Laplace smoothing to avoid 0 or 1
49	double eta = log(p0 / (1.0 - p0));
50	for (int it = 0; it < max_iters; ++it) {
51	double eta_new = newton_update(x, eta);
52	// Dampen if needed to ensure increase in log-likelihood
53	double ll_old = log_likelihood(x, eta);
54	double ll_new = log_likelihood(x, eta_new);
55	double step = 1.0;
56	while (ll_new < ll_old && step > 1e-6) {
57	step *= 0.5;
58	eta_new = eta + step * (eta_new - eta);
59	ll_new = log_likelihood(x, eta_new);
60	}
61	if (fabs(eta_new - eta) < tol) return eta_new;
62	eta = eta_new;
63	}
64	return eta;
65	}
66	};
67
68	int main() {
69	ios::sync_with_stdio(false);
70	cin.tie(nullptr);
71
72	// Example data: 0/1 outcomes (e.g., coin flips)
73	vector<int> x = {1,0,1,1,0,1,0,1,1,0,1,1};
74
75	double eta_hat = BernoulliMLE::fit_eta(x);
76	double p_hat = logistic(eta_hat);
77
78	cout.setf(std::ios::fixed); cout<<setprecision(6);
79	cout << "Estimated natural parameter eta: " << eta_hat << "\n";
80	cout << "Estimated mean parameter p: " << p_hat << "\n";
81	cout << "Log-likelihood at eta_hat: " << BernoulliMLE::log_likelihood(x, eta_hat) << "\n";
82
83	return 0;
84	}
85

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	// Clamp to avoid overflow in exp
5	static inline double safe_exp(double x) {
6	const double MAX_ARG = 709.0; // ~ log(DBL_MAX)
7	if (x > MAX_ARG) return exp(MAX_ARG);
8	if (x < -1000.0) return 0.0; // underflow to 0
9	return exp(x);
10	}
11
12	struct PoissonMLE {
13	static double log_likelihood(const vector<int>& x, double eta) {
14	long long s = 0; for (int xi : x) s += xi; double n = (double)x.size();
15	// l(eta) = eta * sum x - n * exp(eta) - sum log(x!) (last term is const in eta)
16	return eta * (double)s - n * safe_exp(eta);
17	}
18
19	static double newton_update(const vector<int>& x, double eta) {
20	long long s = 0; for (int xi : x) s += xi; double n = (double)x.size();
21	double e = safe_exp(eta);
22	double g = (double)s - n * e; // gradient
23	double H = - n * e; // Hessian (negative)
24	if (fabs(H) < 1e-18) H = (H < 0 ? -1e-18 : 1e-18);
25	return eta - g / H;
26	}
27
28	static double fit_eta(const vector<int>& x, int max_iters = 30, double tol = 1e-12) {
29	// Initialize with log of sample mean (with small epsilon)
30	double mean_x = 0.0; for (int xi : x) mean_x += xi; mean_x /= max<size_t>(1, x.size());
31	double eta = log(max(1e-12, mean_x));
32	for (int it = 0; it < max_iters; ++it) {
33	double eta_new = newton_update(x, eta);
34	double ll_old = log_likelihood(x, eta);
35	double ll_new = log_likelihood(x, eta_new);
36	double step = 1.0;
37	while (ll_new < ll_old && step > 1e-7) {
38	step *= 0.5;
39	eta_new = eta + step * (eta_new - eta);
40	ll_new = log_likelihood(x, eta_new);
41	}
42	if (fabs(eta_new - eta) < tol) return eta_new;
43	eta = eta_new;
44	}
45	return eta;
46	}
47	};
48
49	int main() {
50	ios::sync_with_stdio(false);
51	cin.tie(nullptr);
52
53	// Example count data
54	vector<int> x = {3, 2, 0, 4, 1, 2, 3, 5, 2, 1};
55
56	// Closed-form MLE for rate: lambda_hat = mean(x)
57	double mean_x = 0.0; for (int xi : x) mean_x += xi; mean_x /= (double)x.size();
58	double lambda_hat = mean_x;
59	double eta_closed = log(max(1e-12, lambda_hat));
60
61	// Fit eta via Newton as a cross-check
62	double eta_hat = PoissonMLE::fit_eta(x);
63
64	cout.setf(std::ios::fixed); cout<<setprecision(6);
65	cout << "Closed-form lambda_hat: " << lambda_hat << "\n";
66	cout << "Closed-form eta (log lambda): " << eta_closed << "\n";
67	cout << "Newton-estimated eta: " << eta_hat << "\n";
68	cout << "Log-likelihood at eta_hat: " << PoissonMLE::log_likelihood(x, eta_hat) << "\n";
69
70	return 0;
71	}
72

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	struct BetaBernoulli {
5	// Prior: Beta(alpha, beta) on p
6	// Data: x_i in {0,1}
7	// Posterior: Beta(alpha + s, beta + f), where s = sum x_i, f = n - s
8	static pair<double,double> posterior(double alpha, double beta, const vector<int>& x) {
9	long long s = 0; for (int xi : x) s += xi; long long n = (long long)x.size();
10	return {alpha + (double)s, beta + (double)(n - s)};
11	}
12	static double posterior_mean(double alpha, double beta, const vector<int>& x) {
13	auto [a, b] = posterior(alpha, beta, x);
14	return a / (a + b);
15	}
16	};
17
18	struct GammaPoisson {
19	// Prior: Gamma(a, b) with shape a > 0 and rate b > 0 on lambda
20	// Likelihood: Poisson(x \| lambda)
21	// Posterior: Gamma(a + sum x_i, b + n)
22	static pair<double,double> posterior(double a, double b, const vector<int>& x) {
23	long long s = 0; for (int xi : x) s += xi; long long n = (long long)x.size();
24	return {a + (double)s, b + (double)n};
25	}
26	static double posterior_mean(double a, double b, const vector<int>& x) {
27	auto [ap, bp] = posterior(a, b, x);
28	return ap / bp; // mean of Gamma(shape, rate)
29	}
30	};
31
32	int main() {
33	ios::sync_with_stdio(false);
34	cin.tie(nullptr);
35
36	// Bernoulli example: prior Beta(2,2), data with 7 successes, 3 failures
37	vector<int> bern_x = {1,1,0,1,1,0,1,0,1,1};
38	auto [alpha_post, beta_post] = BetaBernoulli::posterior(2.0, 2.0, bern_x);
39	cout.setf(std::ios::fixed); cout<<setprecision(6);
40	cout << "Beta-Bernoulli posterior: alpha'=" << alpha_post << ", beta'=" << beta_post << "\n";
41	cout << "Posterior mean p: " << BetaBernoulli::posterior_mean(2.0, 2.0, bern_x) << "\n";
42
43	// Poisson example: prior Gamma(3, 1), count data
44	vector<int> pois_x = {3,2,0,4,1,2,3,5,2,1};
45	auto [a_post, b_post] = GammaPoisson::posterior(3.0, 1.0, pois_x);
46	cout << "Gamma-Poisson posterior: a'=" << a_post << ", b'=" << b_post << "\n";
47	cout << "Posterior mean lambda: " << GammaPoisson::posterior_mean(3.0, 1.0, pois_x) << "\n";
48
49	return 0;
50	}
51