Estimates

\[\int_{-\inf}^{\inf} \text{ExpFam}(\theta, g(x)) \mathcal{N} (x | \mu, \sigma^2) \mathrm d x\]

via a very fast and deterministic numerical integration.

Install¶

You can install it via conda

conda install -c conda-forge liknorm

or by building it yourself via the following command:

# DO_CMD=sudo
curl -fsSL https://git.io/JerYI | GITHUB_USER=horta GITHUB_PROJECT=logaddexp bash
curl -fsSL https://git.io/JerYI | GITHUB_USER=limix GITHUB_PROJECT=liknorm bash

The above command should work on Windows, Linux, and MacOS.

Introduction¶

The single-parameter exponential family is a class of distributions that can be expressed as:

f(y; θ, 𝜙) = exp{(yθ - b(θ))/a(𝜙) + c(y,𝜙)}.

for which 𝜙 is assumed to be known. The definition of the functions a(.), b(.), and c(.) determines a probabilistic distribution having the canonical parameter θ. The expectation of y, denoted here by 𝜇, determines the value of θ via the following relation:

b'(θ) = 𝜇

Still, the value 𝜇 is often set indirectly via the natural parameter η, which relates to each other through a link function g(.):

η = g(𝜇)

If g(.) is the so-called canonical function, we have the desirable equality:

θ = η

Usage example¶

Suppose you have the file

/* example.c */
#include "liknorm.h"

#include <stdio.h>

int main()
{
  double log_zeroth, mean, variance;
  double prior_var = 2.5;
  double prior_mean = -2.0;
  double nsuccesses = 2;
  double ntrials = 15;

  struct LikNormMachine *machine = liknorm_create_machine(500);

  liknorm_set_binomial(machine, nsuccesses, ntrials);
  liknorm_set_prior(machine, 1 / prior_var, prior_mean / prior_var);

  liknorm_integrate(machine, &log_zeroth, &mean, &variance);

  printf("%f\n", log_zeroth);
  printf("%f\n", mean);
  printf("%f\n", variance);

  liknorm_destroy_machine(machine);
}

Compiling, linking, and running it via

cc libliknorm.a example.c -o example
./example

should print:

-2.049961
-2.038184
0.524308

Interface¶

Liknorm machine¶

LIKNORM_API void liknorm_integrate(struct LikNormMachine * machine, double * log_zeroth, double * mean, double * variance)

Perform numerical integration.

Parameters

machine: Machine to perform integration.
log_zeroth: Zeroth moment.
log_mean: First moment of the normalized distribution.
log_variance: Variance of the normalized distribution.

LIKNORM_API void liknorm_destroy_machine(struct LikNormMachine * machine)

Destroy a Machine instance.

Parameters

machine: Machine to be destroyed. Always call it before exiting your program, otherwise it will leak memory.

Prior¶

LIKNORM_API void liknorm_set_prior(struct LikNormMachine * machine, double tau, double eta)

Set the natural parameters of Normal prior.

Parameters

machine: Machine to perform integration.
tau: It equals to σ⁻².
eta: It equals to μσ⁻².

Bernoulli¶

LIKNORM_API void liknorm_set_bernoulli(struct LikNormMachine * machine, double k)

Bernoulli distribution.

It is the discrete probability distribution of a random variable which takes the value 1 with probability p and the value 0 with probability 1 − p. (Wikipedia.)

Parameters

machine: Liknorm handler.
k: Number of successes.

Binomial¶

LIKNORM_API void liknorm_set_binomial(struct LikNormMachine * machine, double k, double n)

Binomial distribution.

It is the discrete probability distribution of the number of successes k in a sequence of n independent experiments. (Wikipedia.) The probability mass function is given by:

Binom(k, n) pᵏ (1 - p)ⁿ⁻ᵏ,

for which Binom(m, n) = n! / (m! (n - m)!).

Parameters

machine: Liknorm handler.
k: Number of successes.
n: Number of trials.

Gamma¶

LIKNORM_API void liknorm_set_gamma(struct LikNormMachine * machine, double x, double a)

Gamma distribution.

A gamma distribution is a general type of statistical distribution that is related to the beta distribution and arises naturally in processes for which the waiting times between Poisson distributed events are relevant. (Wolfram.)

Parameters

machine: Liknorm handler.
x: Waiting time.
a: Shape parameter α.

Geometric¶

LIKNORM_API void liknorm_set_geometric(struct LikNormMachine * machine, double x)

Geometric distribution.

Parameters

machine: Liknorm handler.
x: Number of trials to success.

Negative binomial¶

LIKNORM_API void liknorm_set_nbinomial(struct LikNormMachine * machine, double k, double r)

Negative binomial distribution.

It is a discrete probability distribution of the number of successes in a sequence of independent and identically distributed Bernoulli trials before a specified (non-random) number of failures (denoted r) occurs. (Wikipedia.) Let k be the number of successes. The probability mass function is given by:

Binom(k, k + r - 1) pᵏ (1 - p)ʳ,

for which Binom(m, n) = n! / (m! (n - m)!).

Parameters

machine: Liknorm handler.
k: Number of successes.
r: Number of failures.

Poisson¶

LIKNORM_API void liknorm_set_poisson(struct LikNormMachine * machine, double k)

Poisson distribution.

It is a discrete probability distribution that expresses the probability of a given number of events occurring in a fixed interval of time or space if these events occur with a known constant rate and independently of the time since the last event. (Wikipedia.) The probability mass function is given by:

λᵏe^{-λ} / k!,

for which λ is the rate of occurrence and k the number of occurrences.

Parameters

machine: Liknorm handler.
k: Number of occurrences.

Implementation¶

ExpFam¶

struct ExpFam¶

Exponential family of distributions.

We adopt the following representation::

f(y; θ, 𝜙) = exp{(yθ - b(θ))/a(𝜙) + c(y,𝜙)},

for which::

y     : random variable value.
θ     : canonical parameter.
𝜙     : nuisance parameter.
a(𝜙)  :
b(θ)  : log-partition function.
c(y,𝜙): normaliser.

The mean and variance are given by::

E[y]   = b'(θ)
Var[y] = b''(θ)a(𝜙)

In order to define a generalised linear mixed model (GLMM) we use the so-called natural parameter η. Given a link function g(.), the natural parameter relates to the canonical parameter as follows::

η = g(E[y]) = g(b'(θ)).

Every member of the exponential family has a canonical link function, which greatly simplifies the relationship::

η = θ

Public Members

double y¶: Random variable value

double a¶: a(𝜙)

double loga¶: log(a(𝜙))

double c¶: c(y,𝜙)

log_partition *lp¶: b(θ)

log_partition_fderivative *lpfd¶: log(b'(θ))

log_partition_derivatives *lpd¶: log(b''(θ))

Likelihood¶

We assume the canonical link function for every likelihood.

Bernoulli¶

y assumes 1 or 0 for failure. We make use of the Binomial implementation. So, please, refer to the next section for details.

static double bernoulli_log_partition(const double theta)¶

Bernoulli log-partition function.

Please, refer to the binomial_log_partition() function.

static double bernoulli_log_partition_fderivative(const double theta)¶

First derivative of the Bernoulli log-partition function.

Please, refer to the binomial_log_partition_fderivative() function.

static void bernoulli_log_partition_derivatives(const double theta, double *b0, double *logb1, double *logb2)¶

Zeroth, first, and second derivatives of the Bernoulli log-partition function.

Please, refer to the bernoulli_log_partition_fderivative() function.

Binomial¶

The random variable is given by y = k/n. The support is therefore y ϵ {0/n, 1/n, ..., r/n}. The exponential family functions are:

𝜙      = n
a(𝜙)   = 1/𝜙
b(θ)   = log(1 + exp(θ))
c(y,𝜙) = log(binom(n, y𝜙))

Let us define:

𝜇 = E[y] = p.

The canonical link function and its inverse are given by:

canonical(𝜇)     = log(𝜇/(1+𝜇)) = η
canonical_inv(η) = 1/(1 + exp(-η))

double binomial_log_partition(const double theta)¶

Binomial log-partition function.

Definition:

b(𝜃) = log(1 + exp(𝜃)).

double binomial_log_partition_fderivative(const double theta)¶

First derivative of the Binomial log-partition function.

Definition:

log(b'(𝜃)) = 𝜃 - log(1 + exp(𝜃))

void binomial_log_partition_derivatives(const double theta, double *b0, double *logb1, double *logb2)¶

Zeroth, first, and second derivatives of the Binomial log-partition function.

Implements b(𝜃), log(b'(𝜃)), and:

log(b''(𝜃)) = 𝜃 - 2log(1 + exp(𝜃))

Negative Binomial¶

The random variable is given by y = k/r. The support is therefore y ϵ {0/r, 1/r, ..., r/r}. The exponential family functions are:

𝜙 = r
a(𝜙) = 1/𝜙
b(θ) = -log(1 - exp(θ))
c(y,𝜙) = log(binom(y𝜙 + 𝜙 - 1, y𝜙))

Let us define:

𝜇 = E[y] = p / (1 - p)

The canonical link function and its inverse are given by:

canonical(𝜇)     = log(𝜇 / (1 + 𝜇)) = η
canonical_inv(η) = exp(η) / (1 - exp(η))

double nbinomial_log_partition(const double theta)¶

Negative binomial log-partition function.

Definition:

b(𝜃) = -log(1 - exp(𝜃)).

double nbinomial_log_partition_fderivative(const double theta)¶

First derivative of the Negative Binomial log-partition function.

Definition:

log(b'(𝜃)) = 𝜃 - log(1 - exp(𝜃)).

void nbinomial_log_partition_derivatives(const double theta, double *b0, double *logb1, double *logb2)¶

Zeroth, first, and second derivatives of the Negative Binomial log-partition func.

Implements b(𝜃), log(b'(𝜃)), and:

log(b''(𝜃)) = 𝜃 - 2log(1 - exp(𝜃))

Poisson¶

The support is y ϵ {0, 1, ...}. The exponential family functions are:

𝜙      = 1
a(𝜙)   = 𝜙
b(𝜃)   = exp(𝜃)
b'(𝜃)  = exp(𝜃)
b'(𝜃)  = exp(𝜃)
c(y,𝜙) = -log(y!)

Let us define:

𝜇 = E[y] = λ,

for which λ is the Poisson distribution parameter. The canonical link function and its inverse are given by:

canonical(𝜇)     = log(𝜇 / (1 + 𝜇)) = η
canonical_inv(η) = exp(η) / (1 - exp(η))

double poisson_log_partition(const double theta)¶

Poisson log-partition function.

Definition:

b(𝜃) = exp(𝜃)

double poisson_log_partition_fderivative(const double theta)¶

Log of the first derivative of the Poisson log-partition function.

Definition:

log(b'(𝜃)) = 𝜃

void poisson_log_partition_derivatives(const double theta, double *b0, double *logb1, double *logb2)¶

Log of the derivatives of the Poisson log-partition function.

Implements b(𝜃), log(b'(𝜃)), and:

log(b''(𝜃)) = 𝜃

Comments and bugs¶

You can get the source and open issues on Github.