Multiple Multiply Normal Functions

The Cumulative Distribution Unction

Every Which Way Is Normal

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}\)

\[ \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 the corrRnd 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
var m = ak.matrix([[0.7,0.8],[0.3,1.5]]); var s = ak.mul(m, ak.transpose(m)); var rnd = ak.multiNormalRnd(s); var n = 2000; scatter.bounds([[-4, -4], [4, 4]]); ak.asynchLoop( function body() { scatter.plot([rnd().toArray()]); return n-->0; } );

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
var m = ak.matrix([[0.7,0.8],[0.3,1.5]]); var s = ak.mul(m, ak.transpose(m)); var rnd = ak.haltonRnd([2, 3]); var map = ak.multiNormalMap(s); var n = 2000; scatter.bounds([[-4, -4], [4, 4]]); ak.asynchLoop( function body() { var u = ak.vector(rnd()); scatter.plot([map(u).toArray()]); return n-->0; } );

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
var m = ak.matrix([[1.2,0.8],[0.3,1.5]]); var s = ak.mul(m, ak.transpose(m)); var f = ak.multiNormalPDF(s); var n = 2; var pl = ak.ui.level(0.0,'Black'); var pu = ak.ui.level(1.0,'white'); var p = ak.ui.palette(256, [pl,pu]); var bl = [-3,-3, 0]; var bu = [ 3, 3, f(ak.vector(2, 0))]; var ticks = [-3,-2,-1,1,2,3]; var rows = contour.rows(); var cols = contour.cols(); var row = 0; contour.bounds([bl,bu]); contour.pixelSize(n); contour.palette(p); ak.asynchLoop( function body() { var col, xy; for(col=0;col<cols;col+=n) { xy = contour.toXY([col, row]); xy[2] = f(ak.vector(xy)); contour.plot(xy); } row += n; return row<rows; }, function end() { contour.axes([0,0],[ticks,ticks], {strokeStyle:'Gray'}); } );

Program 3: Using ak.multiNormalPDF

var m = ak.matrix([[1.2,0.8],[0.3,1.5]]);
var s = ak.mul(m, ak.transpose(m));
var f = ak.multiNormalPDF(s);

var n = 2;
var pl = ak.ui.level(0.0,'Black');
var pu = ak.ui.level(1.0,'white');
var p = ak.ui.palette(256, [pl,pu]);
var bl = [-3,-3, 0];
var bu = [ 3, 3, f(ak.vector(2, 0))];
var ticks = [-3,-2,-1,1,2,3];
var rows = contour.rows();
var cols = contour.cols();
var row = 0;

contour.bounds([bl,bu]);
contour.pixelSize(n);
contour.palette(p);

ak.asynchLoop(
 function body() {
  var col, xy;

  for(col=0;col<cols;col+=n) {
   xy = contour.toXY([col, row]);
   xy[2] = f(ak.vector(xy));
   contour.plot(xy);
  }
  row += n;
  return row<rows;
 },
 function end() {
  contour.axes([0,0],[ticks,ticks], {strokeStyle:'Gray'});
 }
);

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
var l = ak.matrix([[1,0],[0.5,1]]); var s = ak.mul(l, ak.transpose(l)); var u = ak.vector([-3, 5]); var f = ak.multiNormalCF(u, s); var n = 2; var pl = ak.ui.level(0.0,'Red'); var pm = ak.ui.level(0.5,'Black'); var pu = ak.ui.level(1.0,'Green'); var p = ak.ui.palette(256, [pl,pm,pu]); var bl = [-2,-2,-1]; var bu = [ 2, 2, 1]; var ticks = [-2, -1, 1, 2]; var rows = real.rows(); var cols = real.cols(); var row = 0; real.bounds([bl,bu]); real.pixelSize(n); real.palette(p); imag.bounds([bl,bu]); imag.pixelSize(n); imag.palette(p); ak.asynchLoop( function body() { var col, xy, z; for(col=0;col<cols;col+=n) { xy = real.toXY([col, row]); z = f(ak.vector(xy)); xy[2] = z.re(); real.plot(xy); xy[2] = z.im(); imag.plot(xy); } row += n; return row<rows; }, function end() { real.axes([0,0],[ticks,ticks], {strokeStyle:'Gray'}); imag.axes([0,0],[ticks,ticks], {strokeStyle:'Gray'}); } );

Program 4: Using ak.multiNormalCF

var l = ak.matrix([[1,0],[0.5,1]]);
var s = ak.mul(l, ak.transpose(l));
var u = ak.vector([-3, 5]);
var f = ak.multiNormalCF(u, s);

var n = 2;
var pl = ak.ui.level(0.0,'Red');
var pm = ak.ui.level(0.5,'Black');
var pu = ak.ui.level(1.0,'Green');
var p = ak.ui.palette(256, [pl,pm,pu]);
var bl = [-2,-2,-1];
var bu = [ 2, 2, 1];
var ticks = [-2, -1, 1, 2];
var rows = real.rows();
var cols = real.cols();
var row = 0;

real.bounds([bl,bu]);
real.pixelSize(n);
real.palette(p);

imag.bounds([bl,bu]);
imag.pixelSize(n);
imag.palette(p);

ak.asynchLoop(
 function body() {
  var col, xy, z;

  for(col=0;col<cols;col+=n) {
   xy = real.toXY([col, row]);
   z = f(ak.vector(xy));
   xy[2] = z.re();
   real.plot(xy);
   xy[2] = z.im();
   imag.plot(xy);
  }
  row += n;
  return row<rows;
 },
 function end() {
  real.axes([0,0],[ticks,ticks], {strokeStyle:'Gray'});
  imag.axes([0,0],[ticks,ticks], {strokeStyle:'Gray'});
 }
);

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 n^th Dimension, www.thusspakeak.com, 2016.

[3] Halton, Here's A Pseud, www.thusspakeak.com, 2016.


AKCalc	ECMA


Endarkenment	Turning Sixteen

Search

Multiple Multiply Normal Functions

Next:

Previous:

ak.multiNormalRnd

Multivariate Normal Distribution Constructors

ak.multiNormalMap

ak.multiNormalPDF

ak.multiNormalCF

References

Next:

Previous:

Categories:

Tags:

Leave a comment

Categories

Gallimaufry

Links

	Chrome 26+
	Firefox 20+
	Internet Explorer 9+

Sign In

Search

Multiple Multiply Normal Functions

Next:

Previous:

ak.multiNormalRnd

Multivariate Normal Distribution Constructors

ak.multiNormalMap

ak.multiNormalPDF

ak.multiNormalCF

References

Next:

Previous:

Categories:

Tags:

Leave a comment

Tag Cloud

Categories

Gallimaufry

Links