Last time[1] we took a look at how we could define multivariate normally distributed random variables with linear functions of multiple independent standard univariate normal random variables.
Specifically, given a vector \(\mathbf{Z}\) whose elements are independent standard univariate normal random variables, a constant vector \(\mathbf{\mu}\) and a constant matrix \(\mathbf{L}\)
We got as far as deducing the characteristic function and the probability density function of the multivariate normal distribution, leaving its cumulative distribution function and its complement aside until we'd implemented both them and the random variable itself, which we shall do in this post.
Here
Now its possible, albeit rather unlikely, for a standard normal random variable to return infinite values and this might result in trying to add a positive infinity to a negative infinity or multiply an infinity by zero, yielding NaN. We therefore arbitrarily replace NaNs in the final result with negative infinities. Unfortunately, this means that if
This isn't a very efficient implementation for uncorrelated multivariate normal random variables, however, and so listing 2 provides the
This time
Note that we're adding property access methods to the resulting random number generator
As usual we are dispatching to a
For example, we have already determined that we must populate
For another, we shall often be interested in multivariate normal random variables with means of zero and it would be nice to simply not specify the mean in such cases, having it default to zero.
Now our default
One thing that we should certainly like to be able to do easily is to construct \(n\) dimensional independent standard multivariate normal random variables, like the \(\mathbf{Z}\) that we used at the start of the post. Listing 5 gives constructors which take the number of dimensions and an optional random number generator to do just that.
Note that we require that the number of dimensions is an integer, throwing an error if it is not.
We might also like to create \(n\) dimensional random variables that all have a particular variance, which the constructors given in listing 6 enable with two number arguments and, once again, an optional random number generator.
Here we ensure that the variance is both strictly positive and finite before taking its square root to yield the standard deviation, filling the
To specify a mean as well we pass a third number to invoke the constructors given in listing 7.
This time the first number specifies the number of dimensions, the second the mean and the third the variance. Apart from this they're more or less the same as the previous constructors except that they also ensure that the mean is finite and fill
As noted above it's handy to be able to define an independent multivariate normal random variable with different means and variances for each element by passing them as vectors. The constructors given in listing 8 allow this using
Now we no longer need to explicitly state the number of dimensions since it is implicitly defined by the vector of variances. The main difference between these and the previous constructors is that we must check the validity of each mean and variance, and take the square root of the latter separately.
Note that once again we also allow the means to default to zero if we only pass the variances to the constructor, although we shan't show the implementation here.
With all of the uncorrelated constructors out of the way, we're ready to look at the construction of correlated multivariate normal random variables, although we first need a helper function to check whether a covariance matrix is diagonal and so implies no correlations, as given in listing 9.
Listing 10 gives the constructor that takes an
So first we check that
Then, if
Otherwise we use our
Note that since the
We also have a few last constructors, not listed here, for the covariance matrix without the mean vector and for the Cholesky decomposition of the covariance matrix, either as an
So finally, having spent so much effort on the initialisation of multivariate normal random variables, we are ready to reap our glorious reward of an ever so slightly more convenient way of drawing from them, as demonstrated by program 1.
As we did for the multivariate uniform distribution[2] we shall instead define a function that uniquely maps a vector whose elements all lie between zero and one onto another vector in such a way that if the former's elements are all independent standard uniformly distributed then the latter's will be multivariate normally distributed.
Specifically, we'll just use the elements of the vector as the source of randomness in our definition of multivariate random variables by using the univariate normal inverse CDF to transform them into observations of a standard univariate normal random variable.
Listing 11 show how we do this for an uncorrelated multivariate normal distribution with
After checking that the latter is indeed a vector and has the correct number of dimensions, we have simply substituted each call to
The implementation of the correlated multivariate normal map requires precisely the same changes to the
Finally, we can implement an
One way that we can use
Note that the
For the correlated case we have
This time
Finally, the
Program 3 demonstrates the use of
Note that we're exploiting the identity
The
Finally, the
Program 4 demonstrates the use of
Note that all of these functions can be found in
[2] Into The nth Dimension, www.thusspakeak.com, 2016.
[3] Halton, Here's A Pseud, www.thusspakeak.com, 2016.
Specifically, given a vector \(\mathbf{Z}\) whose elements are independent standard univariate normal random variables, a constant vector \(\mathbf{\mu}\) and a constant matrix \(\mathbf{L}\)
\[
\mathbf{Z}^\prime = \mathbf{L} \times \mathbf{Z} + \mathbf{\mu}
\]
has linearly dependent normally distributed elements, a mean vector of \(\mathbf{\mu}\) and a covariance matrix of
\[
\mathbf{\Sigma} = \mathbf{L}\times\mathbf{L}^\mathrm{T}
\]
where \(\mathbf{L}^\mathrm{T}\) is the transpose of \(\mathbf{L}\) in which the rows and columns are switched.We got as far as deducing the characteristic function and the probability density function of the multivariate normal distribution, leaving its cumulative distribution function and its complement aside until we'd implemented both them and the random variable itself, which we shall do in this post.
ak.multiNormalRnd
Given how simple the definition of a multivariate normal random variable is, it's as good a place as any to start. Listing 1 shows thecorrRnd
function which generates them according to the formula given above.Listing 1: The corrRnd Helper Function
function corrRnd(mu, root, rnd) { var n = mu.length; var z = new Array(n); var i, j, rooti; for(i=0;i<n;++i) z[i] = rnd(); if(z.some(isNaN)) return ak.vector(n, ak.NaN); for(i=n-1;i>=0;--i) { rooti = root[i]; z[i] = mu[i] + rooti[i]*z[i]; for(j=0;j<i;++j) z[i] += rooti[j]*z[j]; if(isNaN(z[i])) z[i] = -ak.INFINITY; } return ak.vector(z); }
Here
rnd
is a standard univariate normal random variable, mu
is an array containing the elements of the vector \(\mathbf{\mu}\) and root
is an array of arrays containing the elements of the matrix \(\mathrm{L}\) which is assumed to be lower triangular having its elements in the top right whose row indices are smaller than their column indices all equal to zero. This means that we can transform the array of standard normal random variables in place if we work backwards since the elements of the multivariate normal random variable only depend upon those of the former whose indices are no greater than theirs.Now its possible, albeit rather unlikely, for a standard normal random variable to return infinite values and this might result in trying to add a positive infinity to a negative infinity or multiply an infinity by zero, yielding NaN. We therefore arbitrarily replace NaNs in the final result with negative infinities. Unfortunately, this means that if
rnd
erroneously returned a NaN we'd translate terms involving it into minus infinity too and so we first check to see whether any of the standard normal values are NaN and, if so, immediately return an ak.vector
filled with NaNs.This isn't a very efficient implementation for uncorrelated multivariate normal random variables, however, and so listing 2 provides the
uncorrRnd
helper function specifically for them.Listing 2: The uncorrRnd Helper Function
function uncorrRnd(mu, root, rnd) { var n = mu.length; var z = new Array(n); var i; for(i=0;i<n;++i) z[i] = mu[i] + root[i]*rnd(); return ak.vector(z); }
This time
root
is an array containing the elements of the vector of standard deviations \(\mathbf{\sigma}\) of the elements of the uncorrelated multivariate normal random variable which we define with
\[
Z^\prime_i = \sigma_i Z_i + \mu_i
\]
The ak.multiNormalRnd
function given in listing 3 creates a random generator that forwards to one or other of these functions depending upon whether the sigma
member of the state
object is an array of numbers or not after initialisation, the former case holding the variances of an uncorrelated random variable and the latter the elements of the covariance matrix of a correlated random variable.Listing 3: ak.multiNormalRnd
ak.multiNormalRnd = function() { var state = {mu: [0], sigma: [1], root: [1], rnd: Math.random}; var arg0 = arguments[0]; var rnd, f; constructors[ak.type(arg0)](state, arg0, arguments); rnd = ak.normalRnd(state.rnd); if(state.sigma.length===0 || ak.nativeType(state.sigma[0])===ak.NUMBER_T) { f = function(){return uncorrRnd(state.mu, state.root, rnd);}; f.sigma = function(){return ak.matrix('diagonal', state.sigma);}; } else { f = function(){return corrRnd(state.mu, state.root, rnd);}; f.sigma = function(){return ak.matrix(state.sigma);}; } f.mu = function(){return ak.vector(state.mu);}; f.rnd = function(){return state.rnd;}; return Object.freeze(f); }; var constructors = {};
Note that we're adding property access methods to the resulting random number generator
f
, with sigma
returning an ak.matrix
representing the covariance matrix, mu
an ak.vector
representing the mean and rnd
returning the underlying random number generator, which we've used to initialise the standard univariate ak.normalRnd
generator rnd
. In the case of uncorrelated random variables, we're converting the array of variances to a diagonal matrix having zeros everywhere except on the leading diagonal where the row and column indices are equal.As usual we are dispatching to a
constructors
object to perform the initialisation but unfortunately its implementation involves a lot more work than the random number generator itself.
Multivariate Normal Distribution Constructors
Thankfully there's nothing particularly difficult to do, its just that there are rather a lot of cases to consider.For example, we have already determined that we must populate
state.sigma
and state.root
with arrays of variances and standard deviations respectively in the uncorrelated case, and with arrays of arrays of the elements of the covariance matrix and the lower triangular \(\mathbf{L}\) in the correlated case. It will consequently be rather convenient to initialise uncorrelated multivariate normal random variables with vectors of variances as well as with diagonal covariance matrices.For another, we shall often be interested in multivariate normal random variables with means of zero and it would be nice to simply not specify the mean in such cases, having it default to zero.
Now our default
state
defines a one dimensional random variable with a mean of zero, a variance of one, a standard deviation of one and uses Math.random
as the underlying random number generator, and the simplest constructors, given in listing 4, are those that leave the mean, variance and standard deviation at their default values.Listing 4: The Trivial Constructors
constructors[ak.UNDEFINED_T] = function() { }; constructors[ak.FUNCTION_T] = function(state, rnd) { state.rnd = rnd; };
One thing that we should certainly like to be able to do easily is to construct \(n\) dimensional independent standard multivariate normal random variables, like the \(\mathbf{Z}\) that we used at the start of the post. Listing 5 gives constructors which take the number of dimensions and an optional random number generator to do just that.
Listing 5: The Uncorrelated Standard Constructors
constructors[ak.NUMBER_T] = function(state, n, args) { var arg1 = args[1]; if(n!==ak.floor(n)) { throw new Error('invalid dimension in ak.multiNormal distribution'); } state.mu.length = n; state.sigma.length = n; state.root.length = n; constructors[ak.NUMBER_T][ak.nativeType(arg1)](state, n, arg1, args); }; constructors[ak.NUMBER_T][ak.UNDEFINED_T] = function(state, n) { while(n--) { state.mu[n] = 0; state.sigma[n] = 1; state.root[n] = 1; } }; constructors[ak.NUMBER_T][ak.FUNCTION_T] = function(state, n, rnd) { while(n--) { state.mu[n] = 0; state.sigma[n] = 1; state.root[n] = 1; } state.rnd = rnd; };
Note that we require that the number of dimensions is an integer, throwing an error if it is not.
We might also like to create \(n\) dimensional random variables that all have a particular variance, which the constructors given in listing 6 enable with two number arguments and, once again, an optional random number generator.
Listing 6: The Identically Distributed Zero Mean Uncorrelated Constructors
constructors[ak.NUMBER_T][ak.NUMBER_T] = function(state, n, arg1, args) { var arg2 = args[2]; constructors[ak.NUMBER_T][ak.NUMBER_T][ak.nativeType(arg2)](state,n,arg1,arg2,args); }; constructors[ak.NUMBER_T][ak.NUMBER_T][ak.UNDEFINED_T] = function(state, n, sigma) { var root; sigma = Number(sigma); if(sigma<=0 || !isFinite(sigma)) { throw new Error('invalid sigma in ak.multiNormal distribution'); } root = Math.sqrt(sigma); while(n--) { state.mu[n] = 0; state.sigma[n] = sigma; state.root[n] = root; } }; constructors[ak.NUMBER_T][ak.NUMBER_T][ak.FUNCTION_T] = function(state, n, sigma, rnd) { var root; sigma = Number(sigma); if(sigma<=0 || !isFinite(sigma)) { throw new Error('invalid sigma in ak.multiNormal distribution'); } root = Math.sqrt(sigma); while(n--) { state.mu[n] = 0; state.sigma[n] = sigma; state.root[n] = root; } state.rnd = rnd; };
Here we ensure that the variance is both strictly positive and finite before taking its square root to yield the standard deviation, filling the
state.sigma
and state.root
with their values respectively and the state.mu
array with zeros as before.To specify a mean as well we pass a third number to invoke the constructors given in listing 7.
Listing 7: The Identically Distributed Non-Zero Mean Uncorrelated Constructors
constructors[ak.NUMBER_T][ak.NUMBER_T][ak.NUMBER_T] = function(state, n, mu, sigma, args) { var arg3 = args[3]; var root; mu = Number(mu); if(!isFinite(mu)) throw new Error('invalid mu in ak.multiNormal distribution'); sigma = Number(sigma); if(sigma<=0 || !isFinite(sigma)) { throw new Error('invalid sigma in ak.multiNormal distribution'); } root = Math.sqrt(sigma); while(n--) { state.mu[n] = mu; state.sigma[n] = sigma; state.root[n] = root; } constructors[ak.NUMBER_T][ak.NUMBER_T][ak.NUMBER_T][ak.nativeType(arg3)](state,arg3); }; constructors[ak.NUMBER_T][ak.NUMBER_T][ak.NUMBER_T][ak.UNDEFINED_T] = function() { }; constructors[ak.NUMBER_T][ak.NUMBER_T][ak.NUMBER_T][ak.FUNCTION_T] = function(state, rnd) { state.rnd = rnd; };
This time the first number specifies the number of dimensions, the second the mean and the third the variance. Apart from this they're more or less the same as the previous constructors except that they also ensure that the mean is finite and fill
state.mu
with it rather than zero.As noted above it's handy to be able to define an independent multivariate normal random variable with different means and variances for each element by passing them as vectors. The constructors given in listing 8 allow this using
ak.vector
objects, firstly for the means and secondly for the variances.Listing 8: The Non-Identically Distributed Non-Zero Mean Uncorrelated Constructors
constructors[ak.VECTOR_T] = function(state, arg0, args) { var arg1 = args[1]; constructors[ak.VECTOR_T][ak.type(arg1)](state, arg0, arg1, args); }; constructors[ak.VECTOR_T][ak.VECTOR_T] = function(state, mu, sigma, args) { var n = sigma.dims(); var arg2 = args[2]; if(mu.dims()!==n) { throw new Error('dimension mismatch in ak.multiNormal distribution'); } state.mu = mu.toArray(); state.sigma = sigma.toArray(); state.root.length = n; while(n--) { if(!isFinite(state.mu[n])) { throw new Error('invalid mu in ak.multiNormal distribution'); } if(state.sigma[n]<=0 || !isFinite(state.sigma[n])) { throw new Error('invalid sigma in ak.multiNormal distribution'); } state.root[n] = Math.sqrt(state.sigma[n]); } constructors[ak.VECTOR_T][ak.VECTOR_T][ak.nativeType(arg2)](state, arg2); }; constructors[ak.VECTOR_T][ak.VECTOR_T][ak.UNDEFINED_T] = function() { }; constructors[ak.VECTOR_T][ak.VECTOR_T][ak.FUNCTION_T] = function(state, rnd) { state.rnd = rnd; };
Now we no longer need to explicitly state the number of dimensions since it is implicitly defined by the vector of variances. The main difference between these and the previous constructors is that we must check the validity of each mean and variance, and take the square root of the latter separately.
Note that once again we also allow the means to default to zero if we only pass the variances to the constructor, although we shan't show the implementation here.
With all of the uncorrelated constructors out of the way, we're ready to look at the construction of correlated multivariate normal random variables, although we first need a helper function to check whether a covariance matrix is diagonal and so implies no correlations, as given in listing 9.
Listing 9: Checking For Diagonal Covariance Matrices
function isDiagonal(m) { var n = m.rows(); var i, j; for(i=0;i<n;++i) { for(j=0;j<i;++j) if(m.at(i, j)!==0) return false; for(j=i+1;j<n;++j) if(m.at(i, j)!==0) return false; } return true; }
Listing 10 gives the constructor that takes an
ak.vector
mean and an ak.matrix
covariance and uses isDiagonal
to switch between uncorrelated and correlated random variables.Listing 10: The Non-Zero Mean Correlated Constructors
constructors[ak.VECTOR_T][ak.MATRIX_T] = function(state, mu, sigma, args) { var n = sigma.cols(); var arg2 = args[2]; var i, chol; if(sigma.rows()!==n) { throw new Error('non-square sigma in ak.multiNormal distribution'); } if(mu.dims()!==n) { throw new Error('dimension mismatch in ak.multiNormal distribution'); } state.mu = mu.toArray(); for(i=0;i<n;++i) if(!isFinite(state.mu[i])) { throw new Error('invalid mu in ak.multiNormal distribution'); } if(isDiagonal(sigma)) { state.sigma.length = n; state.root.length = n; for(i=0;i<n;++i) { state.sigma[i] = sigma.at(i, i); if(state.sigma[i]<=0 || !isFinite(state.sigma[i])) { throw new Error('invalid sigma in ak.multiNormal distribution'); } state.root[i] = Math.sqrt(state.sigma[i]); } } else { try {chol = ak.choleskyDecomposition(sigma);} catch(e) {throw new Error('invalid sigma in ak.multiNormal distribution');} state.root = chol.l().toArray(); state.sigma = chol.toMatrix().toArray(); while(n--) for(i=0;i<=n;++i) if(!isFinite(state.sigma[n][i])) { throw new Error('invalid sigma in ak.multiNormal distribution'); } } constructors[ak.VECTOR_T][ak.MATRIX_T][ak.nativeType(arg2)](state, arg2, args); }; constructors[ak.VECTOR_T][ak.MATRIX_T][ak.UNDEFINED_T] = function() { }; constructors[ak.VECTOR_T][ak.MATRIX_T][ak.FUNCTION_T] = function(state, rnd) { state.rnd = rnd; };
So first we check that
sigma
is a square matrix, that mu
has the same number of dimensions as it has rows and columns and that the elements of mu
are finite.Then, if
sigma
is diagonal we fill state.sigma
with its diagonal elements and state.root
with their square roots, checking that they are positive and finite, much as we did in the vector variance constructor.Otherwise we use our
ak.choleskyDecomposition
to find the lower triangular root \(\mathbf{L}\) from the covariance matrix \(\mathbf{\Sigma}\), storing the former in state.root
and the latter in state.sigma
and checking that all of the covariances are finite, relying upon the ak.choleskyDecomposition
constructor to throw an exception if sigma
isn't positive definite, being the correlated equivalent of having any non-positive variances.Note that since the
ak.choleskyDecomposition
constructor will initialise \(\mathbf{L}\) directly if given a lower triangular matrix, so will the correlated multivariate normal constructor.We also have a few last constructors, not listed here, for the covariance matrix without the mean vector and for the Cholesky decomposition of the covariance matrix, either as an
ak.choleskyDecomposition
object or as an object from which it can be constructed, both with and without the mean vector.So finally, having spent so much effort on the initialisation of multivariate normal random variables, we are ready to reap our glorious reward of an ever so slightly more convenient way of drawing from them, as demonstrated by program 1.
Program 1: Using ak.multiNormalRnd
|
|
---|---|
![]() |
ak.multiNormalMap
Whilst the implementation of the multivariate normal CDF is so tricky that we're going to give it an entire post of its own, the inverse CDF is even worse! Technically speaking it's not even a function since for any given probability there are an infinite number of vectors that a multivariate normal random vector's elements will all be smaller than with that probability.As we did for the multivariate uniform distribution[2] we shall instead define a function that uniquely maps a vector whose elements all lie between zero and one onto another vector in such a way that if the former's elements are all independent standard uniformly distributed then the latter's will be multivariate normally distributed.
Specifically, we'll just use the elements of the vector as the source of randomness in our definition of multivariate random variables by using the univariate normal inverse CDF to transform them into observations of a standard univariate normal random variable.
Listing 11 show how we do this for an uncorrelated multivariate normal distribution with
uncorrMap
by replacing the ak.normalRnd
argument rnd
of the uncorrRnd
function from listing 2 with a ak.normalInvCDF
argument uinv
and passing the vector as an extra argument c
.Listing 11: The Uncorrelated Multivariate Normal Map
function uncorrMap(mu, root, uinv, c) { var n = mu.length; var i; if(ak.type(c)!==ak.VECTOR_T) { throw new Error('invalid argument in ak.multiNormalMap'); } if(c.dims()!==n) { throw new Error('dimension mismatch in ak.multiNormalMap'); } c = c.toArray(); for(i=0;i<n;++i) c[i] = mu[i] + root[i]*uinv(c[i]); return ak.vector(c); }
After checking that the latter is indeed a vector and has the correct number of dimensions, we have simply substituted each call to
rnd
with an inversion of an element of c
.The implementation of the correlated multivariate normal map requires precisely the same changes to the
corrRnd
function from listing 1, as shown by corrMap
in listing 12.Listing 12: The Correlated Multivariate Normal Map
function corrMap(mu, root, uinv, c) { var n = mu.length; var i, j, rooti; if(ak.type(c)!==ak.VECTOR_T) { throw new Error('invalid argument in ak.multiNormalMap'); } if(c.dims()!==n) { throw new Error('dimension mismatch in ak.multiNormalMap'); } c = c.toArray(); if(c.some(isNaN)) return ak.vector(n, ak.NaN); for(i=0;i<n;++i) c[i] = uinv(c[i]); for(i=n-1;i>=0;--i) { rooti = root[i]; c[i] = mu[i] + rooti[i]*c[i]; for(j=0;j<i;++j) c[i] += rooti[j]*c[j]; if(isNaN(c[i])) c[i] = -ak.INFINITY; } return ak.vector(c); }
Finally, we can implement an
ak.multiNormalMap
function that handles both cases by replacing the rnd
random variable in ak.multiNormalRnd
with a uinv
univariate inverse CDF and the calls to the univariate and multivariate random number generators with calls to the univariate and multivariate maps, as illustrated by listing 13.Listing 13: ak.multiNormalMap
ak.multiNormalMap = function() { var uinv = ak.normalInvCDF(); var state = {mu: [0], sigma: [1], root: [1]}; var arg0 = arguments[0]; var f; constructors[ak.type(arg0)](state, arg0, arguments); if(state.sigma.length===0 || ak.nativeType(state.sigma[0])===ak.NUMBER_T) { f = function(c){return uncorrMap(state.mu, state.root, uinv, c);}; f.sigma = function(){return ak.matrix('diagonal', state.sigma);}; } else { f = function(c){return corrMap(state.mu, state.root, uinv, c);}; f.sigma = function(){return ak.matrix(state.sigma);}; } f.mu = function(){return ak.vector(state.mu);}; return Object.freeze(f); };
One way that we can use
ak.multiNormalMap
is to more evenly sample vectors from a multivariate normal distribution than we can with ak.multiNormalRnd
by passing it vectors generated with a quasi random number generator, designed to more evenly fill multidimensional unit cubes than pseudo random number generators do, such as our ak.haltonRnd
[3], as demonstrated by program 2.
Program 2: Using ak.multiNormalMap
|
|
---|---|
![]() |
ak.multiNormalPDF
Next we come to the PDF, which we shall also split into uncorrelated and correlated implementations with the former simply being the product of the PDFs of the independent elements
\[
\mathbf{\phi}_{\mathbf{\mu},\mathbf{\sigma}}
= \prod_i \phi_{\mu_i,\sigma_i}
= \prod_i \frac{1}{\sqrt{2\pi}\sigma_i}e^{-\frac12 \left(\frac{x-\mu_i}{\sigma_i}\right)^2}
= \frac{1}{\sqrt{2\pi}\prod_i\sigma_i}e^{-\frac12 \sum_i \left(\frac{x-\mu_i}{\sigma_i}\right)^2}
\]
where \(\prod\) is the product sign and \(\sum\) is the summation sign, as implemented by the uncorrPDF
function given in listing 14.Listing 14: The Uncorrelated Multivariate Normal PDF
function uncorrPDF(mu, root, scale, x) { var n = mu.length; var z2 = 0; var i, zi; if(ak.type(x)!==ak.VECTOR_T) { throw new Error('invalid argument in ak.multiNormalPDF'); } if(x.dims()!==n) { throw new Error('dimension mismatch in ak.multiNormalPDF'); } x = x.toArray(); for(i=0;i<n;++i) x[i] -= mu[i]; for(i=0;i<n;++i) { zi = x[i]*root[i]; z2 += zi*zi; } return scale*Math.exp(-0.5*z2); }
Note that the
scale
argument is a pre-calculated value for
\[
\frac{1}{\sqrt{2\pi}\prod_i\sigma_i}
\]
and the elements of the root
argument are the reciprocals of the standard deviations.For the correlated case we have
\[
\phi_{\mathbf{\mu},\mathbf{\Sigma}}(\mathbf{z})
= \frac{1}{(2\pi)^{\frac12 n} |\mathbf{\Sigma}|^\frac12} e^{-\frac12 \left(\mathbf{z}-\mathbf{\mu}\right) \times \mathbf{\Sigma}^{-1} \times \left(\mathbf{z}-\mathbf{\mu}\right)}
\]
where \(|\mathbf{\Sigma}|\) is the determinant of the covariance matrix \(\mathbf{\Sigma}\). The corrPDF
function given in listing 15 gives its implementation.Listing 15: The Correlated Multivariate Normal PDF
function corrPDF(mu, root, scale, x) { var n = mu.length; var z2 = 0; var i, j, zi, rooti; if(ak.type(x)!==ak.VECTOR_T) { throw new Error('invalid argument in ak.multiNormalPDF'); } if(x.dims()!==n) { throw new Error('dimension mismatch in ak.multiNormalPDF'); } x = x.toArray(); for(i=0;i<n;++i) x[i] -= mu[i]; if(x.some(isNaN)) return ak.NaN; for(i=0;i<n;++i) { rooti = root[i]; zi = 0; for(j=0;j<i;++j) zi += rooti[j]*x[j]; zi = (x[i]-zi)*rooti[i]; x[i] = zi; z2 += zi*zi; } if(!isFinite(z2)) return 0; return scale*Math.exp(-0.5*z2); }
This time
scale
is a pre-calculated value for
\[
\frac{1}{(2\pi)^{\frac12 n} |\mathbf{\Sigma}|^\frac12}
\]
and root
contains the elements of the lower triangular \(\mathbf{L}\), albeit with the leading diagonal containing their reciprocals instead. Note that rather than explicitly using the inverse of \(\mathbf{\Sigma}\) we're solving for
\[
\mathbf{L} \times \mathbf{z} = \mathbf{x} - \mathbf{\mu}
\]
by noting that
\[
\left(\mathbf{L} \times \mathbf{z}\right)_j = \sum_{j \leqslant i} L_{ij} \times z_j
\]
since \(\mathbf{L}\) is lower triangular, and so
\[
\begin{align*}
L_{ii} \times z_i &= x_j - \mu_j - \sum_{i < j} L_{ij} \times z_j\\
z_i &= \frac{1}{L_{ii}} \times \left(x_j - \mu_j - \sum_{i < j} L_{ij} \times z_j\right)
\end{align*}
\]
In order to do this we're keeping track of the elements of \(\mathbf{z}\) in x
as we calculate them.Finally, the
ak.multiNormalPDF
function given in listing 16 switches between the uncorrelated and correlated helper functions depending on whether state.sigma
is an array of variances or an array of arrays of covariances, computing the appropriate values for the scale
and reciprocating the relevant elements of state.root
in both cases.Listing 16: ak.multiNormalPDF
ak.multiNormalPDF = function() { var state = {mu: [0], sigma: [1], root: [1]}; var arg0 = arguments[0]; var n, scale, i, f; constructors[ak.type(arg0)](state, arg0, arguments); n = state.mu.length; scale = Math.pow(RT2PI, -n); if(state.sigma.length===0 || ak.nativeType(state.sigma[0])===ak.NUMBER_T) { for(i=0;i<n;++i) { scale /= state.root[i]; state.root[i] = 1/state.root[i]; } f = function(x){return uncorrPDF(state.mu, state.root, scale, x);}; f.sigma = function(){return ak.matrix('diagonal', state.sigma);}; } else { for(i=0;i<n;++i) { scale /= state.root[i][i]; state.root[i][i] = 1/state.root[i][i]; } f = function(x){return corrPDF(state.mu, state.root, scale, x);}; f.sigma = function(){return ak.matrix(state.sigma);}; } f.mu = function(){return ak.vector(state.mu);}; return Object.freeze(f); };
Program 3 demonstrates the use of
ak.multiNormalPDF
with a two dimensional contour plot of its values.
Program 3: Using ak.multiNormalPDF
|
|
---|---|
![]() |
ak.multiNormalCF
Next up is the characteristic function, given by
\[
\hat{\phi}_{\mathbf{\mu},\mathbf{\Sigma}}(\mathbf{t})
= e^{i\mathbf{\mu}\mathbf{t} - \frac12 \mathbf{t} \mathbf{\Sigma} \mathbf{t}}
= e^{i\mathbf{\mu}\mathbf{t} - \frac12 \left(\mathbf{t} \mathbf{L}\right) \times \left(\mathbf{L}^\mathrm{T} \mathbf{t}\right)}
\]
where \(i\) is the square root of minus one.
In the uncorrelated case this simplifies to
\[
\hat{\phi}_{\mathbf{\mu},\mathbf{\sigma}}(\mathbf{t})
= e^{\sum_j i \mu_j t_j - \frac12 \left(\sigma_j t_j\right)^2}
\]
as implemented by the uncorrCF
function given in listing 17.Listing 17: The Uncorrelated Multivariate Normal CF
function uncorrCF(mu, root, t) { var n = mu.length; var re = 0; var im = 0; var i, rt; if(ak.type(t)!==ak.VECTOR_T) { throw new Error('invalid argument in ak.multiNormalCF'); } if(t.dims()!==n) { throw new Error('dimension mismatch in ak.multiNormalCF'); } t = t.toArray(); if(t.some(isNaN)) return ak.complex(ak.NaN, ak.NaN); for(i=0;i<n;++i) { rt = root[i]*t[i]; re += rt*rt; im += mu[i]*t[i]; } re = Math.exp(-0.5*re); return re>0 ? ak.complex(re*Math.cos(im), re*Math.sin(im)) : ak.complex(0); }
Note that we're exploiting the identity
\[
e^{x + iy} = e^x \times \left(\cos y + i \sin y\right)
\]
to construct the ak.complex
result. Note that we must be slightly careful to return zero if any elements of t
are infinite, which we do by checking whether re
is greater than zero and returning NaN near the start of the function if any of them are NaN and would therefore also fail that test.The
corrCF
function given in listing 18 implements the correlated case by replacing the simple multiplication for rt
with a matrix-vector multiplication.Listing 18: The Correlated Multivariate Normal CF
function corrCF(mu, root, t) { var n = mu.length; var re = 0; var im = 0; var i, j, rt; if(ak.type(t)!==ak.VECTOR_T) { throw new Error('invalid argument in ak.multiNormalCF'); } if(t.dims()!==n) { throw new Error('dimension mismatch in ak.multiNormalCF'); } t = t.toArray(); if(t.some(isNaN)) return ak.complex(ak.NaN, ak.NaN); for(i=0;i<n;++i) { rt = 0; for(j=i;j<n;++j) rt += root[j][i]*t[j]; re += rt*rt; im += mu[i]*t[i]; } re = Math.exp(-0.5*re); return re>0 ? ak.complex(re*Math.cos(im), re*Math.sin(im)) : ak.complex(0); }
Finally, the
ak.multiNormalCF
selects which of these helper functions to use to compute the CF.Listing 19: ak.multiNormalCF
ak.multiNormalCF = function() { var state = {mu: [0], sigma: [1], root: [1]}; var arg0 = arguments[0]; var f; constructors[ak.type(arg0)](state, arg0, arguments); if(state.sigma.length===0 || ak.nativeType(state.sigma[0])===ak.NUMBER_T) { f = function(t){return uncorrCF(state.mu, state.root, t);}; f.sigma = function(){return ak.matrix('diagonal', state.sigma);}; } else { f = function(t){return corrCF(state.mu, state.root, t);}; f.sigma = function(){return ak.matrix(state.sigma);}; } f.mu = function(){return ak.vector(state.mu);}; return Object.freeze(f); };
Program 4 demonstrates the use of
ak.multiNormalCF
by plotting the real parts of its results in the upper contour chart and their imaginary parts in the lower one with red points representing negative values and green points positive ones in both cases.
Program 4: Using ak.multiNormalCF
|
|
---|---|
![]() |
|
![]() |
Note that all of these functions can be found in
MultiNormalDistribution.js
, although I'd suggest that you ignore the implementation of the CDF until after we've covered it next time.

References
[1] Every Which Way Is Normal, www.thusspakeak.com, 2017.[2] Into The nth Dimension, www.thusspakeak.com, 2016.
[3] Halton, Here's A Pseud, www.thusspakeak.com, 2016.
Leave a comment