Show code cell source Hide code cell source

MAKE_BOOK_FIGURES=True
import numpy as np
import scipy.stats as st

import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib_inline
matplotlib_inline.backend_inline.set_matplotlib_formats('svg')
import seaborn as sns
sns.set_context("paper")
sns.set_style("ticks")

def set_book_style():
    plt.style.use('seaborn-v0_8-white') 
    sns.set_style("ticks")
    sns.set_palette("deep")

    mpl.rcParams.update({
        # Font settings
        'font.family': 'serif',  # For academic publishing
        'font.size': 8,  # As requested, 10pt font
        'axes.labelsize': 8,
        'axes.titlesize': 8,
        'xtick.labelsize': 7,  # Slightly smaller for better readability
        'ytick.labelsize': 7,
        'legend.fontsize': 7,
        
        # Line and marker settings for consistency
        'axes.linewidth': 0.5,
        'grid.linewidth': 0.5,
        'lines.linewidth': 1.0,
        'lines.markersize': 4,
        
        # Layout to prevent clipped labels
        'figure.constrained_layout.use': True,
        
        # Default DPI (will override when saving)
        'figure.dpi': 600,
        'savefig.dpi': 600,
        
        # Despine - remove top and right spines
        'axes.spines.top': False,
        'axes.spines.right': False,
        
        # Remove legend frame
        'legend.frameon': False,
        
        # Additional trim settings
        'figure.autolayout': True,  # Alternative to constrained_layout
        'savefig.bbox': 'tight',    # Trim when saving
        'savefig.pad_inches': 0.1   # Small padding to ensure nothing gets cut off
    })

def set_notebook_style():
    plt.style.use('seaborn-v0_8-white')
    sns.set_style("ticks")
    sns.set_palette("deep")

    mpl.rcParams.update({
        # Font settings - using default sizes
        'font.family': 'serif',
        'axes.labelsize': 10,
        'axes.titlesize': 10,
        'xtick.labelsize': 9,
        'ytick.labelsize': 9,
        'legend.fontsize': 9,
        
        # Line and marker settings
        'axes.linewidth': 0.5,
        'grid.linewidth': 0.5,
        'lines.linewidth': 1.0,
        'lines.markersize': 4,
        
        # Layout settings
        'figure.constrained_layout.use': True,
        
        # Remove only top and right spines
        'axes.spines.top': False,
        'axes.spines.right': False,
        
        # Remove legend frame
        'legend.frameon': False,
        
        # Additional settings
        'figure.autolayout': True,
        'savefig.bbox': 'tight',
        'savefig.pad_inches': 0.1
    })

def save_for_book(fig, filename, is_vector=True, **kwargs):
    """
    Save a figure with book-optimized settings.
    
    Parameters:
    -----------
    fig : matplotlib figure
        The figure to save
    filename : str
        Filename without extension
    is_vector : bool
        If True, saves as vector at 1000 dpi. If False, saves as raster at 600 dpi.
    **kwargs : dict
        Additional kwargs to pass to savefig
    """    
    # Set appropriate DPI and format based on figure type
    if is_vector:
        dpi = 1000
        ext = '.pdf'
    else:
        dpi = 600
        ext = '.tif'
    
    # Save the figure with book settings
    fig.savefig(f"{filename}{ext}", dpi=dpi, **kwargs)

def make_full_width_fig():
    return plt.subplots(figsize=(4.7, 2.9), constrained_layout=True)

def make_half_width_fig():
    return plt.subplots(figsize=(2.35, 1.45), constrained_layout=True)

if MAKE_BOOK_FIGURES:
    set_book_style()
else:
    set_notebook_style()

make_full_width_fig = make_full_width_fig if MAKE_BOOK_FIGURES else lambda: plt.subplots()
make_half_width_fig = make_half_width_fig if MAKE_BOOK_FIGURES else lambda: plt.subplots()

The Principle of Maximum Entropy for Continuous Random Variables

Maximum Entropy Code

Writing generic code for finding maximum entropy distributions can be a lot of work. The compact implementation below is sufficient for our one-dimensional moment examples. It follows the moment-based formulation illustrated by the PyMaxEnt paper and its reference source, while keeping the implementation local and reproducible. It uses fixed Gauss–Legendre quadrature and SciPy’s nonlinear least-squares solver, so the notebook does not download code at runtime.

from scipy.optimize import least_squares


def reconstruct(moments, bnds=(-1.0, 1.0), quadrature_order=256):
    """Reconstruct a 1D maximum-entropy density from raw moments."""
    moments = np.asarray(moments, dtype=float)
    lower, upper = map(float, bnds)
    powers = np.arange(moments.size)

    nodes, weights = np.polynomial.legendre.leggauss(quadrature_order)
    x_quad = 0.5 * (upper - lower) * nodes + 0.5 * (upper + lower)
    w_quad = 0.5 * (upper - lower) * weights
    basis = x_quad[:, None] ** powers[None, :]

    def moment_residual(lambdas):
        density = np.exp(np.clip(basis @ lambdas, -700.0, 700.0))
        fitted = basis.T @ (w_quad * density)
        return fitted - moments

    initial = np.zeros(moments.size)
    initial[0] = np.log(moments[0] / (upper - lower))
    solution = least_squares(moment_residual, initial, max_nfev=5000)
    if not solution.success or np.linalg.norm(solution.fun, ord=np.inf) > 1e-7:
        raise RuntimeError("The requested moments could not be reconstructed.")

    lambdas = solution.x

    def pdf(x):
        x = np.asarray(x, dtype=float)
        polynomial = sum(
            coefficient * x ** power
            for power, coefficient in enumerate(lambdas)
        )
        density = np.exp(np.clip(polynomial, -700.0, 700.0))
        return np.where((x >= lower) & (x <= upper), density, 0.0)

    return pdf, lambdas

The reconstruct function is now defined locally and is ready to use:

# No external module download is required; reconstruct is defined above.

Examples of maximum entropy distributions

We work in a 1D random variable setting. The local reconstruct function requires that you specify the interval support of the distribution, i.e., an interval \([a,b]\) outside of which the probability density function should be zero, and the \(M\) moments of the distribution, i.e.,

\[ \mathbb{E}[X^m] = \mu_m, \]

for \(m=0,\dots,M\). Then, the maximum entropy distribution that satisfies these constraints is given by:

\[ p(x) = 1_{[a,b]}(x)\exp\left\{\sum_{m=1}^M\lambda_mx^m\right\}, \]

where the \(\lambda_0,\dots,\lambda_M\) are fitted so that the constraints are satisfied. Note that there is no need for the normalization constant here because it has been absorbed in \(\lambda_0\). Let’s do some examples to gain some intuition.

No constraints in [-1,1]

The support is \([-1,1]\), and there are no moment constraints. You only have to specify the normalization constraint and the bounds:

mu = [1.0]
pdf, lambdas = reconstruct(mu, bnds=[-1.0, 1.0])

# plot the reconstructed solution
x = np.linspace(-1.0, 1.0, 100)

fig, ax = plt.subplots()
ax.plot(x, pdf(x))
ax.set_xlabel('$x$')
ax.set_ylabel('$p(x)$')
sns.despine(trim=True);

Mean constraint [-1,1]

Same as before, but we are now going to impose a mean constraint:

\[ \mathbb{E}[X] = \mu. \]

mu = [1.0, # The required normalization constraint
      0.0] # The mean constraint

pdf, lambdas = reconstruct(mu, bnds=[-1.0, 1.0])

# plot the reconstructed solution
x = np.linspace(-1.0, 1.0, 100)

fig, ax = plt.subplots()
ax.plot(x, pdf(x))
ax.set_xlabel('$x$')
ax.set_ylabel('$p(x)$');

Questions

• Modify the mean to \(\mu=0.1\) and observe the resulting maximum entropy pdf.

• Modify the mean to \(\mu=-0.1\) and observe the resulting maximum entropy pdf.

• Try \(\mu=0.9\). What happens to the maximum entropy pdf?

• Try \(\mu=1.1\). Why does the code break down?

Variance constraint

In addition to the mean constraint, we now include a variance constraint:

\[ \mathbb{V}[X] = \sigma^2. \]

The local reconstruct function works with raw moment constraints. Therefore, we must connect the variance to the second and first moments. Here is how to do this:

\[ \mathbb{E}[X^2] = \mathbb{V}[X] + \left(\mathbb{E}[X]\right)^2 = \sigma^2 + \mu^2. \]

mu = 0.0
sigma2 = 0.1
mus = [
    1.0, # The required normalization constraint
    mu,  # The mean constraint
    sigma2 + mu ** 2
] # The second moment constraint 

pdf, lambdas = reconstruct(mus, bnds=[-1.0, 1.0])

# plot the reconstructed solution
x = np.linspace(-1.0, 1.0, 100)

fig, ax = plt.subplots()
ax.plot(x, pdf(x))
ax.set_xlabel('$x$')
ax.set_ylabel('$p(x)$')
sns.despine(trim=True);

Questions

• Modify the variance to \(\sigma^2=0.3\) and observe the resulting maximum entropy pdf.

• Modify the variance to \(\sigma^2=0.4\) and observe the resulting maximum entropy pdf. Why did you get this abrupt change?

• Try \(\sigma^2=1\). Why does the code break down?