[pages=-,scale=1.1,offset=0mm -4mm,]imgs_header/front-cover.pdf

Bayesian Matrix Decomposition and Applications

Jun Lu

2023    Jun Lu

Version 0.1

CHAPTER: PREFACE

In 1954, Alston S. Householder published Principles of Numerical Analysis, an early and influential work on matrix decomposition, specifically favoring (block) LU decomposition—a method that factorizes a matrix into the product of lower and upper triangular matrices.
Today, matrix decomposition stands as a foundational technology in machine learning, driven significantly by the backpropagation algorithm's development for neural network fitting. It plays a crucial role in reducing data dimensionality and representing it in a format more amenable to machine learning algorithms.

Bayesian matrix decomposition is a relatively new field within the broader area of matrix decomposition and machine learning.
The concept of Bayesian matrix decomposition is rooted in Bayesian statistics, in which case it integrates the principles of Bayesian statistics with matrix decomposition methods to perform matrix factorization.
The use of Bayesian methods in matrix decomposition was first introduced in the early 2000s, with the aim of addressing the limitations of traditional matrix factorization techniques, e.g., limited explanatory and predictive performance.

The sole aim of this book is to provide a self-contained introduction to concepts and mathematical tools in Bayesian matrix decomposition in order to seamlessly introduce matrix decomposition techniques and their applications in subsequent sections.

However, given the limited scope to provide this discussion, we clearly recognize our inability to cover all the fascinating and useful conclusions regarding Bayesian matrix decomposition, such as the separated analysis of variational inference for carrying out the optimization.
We refer the reader to literature in the field of Bayesian analysis for a more detailed introduction to the related fields.

Primarily, this book serves as a summary of the purpose and significance of key Bayesian matrix decomposition methods, including real-valued decomposition, nonnegative matrix factorization, and Bayesian interpolative decomposition.
It delves into the origin and complexity of these methods, shedding light on their practical applications. The only mathematical prerequisite is a first course in statistics and linear algebra, ensuring accessibility with rigorous proofs provided throughout the development.

CHAPTER: KEYWORDS

Bayesian inference, Gibbs sampling, Conjugate model,
Alternating least squares (ALS), Multiplicative update, Real-valued matrix decomposition, Nonnegative matrix decomposition, Interpolative decomposition, Ordinal matrix decomposition, Poisson matrix decomposition.

We would like to express gratitude towards Ulrich Paquet for providing information on Bayesian ordinal matrix decomposition and towards Gilbert Strang for discussing the proof of the CUR decomposition.
Additionally, we would like to thank Federico Poloni for his comment on the proof of alternating least squares.
The author also wishes to acknowledge the cooperation of Joerg Osterrieder, Christine P. Chai, and Xuanyu Ye towards the Bayesian approach for nonnegative matrix factorization and interpolative decomposition. Their cooperation has helped enlighten the form of many discussions in the book.

CHAPTER: NOTATION

This section provides a concise reference describing notation used throughout this
If you are unfamiliar with any of the corresponding mathematical concepts,
the book describes most of these ideas in Chapter <ref> (p. chapter_introduction).

Numbers and Arrays

$\displaystyle a$       A scalar (integer or real)

$\displaystyle \ba$     A vector

$\displaystyle \bA$     A matrix

$\displaystyle \bI_n$     Identity matrix with $n$ rows and $n$ columns

$\displaystyle \bI$       Identity matrix with dimensionality implied by context

$\displaystyle \ve_i$     Standard basis vector $[0,\dots,0,1,0,\dots,0]$ with a 1 at position $i$

$\displaystyle \text{diag}(\va)$     A square, diagonal matrix with diagonal entries given by $\va$

$\displaystyle \ra$       A scalar random variable

$\displaystyle \rva$      A vector-valued random variable

$\displaystyle \rmA$      A matrix-valued random variable

$\displaystyle \sA$     A set

$\displaystyle \varnothing$     The null set

$\displaystyle \real$     The set of real numbers

$\displaystyle \naturalset$     The set of natural numbers

$\displaystyle \complex$     The set of complex numbers

$\displaystyle \{0, 1\}$     The set containing 0 and 1

$\displaystyle \{0, 1, \dots, n \}$     The set of all integers between $0$ and $n$

$\displaystyle [a, b]$     The real interval including $a$ and $b$

$\displaystyle (a, b]$     The real interval excluding $a$ but including $b$

$\displaystyle \sA \backslash \sB$     Set subtraction, i.e., the set containing the elements of $\sA$ that are not in $\sB$

Matrix indexing

$\displaystyle \eva_i$     Element $i$ of vector $\va$, with indexing starting at 1

$\displaystyle \ba_{-i}$     All elements of vector $\va$ except for element $i$

$\displaystyle  a_{ij}$     Element $i, j$ of matrix $\mA$

$\displaystyle \mA_{i, :}$     Row $i$ of matrix $\mA$

$\displaystyle \mA_{:, i},\, \ba_i$     Column $i$ of matrix $\mA$

Linear Algebra Operations

$\displaystyle \bA^\top$     Transpose of matrix $\mA$

$\displaystyle \bA^+$     Moore-Penrose pseudoinverse of $\mA$

$\displaystyle \bA \circledast \bB $     Element-wise (Hadamard) product of $\mA$ and $\mB$

$\displaystyle \mathrm{det}(\bA)$     Determinant of $\bA$

$\displaystyle \mathrm{rref}(\bA)$     Reduced row echelon form of $\bA$

$\displaystyle \cspace(\bA)$     Column space of $\bA$

$\displaystyle \nspace(\bA)$     Null space of $\bA$

$\displaystyle \mathcalV$     A general subspace

$\displaystyle \rank(\bA)$     Rank of $\bA$

$\displaystyle \trace(\bA)$     Trace of $\bA$

Element-wise product|see Hadamard product
Hadamard product

$\displaystyle\frac{d y} {d x}$     Derivative of $y$ with respect to $x$
$\displaystyle \frac{\partial y} {\partial x} $     Partial derivative of $y$ with respect to $x$

$\displaystyle \nabla_{\bx} y $     Gradient of $y$ with respect to $\bx$

$\displaystyle \nabla_{\bX} y $     Matrix derivatives of $y$ with respect to $\bX$

$\displaystyle \frac{\partial f}{\partial \vx} $     Jacobian matrix $\mJ \in \R^{m\times n}$ of $f: \R^n \rightarrow \R^m$

$\displaystyle \nabla_\vx^2 f(\vx)\text{ or }\mH( f)(\vx)$     The Hessian matrix of $f$ at input point $\vx$

$\displaystyle \int f(\vx) d\vx $     Definite integral over the entire domain of $\vx$

$\displaystyle \int_\sS f(\vx) d\vx$     Definite integral with respect to $\vx$ over the set $\sS$

Jacobian matrix
Hessian matrix

Probability and Information Theory

$\displaystyle \ra \bot \rb$     The random variables $\ra$ and $\rb$ are independent

$\displaystyle \ra \bot \rb \mid \rc $     They are conditionally independent given $\rc$

$\displaystyle P(\ra)$     A probability distribution over a discrete variable

$\displaystyle p(\ra)$     A probability distribution over a continuous variable, or over
a variable whose type has not been specified

$\displaystyle \ra \sim P$     Random variable $\ra$ has distribution $P$

$\displaystyle  \E_{\rx\sim P} [ f(x) ]\text{ or } \E [f(x)]$     Expectation of $f(x)$ with respect to $P(\rx)$

$\displaystyle \Var[f(x)] $      Variance of $f(x)$ under $P(\rx)$

$\displaystyle \Cov[f(x),g(x)] $     Covariance of $f(x)$ and $g(x)$ under $P(\rx)$

$\displaystyle H(\rx) $     Shannon entropy of the random variable $\rx$

$\displaystyle \KL ( P \Vert Q ) $     Kullback-Leibler divergence of P and Q

$\displaystyle \mathcal{N} ( \vx \mid \vmu , \mSigma)$     Gaussian distribution
over $\vx$ with mean $\vmu$ and covariance $\mSigma$

Conditional independence
Kullback-Leibler divergence
Shannon entropy

$\displaystyle f: \sA \rightarrow \sB$     The function $f$ with domain $\sA$ and range $\sB$

$\displaystyle f \circ g $     Composition of the functions $f$ and $g$

$\displaystyle f(\vx ; \vtheta) $     A function of $\vx$ parametrized by $\vtheta$.
(Sometimes we write $f(\vx)$ and omit the argument $\vtheta$ to lighten notation)

$\displaystyle \log(x)$     Natural logarithm of $x$

$\displaystyle \sigma(x),\, \text{Sigmoid}(x)$     Logistic sigmoid, i.e., $\displaystyle \frac{1} {1 + \exp\{-x\}}$

$\displaystyle \zeta(x)$     Softplus, $\log(1 + \exp\{x\})$

$\displaystyle \norm{\bx}_p $     $\ell_p$ norm of $\vx$

$\displaystyle \norm{\bx}=\norm{\bx}_2 $     $\ell_2$ norm of $\vx$

$\displaystyle \norm{\bx}=\norm{\bx}_1 $     $\ell_1$ norm of $\vx$

$\displaystyle \norm{\bx}=\norm{\bx}_\infty $     $\ell_\infty$ norm of $\vx$

$\displaystyle x^+$     Positive part of $x$, i.e., $\max(0,x)$

$\displaystyle u(x)$     Step function with value 1 when $x\geq0$ and value 0 otherwise

$\displaystyle \indicator\{\mathrm{condition}\}$     is 1 if the condition is true, 0 otherwise

Sometimes we use a function $f$ whose argument is a scalar but apply
it to a vector, matrix: $f(\vx)$, $f(\mX)$.
This denotes the application of $f$ to the
array element-wise. For example, if $\bC = \sigma(\bX)$, then $c_{ij} = \sigma(x_{ij})$
for all valid values of $i$ and  $j$.

PD     Positive definite

PSD     Positive semidefinite

MCMC     Markov chain Monte Carlo

i.i.d.     Independently and identically distributed

p.d.f.     Probability density function

p.m.f.     Probability mass function

OLS     Ordinary least squares

NG     Normal-Gamma distribution

NIG     Normal-inverse-Gamma  distribution

NIX     Normal-inverse-Chi-squared distribution

TN     Truncated-normal distribution

GTN     General-truncated-normal distribution

RN     Rectified-normal distribution

IW     Inverse-Wishart distribution

NIW     Normal-inverse-Wishart distribution

ALS     Alternating least squares

GD     Gradient descent

SGD     Stochastic gradient descent

MU     Multiplicative update

MSE     Mean squared error

NMF     Nonnegative matrix factorization

ID     Interpolative decomposition

IID     Intervened interpolative decomposition

BID     Bayesian interpolative decomposition

CHAPTER: MONTE CARLO METHODS

Markov chain Monte Carlo

This book focuses on Markov chain Monte Carlo (MCMC) methods for probabilistic inference, aiming to draw conclusions from a probabilistic model.
This chapter provides a survey of the mathematical details related to probabilistic inference, emphasizing aspects that will establish the foundation for the subsequent chapters.

§ THE BAYESIAN APPROACH

In the  past decade, the Bayesian approach has been used in a wide variety of problems in data analysis, including economic forecasting, medical imaging, and population studies [Besag, 1986, Hill, 1994, Marseille et al., 1996].
In modern statistics, Bayesian approaches
have become increasingly more important and widely used.
The inception of this idea is credited to Thomas Bayes, who came up with this idea but died before publishing it. Fortunately, his friend Richard Price carried on his work and published it in 1764.
The same concept was later independently discovered by Laplace at the end of the 18-th century.
In this section, we describe the basic ideas about the Bayesian approach and use the Beta-Bernoulli model and Bayesian linear model as an appetizer of the pros and prior information of Bayesian models.

Bayes' theorem
Bayesian modeling and statistics are fundamentally driven by  Bayes' theorem. Formally, the theorem is expressed as follows.

Let $\sS$ be a sample space and let $B_1, B_2, \ldots, B_K$ be a partition of $\sA$ such that (1). $\cup_k B_k=\sS$ and (2). $B_i \cap B_j=\varnothing$ for all $i\neq j$.
Let further $A$ be any event. Then it follows that

P(B_k \mid A) = \frac{P(A \mid B_k)P(B_k)}{P(A)} = \frac{P(A\mid B_k)P(B_k)}{\sum_{i=1}^{K}P(A\mid B_i)P(B_i)}.

In Bayesian modeling and statistics, Bayes' theorem provides us with a straightforward rule for updating probabilities when new information, such as observed data, becomes available.
This allows us to adjust our prior beliefs regarding parameters of interest.

To be more specific, let $\mathcalX=\mathcalX (\bx_{1:N})= \{\bx_1, \bx_2, \ldots, \bx_N\}$ be the observations of $N$ data points, and suppose they are independent and identically distributed (i.i.d.) with the probability parameterized by $\btheta$. Note that the parameters $\btheta$ might include  hidden variables, such as   latent variables in a mixture model indicating the cluster to which a data point belongs.

The idea of the Bayesian approach involves assuming a prior probability distribution for $\btheta$ with hyper-parameters $\balpha$ (i.e., $p(\btheta\mid \balpha)$, also known as the probability of the model).
This  distribution represents the plausibility of each possible value of $\btheta$ before observing the data, and it captures our prior uncertainty regarding $\btheta$.
The joint distribution of $\btheta$ and $\mathcalX$ is given by

p(\btheta, \mathcalX) = p(\btheta \mid \balpha) p(\mathcalX \mid \btheta).

And we can integrate out $\btheta$ to obtain the marginal distribution of $\mathcalX$,

p(\mathcalX)= \int_{\btheta} p(\btheta \mid \balpha)p(\mathcalX \mid \btheta) d\btheta.

Then, to make inferences about $\btheta$, one simply considers the conditional distribution of $\btheta$ given the observed data. This is referred to as the posterior distribution, since it represents the plausibility of each possible value of $\btheta$ after seeing the data.
The posterior distribution is the solution space for given problems, since it measures the probability of the present model in light of the data.
Mathematically, this relationship is expressed via Bayes’ theorem,

\begin{equation}\label{equation:posterior_abstract_for_mcmc}
\begin{aligned}
p(\btheta \mid  \mathcalX, \balpha)
&= \frac{p(\mathcalX \mid  \btheta ) p(\btheta \mid  \balpha)}{p(\mathcalX \mid \balpha)} \\
&= \frac{p(\mathcalX \mid  \btheta ) p(\btheta \mid  \balpha)}{\int_{\btheta}  p(\mathcalX, \btheta \mid \balpha) }
= \frac{p(\mathcalX \mid \btheta ) p(\btheta \mid \balpha)}{\int_{\btheta}  p(\mathcalX \mid  \btheta ) p(\btheta \mid  \balpha) }  \propto p(\mathcalX \mid \btheta ) p(\btheta \mid \balpha),
\end{aligned}
\end{equation}

where “$\propto$" means “proportional to" (see Problem <ref>.<ref>), $\mathcalX$ is the observed data set, and $p(\mathcalX \mid \balpha )$ can be disregarded in this case since it acts as a scaling parameter (and we shall see the MCMC algorithm only needs relative probabilities). In other words, we say the posterior is proportional to the product of the likelihood and the prior.
This means that the relative probability at a point in the solution space is determined completely by the likelihood,
which is easily determined by comparing the model to the data, and the prior, which is the probability of the model independent of the data.
The prior encodes any a priori knowledge of the solution irrespective  of the observed data.
For example, a prior for a system reducing over-clustering might assign a higher probability to a larger cluster than to a small cluster [Lu, 2021].

More generally, the Bayesian approach—in a nutshell—is to assume a prior distribution for any unknowns ($\btheta$ in our case), and then just follow the rules of probability to answer any questions of interest. For example, when we find the parameter based on the maximum posterior probability of $\btheta$, we turn to the maximum a posteriori (MAP) estimator.

Frequentists V.S. Bayesian.
The frequentist approach to statistics, developed by Neyman, evaluates statistical procedures based on a probability distribution over all possible data sets.
To be more specific, frequentists consider the parameter vector $\btheta$ to be fixed (albeit unknown), while introducing uncertainty over possible data sets $\mathcalX$. Frequentist methods are often considered more objective as they avoid incorporating subjective prior information.
In contrast, Bayesian methods allow for the incorporation of prior beliefs. The Bayesian approach treats the data set $\mathcalX$ as given, while introducing uncertainty over $\btheta$.
However, statisticians nowadays tend to move comfortably between these approaches and popular statistical procedures often combine both of them,  incorporating Bayesian methods for certain aspects of the analysis while using frequentist methods for others.
For instance, empirical Bayesian methods  have a Bayesian spirit but are not strictly Bayesian; their analysis is frequently frequentist [Haugh, 2021].

§ APPROXIMATE INFERENCE

In this book, our emphasis is on approximate probabilistic inference methods.
In certain cases, it is computationally feasible to compute the posterior exactly. For example, exponential families with conjugate priors often enable analytical solutions.
Although exact inference methods exist and they are precise and useful for certain classes of problems, exact inference methods in complicated models are usually intractable, because these methods typically depend on integrals, summations, or intermediate representations that grow large and become impractical as the state space grows too large so as to make the computation inefficient. For instance, we may use conjugate priors in a Gaussian mixture model. However, the model is hierarchical and is too complicated to compute the exact posterior.
Consequently, in such cases,  approximate probabilistic inference methods are rather useful and necessary.

Generally, variational methods and Monte Carlo methods  are two primary classes of approximate inference. We here give a brief comparison of the two methods [Bonawitz, 2008].
In variational inference methods, we first approximate the full model with a simpler model, in which the inference questions are tractable. Then, the parameters of this simplified model are calculated by some methods (e.g., using optimization methods) to minimize a measure of the dissimilarity between the original model and the simplified version; this calculation usually performs deterministically because   optimization methods are used. Finally, certain queries can be calculated and executed in the simplified model. In other words, the main idea behind variational methods is to pick a family of distributions over the parameters with its own variational parameters—$q(\btheta \mid  \boldsymbol\nu)$, where $\boldsymbol\nu$ represents the variable containing  variational parameters. Then, find the setting of the parameters that makes $q$ close to the posterior of interest. For a detailed example,  refer to Ma et al., 2014. The main advantage of variational methods is deterministic; however, the corresponding results are in the form of a lower bound of the desired quantity, and the tightness of this bound depends on the degree to which the simplified distribution can model the original posterior distribution. The variational inference is an important tool for Bayesian deep learning [Jordan et al., 1999, Graves, 2011, Hoffman et al., 2013, Ranganath et al., 2014, Mandt and Blei, 2014].

On the contrary, in Monte Carlo methods, we first draw a sequence of samples from the true target posterior distribution. Then certain inference questions are then answered by using this set of samples as an approximation of the target distribution itself. Monte Carlo methods are guaranteed to converge—if you want a more accurate answer, you just need to run the inference for longer; in the limit of running the Monte Carlo algorithm forever, the approximation results from the samples converge to the target distribution (see Section <ref>).

§ MONTE CARLO (MC) METHODS

In Monte Carlo methods, we first draw $N$ samples $\btheta_1, \btheta_2, \ldots, \btheta_N$ from the posterior distribution $p(\btheta \mid  \mathcalX, \balpha)$ in Equation (<ref>).
Subsequently, we approximate the distribution of interest by

\begin{equation}
p(\btheta \mid  -) \approx \overset{\sim}{p}(\btheta \mid  -) = \frac{1}{N} \sum_{n=1}^N \delta_{\btheta_n}(\btheta),
\end{equation}

where $\delta_{\btheta_i}(\btheta)$ is the Dirac delta function [The Dirac delta function $\delta_{\bx_0}
(\bx)$ has the property that it is nonzero and equals  1 only at $\bx = \bx_0$.]. As the number of samples increases, the approximation
(almost surely) converges to the true target distribution, i.e., $\overset{\sim}{p}(\btheta) \overset{\overset{a.s.}{N\rightarrow \infty} }{\longrightarrow} p(\btheta)$.

These kinds of sampling-based methods are extensively used in modern statistics due to their ease of use
and the generality with which they can be applied. The fundamental problem solved by
these methods is the approximation of expectations such as

\begin{equation}
\E {h(\bTheta)} = \int_{\btheta} h(\btheta) p(\btheta) d \btheta,
\end{equation}

in the case of a continuous random variable $\bTheta$ with probability density function (p.d.f.) $p$. Or

\begin{equation}
\E {h(\bTheta)} = \sum_{\btheta} h(\btheta) p(\btheta),
\end{equation}

in the case of a discrete random variable $\bTheta$ with probability mass function (p.m.f.) $p$. The general principle at work is that such expectations can be approximated by

\begin{equation}
\E {h(\bTheta)} \approx \sum_{n=1}^N h(\btheta_n).
\end{equation}

If it were generally easy to draw samples directly from $p(\btheta \mid  \mathcalX, \balpha)$, the Monte Carlo
story would end here. Unfortunately, this is usually intractable. We can consider the posterior form $p(\btheta \mid  \mathcalX, \balpha) = \frac{p(\mathcalX \mid  \btheta ) p(\btheta \mid  \balpha)}{p(\mathcalX \mid  \balpha)}$, where in many problems $p(\mathcalX \mid  \btheta ) p(\btheta \mid  \balpha)$ can be computed easily, but computing $p(\mathcalX \mid  \balpha)$ is challenging due to integrals, summations, etc. In such cases, Markov chain Monte Carlo is especially useful.

§.§ Markov Chain Monte Carlo (MCMC)

Markov chain Monte Carlo (MCMC) algorithms, also called MCMC samplers, are numerical approximation algorithms.
Since MCMC algorithms directly sample the solution space, uncertainty estimates
are determined simultaneously with a “best" solution.
Additionally, provided that the data support them, multiple solutions are feasible.
Intuitively, it is a stochastic hill-climbing approach to inference, operating over the complete data set. This inference method is designed to spend most of the computational efforts to sample points from the high probability regions of true target posterior distribution $p(\btheta \mid  \mathcalX, \balpha)$ [Andrieu et al., 2003, Bonawitz, 2008, Hoff, 2009, Geyer, 2011].

In this sampler, a Markov chain stochastic walk is taken through the state space $\bTheta$ such that the probability of being in a particular state $\btheta_t$ at any point in the walk is $p(\btheta_t \mid  \mathcalX, \balpha)$. Therefore, samples from the true posterior distribution $p(\btheta \mid \mathcalX, \balpha)$ can be approximated by recording the samples (states) visited by the stochastic walk and some other post-processing methods such as thinning. The stochastic walk is a Markov chain, i.e. the choice of state at time $t + 1$ depends only on its previous state—the state at time $t$.  Formally, if $\btheta_t$ is the state of the chain at time $t$, then it follows that $p(\btheta_{t+1}\mid \btheta_1, \btheta_2,\ldots, \btheta_t) = p(\btheta_{t+1}\mid  \btheta_t)$. That is, Markov chains are history-free. This history-free property of Markov chains provides two main advantages:

* Due to this history-free property, the Markov chain Monte Carlo methods  can be run for an unlimited number of iterations without consuming additional memory space.

* The history-free property also indicates that the MCMC stochastic walk can be completely characterized by $p(\btheta_{t+1}\mid  \btheta_t)$, known as the transition kernel.

We then focus on the discussion of the transition kernel. The transition kernel $\bK$ can also be expressed as a linear transform.
If $p_t = p_t (\btheta)$ is a row vector that encodes the probability of the walk being in state $\btheta$ at time $t$, then $p_{t+1} = p_t \bK$. If the stochastic walk starts from state $\btheta_0$, then the distribution from this initial state is the delta distribution $p_0 = \delta_{\btheta_0} (\btheta)$, and the state distribution for the chain after step $t$ is $p_t = p_0\bK^t$. We can easily find that the key to Markov chain Monte Carlo lies in choosing a kernel $\bK$ such that $\underset{t\rightarrow \infty}{\mathrm{lim}} p_t = p(\btheta \mid  \mathcalX, \balpha)$, independent of the choice of $\btheta_0$. Kernels exhibiting this property are said to converge to an equilibrium distribution $p_{eq} = p(\btheta \mid  \mathcalX)$. Convergence is guaranteed if both of the following criteria are met (see Bonawitz, 2008):

* $p_{eq}$ is an invariant (or stationary) distribution for $\bK$. A distribution $p_{inv}$ is considered an invariant distribution for $\bK$ if $p_{inv} = p_{inv} \bK$;

* $\bK$ is ergodic. A kernel is called ergodic if it is irreducible (meaning that any state can be reached from any other state) and aperiodic (indicating that the stochastic walk never gets stuck in cycles).

There are a large number of MCMC algorithms, too extensive for comprehensive coverage in this context. Popular families include Gibbs sampling, Metropolis-Hastings (MH), slice sampling, Hamiltonian Monte Carlo, adaptive rejection sampling, and many others. Though the name is potentially  misleading, Metropolis-within-Gibbs (MWG) was initially developed by [Metropolis et al., 1953], and MH subsequently emerged as a generalization of MWG [Hastings, 1970]. All MCMC algorithms are recognized as special instances of the MH algorithm. Regardless of the algorithm, the goal of Bayesian inference is to maximize the unnormalized joint posterior distribution and collect samples of the target distributions, which are marginal posterior distributions, later to be used for inference queries.

The most generalizable MCMC algorithm is the Metropolis-Hastings (MH) generalization [Metropolis et al., 1953, Hastings, 1970] of the MWG algorithm. The MH algorithm extended MWG to include asymmetric proposal distributions.

In this method, it converts an arbitrary proposal kernel $q(\btheta_{\star} \mid  \btheta_t  )$ into a transition kernel with the desired invariant distribution $p_{eq}(\btheta)$.
In order to generate a sample from a MH transition kernel, the process involves drawing a proposal $\btheta_{\star} \sim q( \btheta_{\star} \mid  \btheta_t )$ and subsequently evaluating the MH acceptance probability by

\begin{equation}
P[A(\btheta_{\star}\mid  \btheta_{t} )]  = \min \left(1, \frac{p(\btheta_{\star} \mid  \balpha)q(\btheta_{t} \mid  \btheta_{\star}  )}{p(\btheta_{t} \mid  \balpha)q(\btheta_{\star} \mid  \btheta_{t}  )} \right),
\end{equation}

with probability $P[A(\btheta_{\star} \mid  \btheta_{t}  )] $ the proposal is accepted and we set $\btheta_{t+1} =  \btheta_{\star}$; otherwise the
proposal is rejected and we set $\btheta_{t+1} =  \btheta_{t}$. That is

\begin{equation}
\btheta_{t+1} =\left\{
\begin{array}{ll}
\btheta_\star,  \text{ with probability } P[A(\btheta_{\star} \mid  \btheta_{t}  )]; \\
\btheta_{t},  \text{ with probability } 1 - P[A( \btheta_{\star} \mid  \btheta_{t}  )] .
\end{array}
\right.
\end{equation}

Intuitively, we may find that the $\frac{p(\btheta_{\star} \mid  \balpha)}{p(\btheta_{t} \mid  \balpha)}$ term tends to accept moves that lead to higher probability regions of the state space, while  the $\frac{q( \btheta_{t} \mid  \btheta_{\star} )}{q(\btheta_{\star} \mid  \btheta_{t}  )}$ term tends to accept moves that are easy to undo.

Since in MH, we only evaluate $p(\btheta)$ as a part of the ratio $\frac{p(\btheta_{\star} \mid  \balpha)}{p(\btheta_{t} \mid  \balpha)}$, we do not need to compute $p(\mathcalX \mid  \balpha)$ as mentioned in Section <ref>.

The key in MH is the proposal kernel $q(\btheta_{\star} \mid  \btheta_t  )$. However, the transition kernel is not $q(\btheta_{\star} \mid  \btheta_t  )$. Informally, the kernel $K(\btheta_{t+1} \mid  \btheta_t  )$ in MH is

p (\btheta_{t+1} \mid  \mathrm{accept}) P[\mathrm{accept}] + p (\btheta_{t+1}\mid \mathrm{reject}) P[\mathrm{reject}].

While Tierney, 1998 introduced that the precise transition kernel is

\begin{equation}
\begin{aligned}
K(\btheta_t \rightarrow \btheta_{t+1})
&= p(\btheta_{t+1} \mid  \btheta_t) \\
&= q(\btheta_{t+1} \mid  \btheta_t  ) A(\btheta_{t+1} \mid  \btheta_t  ) + \delta_{\btheta_t}(\btheta_{t+1}) \int_{\btheta_\star} q(\btheta_{\star} \mid  \btheta_t  ) (1 - A(\btheta_{\star} \mid  \btheta_{t}  )) .
\end{aligned}
\end{equation}

§.§ MC V.S. MCMC

As shown in previous sections, the purpose of Monte Carlo or Markov chain Monte Carlo approximation is to obtain a sequence of parameter values $\{\btheta^{(1)}, \ldots, \btheta^{(N)}\}$ such that

\begin{equation}
\frac{1}{N} \sum_{n=1}^N h(\btheta^{(n)}) \approx \int_{\btheta} h(\btheta) p(\btheta) d\btheta,
\end{equation}

for any functions $h$ of interest in the case of continuous random variables. In other words, we want the empirical average of $\{h(\btheta^{(1)}), \ldots,$ $h(\btheta^{(N)})\}$ to approximate the expected value of $h(\btheta)$ under a target probability distribution $p(\btheta)$. In order for this to be a good approximation for a wide range of functions $h$, we require the empirical distribution of the simulated sequence $\{\btheta^{(1)}, \ldots, \btheta^{(N)}\}$ to look like the target distribution $p(\btheta)$. MC and MCMC are two ways of generating such a sequence. MC simulation, in which we generate independent samples from the target distribution, is in some sense the “true situation." Independent MC samples automatically create a sequence that is representative of $p(\btheta)$, which means the probability that $\btheta^{(n)} \in A$ for any set $A$ is

\begin{equation}
\int_{A} p(\btheta) d\btheta,
\end{equation}

where $n \in \{1, \ldots, N\}$. However, this is not true for MCMC samples, in which case all we are sure of is that

\begin{equation}
\lim_{n \rightarrow \infty} Pr(\theta^{(n)} \in A) = \int_A p(\btheta) d\btheta.
\end{equation}

§.§ Gibbs Sampler

Gibbs sampling was introduced by Turchin [Turchin, 1971] and later by brothers Geman and Geman [Geman and Geman, 1984] in the context of image restoration. The Geman brothers named the algorithm after the physicist J. W. Gibbs, some eight decades after his death, in reference to an analogy between the sampling algorithm and statistical physics.

Gibbs sampling is applicable when the joint distribution is not explicitly known  or is challenging  to sample from directly.
Instead, it relies on the known and easily samplable conditional distribution of each variable.
A Gibbs sampler sequentially generates a draw from the distribution of each parameter or variable in turn, conditional on the current values of the other parameters or variables. Therefore, a Gibbs sampler is a componentwise algorithm. In our example, given some data $\mathcalX$ and a probability distribution $p(\btheta \mid  \mathcalX, \balpha)$ parameterized by $\btheta = \{\theta_1, \theta_2, \ldots, \theta_p\}$. We can successively draw samples from the distribution by sampling from

\begin{equation}\label{equation:gibbs_thetai_t}
\theta_i^{(t)} \sim p(\theta_i \mid  \btheta_{-i}^{(t-1)}, \mathcalX, \balpha),
\end{equation}

where $\btheta_{-i}^{(t-1)}$ is all current values of $\btheta$ in the $(t-1)$-th iteration except for $\theta_i$. If we sample sufficiently long, these $\theta_i$ values will be random samples from the distribution $p$.
If we sample new values in turn   for each parameter $\theta_i$ from Equation (<ref>), we will eventually converge to draws from the posterior $p(\btheta \mid  \mathcalX, \balpha)$.
When doing this Gibbs sampler, we also have to discard the first $k$ draws since it takes a while to converge (i.e., burn-in).
And because the consecutive draws are correlated, we only use  every $j$-th sample (i.e., thinning).

In deriving a Gibbs sampler, it is often helpful to observe that

\begin{equation}
p(\theta_i  \mid   \btheta_{- i}, \mathcalX)
= \frac{
p(\theta_1, \theta_2, \ldots,\theta_p, \mathcalX)
p(\btheta_{- i}, \mathcalX)
} \propto p(\theta_1, \theta_2, \ldots,\theta_p, \mathcalX).
\end{equation}

That is, the conditional distribution is proportional to the joint distribution. We will get a lot of benefits from this simple observation by dropping constant terms from the joint distribution (relative to the parameters we are conditioned on).

Shortly, as a simplified example, given a joint probability distribution $p(\theta_1,\theta_2\mid \mathcalX)$, a Gibbs sampler would draw $p(\theta_1\mid \theta_2,\mathcalX)$ , then $p(\theta_2\mid \theta_1,\mathcalX)$ iteratively.
The procedure generates a sequence of realization for random variables $\theta_1$ and $\theta_2$:

\begin{equation}
(\theta_1^0, \theta_2^0),\,\,\,
(\theta_1^1, \theta_2^1),
\,\,\,(\theta_1^2, \theta_2^2),\,\,\, \cdots, \nonumber
\end{equation}

which converges to the joint distribution $p(\theta_1, \theta_2)$. Further details about Gibbs sampling can be found in Turchin, 1971, Geman and Geman, 1984, Hoff, 2009, Gelman et al., 2013.

§.§ Adaptive Rejection Sampling (ARS)

The purpose of the adaptive rejection sampling algorithm is to provide a relatively efficient way to sample from a distribution from the large class of log-concave densities [Gilks and Wild, 1992, Wild and Gilks, 1993]. We provide a brief overview of the algorithm here, with additional details available in Gilks and Wild, 1992 and Wild and Gilks, 1993.

Rejection sampling. Figure is  from Michael I. Jordan's lecture notes.

§.§.§ Rejection Sampling

In rejection sampling, the objective is to sample from a target probability density function $p(x)$, given that we can sample from a probability density function $q(x)$ easily.
Although the target density $p(x)$ is unknown, the approach relies on establishing an envelop by considering $M \times q(x)$ such that it covers $p(x)$ for some $M>1$, as illustrated in Figure <ref>.
This is expressed as:

\begin{equation}
\frac{p(x)}{q(x)} < M, \text{ for all $x$.}
\end{equation}

Subsequently, when sampling $x_i$ from $q(x)$, and if $y_i=u \times M\times q(x_i)$ lies below the region under $p(x)$ for some $u \sim \mathrm{Uniform}(0,1)$, then we accept $x_i$; otherwise, it is rejected.

In essence, the method involves sampling $x_i$ from a distribution and making an acceptance or rejection decision based on the comparison with the envelope.

§.§.§ Adaptive Rejection Sampling

Adaptive rejection sampling. Figure is from Michael I. Jordan's lecture notes.

The adaptive rejection sampling method, tailored for log-concave densities, is an extension of rejection sampling.
The basic idea involves dynamically constructing an upper envelope (the upper bound on $p(x)$), serving as an adaptive replacement for  $M\times q(x)$ in rejection sampling.

As shown in Figure <ref>, the logarithm of the density, $\log p(x)$, is considered. We then sample $x_i$ from the upper envelope, and the sample is  either accepted or rejected akin to rejection sampling.
In case of rejection, a tangent is drawn passing through $x = x_i$ and $y = \log(p)$; and the tangent is used to reduce the upper envelope to decrease the number of rejected samples. The intersections of these tangent planes enable the formation of an envelope adaptively. To sample from the upper envelope, we need to transform from log space by exponentiating and using properties of the exponential distribution.

§ BAYESIAN APPETIZERS

This section delves into semi-conjugate priors with the Gibbs sampler and fully conjugate priors without approximate inference, offering an in-depth exploration  of Bayesian approaches.
Readers who already have a basic knowledge of Bayesian inference may choose to skip this section.

§.§ Beta-Bernoulli Model

We formally introduce a Beta-Bernoulli model to illustrate how the Bayesian approach works.
The Bernoulli distribution serves to model binary outcomes, i.e., outputting two possible values. The likelihood under this model corresponds to the probability mass function of the Bernoulli distribution with parameter $\theta$:

\begin{equation}\label{equation:bernoulli_distribution}
\bernoulli(x\mid \theta) = p(x\mid \theta) = \theta^x (1-\theta)^{1-x} \indicator(x\in \{0,1\}). \nonumber
\end{equation}

That is,

\bernoulli(x\mid \theta)=p(x\mid \theta)=\left\{
\begin{aligned}
&1-\theta ,& \mathrm{\,\,if\,\,} x = 0;  \\
&\theta , &\mathrm{\,\,if\,\,} x =1,
\end{aligned}
\right.

where $\theta$ is the probability of outputting 1, and  $1-\theta$ is the probability of outputting 0.
The mean of the Bernoulli distribution is $\theta$.
Suppose $\mathcalX=\{x_1, x_2, \ldots , x_n\}$ are drawn i.i.d. from the Bernoulli distribution $\bernoulli( \theta)$. Then the likelihood under the  Bernoulli distribution with parameter $\theta$ is given by

\begin{aligned}
\text{likelihood} =  p(\mathcalX \mid \theta) &= \theta^{\sum x_i} (1-\theta)^{n-\sum x_i},
\end{aligned}

which is a distribution on $\mathcalX$ and is called the likelihood function on $\mathcalX$.

And we will see the prior under this model is the probability density function of a Beta distribution:

\begin{equation}
\mathrm{prior} = \betadist(\theta\mid a, b)=	p(\theta\mid a, b) =\frac{1}{B(a,b)} \theta^{a-1}(1-\theta)^{b-1} \indicator(0\leq \theta\leq 1), \nonumber
\end{equation}

where $B(a,b)$ denotes the Euler's beta function and it can be simply regarded as a constant normalization term, and $\indicator(a\leq x\leq b)$ is a step function that has a value of 1 when $a\leq x\leq b$  and 0 otherwise (when $x<a$ or $a>b$). Figure <ref> compares different parameters for the Beta distribution.
Specifically, when $a=b=1$, the beta distribution reduces to a uniform distribution in the support of $[0,1]$.

Beta probability density functions for different values of the parameters $a$ and $b$.
When $a=b=1$, the beta distribution reduces to a uniform distribution in the support of $[0,1]$.
The mean, variance, and mode  of the Beta distribution are $\Exp[\rx]=\frac{a}{a+b}$,
$\Var[\rx]=\frac{ab}{(a+b+1)(a+b)^2}$, and $\mathrm{mode}[\rx]=\frac{a-1}{(a-1)+(b-1)}$ if $a>1, b>1$, respectively.

A Beta prior is assigned to the parameter  $\theta$ of the Bernoulli distribution. The posterior is proportional to the product of the likelihood and prior densities, given by

\begin{equation}
\begin{aligned}
\mathrm{posterior} = p(\theta\mid \mathcalX)
&\propto p(\mathcalX \mid \theta) p(\theta\mid a,b) \\
&=\theta^{\sum x_i} (1-\theta)^{n-\sum x_i} \times \frac{1}{B(a,b)} \theta^{a-1}(1-\theta)^{b-1}\cdot \indicator(0\leq\theta\leq1) \\
&\propto \theta^{a+\sum x_i-1}(1-\theta)^{b+n-\sum x_i-1}\cdot \indicator(0\leq\theta\leq 1) \\
&\propto \betadist\left(\theta \,\bigg|\, a+\sum_{i=1}^n x_i, b+n-\sum_{i=1}^n x_i\right). \nonumber
\end{aligned}
\end{equation}

We find that the posterior distribution shares the same form as the prior distribution; both stem from the Beta distribution, albeit with distinct parameters.  When this happens, we call the prior a conjugate prior. A conjugate prior exhibits an advantageous structure, facilitating straightforward computations of the posterior probability density function, its derivatives, and the sampling process from the posterior.
The advantage of using conjugate priors lies in their ability to preserve the mathematical form of the prior throughout the updating process. This results in closed-form expressions for the posterior distribution, avoiding the need for complex numerical methods.

Upon comparing the prior and posterior formulations, it becomes evident that the hyper-parameter $a$ represents the prior count of occurrences of $1$'s in the output, while $b$ corresponds to the prior count of 0's in the output.
And the sum $a+b$ aggregates  the prior information about the sample size. Uninformative prior parameters can  then be characterized by $a=b=1$, i.e., a uniform distribution.

From this Beta-Bernoulli model example, like the maximum likelihood estimator and method of moment (MoM, i.e., using the moment information to calculate the model parameter), Bayesian modeling also serves as a type of point estimator. However, Bayesian models output a probability of the parameter of interest, denoted as $p(\theta \mid \mathcalX)$ in the example.

When predicting outcomes for new incoming data, instead of providing a direct model $p(x_{n+1} \mid  \theta)$, the approach involves integration:

\begin{equation}
p(x_{n+1} \mid  \mathcalX) = \int p(x_{n+1} \mid  \theta) p(\theta \mid  \mathcalX) d\theta.\nonumber
\end{equation}

In other words, $x_{n+1}$ depends on $\mathcalX$. The observed data $\mathcalX$ provide insight into $\theta$, which in turn provides information about $x_{n+1}$ (i.e., $\mathcalX \rightarrow \theta \rightarrow x_{n+1}$).

[Amount of Data Matters]
Suppose we have three observations for the success in the Bernoulli distribution:

* 10 out of 10 are observed to be success (1's);

* 48 out of 50 are observed to be success (1's);

* 186 out of 200 are observed to be success (1's).

So, what is the probability of success in the Bernoulli model? The normal answer to case 1, 2, and 3 are 100%, 96%, and 93%, respectively. But an observation of 10 inputs is rather a small amount of data, and the statistical significance may be limited due to potential noise.

To address this, we introduce a $Beta(1,1)$ prior over the Bernoulli distribution. The posterior probability of success for each case would be $\frac{11}{12}=91.6\%$, $\frac{49}{52}=94.2\%$, and $\frac{187}{202}=92.6\%$, respectively. Notably, case 1 now exhibits a lower probability of success compared to case 2.

A Bayesian perspective on the problem naturally incorporates the amount of data as well as its average. This special case shown here is also called  Laplace's rate of succession [Ollivier, 2015]. Laplace's “add-one" rule (i.e., employing a $Beta(1,1)$ prior) of succession modifies the observed frequencies in a sequence of successes and failures by adding one to the observed counts.
This modification improves prediction by avoiding zero probabilities and corresponds to a uniform Bayesian prior on the parameter.
Suppose further the prior parameters are $a=b=2$, Figure <ref> compares the prior distribution and the posterior distributions for the three cases.

The prior distribution is $\betadist(x\mid 2,2)$. The posterior distributions for the three cases in Example <ref> are $\betadist(x\mid 12,2)$, $\betadist(x\mid 50,4)$, and $\betadist(x\mid 188,16)$, respectively.

[hidealllines=,backgroundcolor=,frametitle=Why Bayes?]
The presented example underscores that Bayesian models incorporate prior information on model parameters, making them particularly useful for regularizing regression problems with limited data information. This characteristic has propelled the Bayesian approach to global recognition over the decades.

In this framework, the prior information $p(\theta)$ and likelihood function $p(x\mid \theta)$ collectively represent a rational person's belief. And then the Bayes' rule emerges as the optimal method for updating this person's beliefs about parameter $\theta$ given new information from the data [Fahrmeir et al., 2007, Hoff, 2009].

The prior information given by $p(\theta)$ might be wrong if it does not accurately represent our prior beliefs.
Nevertheless, this does not mean that the posterior $p(\theta \mid  x)$ is not useful.
As encapsulated by the famous quote “all models are wrong, but some are useful" [Box and Draper, 1987]. If the prior $p(\theta)$ approximates our beliefs, then the posterior $p(\theta \mid  x)$ serves as a valuable approximation to posterior beliefs.

§.§ Bayesian Linear Model with Zero-Mean Prior

In the linear model, given the input data matrix $\bX\in\real^{n\times p}$ and the observation vector $\by\in\real^n$, the model considers the overdetermined system $\by=\bX\bbeta$ when $n>p$, where the vector $\bbeta\in\real^p$ represents a vector of weights.
It often happens that $\by=\bX\bbeta$ has no solution since there are too many equations, i.e., the matrix $\bX$ has more rows than columns ($n>p$). Define the column space of $\bX$ by $\{\bX\bgamma: \forall \bgamma \in\real^p \}$ and denoted by $\cspace(\bX)$. Thus the absence of a solution for the linear system $\by=\bX\bbeta$ implies  that the observation vector $\by$ is outside the column space of $\bX$.  Addressing this issue involves minimizing  the  mean squared error (MSE).

Instead of directly seeking the least MSE, we extend the analysis by introducing a Gaussian noise vector  $\bepsilon\in\real^n$ such that  $\boldsymbol{y} = \boldsymbol{X}\bbeta + \bepsilon$, where $\bepsilon \sim \normal(\bzero, \sigma^2 \boldsymbol{I})$ and $\sigma^2$ is fixed (where $\normal(\ba, \bB)$ denotes a multivariate Gaussian distribution [We delay the definition in Chapter <ref>, p. chapter:conjugate_models_bmf when we discuss regular conjugate models, see Definition <ref>, p. definition:multivariate_gaussian.] with mean $\ba$ and covariance $\bB$. And a detailed analysis of this model can be found in Rasmussen, 2003, Hoff, 2009, Lu, 2022), this additive Gaussian noise assumption gives rise to the likelihood. Let $\mathcalX (\bx_{1:n})= \{\bx_1, \bx_2, \ldots, \bx_n\}$ be the observations of $n$ data points, the likelihood function under this Gaussian additive noise model is

\begin{equation}\label{equation:linear_gaussian_addi_likelihood}
\mathrm{likelihood} = \by \mid  \bX, \bbeta, \sigma^2 \sim \normal(\bX\bbeta, \sigma^2\bI).
\end{equation}

Suppose we specify a multivariate Gaussian prior  with a zero-mean and a covariance matrix $\bSigma_0$ over the weight parameter $\bbeta$,

\begin{equation}
\mathrm{prior} = 	\bbeta \sim \normal(\bzero, \bSigma_0). \nonumber
\end{equation}

Applying Bayes' theorem, “$\mathrm{posterior} \propto \mathrm{likelihood} \times \mathrm{prior} $," we obtain the posterior

\begin{equation}
\begin{aligned}
&\gap \mathrm{posterior}
= p(\bbeta\mid \by,\bX, \sigma^2) \\
&\propto p(\by\mid \bX, \bbeta, \sigma^2) \cdot  p(\bbeta \mid  \bSigma_0) \\
&=  \frac{1}{(2\pi \sigma^2)^{n/2}} \exp\left\{-\frac{1}{2\sigma^2} (\by-\bX\bbeta)^\top(\by-\bX\bbeta)\right\}
\times \frac{1}{(2\pi)^{n/2}\abs{\bSigma_0}^{1/2}}\exp\left(-\frac{1}{2} \bbeta^\top\bSigma_0^{-1}\bbeta\right) \\
&\propto \exp\left\{-\frac{1}{2} (\bbeta - \bbeta_1)^\top \bSigma_1^{-1} (\bbeta - \bbeta_1)\right\}
\propto \normal(\bbeta_1, \bSigma_1)
, \nonumber
\end{aligned}
\end{equation}

\bSigma_1 = \left(\frac{1}{\sigma^2} \bX^\top\bX + \bSigma_0^{-1}\right)^{-1}, \qquad
\bbeta_1 = \left(\frac{1}{\sigma^2}\bX^\top\bX + \bSigma_0^{-1}\right)^{-1}
\left(\frac{1}{\sigma^2}\bX^\top\by\right).

Therefore, the posterior distribution is also a multivariate Gaussian distribution (same form as the prior distribution, i.e., a conjugate prior):

\begin{equation}
\mathrm{posterior} = \bbeta\mid \by,\bX, \sigma^2 \sim \normal(\bbeta_1, \bSigma_1). \nonumber
\end{equation}

A word on the notation. Note that we use $\{\bbeta_1,\bSigma_1\}$ to denote the posterior mean vector and posterior covariance matrix in the zero-mean prior model. Similarly, the posterior mean vector and posterior covariance matrix in semi-conjugate prior and fully conjugate prior models will be denoted as $\{\bbeta_2,\bSigma_2\}$ and $\{\bbeta_3,\bSigma_3\}$, respectively, for clarity (refer to the sections below).

Connection to ordinary least squares (OLS).
In the Bayesian linear model, there is no general requirement for  $\bX$ to have full rank.
Note further that if we assume full rank for $\bX$ (i.e., $\bX^\top\bX$ is invertible when $n>p$), in the limit, when $\bSigma_0^{-1} \rightarrow \bzero$, $\bbeta_1 \rightarrow \hat{\bbeta} = (\bX^\top\bX)^{-1}\bX\by$, in which case, the maximum a posteriori (MAP) estimator from the Bayesian model reduces to the ordinary least squares (OLS) estimator. And the posterior is $\bbeta\mid \by,\bX, \sigma^2 \sim \normal (\hat{\bbeta}, \sigma^2(\boldsymbol{X}^\top \boldsymbol{X})^{-1})$, which shares a similar form with the OLS estimator $\hat{\bbeta} \sim \normal(\bbeta, \sigma^2(\boldsymbol{X}^\top \boldsymbol{X})^{-1})$ under the Gaussian disturbance (see Lu, 2022).

In the least squares approximation problem, we use $\bX\bbeta$ to approximate $\by$.
However, two issues may arise: the model can potentially overfit and $\bX$ might not have full rank.
To address these concerns, in a ridge regression model, we regularize large values of $\bbeta$, prompting  simpler models. Instead of minimizing $\norm{\by-\bX\bbeta}^2$, we minimize $\norm{\by-\bX\bbeta}^2+\lambda\norm{\bbeta}^2$, where $\lambda$ is a hyper-parameter that can be tuned accordingly, e.g., via cross-validation (CV):

\begin{equation}
\mathop{\arg\min}_{\bbeta}{(\by-\bX\bbeta)^\top(\by-\bX\bbeta) + \lambda \bbeta^\top\bbeta}. \nonumber
\end{equation}

By differentiating and setting the gradient to zero, we obtain

\begin{equation}
\hat{\bbeta}_{ridge} = \left(\bX^\top\bX + \lambda \bI\right)^{-1} \bX^\top\by, \nonumber
\end{equation}

in which case, $(\bX^\top\bX + \lambda \bI)$ is invertible even when $\bX$ does not have full rank. Further details on ridge regression will be left to the readers.

Connection to ridge regression.
We realize that when we set $\bSigma_0 = \bI$, we obtain $\bbeta_1 = \left(\bX^\top\bX + \sigma^2 \bI\right)^{-1} \bX^\top\by$ and $\bSigma_1 = \left(\frac{1}{\sigma^2}\bX^\top\bX+ \bI\right)^{-1}$. Since the posterior is $\bbeta\mid \by,\bX, \sigma^2 \sim \normal(\bbeta_1, \bSigma_1)$. The MAP estimator of $\bbeta = \bbeta_1 =  (\bX^\top\bX + \sigma^2 \bI)^{-1} \bX^\top\by$, which shares the same form as the ridge regression by letting $\sigma^2 = \lambda$.
This implies that the ridge regression estimator is a special instance of the Bayesian linear model with zero-mean prior. And the ridge regression has a nice interpretation from the Bayesian approach—finding the mode of the posterior.

An illustrative example of this Bayesian linear model is presented in Rasmussen, 2003, where the “well determined" (i.e., the distribution around the slope is more compact) slope of $\bbeta$ is almost unchanged after the posterior process while the intercept, which is more dispersed, is shrunk towards zero. This is actually a regularization effect on the parameter as the ridge regression does.

§.§ Bayesian Linear Model with Semi-Conjugate Prior

We will employ the Gamma distribution as the prior density for the inverse variance (precision) parameter of a Gaussian distribution. The rigorous definition of the Gamma distribution can be found in Chapter <ref> (Definition <ref>, p. definition:gamma-distribution) when we discuss conjugate models.

As for the reason of using the Gamma distribution as the prior for precision, we quote the description from Kruschke, 2014:
Because of its role in conjugate priors for Gaussian likelihood functions, the Gamma distribution is routinely used as a prior for precision (i.e., inverse variance). But there is no logical necessity to do so, and modern MCMC methods permit more flexible specification of priors. Indeed, because precision is less intuitive than standard deviation, it can be more useful to give standard deviation a uniform prior that spans a wide range.

Same setting as Section <ref>, we now consider the case where the variance $\sigma^2$ of the Gaussian likelihood is not fixed. Again, The likelihood function is given by:

\begin{equation}
\mathrm{likelihood} = \by \mid  \bX, \bbeta, \sigma^2 \sim \normal(\bX\bbeta, \sigma^2\bI). \nonumber
\end{equation}

We specify a non-zero-mean Gaussian prior over the weight parameter $\bbeta$,

\begin{equation}
\begin{aligned}
{\color{mylightbluetext}\mathrm{prior:\,}} &\bbeta \sim \normal({\color{mylightbluetext}\bbeta_0}, \bSigma_0) \\
&{\color{mylightbluetext}\gamma = 1/\sigma^2 \sim \gammadist(a_0, b_0)}, \nonumber
\end{aligned}
\end{equation}

where we differentiate from previous descriptions by blue text, and $\gammadist(a, b)=\frac{b^a}{\Gamma(a)} x^{a-1}\exp(-bx)$ denotes a Gamma distribution with parameters $a,b$. And the function $\Gamma(a)=\int_{0}^{\infty} a^{t-1} \exp(-a)dt$ is the Gamma function.

Step 1, conditioned on $\sigma^2$.
Then, given $\sigma^2$, by the Bayes' theorem “$\mathrm{posterior} \propto \mathrm{likelihood} \times \mathrm{prior} $," we get the conditional posterior density of $\bbeta$:

\begin{equation}
\begin{aligned}
\mathrm{posterior}&= p(\bbeta\mid \by,\bX, \sigma^2)
\propto p(\by\mid \bX, \bbeta, \sigma^2) \cdot p(\bbeta \mid  \bbeta_0, \bSigma_0) \\
&=  \frac{1}{(2\pi \sigma^2)^{n/2}} \exp\left\{-\frac{1}{2\sigma^2} (\by-\bX\bbeta)^\top(\by-\bX\bbeta)\right\} \\
&\gap \times \frac{1}{(2\pi)^{n/2}\abs{\bSigma_0}^{1/2}}\exp\left\{-\frac{1}{2} (\bbeta-\bbeta_0)^\top\bSigma_0^{-1}(\bbeta-\bbeta_0)\right\} \\
&\propto \exp\left\{-\frac{1}{2} (\bbeta - \bbeta_2)^\top \bSigma_2^{-1} (\bbeta - \bbeta_2)\right\}
\propto \normal(\bbeta_2, \bSigma_2)
, \nonumber
\end{aligned}
\end{equation}

where the parameters are

\begin{aligned}
\bSigma_2 &= \left(\frac{1}{\sigma^2} \bX^\top\bX + \bSigma_0^{-1}\right)^{-1},\\
\bbeta_2 &= \bSigma_2 (\bSigma_0^{-1}\bbeta_0+\frac{1}{\sigma^2}\bX^\top\by)
= \left(\frac{1}{\sigma^2}\bX^\top\bX + \bSigma_0^{-1}\right)^{-1}
\left(\textcolor{mylightbluetext}{\bSigma_0^{-1}\bbeta_0}+\frac{1}{\sigma^2}\bX^\top\by\right).
\end{aligned}

Therefore, the conditional posterior follows from a Gaussian distribution:

\begin{equation}
\mathrm{posterior} = \bbeta\mid \by,\bX, \sigma^2 \sim \normal(\bbeta_2, \bSigma_2). \nonumber
\end{equation}

Connection to the zero-mean prior model.
We highlight the connection between the zero-mean prior model and the semi-conjugate prior model as follows:

* We note that $\bbeta_1$ in Section <ref> is a special case of $\bbeta_2$ when $\bbeta_0=\bzero$.

* And if we assume further $\bX$ has full rank, when $\bSigma_0^{-1} \rightarrow \bzero$, we have $\bbeta_2 \rightarrow \hat{\bbeta} = (\bX^\top\bX)^{-1}\bX\by$, which reduces to  the OLS estimator.

* When $\sigma^2 \rightarrow \infty$, $\bbeta_2$ is approximately approaching  $\bbeta_0$, the prior expectation of parameter.
However, in the zero-mean prior model, $\sigma^2 \rightarrow \infty$ will make $\bbeta_1$ approach  $\bzero$.

* Weighted average: we reformulate $\bbeta_2$ by

\begin{equation}
\begin{aligned}
\bbeta_2 &=  \left(\frac{1}{\sigma^2}\bX^\top\bX + \bSigma_0^{-1}\right)^{-1}
\left(\bSigma_0^{-1}\bbeta_0+\frac{1}{\sigma^2}\bX^\top\by\right) \\
&= \left(\frac{1}{\sigma^2}\bX^\top\bX + \bSigma_0^{-1}\right)^{-1} \bSigma_0^{-1}\bbeta_0 + \left(\frac{1}{\sigma^2}\bX^\top\bX + \bSigma_0^{-1}\right)^{-1} \frac{\bX^\top\bX}{\sigma^2} (\bX^\top\bX)^{-1}\bX^\top\by \\
&=(\bI-\bA)\bbeta_0 + \bA \hat{\bbeta}, \nonumber
\end{aligned}
\end{equation}

where $\hat{\bbeta}=(\bX^\top\bX)^{-1}\bX^\top\by$ is the OLS estimator of $\bbeta$, and $\bA=(\frac{1}{\sigma^2}\bX^\top\bX + \bSigma_0^{-1})^{-1} \frac{\bX^\top\bX}{\sigma^2}$.
We observe that the posterior mean of $\bbeta$ is a weighted average of the prior mean and the OLS estimator of $\bbeta$. Thus, if we set the prior parameter $\bbeta_0 = \hat{\bbeta}$, the posterior mean of $\bbeta$ becomes precisely $\hat{\bbeta}$.

Step 2, conditioned on $\bbeta$.
Given $\bbeta$, again, by Bayes' theorem, we obtain the conditional posterior density of $\gamma = \frac{1}{\sigma^2}$:

\begin{equation}
\begin{aligned}
\mathrm{posterior}&= p(\gamma=\frac{1}{\sigma^2}\mid \by,\bX, \bbeta)
\propto p(\by\mid \bX, \bbeta, \gamma) \cdot p(\gamma \mid  a_0, b_0) \\
&=  \frac{\gamma^{n/2}}{(2\pi )^{n/2}}
\exp\left\{-\frac{\gamma}{2} (\by-\bX\bbeta)^\top(\by-\bX\bbeta)\right\} \\
&\gap \times \frac{{b_0}^{a_0}}{\Gamma(a_0)} \gamma^{a_0-1} \exp(-b_0 \gamma) \\
&\propto \gamma^{(a_0+\frac{n}{2}-1)} \exp\left\{-\gamma\left[b_0+\frac{1}{2}(\by-\bX\bbeta)^\top(\by-\bX\bbeta)\right]\right\}, \nonumber
\end{aligned}
\end{equation}

and the conditional posterior follows from a Gamma distribution:

\begin{equation}
\mathrm{posterior\,\, of\,\,} \gamma \mathrm{\,\,given\,\,} \bbeta  = \gamma\mid \by,\bX, \bbeta \sim \gammadist\left(a_0+\frac{n}{2}, \left[b_0+\frac{1}{2}(\by-\bX\bbeta)^\top(\by-\bX\bbeta)\right]\right). \nonumber
\end{equation}

Prior information on the noise/precision.
We can find an intuitive prior interpretation as follows:

* We notice that the prior mean and posterior mean of $\gamma$ are $\Exp[\gamma]=\frac{a_0}{b_0}$ and $\Exp[\gamma \mid \bbeta]=\frac{a_0 + \frac{n}{2}}{b_0 +\frac{1}{2}(\by-\bX\bbeta)^\top(\by-\bX\bbeta)}$, respectively. So the latent meaning of $2 a_0$ is the prior sample size for the noise $\sigma^2 = \frac{1}{\gamma}$.

* As we assume $\by=\bX\bbeta +\bepsilon$ where $\bepsilon \sim \normal(\bzero, \sigma^2\bI)$, then $\frac{(\by-\bX\bbeta)^\top(\by-\bX\bbeta)}{\sigma^2} \sim \chi^2(n)$ and $\Exp\left[\frac{1}{2}(\by-\bX\bbeta)^\top(\by-\bX\bbeta)\right] = \frac{n}{2}\sigma^2$ [$\chi^2(n)$ is a Chi-squared distribution with $n$ degrees of freedom. See Definition <ref> (p. definition:chisquare_distribution).]. So the latent meaning of $\frac{b_0}{a_0}$ is the prior variance of the noise.

* Some textbooks would write $\gamma \sim  \gammadist(n_0/2, n_0\sigma_0^2/2)$ to make this explicit (in which case, $n_0$ is the prior sample size, and $\sigma_0^2$ is the prior variance). But a prior in this form seems coming from nowhere at first glance.

Gibbs sampler.
By this Gibbs sampling method introduced in Section <ref>, we can construct a Gibbs sampler for the Bayesian linear model with semi-conjugate prior:

0. Set initial values to $\bbeta$ and $\gamma = \frac{1}{\sigma^2}$;

1. update $\bbeta$: $\mathrm{posterior} = \bbeta\mid \by,\bX, \gamma \sim \normal(\bbeta_2, \bSigma_2)$;

2. update $\gamma$: $\mathrm{posterior}  = \gamma\mid \by,\bX, \bbeta \sim \gammadist\left(a_0+\frac{n}{2}, [b_0+\frac{1}{2}(\by-\bX\bbeta)^\top(\by-\bX\bbeta)]\right)$.

§.§ Bayesian Linear Model with Full Conjugate Prior

Introducing a Gamma prior over the inverse variance is equivalent to imposing an inverse-Gamma prior [We again delay the definition in Definition <ref> (p. definition:inverse_gamma_distribution) when we discuss regular conjugate models.] on the variance.

Same setting as the semi-conjugate prior distribution in Section <ref>. We have the likelihood function:

\begin{equation}
\mathrm{likelihood} = \by \mid  \bX, \bbeta, \sigma^2 \sim \normal(\bX\bbeta, \sigma^2\bI). \nonumber
\end{equation}

But now we specify a joint Gaussian and inverse-Gamma prior over the weight and variance parameters by

\begin{equation}
\begin{aligned}
{\color{mylightbluetext}\mathrm{prior:\,}} &\bbeta\mid \sigma^2 \sim \normal(\bbeta_0, {\color{mylightbluetext}\sigma^2} \bSigma_0) \\
&{\color{mylightbluetext}\sigma^2 \sim \inversegammadist(a_0, b_0)}, \nonumber
\end{aligned}
\end{equation}

where again we differentiate  this description from the previous ones with blue text.
Equivalently, we can formulate the prior into a joint one, which is called the normal-inverse-Gamma (NIG) distribution:

\begin{equation}
\begin{aligned}
\mathrm{prior:\,} &\bbeta,\sigma^2 \sim \nig(\bbeta_0, \bSigma_0, a_0, b_0) =  \normal(\bbeta_0, \sigma^2 \bSigma_0)\cdot \inversegammadist(a_0, b_0) . \nonumber
\end{aligned}
\end{equation}

Once again, applying  Bayes' theorem, “$\mathrm{posterior} \propto \mathrm{likelihood} \times \mathrm{prior} $," we obtain the posterior:

\begin{equation}
\begin{aligned}
\mathrm{posterior}&= p(\bbeta,\sigma^2\mid \by,\bX)
\propto p(\by\mid \bX, \bbeta, \sigma^2)\cdot  p(\bbeta, \sigma^2 \mid  \bbeta_0, \bSigma_0, a_0, b_0) \\
&=  \frac{1}{(2\pi \sigma^2)^{n/2}} \exp\left\{-\frac{1}{2\sigma^2} (\by-\bX\bbeta)^\top(\by-\bX\bbeta)\right\} \\
&\gap \times \frac{1}{(2\pi \sigma^2)^{p/2} \abs{\bSigma_0}^{1/2}} \exp\left\{-\frac{1}{2\sigma^2} (\bbeta - \bbeta_0)^\top\bSigma_0^{-1} (\bbeta - \bbeta_0)\right\} \\
&\gap \times  \frac{{b_0}^{a_0}}{\Gamma(a_0)} \frac{1}{(\sigma^2)^{a_0+1}} \exp(-\frac{b_0}{\sigma^2}) \\
&\propto \frac{1}{(2\pi \sigma^2)^{p/2} }  \exp\left\{  \frac{1}{2\sigma^2} (\bbeta -\bbeta_3)^\top\bSigma_3^{-1}(\bbeta -\bbeta_3) \right\} \\
&\gap \times \frac{1}{(\sigma^2)^{a_0 +\frac{n}{2}+1}} \exp\left\{-\frac{1}{\sigma^2} \left[b_0+\frac{1}{2} (\by^\top\by +\bbeta_0^\top\bSigma_0^{-1}\bbeta_0 -\bbeta_3^\top\bSigma_3^{-1}\bbeta_3) \right]\right\}, \nonumber
\end{aligned}
\end{equation}

where the parameters are

\begin{aligned}
\bSigma_3 &= \left( \bX^\top\bX + \bSigma_0^{-1}\right)^{-1}, \\
\bbeta_3 &= \bSigma_3(\bX^\top\by + \bSigma_0^{-1}\bbeta_0) =
\left( \bX^\top\bX + \bSigma_0^{-1}\right)^{-1}(\bSigma_0^{-1}\bbeta_0 + \bX^\top\by).
\end{aligned}

Let $a_n = a_0 +\frac{n}{2}+1$ and $b_n=b_0+\frac{1}{2} (\by^\top\by +\bbeta_0^\top\bSigma_0^{-1}\bbeta_0 -\bbeta_3^\top\bSigma_3^{-1}\bbeta_3) $. The posterior is thus a NIG distribution:

\begin{equation}
\begin{aligned}
\mathrm{posterior}&=  \bbeta, \sigma^2 \mid  \by, \bX \sim \nig(\bbeta_3, \bSigma_3, a_n, b_n). \nonumber
\end{aligned}
\end{equation}

Connection to zero-mean prior and semi-conjugate prior models.
We highlight the connection of the fully conjugate model to the zero-mean prior and semi-conjugate prior models as follows:

* If we assume further $\bX$ has full rank, when $\bSigma_0^{-1} \rightarrow \bzero$, we have $\bbeta_3 \rightarrow \hat{\bbeta} = (\bX^\top\bX)^{-1}\bX\by$, which reduces to the OLS estimator.

* When $b_0 \rightarrow \infty$, then $\sigma^2 \rightarrow \infty$ and $\bbeta_3$ is approximately approaching $\bbeta_0$, the prior expectation of parameter. Compared to $\bbeta_2$ in Section <ref>, $\sigma^2 \rightarrow \infty$ will make $\bbeta_2$ approach to $\bbeta_0$, where $\sigma^2$ is a fixed hyper-parameter.

* Weighted average: we reformulate

\begin{equation}
\begin{aligned}
\bbeta_3 &=  \left(\bX^\top\bX + \bSigma_0^{-1}\right)^{-1}(\bSigma_0^{-1}\bbeta_0+\bX^\top\by) \\
&= \left(\bX^\top\bX + \bSigma_0^{-1}\right)^{-1} \bSigma_0^{-1}\bbeta_0 +
\left(\bX^\top\bX + \bSigma_0^{-1}\right)^{-1} (\bX^\top\bX) (\bX^\top\bX)^{-1}\bX^\top\by \\
&=(\bI-\bC)\bbeta_0 + \bC \hat{\bbeta}, \nonumber
\end{aligned}
\end{equation}

where $\hat{\bbeta}=(\bX^\top\bX)^{-1}\bX^\top\by$ is the OLS estimator of $\bbeta$, and $\bC=(\bX^\top\bX + \bSigma_0^{-1})^{-1} (\bX^\top\bX)$. We observe that the posterior mean of $\bbeta$ is a weighted average of the prior mean and the OLS estimator of $\bbeta$. Thus, if we set $\bbeta_0 = \hat{\bbeta}$, the posterior mean of $\bbeta$ will be exactly $\hat{\bbeta}$.

* From the relationship of $a_n = a_0 +\frac{n}{2}+1$, we know that $2a_0$ is the prior sample size for $\sigma^2$.

* $\bSigma_3^{-1} = \bX^\top\bX + \bSigma_0^{-1}$: the posterior precision matrix (inverse covariance matrix) is equal to data precision $\bX^\top\bX$ + prior precision.

*  Suppose both $p(x)$ and $q(x)$ are probability density functions with

p(x) \propto q(x).

Show that $p(x) = q(x)$ for all $x$.

CHAPTER: REGULAR PROBABILITY MODELS AND CONJUGACY

§ CONJUGATE PRIORS

In Section <ref> (p. sec:beta-bernoulli), we briefly discussed  conjugate priors.
Conjugate priors play a crucial role in Bayesian statistics, providing a convenient mathematical property that simplifies the computation of posterior distributions.
We now present the formal definition as follows.

Conjugate prior
Given a family $\{p(\mathcalX \mid \btheta): \btheta \in \bTheta\}$ of generating distributions, a collection of priors $p_\omega (\btheta)$ indexed by $\bomega \in \bOmega$ is called a conjugate prior family if for any $\bomega$ and any data, the resulting posterior equals to $p_{\bomega^\prime} (\btheta \mid \mathcalX)$ for some $\bomega^\prime \in \bOmega$.

In Bayesian inference, a prior distribution represents our beliefs about the parameters of a statistical model before observing any data. A prior is considered conjugate to a likelihood function if their combination results in a posterior distribution that belongs to the same family of distributions as the prior. This property facilitates analytical solutions, making the computation of the posterior more tractable.
A toy example of the Beta-Bernoulli model can  enhance our comprehension of conjugate priors.
Suppose $\mathcalX=\{x_1, x_2, ..., x_N\}$ are drawn independently and identically distributed (i.i.d.) from a Bernoulli distribution with parameter $\theta$, i.e., $Bernoulli(x\mid \theta)$.
$Beta(\theta \mid a,b)$ distribution, with $a, b >0$, is conjugate to $Bernoulli(x\mid\theta)$, since the posterior density is $p(\theta \mid \mathcalX) = Beta(\theta \mid a+\sum x_i, b+N-\sum x_i)$.

Conjugate priors make it possible to do Bayesian reasoning in a computationally efficient manner, as well as having the philosophically satisfying interpretation of representing real or imaginary prior data. Generally, there are basically two reasons why models with conjugate priors are popular [Robert et al., 2007, Bernardo and Smith, 2009, Hoff, 2009, Gelman et al., 2013]:

* They usually allow us to derive a closed-form expression for the posterior distribution.

* They are easy to interpret, as we can easily observe how the parameters of the prior change after the Bayesian updates.

§ REGULAR UNIVARIATE MODELS AND CONJUGACY

[definition:gaussian_distribution]Gaussian, p. definition:gaussian_distribution
[definition:gamma-distribution]Gamma, p. definition:gamma-distribution
[definition:gaussian_distribution]Student's $t$, p. equation:student_t_dist

[definition:inverse_gamma_distribution]Inverse-Gamma, p. definition:inverse_gamma_distribution     [definition:truncated_normal]Truncated-Normal, p. definition:truncated_normal
[definition:gaussian_distribution]Inverse-Gaussian, p. definition:inverse_gaussian_distribution

[definition:chisquare_distribution]Chi-Square, p. definition:chisquare_distribution        [definition:gaussian_distribution]Normal-Inv-Gamma, p. definition:normal_inverse_gamma
[definition:gaussian_distribution]Inverse-Chi-Squared, p. definition:inverse-chi-square

[definition:gaussian_distribution]Nor-Inv-Chi-Squared, p. definition:normal_inverse_chi_square       [definition:general_truncated_normal]General-Truncated-Nor, p. definition:general_truncated_normal     [definition:gaussian_distribution]Half-Normal, p. definition:half_normal
[definition:gaussian_distribution]Laplace, p. definition:laplace_distribution
[definition:gaussian_distribution]Skew-Laplace, p. definition:skew_laplace_distribution      [definition:gaussian_distribution]Rectified-Normal, p. definition:reftified_normal_distribution
[definition:multinomial_dist]Multinomial, p. definition:multinomial_dist

[definition:dirichlet_dist]Dirichlet, p. definition:dirichlet_dist

[definition:poisson_distribution]Poisson, p. definition:poisson_distribution
[definition:exponential_distribution]Exponential, p. definition:exponential_distribution

[definition:multivariate_gaussian]Multi Gaussian, p. definition:multivariate_gaussian

[definition:multivariate-stu-t]Multi Student's $t$, p. definition:multivariate-stu-t

[definition:wishart_dist]Wishart, p. definition:wishart_dist

[definition:multi_inverse_wishart]Inverse-Wishart, p. definition:multi_inverse_wishart

[definition:normal_inverse_wishart]Normal-Inv-Wishart, p. definition:normal_inverse_wishart

Links for common distributions.

In most of our Bayesian matrix decomposition developments, we express the models with univariate distributions.
While in the Gaussian model, we also apply multivariate distributions.
In this section, we provide rigorous definitions for common univariate distributions and their conjugate priors.
Table <ref> provides an overview of the topics covered in this section.
Additionally, with special considerations, we also discuss the multivariate Gaussian distribution and its conjugacy in the next section.

Gaussian distribution

A random variable $\rx$ is said to follow the Gaussian distribution (a.k.a., a normal distribution) with mean and variance parameters $\mu$ and $\sigma^2>0$, denoted by $\rx \sim \normal(\mu,\sigma^2)$ [Note if two random variables $\ra$ and $\rb$ have the same distribution, then we write $\ra \sim \rb$.], if

f(x; \mu,\sigma^2)=\frac{1}{\sqrt{2\pi\sigma^2}} \exp \left\{-\frac{1}{2\sigma^2 }(x-\mu)^2 \right\}
=\sqrt{\frac{\tau}{2\pi}}\exp \left\{ -\frac{\tau}{2}(x-\mu)^2 \right\} .

The mean and variance of $\rx \sim \normal( \mu,\sigma^2)$ are given by

\Exp[\rx] = \mu, \qquad \Var[\rx] =\sigma^2=\tau^{-1},

where $\tau$ is also known as the precision of the Gaussian distribution.
Figure <ref> illustrates the impact of different parameters $\mu, \sigma^2$ for the Gaussian distribution.

Gaussian probability density functions for different values of the mean and variance parameters $\mu$ and $\sigma^2$.

Suppose $\mathcalX=\{x_1, x_2, ..., x_N\}$ are drawn i.i.d. from a Gaussian distribution of $\normal(x\mid \mu, \sigma^2)$.
For conjugate Bayesian analysis, we may rewrite the Gaussian probability density function as follows,

\begin{equation}\label{equation:uni_gaussian_likelihood}
\begin{aligned}
p(\mathcal{X} \mid \smu, \ssigma^2) &= \prod^N_{i=1} \mathcal{N} (x_i\mid\smu, \ssigma^2) \\
&= (2\pi)^{-N/2}  (\ssigma^2)^{-N/2} \exp\left\{-\frac{1}{2 \ssigma^2}  \left[ N(\overline{x} - \smu)^2 + N \sum_{n=1}^N(x_n - \overline{x})^2   \right] \right\}  \\
&= (2\pi)^{-N/2}  (\ssigma^2)^{-N/2} \exp\left\{-\frac{1}{2 \ssigma^2}  \left[  N(\overline{x} - \smu)^2 +  N S_{\overline{x}} \right] \right\},
\end{aligned}
\end{equation}

where $S_{\overline{x}}=\sum_{n=1}^N(x_n - \overline{x})^2$ and $\overline{x} = \frac{1}{N} \sum_{i=1}^{N}x_i$.
This form can help find the conditional posterior of Gaussian likelihood under normal-inverse-Gamma prior (Equation (<ref>)).

Given fixed mean $\mu$ and variance $\sigma^2$ parameters, we have

\begin{equation}\label{equation:gaussian_form_conform}
\begin{aligned}
p(x\mid \mu, \sigma^2) &=\normal(x \mid \mu, \sigma^2)
\propto \exp\left\{ -\frac{1}{2\sigma^2} x^2 + \frac{\mu}{\sigma^2} x \right\},
\end{aligned}
\end{equation}

where “$\propto$" means “proportional to."
Therefore, if we find a form conforming to the above equation, we can say the random variable $\rx$ follows the Gaussian distribution $\rx\sim \normal(\mu, \sigma^2)$.
Refer to the example of a Bayesian GGG matrix decomposition model in Equation (<ref>) (p. equation:ggg_poster_wmk1).

While the product of two Gaussian variables remains an open problem, the sum of Gaussian variables follows from a new Gaussian distribution.

Let $\rx$ and $\ry$ be two Gaussian distributed variables with means $\mu_x, \mu_y$ and variance $\sigma_x^2, \sigma_y^2$, respectively.

* When there is no correlation between the two variables, then it follows that

\rx+\ry \sim \normal(\mu_x+\mu_y, \sigma_x^2+\sigma_y^2).

* When there exists a correlation of $\rho$ between the two variables, then  it follows that

\rx+\ry \sim \normal(\mu_x+\mu_y, \sigma_x^2+\sigma_y^2+2\rho\sigma_x\sigma_y).

Conjugate prior for mean of a Gaussian distribution and Normal-Normal model.
The Gaussian distribution serves as a conjugate prior for the mean parameter of a Gaussian distribution when the variance is fixed. To see this, suppose  $\mathcalX=\{x_1, x_2, \ldots, x_N\}$ are i.i.d. normal with mean $\theta$ and precision $\lambda$, i.e., the likelihood is $\normal(x_i \mid \theta, \lambda^{-1})$, where the variance $\sigma^2=\lambda^{-1}$ is fixed, and $\theta$ is given a $\normal(\mu_0, \lambda^{-1}_0)$ prior: $\theta \sim \normal( \mu_0, \lambda^{-1})$. Using Bayes' theorem, “posterior $\propto$ likelihood $\times$ prior,"  the posterior density is

p(\theta \mid \mathcalX ) \propto \prod_{i=1}^{N} \normal(x_i \mid \theta, \lambda^{-1})
\times \normal(\theta \mid \mu_0, \lambda_0^{-1})\propto \normal(\theta \mid \widetilde{\mu}, \widetilde{\lambda}^{-1}),

where, given $\overline{x} = \frac{1}{N}\sum_{i=1}^{N} x_i$,

\begin{equation}\label{equation:posterior-param-normal-normal}
\begin{aligned}
\widetilde{\mu}
&= \frac{\lambda_0 \mu_0 + \lambda \sum_{i=1}^{N} x_i}{ \lambda_0 + N\lambda}
= \frac{\lambda_0}{\lambda_0 + N\lambda} \mu_0 + \frac{N\lambda}{\lambda_0 + N\lambda} \overline{x}, \\
\widetilde{\lambda} &= \lambda_0 + N \lambda.
\end{aligned}
\end{equation}

Thus, the posterior mean is a weighted mean of the prior mean $\mu_0$ and the sample average $\overline{x}$; the posterior precision is the sum of the prior precision and sample precision $N\lambda$.
We show that  the Gaussian distribution is, itself, a conjugate prior for the mean parameter of a Gaussian distribution when fixing variance. This is often referred to as the Normal-Normal model, and this model can be used to provide evidence for the bimodal property of human heights [Schilling et al., 2002].

Student's $t$ distribution

A random variable $\rx$ is said to follow the Student's $t$ distribution with parameters $\mu$, $\sigma^2>0$, and $\nu$, denoted by $\rx \sim \studentt(\mu,\sigma^2, \nu)$, if

\begin{equation}\label{equation:stut_defini}
\begin{aligned}
f(x; \mu, \sigma^2, \nu)&= \frac{\Gamma(\frac{\nu+1}{2})}{\Gamma(\frac{\nu}{2})} \frac{1}{\sigma\sqrt{\nu\pi}} \times \left[ 1+ \frac{(x-\mu)^2}{\nu \sigma^2}   \right]^{-(\frac{\nu+1}{2})},
%&= \frac{\Gamma(\nu/2 + D/2)}{\Gamma(\nu/2)} |\pi\bV|^{-1/2} \times \left[ 1+ \frac{1}{\nu} (\bx-\bmu)^\top \bV^{-1} (\bx-\bmu)  \right]^{-(\frac{\nu+1}{2})},
\end{aligned}
\end{equation}

where $\sigma^2$ is called the scale parameter, and $\nu$ is the degree of freedom.
The distribution has fatter tails than a Gaussian distribution. The smaller $\nu$ is, the fatter the tail.
As $\nu\rightarrow \infty$, the distribution converges towards a Gaussian.
A particular case of the Student's $t$ distribution is the Cauchy distribution,
$\rx\sim \cauchydist(\mu, \sigma^2)$ if $\rx\sim \studentt(\mu, \sigma^2, 1)$, i.e., $\nu=1$.
The mean and variance of $\rx \sim \studentt(\mu,\sigma^2, \nu)$ are given by

$$ \Exp[\rx]=\left\{
\begin{aligned}
&\mu , \, &\mathrm{if\,} \nu >1; \\
&\text{undefined}, \, &\mathrm{if\,} \nu \leq 1.
\end{aligned}
\right.\qquad
\Var[\rx]=\left\{
\begin{aligned}
&\frac{\nu}{\nu-2}\sigma^2, \, &\mathrm{if\,} \nu >2; \\
&\infty, \, &\mathrm{if\,} 1<\nu\leq 2.
\end{aligned}
\right.

Figure <ref> compares different parameters $\mu,\sigma^2, \nu$ for the Student's $t$ distribution.

In Figure <ref>, we vary the $\nu$ parameter for the Student's $t$ distribution.
As $\nu$ decreases, the distribution becomes more spread out, leading to fatter tails compared to a Gaussian distribution. This allows for more flexibility in modeling data with greater uncertainty or outliers since the Student's $t$ distribution has a greater probability of observing extreme values. In Bayesian modeling, the Student's $t$ distribution is often used as a prior for the mean parameter of a Gaussian likelihood, allowing for estimation of both the mean and precision of the data. This results in a Student's $t$-Normal model.

[Student's $t$ distribution by varying parameter $\nu$. When $\nu=100$, the distribution is very close to a Gaussian distribution.]

[Student's $t$ distribution  by varying parameter $\sigma^2$.]

Student's $t$ distribution for different values of the parameters $\nu$ and $\sigma^2$.

Gamma distribution

A random variable $\rx$ is said to follow the Gamma distribution with shape parameter $r>0$ and rate parameter $\lambda>0$ [Note the inverse rate parameter $1/\lambda$ is called the scale parameter.
In probability theory and statistics,
the location parameter shifts the entire distribution left or right, e.g., the mean parameter of a Gaussian distribution;
the shape parameter compresses or stretches the entire distribution;
the scale parameter changes the shape of the distribution in some manner.
], denoted by $\rx \sim \gammadist(r, \lambda)$, if

$$ f(x; r, \lambda)=\left\{
\begin{aligned}
&\frac{\lambda^r}{\Gamma(r)} x^{r-1} \exp(-\lambda x) ,& \mathrm{\,\,if\,\,} x \geq 0; \\
&0 , &\mathrm{\,\,if\,\,} x <0,
\end{aligned}
\right.

where $\Gamma(x)=\int_{0}^{\infty} x^{t-1}\exp(-x)dt$ is the Gamma function,  and we can just take it as a function to normalize the distribution into sum to 1. In special cases when $y$ is a positive integer, $\Gamma(y) = (y-1)!$.

The mean and variance of $\rx \sim \gammadist(r, \lambda)$ are given by

\begin{equation}
\Exp[\rx] = \frac{r}{\lambda}, \qquad \Var[\rx] = \frac{r}{\lambda^2}. \nonumber
\end{equation}

Specially, let $\rx_1, \rx_2, \ldots, \rx_n$ be i.i.d. random variables drawn from $\gammadist(r_i, \lambda)$ for each $i \in \{1, 2, \ldots, n\}$. Then $\ry = \sum_{i=1}^{n} \rx_i$ is a random variable following from $\gammadist(\sum_{i=1}^{n}r_i, \lambda)$.
Figure <ref> compares different parameters $r, \lambda$ for the Gamma distribution.

It's crucial to keep in mind that the definition of the Gamma distribution does not restrict $r$ to be a natural number and it allows for $r$ to be any positive number.
However, when $r$ is a positive integer, the Gamma distribution can be interpreted as a sum of $r$ exponentials of rate $\lambda$ (see Definition <ref>).
The summation property holds true more generally for Gamma variables with the same rate parameter. If $\rx_1$ and $\rx_2$ are random variables drawn from $\gammadist(r_1, \lambda)$ and $\gammadist(r_2, \lambda)$, respectively, then their sum $\rx_1+\rx_2$ is a Gamma random variable from $\gammadist(r_1+r_2, \lambda)$.

[Gamma distribution.]

[Inverse-Gamma distribution.]

Gamma and inverse-Gamma probability density functions for different values of the parameters $r$ and $\lambda$.

Conjugate prior for rate  of a Gamma distribution and Gamma-Gamma model.
The Gamma distribution is a conjugate prior of the rate parameter of another Gamma distribution.
Suppose data $\mathcalX =\{x_1,x_2,\ldots,x_N\}$ follows i.i.d. from a Gamma distribution $\rx_i\sim \gammadist(r, \lambda)$. Suppose further the rate parameter is given a Gamma prior $\lambda \sim \gammadist(a, \frac{a}{b})$. Using Bayes' theorem, the posterior density is

\begin{aligned}
p(\lambda \mid \mathcalX) &\propto
\prod_{i=1}^{N}\gammadist(x_i \mid r, \lambda) \times
\gammadist(\lambda \mid a, \frac{a}{b})\\
\prod_{i=1}^{N} \frac{\lambda^r}{\Gamma(r)} x_i^{r-1} \exp(-\lambda x_i)
\times
\frac{(\frac{a}{b})^a}{\Gamma(a)}
\lambda^{a-1} \exp(-\frac{a}{b} \lambda) \\
&\propto \lambda^{Nr+a-1} \exp\left\{-\bigg(\sum_{i=1}^{N}x_i +\frac{a}{b}\bigg) \lambda\right\}
\propto \gammadist(\lambda \mid\widetilde{\alpha}, \widetilde{\beta}),
\end{aligned}

\widetilde{\alpha} =Nr+a,
\qquad
\widetilde{\beta}=\sum_{i=1}^{N}x_i +\frac{a}{b}.

That is, the posterior density of rate $\lambda$ follows from a Gamma distribution.
We show that  the Gamma distribution is, itself, a conjugate prior for the rate parameter of a Gamma distribution when fixing the shape parameter. This is often referred to as the Gamma-Gamma model.

Conjugate prior for precision of a Gaussian  distribution.

The Gamma distribution  also serves as a conjugate prior for the precision parameter of a Gaussian distribution. To see this, suppose each entry $a_{mn}$ of matrix $\bA$ is i.i.d. normal model with mean $b_{mn}$ and precision $\tau$, i.e., the likelihood is $p(\bA \mid \bB, \tau^{-1})=\normal(\bA\mid\bB, \tau^{-1})$, the prior of $\tau$ is $p(\tau)=\gammadist(\tau \mid \alpha, \beta)$, where $\bA, \bB\in \real^{M\times N}$ are two matrices containing $a_{mn}$ and $b_{mn}$, respectively (the result can be applied to vector or scalar cases). Using Bayes' theorem, it can be shown that

\begin{equation}\label{equation:gamma_conjugacy_general}
\begin{aligned}
p(\tau \mid \bA, \bB, \alpha, \beta)
&\propto \normal(\bA\mid\bB, \tau^{-1})\times \gammadist(\tau\mid \alpha, \beta)\\
&=\prod_{i,j=1}^{M,N} \normal(a_{ij}\mid b_{ij}, (\tau)^{-1})
\times \frac{\beta^\alpha}{\Gamma(\alpha)} \tau^{\alpha-1} \exp(-\beta \tau)\\
&\propto \tau^{\frac{MN}{2}}\exp\left\{   -\frac{\tau}{2}  \sum_{i,j=1}^{M,N}(a_{ij} - b_{ij}  )^2\right\}
\cdot \tau^{\alpha-1}\exp(-\beta \tau)\\
&=\tau^{\frac{MN}{2}+\alpha-1} \exp\left\{   -\tau  \left(\sum_{i,j=1}^{M,N}\frac{1}{2}(a_{ij} - b_{ij}  )^2 +\beta\right)\right\}\\
&\propto \gammadist(\tau \mid \widetilde{\alpha}, \widetilde{\beta}),\\
\end{aligned}
\end{equation}

\begin{equation}\label{equation:gamma_conju_posterior}
\widetilde{\alpha}=\frac{MN}{2}+\alpha,
\qquad
\widetilde{\beta}=  \sum_{i,j=1}^{M,N}\frac{1}{2}(a_{ij} - b_{ij}  )^2 +\beta.
\end{equation}

That is, the posterior density of precision $\tau$ follows from a  Gamma distribution.

Joint conjugate prior for  Gaussian mean and precision.
Going further, when the variance/precision parameter of the Gaussian distribution is not fixed with $x_1, x_2, \ldots, x_N$ drawn i.i.d. from a normal distribution with mean $\theta$ and precision $\lambda$. The normal-Gamma distribution $\normalgamma(\alpha, \beta, \mu, c)$, with $\mu\in\real$  and $\alpha, \beta, c\in \real_+$  is a joint distribution on $(\theta, \lambda)$ by letting

\begin{aligned}
\lambda &\sim \gammadist(\alpha, \beta); \\
\theta \mid \lambda &\sim \normal(\mu, (c\lambda)^{-1}).
\end{aligned}

That is, the joint p.d.f. is

p(\theta, \lambda ) = \normal(\theta \mid \mu, (c\lambda)^{-1})\cdot \gammadist(\lambda \mid \alpha, \beta)
=\normalgamma(\theta, \lambda\mid \alpha, \beta, \mu, c).

It turns out the posterior density is again a normal-Gamma distribution with

p(\theta, \lambda ) \propto
\prod_{i=1}^{N} \normal(x_i \mid \theta, \lambda^{-1})
\cdot
\normalgamma(\theta, \lambda\mid \alpha, \beta, \mu, c)
\propto \normalgamma(\theta, \lambda \mid \widetilde{\alpha},\widetilde{\beta}, \widetilde{\mu}, \widetilde{c}),

\begin{aligned}
\widetilde{\mu} &= \frac{c\mu +\sum_{i=1}^{N}x_i}{c+N} , \gap
&\widetilde{c} &= c+N, \\
\widetilde{\alpha}&= \alpha+ \frac{N}{2}, \gap &\widetilde{\beta} &= \beta + \frac{1}{2} \bigg(c\mu^2 -\widetilde{c}\widetilde{\mu}^2  + \sum_{i=1}^{N}x_i\bigg).
\end{aligned}

In contrast to the Normal-Normal, this model is often referred to as the NormalGamma-Normal model.
The posterior mean for $\theta$ is a weighted average of the prior mean and the sample mean,

\widetilde{\mu} = \frac{c\mu +\sum_{i=1}^{N}x_i}{c+N}= \frac{c}{c+N} \mu + \frac{N}{c+N} \overline{x},

where $\overline{x} = \frac{1}{N}\sum_{i=1}^{N}x_i$.
From the posterior form of $\widetilde{c}$, the prior interpretation of $c$ can be described as the prior sample size for estimating mean parameter $\theta$. The posterior shape parameter $\widetilde{\alpha}$ grows linearly with the sample size. And the posterior rate parameter $\widetilde{\beta}$ can be written as

\widetilde{\beta}
= \beta + \frac{1}{2} \bigg(c\mu^2 -\widetilde{c}\widetilde{\mu}^2  + \sum_{i=1}^{N}x_i\bigg)
= \beta + \frac{1}{2}\sum_{i=1}^{N}(x_i - \overline{x})^2 + \frac{1}{2} \frac{cN}{c+N} (\overline{x}-\mu)^2.

In other words, it is decomposed into the sum of a prior variation, the observed variation (sample variance), and the variation between the prior mean and sample mean:

\widetilde{\beta}=\text{(prior variation)} + \frac{1}{2}N \text{(observed variation)} + \frac{1}{2}\frac{cN}{c+N} \text{(variation between means)}.

Inverse-Gamma distribution
Putting a Gamma prior over the inverse variance of a Gaussian distribution is equivalent to placing an inverse-Gamma prior on the variance.
We will now provide the formal definition of the inverse-Gamma distribution.

A random variable $\rx$ is said to follow the inverse-Gamma distribution with shape parameter $r>0$ and scale parameter $\lambda>0$, denoted by $\rx\sim \inversegammadist(r, \lambda)$, if

$$ f(x; r, \lambda)=\left\{
\begin{aligned}
&\frac{\lambda^r}{\Gamma(r)} x^{-r-1} \exp(- \frac{\lambda}{x} ) ,& \mathrm{\,\,if\,\,} x > 0;  \\
&0 , &\mathrm{\,\,if\,\,} x \leq 0.
\end{aligned}
\right.

And it is denoted by $\rx \sim \inversegammadist(r, \lambda)$.
The mean and variance of inverse-Gamma distribution are given by

$$ \Exp[\rx]=\left\{
\begin{aligned}
&\frac{\lambda}{r-1}, \, &\mathrm{if\,} r\geq 1; \\
&\infty, \, &\mathrm{if\,} 0<r<1.
\end{aligned}
\right.\qquad
\Var[\rx]=\left\{
\begin{aligned}
&\frac{\lambda^2}{(r-1)^2(r-2)}, \, &\mathrm{if\,} r> 2; \\
&\infty, \, &\mathrm{if\,} 0<r\leq 2.
\end{aligned}
\right.

Figure <ref> illustrates the impact of different parameters $r, \lambda$ for the inverse-Gamma distribution.

If $\rx$ is Gamma distributed, then $\ry=1/\rx$ is inverse-Gamma distributed.
Note that the inverse-Gamma density is not simply the Gamma density with
$x$ replaced by $\frac{1}{y}$. There is an additional factor of $y^{-2}$. [Which is from the Jacobian in the change-of-variables formula. A short proof is provided here. Let $y=\frac{1}{x}$ where $y\sim \inversegammadist(r, \lambda)$ and $x\sim \gammadist(r, \lambda)$. Then, $f(y) |dy| = f(x) |dx|$, which results in $f(y) = f(x) \abs{\frac{dx}{dy}} = f(x)x^2 \xlongequal{ \mathrm{y}=\frac{1}{x}} \frac{\lambda^r}{\Gamma(r)} y^{-r-1} \exp(- \frac{\lambda}{y})$ for $y>0$. ]  The inverse-Gamma distribution is useful as a prior for positive parameters. It imparts a quite heavy tail and keeps probability further from zero than the Gamma distribution (see examples in Figure <ref>).

Conjugate prior for variance of a  Gaussian distribution.
The inverse-Gamma distribution is a conjugate prior for the variance parameter of a Gaussian distribution with a fixed mean parameter. To see this, let the likelihood be $p(\bA \mid \bB, \sigma^2)=\normal(\bA\mid\bB, \sigma^2)$, where $\bA,\bB\in\real^{M\times N}$ are two matrices containing elements of $a_{mn}$ and $b_{mn}$, respectively (again the result can be applied to vector or scalar cases), and let the prior of $\sigma^2$ be $p(\sigma^2)=\inversegammadist(\sigma^2 \mid \alpha, \beta)$. Using Bayes' theorem, it can be shown that

\begin{equation}\label{equation:inverse_gamma_conjugacy_general}
\begin{aligned}
p(\sigma^2 \mid \bA, \bB, \alpha, \beta)
&\propto \normal(\bA\mid \bB, \sigma^2)\times \inversegammadist(\sigma^2 \mid \alpha, \beta)\\
&=\prod_{i,j=1}^{M,N} \normal(a_{ij}\mid b_{ij}, \sigma^2)
\times \frac{\beta^\alpha}{\Gamma(\alpha)} (\sigma^2)^{-\alpha-1} \exp(-\frac{\beta}{\sigma^2})\\
&\propto \frac{1}{\sigma^{MN}}\exp\left\{   -\frac{1}{2\sigma^2}  \sum_{i,j=1}^{M,N}(a_{ij} - b_{ij}  )^2\right\}
\cdot (\sigma^2)^{-\alpha-1}\exp(-\frac{\beta}{\sigma^2})\\
&=(\sigma^2)^{-\frac{MN}{2}-\alpha-1} \exp\left\{   -\frac{1}{\sigma^2}  \left(\sum_{i,j=1}^{M,N}\frac{1}{2}(a_{ij} - b_{ij}  )^2 +\beta\right)\right\}\\
&\propto \inversegammadist(\sigma^2 \mid \widetilde{\alpha}, \widetilde{\beta}),\\
\end{aligned}
\end{equation}

\begin{equation}\label{equation:inversegamma_conjugate_posterior}
\widetilde{\alpha}=\frac{MN}{2}+\alpha,
\qquad
\widetilde{\beta}=  \sum_{i,j=1}^{M,N}\frac{1}{2}(a_{ij} - b_{ij}  )^2 +\beta.
\end{equation}

That is, the posterior density of the variance $\sigma^2$ is also an inverse-Gamma distribution. And as claimed, we find the posterior parameters in Equation (<ref>) are exactly the same as that in Equation (<ref>) from a Gamma prior.

Normal-inverse-Gamma distribution
As we have seen that the normal-Gamma density is a joint conjugate prior for the mean and precision parameters of a Gaussian distribution. The normal-inverse-Gamma (NIG) distribution defined as follows is a joint conjugate prior for the mean and variance parameters of a Gaussian distribution.

The joint density of  normal-inverse-Gamma distribution is a  density defined as

\begin{equation}
\begin{aligned}
&\gap \nig (\smu, \sigma^2 \mid m, \kappa, r, \lambda)
= \mathcal{N} (\mu\mid m, \frac{\ssigma^2}{\kappa})  \cdot \inversegammadist (\ssigma^2 \mid r,  \lambda) \\
&=\frac{1}{Z_{\nig}(\kappa, r, \lambda)}  (\ssigma^2)^{-\frac{2r +3}{2}} \exp\left\{-\frac{1}{2 \ssigma^2}\left[\kappa(m-\mu)^2 + 2\lambda \right] \right\},  \\
\end{aligned}
\label{equation:uni_gaussian_prior-nig}
\end{equation}

where $\sigma^2, r, \lambda>0$, and $Z_{\nig}(\kappa, r, \lambda)$ is a normalizing constant:

\begin{equation}\label{equation:uni_gaussian_giw_constant-nig}
Z_{\nig}(\kappa, r, \lambda) = \frac{\Gamma(r)}{\lambda^{r}} \sqrt{\frac{2\pi}{\kappa}}.
\end{equation}

Figure <ref> shows some normal-inverse-Gamma probability density functions by varying different parameters.

Joint conjugate prior for the Gaussian mean and variance (under NIG).
The normal-inverse-Gamma defines an equivalent prior over the mean and variance parameters of a Gaussian distribution as the normal-Gamma prior, but is sometimes more convenient than the latter one.
Similar to the normal-Gamma prior, when the variance and mean parameters of the Gaussian distribution are not fixed with $N$ data points $\mathcalX=\{x_1, x_2, \ldots, x_N\}$ drawn i.i.d. from a normal distribution with mean $\mu$ and variance $\sigma^2$.
The normal-inverse-Gamma $\normalinversegamma(m_0, \kappa_0, r_0, \lambda_0)$ with $m_0\in\real$ and $r_0, \lambda_0, \kappa_0\in\real_+$ is a joint distribution on $\mu, \sigma^2$ by letting

\begin{aligned}
\sigma^2 &\sim \inversegammadist(r_0, \lambda_0);\\
\mu \mid \sigma^2 &\sim \normal(m_0, \frac{\sigma^2}{\kappa_0}).
\end{aligned}

With this prior, $\mu$ and $\sigma^2$ decouple, and the posterior conditional densities of $\mu$ and $\sigma^2$ are Gaussian and inverse-Gamma, respectively.
The joint p.d.f of NIG prior can be expressed as

p(\mu, \sigma^2) = \normal(m_0, \frac{\sigma^2}{\kappa_0}) \cdot \inversegammadist(r_0, \lambda_0)
= \normalinversegamma(\mu, \sigma^2 \mid m_0, \kappa_0, r_0, \lambda_0).

Again, by  Bayes' theorem “$\mathrm{posterior} \propto \mathrm{likelihood} \times \mathrm{prior} $," the  posterior of the $\mu$ and $\sigma^2$ parameters under the NIG prior is

\begin{equation}\label{equation:conjugate_nigamma_general}
\begin{aligned}
&\gap p(\mu, \sigma^2\mid \mathcalX, \bbeta ) \\
&\propto \normal(\mathcalX \mid \mu, \sigma^2)
\cdot  \normalinversegamma(\mu, \sigma^2 \mid \bbeta)\\
&\propto \prod_{i=1}^{N}\normal(x_i\mid \mu, \sigma^2)
\cdot \normalinversegamma(\mu, \sigma^2 \mid  m_0, \kappa_0, r_0, \lambda_0)\\
%  = p(\mathcalX, \mu, \sigma^2 | m_0, \kappa_0, \alpha_0, \beta_0)\\
%&= (2\pi)^{-N/2}  (\ssigma^2)^{-N/2} \exp\left(-\frac{1}{2 \ssigma^2}  \left[  N \sum_{n=1}^N(x_n - \overline{x})^2 + N(\overline{x} - \smu)^2  \right] \right)\\
%&\gap \times \frac{1}{Z_{\nix}(\kappa_0, \nu_0, S_0)} (\ssigma^2)^{-(\nu_0/2 + 3/2)} \exp\left( -\frac{1}{2 \ssigma^2} \left[\nu_0 S_0 + \kappa_0(m_0-\mu)^2\right] \right)\\
&\stackrel{\star}{=} \frac{C}{(\sigma^2)^{\frac{2r_0 + 3+N}{2}}} \exp\left\{-\frac{1}{2 \ssigma^2}  \left[  N(\overline{x} - \smu)^2 +  N S_{\overline{x}} \right] \right\}
\exp\left\{ -\frac{1}{2 \ssigma^2} \left[2\lambda_0 + \kappa_0(m_0-\mu)^2\right] \right\}\\
&\propto (\sigma^2)^{-\frac{2r_N + 3}{2}}\exp\left\{ -\frac{1}{2 \ssigma^2} \left[ \lambda_N + \kappa_N(m_N-\mu)^2\right] \right\}\\
&\propto\nig(\mu, \sigma^2 \mid m_N, \kappa_{N}, r_N, \lambda_N),
\end{aligned}
\end{equation}

where $\bbeta=\{m_0, \kappa_0, r_0, \lambda_0\}$, $C=\frac{(2\pi)^{-N/2}}{Z_{\nig}(\kappa_0, r_0, \lambda_0)}$, equality $(\star)$ is from Equation (<ref>), and

\begin{aligned}
m_N &= \frac{\kappa_0 m_0 + N\overline{x}}{\kappa_{N}} = \frac{\kappa_0 }{\kappa_{N}}m_0 + \frac{N}{\kappa_{N}}\overline{x},\\
\kappa_{N}&= \kappa_{0} +N,\\
r_N &= r_0 +\frac{N}{2},\\
\lambda_N &=\lambda_0 +\frac{1}{2}(NS_{\overline{x}} + N\overline{x}^2 + \kappa_{0} m_0^2 -\kappa_{N}m_N^2)\\
&= \lambda_0+\frac{1}{2}\left(NS_{\overline{x}} + \frac{\kappa_0 N }{\kappa_{0}+N} (\overline{x} - m_0)^2\right),
\end{aligned}

where $S_{\overline{x}}=\sum_{n=1}^N(x_n - \overline{x})^2$ and $\overline{x} = \frac{1}{N} \sum_{i=1}^{N}x_i$.
Note in the above derivation, we use the  fact about the likelihood under Gaussian in Equation (<ref>).
We will discuss the posterior marginal likelihood in the normal-inverse-Chi-squared (NIX) case. Further discussion on the posterior marginal likelihood for the NIG prior can be found in Murphy, 2007. We will leave this to the readers as it is rather similar as that in the NIX prior.

[Contour plot of normal-inverse-Gamma density by varying parameter $r$ (purple=low, yellow=high). ]

[Contour plot of normal-inverse-Gamma density by varying parameter $\lambda$ (purple=low, yellow=high).]

[Contour plot of normal-inverse-Gamma density by varying parameter $\kappa$ (purple=low, yellow=high).]

[Contour plot of normal-inverse-Gamma density by varying parameter $m$ (purple=low, yellow=high).]

Normal-inverse-Gamma probability density functions by varying different parameters.

Chi-squared distribution
Another distribution that is closely related to the Gamma distribution is called the Chi-squared distribution, which is extensively used in the distribution theory of linear models [Lu, 2022].
The rigorous definition is given as follows.

Let $\bA \sim \normal(\bzero, \bI_{p})$, where $\bI_p$ is a $p\times p$ identity matrix. Then $\rx=\sum_i^p a_{ii}$ follows the Chi-squared distribution with $p$ degrees of freedom. We write $\rx \sim \chi^2(p)$, and we can see this is equivalent to $\rx\sim \gammadist(p/2, 1/2)$:

f(x; p)=\left\{
\begin{aligned}
&\frac{1}{2^{p/2}\Gamma(\frac{p}{2})} x^{\frac{p}{2}-1} \exp(-\frac{x}{2}) ,& \mathrm{\,\,if\,\,} x \geq 0;  \\
&0 , &\mathrm{\,\,if\,\,} x <0.
\end{aligned}
\right.

The mean, variance of $\rx\sim \chi^2(p)$ are given by

\Exp[\rx]=p, \qquad \Var[\rx]=2p.

Figure <ref> compares different degrees of freedom $p$ for the Chi-squared distribution.

[Chi-squared distribution.]

[Inverse-Chi-squared distribution.]

Chi-squared and inverse-Chi-squared probability density functions for different values of  parameters.

An counterpart of the inverse-Gamma distribution is known as the inverse-Chi-squared distribution.
Following  the definition of the inverse-Gamma distribution in Definition <ref>, we provide the rigorous definition of the inverse-Chi-squared distribution as follows.

A random variable $\rx$ is said to follow the inverse-Chi-squared distribution with parameter $\nu>0$ and $s^2>0$, denoted by $\rx\sim \inversegammadist(\frac{\nu}{2}, \frac{\nu s^2}{2})$ if

$$ f(x; \nu, s^2)=\left\{
\begin{aligned}
&\frac{{(\frac{\nu s^2}{2})}^{\frac{\nu}{2}}}{\Gamma(\frac{\nu}{2})} x^{-\frac{\nu}{2}-1} \exp(- \frac{\nu s^2}{2x} ) ,& \mathrm{\,\,if\,\,} x > 0;  \\
&0 , &\mathrm{\,\,if\,\,} x \leq 0.
\end{aligned}
\right.

And it is also compactly denoted by $\rx \sim \inversechidist(\nu, s^2)$. The parameter $\nu >0$ is called the degrees of freedom, and $s^2 > 0$ is the scale parameter. And it is also known as the scaled inverse-Chi-squared distribution.
The mean and variance of the inverse-Chi-squared distribution are given by

$$ \Exp[\rx]=\left\{
\begin{aligned}
&\frac{\nu s^2}{\nu-2}, \, &\mathrm{if\,\,} \nu\geq 2; \\
&\infty, \, &\mathrm{if\,\,} 0<\nu<2.
\end{aligned}
\right.\qquad
\Var[\rx]=\left\{
\begin{aligned}
&\frac{2\nu^2 s^4}{(\nu-2)^2(\nu-4)}, \, &\mathrm{if\,\,} \nu\geq 4; \\
&\infty, \, &\mathrm{if\,\,} 0<\nu<4.
\end{aligned}
\right.

To establish a connection with the inverse-Gamma distribution, we can set $S=\nu s^2$. Then the inverse-Chi-squared distribution can also be denoted by $\rx\sim \inversegammadist(\frac{\nu}{2}, \frac{S}{2})$ if $\rx \sim \inversechidist(\nu, s^2)$,  the form of which conforms to the univariate case of the inverse-Wishart distribution (Definition <ref>).
And we will observe the similarities in the posterior parameters as well.
Figure <ref> illustrates the impact of different parameters $\nu, s^2$ for the inverse-Chi-squared distribution.

[Conjugate prior for the variance of a Gaussian distribution]
Show that the inverse-Chi-squared distribution is a conjugate prior for the Gaussian variance parameter when the mean parameter is fixed. Hint: the derivation follows the same procedure as in the inverse-Gamma case.

Joint conjugate prior
Normal-inverse-Chi-squared distribution
As we have seen that the normal-inverse-Gamma distribution is a joint conjugate prior for Gaussian mean and variance parameters. The normal-inverse-Chi-squared (NIX) distribution defined as follows is an alternative joint conjugate prior.

Similar to the normal-inverse-Gamma distribution, the normal-inverse-Chi-squared (NIX) distribution is defined as (where again we set $S=\nu s^2$ as that in the inverse-Chi-square distribution to establish a connection with the normal-inverse-Gamma density)

\begin{equation}
\begin{aligned}
&\gap \nix (\smu, \sigma^2 \mid m, \kappa, \nu, S)
= \mathcal{N} (\mu\mid m, \frac{\ssigma^2}{\kappa})
\cdot \inversechidist (\ssigma^2 \mid \nu,  s^2) \\
&=\frac{1}{Z_{\nix}(\kappa, \nu, s^2)} (\ssigma^2)^{-(\nu/2 + 3/2)} \exp\left\{ -\frac{1}{2 \ssigma^2} \left[\nu s^2 + \kappa(m-\mu)^2\right] \right\} \\
&\xlongequal{S = \nu s^2}
\frac{1}{Z_{\nix}(\kappa, \nu, s^2)} (\ssigma^2)^{-(\nu/2 + 3/2)} \exp\left\{ -\frac{1}{2 \ssigma^2} \left[S + \kappa(m-\mu)^2\right] \right\},
\end{aligned}
\label{equation:uni_gaussian_prior}
\end{equation}

where $\sigma^2, \nu, s^2>0$, and $Z_{\nix}(\kappa, \nu, s^2)$ is a normalizing constant:

\begin{equation}\label{equation:uni_gaussian_giw_constant}
Z_{\nix}(\kappa, \nu, s^2) = \Gamma\big(\frac{\nu}{2}\big)
\big(\frac{2}{\nu s^2}\big)^{\nu/2} \sqrt{\frac{2\pi}{\kappa}}
=  \Gamma\big(\frac{\nu}{2}\big) \big(\frac{2}{S}\big)^{\nu/2} \sqrt{\frac{2\pi}{\kappa}}.
\end{equation}

The  normal-inverse-Chi-squared distribution can also be denoted by
$\rx\sim \normalinversegamma(m, \kappa, \frac{\nu}{2}, \frac{S}{2})$ if
$\rx \sim \nix(m, \kappa, \nu, s^2)$, the form of which  conforms to the univariate case of the normal-inverse-Wishart distribution (Equation (<ref>)). And we will see the similarities in the posterior parameters as well.

Joint conjugate prior for the Gaussian mean and variance (under NIX).
Similar to the normal-inverse-Gamma prior, when the variance and mean parameters of the Gaussian distribution are not fixed with $N$ data points $\mathcalX=\{x_1, x_2, \ldots, x_N\}$ drawn i.i.d. from a normal distribution with mean $\mu$ and variance $\sigma^2$.
The normal-inverse-Chi-squared $\nix(m_0, \kappa_0, \nu_0, S_0=\nu_0\sigma_0^2)$ with $m_0\in\real$ and $\kappa_0, \mu_0, S_0\in\real_+$ is a joint distribution on $\mu, \sigma^2$ by letting

\begin{aligned}
\sigma^2 &\sim \inversechidist(\nu_0, \sigma_0^2);\\
\mu \mid \sigma^2 &\sim \normal(m_0, \frac{\sigma^2}{\kappa_0}).
\end{aligned}

Again, by Bayes' theorem “$\mathrm{posterior} \propto \mathrm{likelihood} \times \mathrm{prior} $," the conditional posterior of the $\mu$ and $\sigma^2$ parameters under the NIX prior is

\begin{equation}\label{equation:nix-posterior}
\begin{aligned}
&\gap p(\mu, \sigma^2\mid \mathcalX, \bbeta ) \\
&\propto p(\mathcalX \mid \mu, \sigma^2) p(\mu, \sigma^2 \mid \bbeta) = p(\mathcalX, \mu, \sigma^2 \mid \bbeta)\\
%&= (2\pi)^{-N/2}  (\ssigma^2)^{-N/2} \exp\left(-\frac{1}{2 \ssigma^2}  \left[  N \sum_{n=1}^N(x_n - \overline{x})^2 + N(\overline{x} - \smu)^2  \right] \right)\\
%&\gap \times \frac{1}{Z_{\nix}(\kappa_0, \nu_0, S_0)} (\ssigma^2)^{-(\nu_0/2 + 3/2)} \exp\left( -\frac{1}{2 \ssigma^2} \left[\nu_0 S_0 + \kappa_0(m_0-\mu)^2\right] \right)\\
&= \frac{C}{(\sigma^2)^{\frac{\nu_0 + 3+N}{2}}} \exp\left\{ -\frac{1}{2 \ssigma^2}  \left[  N(\overline{x} - \smu)^2 +  N S_{\overline{x}} \right] \right\}
\exp\left\{ -\frac{1}{2 \ssigma^2} \left[S_0 + \kappa_0(m_0-\mu)^2\right] \right\}\\
&= C\times  (\sigma^2)^{-\frac{\nu_N + 3}{2}}\exp\left\{ -\frac{1}{2 \ssigma^2}
\left[ S_N + \kappa_N(m_N-\mu)^2\right] \right\}\\
&\propto\nix(\mu, \sigma^2\mid m_N, \kappa_{N}, \nu_N, \textcolor{mylightbluetext}{S_N})
= \normal (\mu\mid m_N, \frac{\ssigma^2}{\kappa_N})  \cdot \inversechidist (\ssigma^2 \mid \nu_N,  \textcolor{mylightbluetext}{\sigma^2_N}),
\end{aligned}
\end{equation}

where $\bbeta=\{m_0, \kappa_0, \nu_0, S_0=\nu_0\sigma_0^2\}$, $C=\frac{(2\pi)^{-N/2}}{Z_{\nix}(\kappa_0, \nu_0, \sigma^2_0)}$, and

\begin{aligned}
m_N &= \frac{\kappa_0 m_0 + N\overline{x}}{\kappa_{N}} = \frac{\kappa_0 }{\kappa_{N}}m_0 + \frac{N}{\kappa_{N}}\overline{x},\\
\kappa_{N}&= \kappa_{0} +N,\\
\nu_N &= \nu_0 +N,\\
S_N &= S_0 +NS_{\overline{x}} + N\overline{x}^2 + \kappa_{0} m_0^2 -\kappa_{N}m_N^2\\
&=S_0 +NS_{\overline{x}} + \frac{\kappa_0 N }{\kappa_{0}+N} (\overline{x} - m_0)^2,\\
\nu_N \sigma_N^2 &= S_N \leadto     \sigma_N^2 = \frac{S_N}{\nu_N} ,
\end{aligned}

Then the posterior density is a normal-inverse-Chi-squared density. [This posterior shares the  same form as that in the multivariate case from Equation (<ref>) except the $N$ in $NS_{\overline{x}}$, which results from the difference between the multivariate Gaussian distribution and the univariate Gaussian distribution. Similarly, in the inverse-Chi-squared language, we can show that $\nu_N \sigma_N^2 = S_N$.]

Suppose $\nu_0\geq 2$, or $N\geq 2$ such that $\nu_N\geq 2$, the posterior expectations are given by

\Exp[\mu \mid \mathcalX, \bbeta] = m_N, \qquad
\Exp[\sigma^2 \mid \mathcalX, \bbeta] = \frac{S_N}{\nu_N-2}.

Marginal posterior of $\sigma^2$.
Integrating out $\mu$ in the posterior, we have

\begin{aligned}
p(\sigma^2 \mid \mathcalX, \bbeta)
&= \int_{\mu} p(\mu, \sigma^2 \mid \mathcalX, \bbeta) d \mu \\
&= \int_{\mu }  \mathcal{N} (\mu\mid m_N, \frac{\ssigma^2}{\kappa_N})
\cdot \inversechidist (\ssigma^2 \mid \nu_N,  \sigma^2_N)d\mu \\
&= \inversechidist (\ssigma^2\mid \nu_N,  \sigma^2_N),
\end{aligned}

which is just an integral over a Gaussian distribution.

Marginal posterior of $\mu$.
Integrating out $\sigma^2$ in the posterior, we have

\begin{aligned}
p(\mu \mid \mathcalX, \bbeta)
&= \int_{\sigma^2} p(\mu, \sigma^2 \mid \mathcalX, \bbeta) d \sigma^2 \\
&= \int_{\sigma^2 }  \mathcal{N} (\mu\mid m_N, \frac{\ssigma^2}{\kappa_N})
\cdot \inversechidist (\ssigma^2 \mid \nu_N,  \sigma^2_N)d\sigma^2 \\
&= \int_{\sigma^2}C(\sigma^2)^{-\frac{\nu_N + 3}{2}}
\exp\left \{-\frac{1}{2 \ssigma^2} \left[ S_N + \kappa_N(m_N-\mu)^2\right] \right\}d\sigma^2. %\qquad \text{$C$ given in Equation~\eqref{equation:nix-posterior}}\\
\end{aligned}

Let $\phi = \sigma^2$ and $\alpha = (\nu_N+1)/2$, $A =  S_N + \kappa_N(m_N-\mu)^2$, and $x = \frac{A}{2\phi}$, we have

\frac{d \phi}{d x} = -\frac{A}{2}x^{-2},

where $A$ can be easily verified to be positive and $\phi=\sigma^2>0$.
It follows that

\begin{equation*}
\begin{aligned}
p(\mu \mid \mathcalX, \bbeta) &=\int_{0}^{\infty} C(\phi)^{-\alpha-1} \exp\left(-\frac{A}{2 \phi} \right)d\phi\\
&=\int_{\textcolor{black}{\infty}}^{\textcolor{black}{0}}  C(\frac{A}{2x})^{-\alpha-1} \exp\left( -x\right)   ( \textcolor{black}{-}\frac{A}{2}x^{-2}) dx \qquad &\text{(since $x=\frac{A}{2\phi}$)}\\
&=\int_{\textcolor{black}{0}}^{\textcolor{black}{\infty}}  C(\frac{A}{2x})^{-\alpha-1} \exp\left( -x\right)   ( \frac{A}{2}x^{-2}) dx\\
&= (\frac{A}{2})^{-\alpha} \int_{x}  Cx^{\alpha-1} \exp\left( -x\right)   dx \\
&= (\frac{A}{2})^{-\alpha}  (C\cdot \Gamma(1)) \int_{x}\gammadist(x\mid \alpha, 1) dx\qquad &\text{(see Definition~\ref{definition:gamma-distribution})}\\
&= (C\cdot \Gamma(1))\left[ \nu_N\sigma_N^2 +\kappa_{N}(m_N-\mu)^2 \right]^{-\frac{\nu_N+1}{2}}\\
&\overset{(a)}{=} (C\cdot \Gamma(1)) (\nu_N\sigma_N^2)^{-\frac{\nu_N+1}{2}} \left[ 1 +\frac{\kappa_{N}}{\nu_N\sigma_N^2}(m_N-\mu)^2 \right]^{-\frac{\nu_N+1}{2}}
\end{aligned}
\end{equation*}

We notice that $C$ is defined in Equation (<ref>) (in terms of $\{\kappa_N, \nu_N, \sigma^2_N\}$) with

C\overset{(b)}{=}\frac{(2\pi)^{-N/2}}{Z_{\nix}(\kappa_N, \nu_N, \sigma^2_N)} = \frac{(2\pi)^{-N/2}}{\frac{\sqrt{(2\pi)}}{\sqrt{\kappa_N}} \Gamma(\frac{\nu_N}{2}) (\frac{2}{\nu_N \sigma^2_N})^{\nu_N/2}} \propto (\nu_N \sigma^2_N)^{\nu_N/2}.

Combining equalities (a) and (b) above, we obtain

p(\mu\mid \mathcalX, \bbeta)  \propto \frac{1}{\sigma_N/\sqrt{\kappa_N}}
\left[ 1 +\frac{\kappa_{N}}{\nu_N\sigma_N^2}(\mu-m_N)^2 \right]^{-\frac{\nu_N+1}{2}}
\propto \tau(\mu\mid m_N, \sigma_N^2/\kappa_N, \nu_N),

which is a univariate Student's $t$ distribution (Definition <ref>).

Marginal likelihood of data.
By Equation (<ref>), we can obtain the marginal likelihood of data under hyper-parameter $\bbeta=( m_0, \kappa_0, \nu_0, S_0=\nu_0\sigma_0^2)$

\begin{aligned}
p(\mathcalX \mid \bbeta)
&= \int_{\mu}
\int_{\sigma^2}  p(\mathcalX, \mu, \sigma^2 \mid\bbeta)
d\mu d\sigma^2\\
%&= \int_{\mu} \int_{\sigma^2} \normal (\mu| m_N, \frac{\ssigma^2}{\kappa_N})  \cdot \inversechidist (\ssigma^2 | \nu_N,  \sigma^2_N)d\mu d\sigma^2 \qquad \text{(see Equation~\ref{equation:nix-posterior})}\\
&=\frac{(2\pi)^{-N/2}}{Z_{\nix}(\kappa_0, \nu_0, \sigma^2_0)}  \int_{\mu} \int_{\sigma^2}
(\sigma^2)^{-\frac{\nu_N + 3}{2}}\exp\left\{ -\frac{1}{2 \ssigma^2} \left[ S_N + \kappa_N(m_N-\mu)^2\right] \right\}
d\mu d\sigma^2\\
&= (2\pi)^{-N/2}\frac{Z_{\nix}(\kappa_N, \nu_N, \sigma^2_N)}{Z_{\nix}(\kappa_0, \nu_0, \sigma^2_0)} \\
&= (\pi)^{-N/2} \frac{\Gamma(\nu_N/2)}{\Gamma(\nu_0/2)} \sqrt{\frac{\kappa_0}{\kappa_N}} \frac{(\nu_0\sigma^2_0)^{\nu_0/2}}{(\nu_N\sigma^2_N)^{\nu_N/2}}.
\end{aligned}

Posterior predictive for new data with observations.
Let the number of samples for data set $\{\xstar, \mathcalX\}$ be $\Nstar = N+1$, we have

\begin{equation}\label{equation:nix-posterior-new-withobser}
\begin{aligned}
&\gap p(\xstar \mid\mathcalX, \bbeta)
= \frac{p(\xstar, \mathcalX \mid \bbeta)}{p(\mathcalX\mid \bbeta)}\\
&=\left\{(2\pi)^{-\Nstar/2}\frac{Z_{\nix}(\kappa_\Nstar, \nu_\Nstar, \sigma^2_\Nstar)}{Z_{\nix}(\kappa_0, \nu_0, \sigma^2_0)}\right\}
\bigg/\left\{(2\pi)^{-N/2}\frac{Z_{\nix}(\kappa_N, \nu_N, \sigma^2_N)}{Z_{\nix}(\kappa_0, \nu_0, \sigma^2_0)}\right\}\\
&=(2\pi)^{-1/2} \frac{Z_{\nix}(\kappa_\Nstar, \nu_\Nstar, \sigma^2_\Nstar)}{Z_{\nix}(\kappa_N, \nu_N, \sigma^2_N)}\\
\sqrt{\frac{\kappa_N}{\kappa_{N^{\star}}}}    \frac{\Gamma(\frac{\nu_{N^{\star}}}{2})}{\Gamma(\frac{\nu_{N}}{2})}
\frac{(\nu_N \sigma_N^2)^{\frac{\nu_N}{2}}}{(\nu_{\Nstar}\sigma_{\Nstar}^2)^{\frac{\nu_{\Nstar}}{2}}}\\
\sqrt{\frac{\kappa_N}{(\kappa_{N}+1)}   \frac{1}{(\pi\nu_{N}\sigma_{N}^2)}}
%\frac{(\nu_{\Nstar}\sigma_{\Nstar}^2)^{-\frac{\nu_{N}+1}{2}}}{(\nu_N \sigma_N^2)^{-\frac{\nu_N +1}{2}}}
\left[\frac{(\nu_{\Nstar}\sigma_{\Nstar}^2)}{(\nu_N \sigma_N^2)}\right]^{-\frac{\nu_{N}+1}{2}}.
\end{aligned}
\end{equation}

We realize that

\begin{aligned}
m_N &= \frac{\kappa_{N^{\star}}m_{N^{\star}} - x^\star}{\kappa_N}=\frac{(\kappa_0 + N + 1)m_{N^{\star}} - x^\star}{\kappa_0 + N} , \\
m_{N^{\star}} &= \frac{\kappa_{N} m_N +x^{\star}}{\kappa_{N^{\star}}}
= \frac{(\kappa_{0}+N) m_N +x^{\star}}{\kappa_{0}+N+1},\\
S_{N^{\star}} &= S_N + x^{\star} x^{\star T} - \kappa_{N^{\star}} m_{N^{\star}}^2  + \kappa_N m_N^2 \\
&= S_N + \frac{\kappa_N + 1}{\kappa_N}(m_{N^{\star}} - x^\star)^2\\
&=S_N + \frac{\kappa_N }{\kappa_N+ 1}(m_{N} - x^\star)^2,
\end{aligned}

Thus, we have

\begin{equation}\label{equation:nix-substitute-posterior-withobser}
\begin{aligned}
\left[\frac{(\nu_{\Nstar}\sigma_{\Nstar}^2)}{(\nu_N \sigma_N^2)}\right]^{-\frac{\nu_{N}+1}{2}}&=
\left(\frac{S_\Nstar}{S_N}\right)^{-\frac{\nu_{N}+1}{2}} =1 +  \frac{\kappa_N(m_{N} - x^\star)^2 }{(\kappa_N+ 1)\nu_N\sigma_N^2}.
\end{aligned}
\end{equation}

Substituting Equation (<ref>) into Equation (<ref>), it follows that

\begin{aligned}
p(\bxstar \mid\mathcalX, \bbeta)
\sqrt{\frac{\kappa_N}{(\kappa_{N}+1)}   \frac{1}{(\pi\nu_{N}\sigma_{N}^2)}}
\left[1 +  \frac{\kappa_N(m_{N} - x^\star)^2 }{(\kappa_N+ 1)\nu_N\sigma_N^2}\right]^{-\frac{\nu_{N}+1}{2}}\\
&= \tau(x^\star \mid m_N, \frac{\kappa_{N}+1}{\kappa_{N}}\sigma^2_N, \nu_N  ).
\end{aligned}

Posterior predictive for new data without observations.
Similarly, considering the case without observations, we have

\begin{aligned}
p(\xstar \mid \bbeta)
&=  \int_{\mu} \int_{\sigma^2} p(\xstar, \mu, \sigma^2 \mid \bbeta) d\mu d\sigma^2 \\
&=(2\pi)^{-1/2}\frac{Z_{\nix}(\kappa_1, \nu_1, \sigma^2_1)}{Z_{\nix}(\kappa_0, \nu_0, \sigma^2_0)}\\
\sqrt{\frac{\kappa_0}{\kappa_{1}}}    \frac{\Gamma(\frac{\nu_{1}}{2})}{\Gamma(\frac{\nu_{0}}{2})}
\frac{(\nu_0 \sigma_0^2)^{\frac{\nu_0}{2}}}{(\nu_{1}\sigma_{1}^2)^{\frac{\nu_{1}}{2}}}\\
\sqrt{\frac{\kappa_0}{(\kappa_{0}+1)}   \frac{1}{(\pi\nu_{0}\sigma_{0}^2)}}
\left[\frac{(\nu_{1}\sigma_{1}^2)}{(\nu_0 \sigma_0^2)}\right]^{-\frac{\nu_{0}+1}{2}}\\
&=\tau\big(x^\star \mid m_0, \frac{\kappa_{0}+1}{\kappa_{0}}\sigma^2_0, \nu_0  \big).
\end{aligned}

§ EXPONENTIAL AND CONJUGACY

Exponential distribution

The exponential distribution is a  probability distribution commonly used in modeling events occurring  randomly over time, such as  the time elapsed until the occurrence of a certain event, or the time between two consecutive events.
It is a special Gamma distribution with support on nonnegative real values.

A random variable $\rx$ is said to follow the exponential distribution with rate parameter $\lambda>0$, denoted by $\rx \sim \exponential(\lambda)$, if

$$ f(x; \lambda)=\left\{
\begin{aligned}
& \lambda \exp(-\lambda x)
,& \mathrm{\,\,if\,\,} x \geq 0;  \\
&0 , &\mathrm{\,\,if\,\,} x <0.
\end{aligned}
\right.

We can see this is equivalent to $\rx\sim \gammadist(1, \lambda)$.
The mean and variance of $\rx \sim \exponential(\lambda)$ are given by

\begin{equation}
\Exp[\rx] = \lambda^{-1}, \qquad \Var[\rx] =\lambda^{-2}. \nonumber
\end{equation}

The support of an exponential distribution is on $(0,\infty)$.
Figure <ref> compares different parameters $\lambda$ for the exponential distribution.

Note that the average $\lambda^{-1}$ is the average time until the occurrence of the event of interest, interpreting $\lambda$   as a rate parameter.
An important property of the exponential distribution is that it is “memoryless," meaning that  the probability of waiting for an additional amount of time $x$  depends only on $x$, not on the past waiting time.

Let $\rx\sim \exponential(\lambda)$. Then we have $p(\rx\geq x + s \mid \rx \geq s) = p(\rx\geq x)$.

Exponential probability density functions for different values of the rate parameter $\lambda$.

Conjugate prior for the exponential rate parameter.
The Gamma distribution is a conjugate prior for the rate parameter of an exponential distribution. To see this, suppose  $\mathcalX=\{x_1, x_2, \ldots, x_N\}$ are drawn i.i.d. from an exponential distribution with rate $\lambda$, i.e., the likelihood is $\exponential(x \mid \lambda)$, and $\lambda$ is given a $\gammadist(\alpha_0, \beta_0)$ prior: $\lambda \sim \gammadist(\alpha_0, \beta_0)$. Using Bayes' theorem, the posterior is

p(\lambda \mid \mathcalX) \propto \prod_{i=1}^{N} \exponential(x_i \mid \lambda)
\times \gammadist(\lambda\mid \alpha_0, \beta_0)
\propto \gammadist(\theta \mid \widetilde{\alpha}, \widetilde{\beta}),

\begin{equation}\label{equation:posterior-param-exponenral-gamma}
\widetilde{\alpha}= \alpha_0+N,
\qquad
\widetilde{\beta} =\beta_0 + \sum_{i=1}^{N} x_i.
\end{equation}

From this posterior form, the prior parameter $\alpha_0$ can be interpreted as the number of prior observations, and $\beta_0$ as the sum of the prior observations. The posterior mean  is

\frac{\widetilde{\alpha}}{\widetilde{\beta}} = \frac{\alpha_0+N}{\beta_0 + \sum_{i=1}^{N} x_i}.

§ UNIVARIATE GAUSSIAN-RELATED MODELS

The truncated-normal (TN) distribution is a variant of the normal distribution,excluding values smaller than zero.
In other words, it is a normal distribution that is “cut off" at zero. The support of the distribution is nonnegative real such that it can be applied for application in a nonnegative matrix factorization context.

Truncated-normal distribution

A random variable $\rx$ is said to follow the truncated-normal distribution with “parent" mean $\mu$ and “parent" precision $\tau>0$, denoted by $\rx \sim \truncatednormal(\mu, \tau^{-1})$, if

$$ f(x; \mu, \tau^{-1})=\left\{
\begin{aligned}
&\frac{\sqrt{\frac{\tau}{2\pi}} \exp \{-\frac{\tau}{2}(x-\mu)^2  \}  }{1-\Phi(-\mu\sqrt{\tau})}
,& \mathrm{\,\,if\,\,} x \geq 0;  \\
&0 , &\mathrm{\,\,if\,\,} x <0,
\end{aligned}
\right.

where $\Phi(y) = \int_{-\infty}^{y} \normal(u\mid 0,1)du= \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{y} \exp(-\frac{u^2}{2}) du $ is the cumulative distribution function (c.d.f.) of $\normal(0,1)$, the standard normal distribution.
Generally, the cumulative density function of $\ry\sim \normal(\mu, \sigma^2)$ can be written as

F(y) = p(\ry\leq y)
= \Phi\left(\frac{y-\mu}{\sigma}\right)
= \Phi\left((y-\mu)\cdot \sqrt{\tau}\right).~
\footnote{
Or equivalently, for general Gaussian distribution $\ry\sim\normal(\mu,\sigma^2)$, the c.d.f. is $F(y)=\frac{1}{2}\left\{
\right\}$, where the error function is $\text{erf}(t)=\frac{2}{\sqrt{\pi}} \int_{0}^{t}\exp(-y^2) dy$.}

The mean and variance of $\rx \sim \truncatednormal(\mu, \tau^{-1})$ are given by

\begin{aligned}
\Exp[\rx] &= \mu - \frac{1}{\sqrt{\tau}}\cdot \frac{ - \phi(\alpha)}{1 - \Phi(\alpha)}, \\
\Var[\rx] &= \frac{1}{\tau}
\left( 1 +
\frac{  \alpha\phi(\alpha)}{1 - \Phi(\alpha)}+
\left(\frac{ \alpha\phi(\alpha)}{1 - \Phi(\alpha)}\right)^2
\right),
\end{aligned}

where $\phi(y)=\frac{1}{\sqrt{2\pi} } \exp(-\frac{y^2}{2})$ is the p.d.f. of the standard normal distribution, and
$\alpha = -\mu\cdot \sqrt{\tau}$ [Burkardt, 2014].
Figure <ref> compares different parameters $\mu, \tau$ for the TN distribution. Figure <ref> shows the mean value of the TN distribution by varying $\mu$ given fixed $\tau$; we can find when $\mu\rightarrow -\infty$, the mean is approaching zero.

Conjugate prior for the nonnegative mean parameter of a Gaussian.
We have previously shown that a Gaussian distribution can serve as  a conjugate prior for the mean parameter of another Gaussian distribution when the variance is fixed.
The truncated-normal distribution  also acts as a conjugate prior for the nonnegative mean parameter of a Gaussian distribution when the variance is fixed. To see this, suppose  $\mathcalX=\{x_1, x_2, \ldots, 0x_N\}$ are drawn i.i.d. from a normal distribution with mean $\theta$ and precision $\tau$, i.e., the likelihood is $\normal(x \mid \theta, \tau^{-1})$ with fixing the variance $\sigma^2=\tau^{-1}$, and $\theta$ is given a $\truncatednormal(\mu_0, \tau^{-1}_0)$ prior: $\theta \sim \truncatednormal(\mu_0, \tau_0^{-1})$. Using Bayes' theorem, the posterior is

\begin{equation}\label{equation:conjugate_truncated_nonnegative_mean}
\begin{aligned}
&\gap p(\theta \mid \mathcalX) \propto \prod_{i=1}^{N} \normal(x_i \mid \theta, \tau^{-1})
\times \truncatednormal(\theta \mid \mu_0, \tau_0^{-1})\\
\exp\left\{ -\frac{\tau_0+ N\tau}{2} \theta^2+ \big(\tau \sum_{i=1}^{N}x_i + \tau_0 \mu_0\big)\theta  \right\}\cdot u(\theta)\\
&\propto \truncatednormal(\theta \mid \widetilde{\mu}, \widetilde{\tau}^{-1}),
\end{aligned}
\end{equation}

where $u(y)$ is the step function with value 1 if $y\geq 0$ and value 0 if $y<0$, and

\widetilde{\mu}= \frac{\tau_0 \mu_0 + \tau \sum_{i=1}^{N} x_i}{ \tau_0 + N\tau},
\qquad
\widetilde{\tau} = \tau_0 + N \tau.

The posterior parameters are exactly the same as those in the Normal-Normal model (Equation (<ref>)), and the posterior “parent" mean can also be expressed as a weighted mean of $\mu_0$ and $\overline{x}$.

Truncated-normal  and general-truncated-normal  probability density functions for different values of the parameters $\mu$ and $\tau$.

Mean of truncated-normal and general-truncated-normal distribution by varying $\mu, \tau, a$, and $b$ parameters.

Expanding beyond   the truncated-normal distribution,
the general-truncated-normal (GTN) distribution is also a variant of the normal distribution, excluding values outside a specified range. In other words, it is a normal distribution that is “cut off" at some specified lower and/or upper bound. The range between the lower and upper bound is called the support of the distribution. [In some contexts, the general-truncated-normal distribution is named the truncated-normal directly. We here differentiate the two as the definitions describe.]

General-truncated-normal distribution

A random variable $\rx$ is said to follow the general-truncated-normal distribution with “parent" mean $\mu$ and “parent" precision $\tau>0$, denoted by $\rx \sim \generaltruncatednormal(\mu, \tau^{-1}, \textcolor{mylightbluetext}{a, b})$, if

$$ f(x; \mu, \tau^{-1}, a, b)=\left\{
\begin{aligned}
&0,  & \mathrm{\,\,if\,\,} x < a;  \\
&\frac{\sqrt{\frac{\tau}{2\pi}} \exp \{-\frac{\tau}{2}(x-\mu)^2  \}  }{\Phi((b-\mu)\cdot \sqrt{\tau})-\Phi((a-\mu)\cdot \sqrt{\tau})}
,& \mathrm{\,\,if\,\,} a\leq x \leq b;  \\
&1 , &\mathrm{\,\,if\,\,} 0 >b,
\end{aligned}
\right.

where $\Phi(\cdot)$ is the c.d.f. of $\normal(0,1)$, the standard normal distribution.
The mean and variance of $\rx \sim \generaltruncatednormal(\mu, \tau^{-1}, a, b)$ are given by

\begin{aligned}
\Exp[\rx] &= \mu - \frac{1}{\sqrt{\tau}}\cdot \frac{\phi(\beta) - \phi(\alpha)}{\Phi(\beta) - \Phi(\alpha)}, \\
\Var[\rx] &= \frac{1}{\tau}
\left( 1 -
\frac{\beta \phi(\beta) - \alpha\phi(\alpha)}{\Phi(\beta) - \Phi(\alpha)}-
\left(\frac{\beta \phi(\beta) - \alpha\phi(\alpha)}{\Phi(\beta) - \Phi(\alpha)}\right)^2
\right),
\end{aligned}

where $\phi(\cdot)$ is the p.d.f. of the standard normal distribution, and

\alpha = (a-\mu)\cdot \sqrt{\tau}, \qquad
\beta = (b-\mu)\cdot \sqrt{\tau}.

Note that, the truncated-normal distribution is a special general-truncated-normal with $a=0$ and $b=\infty$ [Burkardt, 2014].
Figure <ref> compares different parameters $\mu, \tau$ for the GTN distribution.
Figure <ref> shows the mean value of the GTN distribution by varying $\mu$ given fixed $\tau, a, b$; we again find when $\mu\rightarrow -\infty$, the mean is approaching zero.

Conjugate prior for the constrained mean parameter of a Gaussian.
We have previously shown that a Gaussian distribution can serve as a conjugate prior for the mean parameter of another Gaussian distribution when the variance is fixed.
The truncated-normal distribution is also a conjugate prior for the nonnegative mean parameter of a Gaussian distribution when the variance is fixed. To see this, suppose  $\mathcalX=\{x_1, x_2, \ldots, x_N\}$ are drawn i.i.d. from a  normal distribution with mean $\theta$ and precision $\tau$, i.e., the likelihood is $\normal(x \mid \theta, \tau^{-1})$, where the variance $\sigma^2=\tau^{-1}$ is fixed, and $\theta$ is given a $\generaltruncatednormal(\mu_0, \tau^{-1}_0, a, b)$ prior: $\theta \sim \generaltruncatednormal( \mu_0, \tau_0^{-1}, a, b)$, with $a>0$ for the nonnegativity constrain. Using Bayes' theorem, the posterior is

\begin{equation}\label{equation:conjugate_truncated_constrained_mean}
\begin{aligned}
&\gap p(\theta \mid \mathcalX) \propto \prod_{i=1}^{N} \normal(x_i \mid \theta, \tau^{-1})
\times \generaltruncatednormal(\theta \mid \mu_0, \tau_0^{-1}, a, b)\\
&\propto \exp\left\{ -\frac{\tau_0+ N\tau}{2} \theta^2
+ \big(\tau \sum_{i=1}^{N}x_i + \tau_0 \mu_0\big)\theta  \right\}
\cdot \indicator(a\leq \theta\leq b)\\
&\propto \generaltruncatednormal(\theta \mid \widetilde{\mu}, \widetilde{\tau}^{-1}, a, b),
\end{aligned}
\end{equation}

where $\indicator(a\leq y\leq b)$ is the step function with value 1 if $ a\leq y\leq b$ and value 0 otherwise, and

\widetilde{\mu}= \frac{\tau_0 \mu_0 + \tau \sum_{i=1}^{N} x_i}{ \tau_0 + N\tau},
\qquad
\widetilde{\tau} = \tau_0 + N \tau.

The posterior parameters are again exactly the same as those in the Normal-Normal model (Equation (<ref>)).

[Normal, $\tau$ is the precision parameter.]

Half-normal and normal probability density functions for different values of the parameters $\mu$ and $\tau$.
The probability density of any value $x\geq \mu$ in the half-normal distribution is twice as that in the normal distribution with the same parameters $\mu, \tau$.

The half-normal distribution is a special case of the normal distribution, where the support is on the range with values larger than the mean parameter $\mu$ of a normal distribution (known as the “parent" mean parameter of the normal distribution) and the distribution is unsymmetrical around $\mu$.

This distribution is often used in modeling the scale or standard deviation of a process where the values cannot be smaller than $\mu$.
Half-normal distribution

A random variable $\rx$ is said to follow the half-normal distribution with “parent" mean $\mu$ and “parent" precision $\tau>0$, denoted by $\rx \sim \halfnormal(\mu, \tau^{-1})$, if

$$ f(x; \mu, \tau^{-1})=\left\{
\begin{aligned}
&\sqrt{\frac{2\tau}{\pi}} \exp\left\{- \frac{\tau}{2}(x-\mu)^2 \right\}
,& \mathrm{\,\,if\,\,} x \geq \mu;  \\
&0 , &\mathrm{\,\,if\,\,} x <\mu.
\end{aligned}
\right.

The mean and variance of $\rx \sim \halfnormal(\mu, \tau^{-1})$ are given by

\Exp[\rx] = \mu + \sqrt{\frac{2}{\pi \tau}} , \qquad
\Var[\rx] = \frac{1}{\tau}\big(1-\frac{2}{\pi} \big).

Figure <ref> compares different parameters $\mu, \tau$ for the half-normal distribution.

Rectified-normal distribution
Exponentially rectified-normal distribution

In this text, the rectified-normal (RN) distribution is defined as proportional to the product of a Gaussian distribution and an exponential distribution. And it can also be called the exponentially rectified-normal distribution.

A random variable $\rx$ is said to follow the rectified-normal distribution (or exponentially rectified-normal distribution) with “parent" mean $\mu$, “parent" precision $\tau>0$, and “parent" rate $\lambda >0$, denoted by $\rx \sim \rectifieddist(\mu, \tau^{-1}, \lambda)$, if

\begin{aligned}
f(x; \mu, \tau^{-1}, \lambda)
&=\frac{1}{C}\cdot \normal(x\mid  \mu, \tau^{-1}) \cdot \exponential(x\mid \lambda)\\
&\propto \exp\left\{ -\frac{\tau}{2} \left( x - \frac{\tau\mu-\lambda}{\tau} \right)^2 \right\} \cdot u(x)\\
&\propto \truncatednormal(x\mid\textcolor{mylightbluetext}{\frac{\tau\mu-\lambda}{\tau}}, \tau^{-1}),
\end{aligned}

where $\truncatednormal(\cdot)$ is the density function of a truncated-normal distribution, and $C$ is a constant  value,

\begin{equation}\label{equation:rf_constant}
C =C^{RN}(\mu, \tau, \lambda)
= \lambda \left\{1 - \Phi\bigg(-\frac{\tau\mu-\lambda}{\sqrt{\tau}} \bigg)\right\}
\cdot \exp\left( -\mu\lambda + \frac{\lambda^2}{2\tau}\right).
\end{equation}

That is, the rectified-normal distribution is a special truncated-normal distribution with more flexibility.
The mean and variance of $\rx \sim \rectifieddist(\mu, \tau^{-1}, \lambda)$ are given by

\begin{aligned}
\Exp[\rx] &= \frac{\tau\mu-\lambda}{\tau} - \frac{1}{\sqrt{\tau}}\cdot \frac{ - \phi(\alpha)}{1 - \Phi(\alpha)}, \\
\Var[\rx] &= \frac{1}{\tau}
\left\{ 1 +
\frac{  \alpha\phi(\alpha)}{1 - \Phi(\alpha)}+
\left(\frac{ \alpha\phi(\alpha)}{1 - \Phi(\alpha)}\right)^2
\right\},
\end{aligned}

\alpha = -\frac{\tau\mu-\lambda}{\tau}\cdot \sqrt{\tau}.
Figure <ref> compares different parameters $\mu, \tau$ for the RN distribution.

The comparison between the truncated-normal and rectified-normal distributions is presented in Figure <ref>.

Conjugate prior for the nonnegative mean parameter of a Gaussian by RN.
Similar to the TN distribution, the RN distribution also serves to enforce nonnegative constraint, and is conjugate to the nonnegative mean parameter of a Gaussian likelihood. However, due to the extra parameter $\lambda$, the RN distribution is more flexible in this sense. And the derivation follows from Equation (<ref>).

To see this, suppose  $\mathcalX=\{x_1, x_2, \ldots, x_N\}$ are drawn i.i.d. from a  normal distribution with mean $\theta$ and precision $\tau$, i.e., the likelihood is $\normal(x \mid \theta, \tau^{-1})$ where the variance $\sigma^2=\tau^{-1}$ is fixed, and $\theta$ is given a $\rectifieddist(\mu_0, \tau^{-1}_0, \lambda_0)$ prior: $\theta \sim \rectifieddist( \mu_0, \tau_0^{-1}, \lambda_0)$. Using Bayes' theorem, the posterior is

\begin{equation}\label{equation:conjugate_rectified_nonnegative_mean}
\begin{aligned}
&\gap p(\theta \mid \mathcalX) \propto \prod_{i=1}^{N} \normal(x_i \mid \theta, \tau^{-1})
\times \rectifieddist(\theta\mid \mu_0, \tau_0^{-1}, \lambda_0) \\
&= \prod_{i=1}^{N} \normal(x_i \mid \theta, \tau^{-1})
\times \truncatednormal(\theta\mid m_0, \tau_0^{-1})
% &\propto \exp\left( -\frac{\tau_0+ N\tau}{2} \theta^2+ (\tau \sum_{i=1}^{N}x_i + \tau_0 \mu_0)\theta  \right)\cdot u(\theta)\\
&\propto \truncatednormal(\theta \mid \widetilde{\mu}, \widetilde{\tau}^{-1}),
\end{aligned}
\end{equation}

\widetilde{\mu}= \frac{\tau_0 m_0 + \tau \sum_{i=1}^{N} x_i}{ \tau_0 + N\tau},
\qquad
\widetilde{\tau} = \tau_0 + N \tau.

That is, the posterior density is a special RN or TN distribution.

[Truncated-normal. Same as Figure <ref>.]

Truncated-normal and rectified-normal probability density functions for different values of the parameters $\mu$, $\tau$, and $\lambda$.

Inverse-Gaussian distribution
The inverse-Gaussian distribution, also known as the Wald distribution, is a continuous probability distribution with two parameters, $\mu > 0$ and $b > 0$.
It is a versatile distribution that is used in various applications including modeling waiting times, stock prices, and lifetimes of mechanical systems.
An important property of the distribution is that it is well-suited for modeling nonnegative, continuous, and positively skewed data with finite mean and variance.

A random variable $\rx$ is said to follow the inverse-Gaussian distribution with parameters $\mu>0$ and $b>0$, denoted by $\rx \sim \inversenormaldist(\mu,\lambda)$ [Note we use $\inversenormaldist$ to denote the inverse-Gaussian distribution and use $\inversegammadist$ to denote the inverse-Gamma distribution.], if

f(x; \mu,\lambda)=\left\{
\begin{aligned}
& \sqrt{\frac{\lambda}{2\pi x^3}} \exp \left( -\frac{\lambda(x-\mu)^2}{2\mu^2 x} \right)
,& \mathrm{\,\,if\,\,} x \geq 0;  \\
&0 , &\mathrm{\,\,if\,\,} x <0.
\end{aligned}
\right.

The mean and variance of $\rx \sim \inversenormaldist( \mu,\lambda)$ are given by

\begin{equation}
\Exp[\rx] = \mu, \qquad \Var[\rx] = \frac{\mu^3}{\lambda}. \nonumber
\end{equation}

The support of an inverse-Gaussian distribution is on $(0,\infty)$.
Figure <ref> compares different parameters $\mu, \lambda$ for the inverse-Gaussian distribution.

Inverse-Gaussian probability density functions for different values of the parameters $\mu$ and $\lambda$.

The Laplace distribution, also known as  the double exponential distribution, is named after Pierre-Simon Laplace (1749–1827), who obtained the distribution in 1774 [Kotz et al., 2001, Härdle and Simar, 2007].
The Laplace distribution is useful in modeling heavy-tailed data since it has heavier tails than the normal distribution, and it is used extensively in sparse-favoring models since it expresses a high peak with heavy tails (same as the $\ell_1$ regularization term in non-probabilistic or non-Bayesian optimization methods).
When we have a prior  belief that the parameter of interest is likely to be close to the mean  with the potential for large deviations, the Laplace distribution is then used in Bayesian modeling as a prior distribution for this context.
Laplace distribution

A random variable $\rx$ is said to follow the Laplace distribution with location and scale parameters $\mu$ and $b>0$, respectively, denoted by $\rx \sim \laplacedist(\mu,b)$, if

f(x; \mu,b)=\frac{1}{2b} \exp \left( -\frac{\abs{x-\mu}}{b } \right).

The mean and variance of $\rx \sim \laplacedist( \mu,b)$ are given by

\begin{equation}
\Exp[\rx] = \mu, \qquad \Var[\rx] =2b^2. \nonumber
\end{equation}

Figure <ref> compares different parameters $\mu$ and $b$ for the Laplace distribution.

[Laplace distribution.]

[Skew-Laplace distribution.]

Laplace and skew-Laplace probability density functions for different values of the parameters.

Laplace as a mixture of normal distributions.