title: Implementing a multivariate gaussian probability density function for >2 dimensions in C++
body: <p>I'm working on implementing a probability density function of a multivariate Gaussian in C++, and I'm stuck on how to best handle the cases where dimension > 2.</p>
  
  <p>The pdf of a gaussian can be written as</p>
  
  <p><img src="https://i.stack.imgur.com/Z97fh.png" alt="multivariate gaussian pdf"></p>
  
  <p>where (A)' or A' represents the transpose of the 'matrix' created by subtracting the mean from all elements of x. In this equation, k is the number of dimensions we have, and sigma represents the covariance matrix, which is a k x k matrix. Finally, |X| means the determinant of the matrix X.</p>
  
  <p>In the univariate case, implementing the pdf is trivial. Even in the bivariate (k = 2) case, it's trivial. However, when we go further than two dimensions, the implementation is much harder.</p>
  
  <p>In the bivariate case, we'd have</p>
  
  <p><img src="https://i.stack.imgur.com/oKBPv.png" alt="bivariate gaussian pdf"></p>
  
  <p>where rho is the correlation between x and y, with correlation equal to</p>
  
  <p><img src="https://i.stack.imgur.com/fj4Px.png" alt="correlation between two random variables X and Y"></p>
  
  <p>In this case, I could use <code>Eigen::Matrix&lt;double, Eigen::Dynamic, Eigen::Dynamic&gt;</code> to implement the first equation, or just calculate everything myself using the second equation, without benefiting from Eigen's simplified linear algebra interface.</p>
  
  <p>My thoughts for an attempt at the multivariate case would probably begin by extending the above equations to the multivariate case</p>
  
  <p><img src="https://i.stack.imgur.com/rODv5.png" alt="multivariate pdf"></p>
  
  <p>with</p>
  
  <p><img src="https://i.stack.imgur.com/Gvmu7.png" alt="multivariate pdf"></p>
  
  <p>My questions are:</p>
  
  <ol>
  <li>Would it be appropriate/advised to use a <code>boost::multi_array</code> for the
  n-dimensional array, or should I instead try to leverage Eigen? </li>
  <li>Should I have separate functions for the univariate/bivariate cases,
  or should I just abstract it all to the multivariate case using
  boost::multi_array (or an appropriate alternative)?</li>
  </ol>