This is a attempt to start discussing the interface of multivariate distributions. In this example I've used the bivariate normal distribution, but it the idea's should extend to multivariate distributions. An instance of the bivariate normal distribution is uniquely defined by its mean (2 elements) and its covariance matrix. The covariance matrix is symmetric by definition. The covariance matrix should also be "semi-positive definite", a check for that might end up in the "detail/common_error_handling.hpp" Below are three version of construction a (bi)(multi)variate normal distribution. The first version uses vectors and matrices as argument to present the info. The other versions are variadic constructors, presents all the necessary *elements* of the vector and matrix, circumventing the need to "use" and "fill" specific containers for vectors and matrices. The dimension should be deducible from the number of arguments. // Constructor #1: // pass the mean vector, as an object of vector type // pass the covariance matrix as an object of matrix type // // traits: real the scalar type used for precision, return type of functions... // vector should have an operator(i) // matrix should have an operator(i,j) // VectorType mean; MatrixType cov; ... bivariate_normal_distribution<RealType,VectorType,MatrixType> BN1(mean,cov); // Constructor #2: // pass the mean vector, element-by-element: x,y // pass the upper covariance matrix: COV(x,x), COV(x,y), COV(y,y) // // the generic N-variate case will have N means and N*(N+1)/2 COV elements // N | Number of argument // 2 | 2+3 = 5 // 3 | 3+6 = 9 // 4 | 4+10 = 14 // ... // N | N*(N+3)/2 bivariate_normal_distribution<> BN1(0.1, 0.1, 1.0, 0.9, 1.0); // Constructor #3: // pass the mean vector, element-by-element: x,y // pass the full covariance matrix: COV(x,x), COV(x,y),COV(y,x), COV(y,y) // // the generic N-variate case will have N means and N*N COV elements // // the generic N-variate case will have N means and N*(N+1)/2 COV elements // N | Number of argument // 2 | 2+4 = 6 // 3 | 3+9 = 12 // 4 | 4+16 = 20 // ... // N | N*(N+1) bivariate_normal_distribution<> BN1(0.1, 0.1, 1.0, 0.9, 0.9, 1.0);