1.1. Model evidence

Calculation of the marginalized log likelihood for models of different polynomial order.

1.1.1. Import requirments

import numpy as np
import plotly.graph_objects as go

import bayesianLinearRegression

1.1.2. Define functions to generate the design matrix sinusoidal regression data

def getPolynomialBasisFunctions(M):
    basis_functions = [None for m in range(M+1)]
    for m in range(M+1):
        basis_functions[m] = lambda x, m=m: x**m
    return basis_functions


def buildDesignMatrixRow(x, basis_functions):
    M = len(basis_functions)
    design_matrix_row = np.empty(shape=M, dtype=np.double)
    for m in range(M):
        design_matrix_row[m] = basis_functions[m](x)
    return design_matrix_row


def buildDesignMatrix(x, basis_functions):
    M = len(basis_functions)
    N = len(x)
    design_matrix = np.empty(shape=(N, M), dtype=np.double)
    for n in range(N):
        design_matrix[n, :] = buildDesignMatrixRow(
            x=x[n], basis_functions=basis_functions)
    return design_matrix

1.1.3. Define a function to generate polynomial data

def generateData(x, sigma, coefs):
    basis_functions = getPolynomialBasisFunctions(M=len(coefs)-1)
    Phi = buildDesignMatrix(x=x, basis_functions=basis_functions)
    y = Phi @ coefs
    noise = np.random.normal(loc=0, scale=sigma, size=len(y))
    t = y + noise
    return y, t

1.1.4. Set estimation parameters

prior_precision = 10.0
likelihood_precision = 10.0

1.1.5. Generate data

N = 50
x = 1.0 + np.random.uniform(size=N)

# we generate data with M+1=5 coefficients, so that the marginalized log
# likelihood should attain its maximum at M=4 (see Figure at the bottom).
_, y = generateData(x=x, sigma=1.0/likelihood_precision,
                    coefs=np.array([-0.5, 0.5, -0.5, 0.5, -0.5]))

1.1.6. Calculate model evindences

Ms = np.arange(10)
log_evidences = [None for m in Ms]
for M in Ms:
    basis_functions = getPolynomialBasisFunctions(M=M)
    Phi = buildDesignMatrix(x=x, basis_functions=basis_functions)
    mN, SN = bayesianLinearRegression.batchWithSimplePrior(
        Phi=Phi, y=y, alpha=prior_precision, beta=likelihood_precision)
    log_evidences[M] = bayesianLinearRegression.computeLogEvidence(
        Phi=Phi, y=y, mN=mN, SN=SN,
        alpha=prior_precision, beta=likelihood_precision)

Plot models’ log evidences (log evidence should maximize for M=4; see Generate data above) —————————————————————–

fig = go.Figure()
trace = go.Scatter(x=Ms, y=log_evidences, mode="lines+markers",
                   line=dict(color="blue"))
fig.add_trace(trace)
fig.update_layout(xaxis_title="M",
                  yaxis_title=r"$\log p(\mathbf{y}|\alpha,\beta)$")
fig


# sphinx_gallery_thumbnail_path = '_static/model_evidence.png'