Многомерное нормальное распределение

Математический алгоритм:

1. Перейти в новуй систему координат (x->y), сдвинув начало координат, в точку математического ожидания. Пепесчитать матрицу ковариаций. Ay[i,j]=Ax[i,j]-m[i]*m[j]
2. Просверим что матрица Ay не диагональная.
3. Найдем собственные вектора матрицы Ay E=EigenVectors(Ay).
4. Перейдем в новую систему координат (вращение), так чтобы оси новой систему соотвесвовали собственным векторам (y->z). Az = E * Ay * E^-1. В этой системе координат, матрица Az должна быть диагональной.
5. Нахидом случайное число z.
6. Переходим в исходную систему координат z->y->x

В случае, если нужно сгенерировать много случайных чисел, то шаги 1-4 можно выполнить один раз.

Критическим здесь является являестя этап нахиждения собственных векторов. К сожелению, в boost.uBLAS, нет подходящей функции, но можно воспользоваться eigen или MKL. Кроме того, задача нахождения собственных векторов - очень трудоемкая (сложенее нахождения обратной матрицы), и в случае плохо обусловленных матриц большого размера, удовлетворительного алгоритма просто не сущесвует. Даже размерность 10 может оказаться критической, с точки зрения потери точности.