在 C++ 中从零实现 Wishart 分布随机矩阵生成的方法。基于 Bartlett 分解，通过 Cholesky 分解尺度矩阵 Sigma 和构造包含标准正态及卡方分布元素的矩阵 A，最终计算 W = B * B^T 得到样本。代码不依赖 BLAS/LAPACK，数值稳定且注释完整，适用于多元统计、贝叶斯建模等场景。
在多元统计分析、随机矩阵理论、贝叶斯统计、机器学习以及信号处理等领域中，Wishart 分布（Wishart Distribution） 是一个极其重要的矩阵值概率分布。
它可以被看作是：
多元正态分布样本协方差矩阵的分布
换句话说：
Wishart 分布广泛应用于：
📌 在工程和科研中，能够正确、高效地从 Wishart 分布采样是很多高级算法的基础模块。
尽管很多统计库提供现成函数，但手写 Wishart 采样具有重要价值：
本项目目标是：
👉 使用 C++ 从零实现 Wishart 分布随机矩阵生成
核心功能包括：
nSigma (p × p)W (p × p)Wishart 分布的采样可转化为两个关键步骤：
卡方分布可以通过对标准正态分布变量的平方求和来生成。
1. 对 Sigma 进行 Cholesky 分解 2. 构造 Bartlett 矩阵 A 3. 计算 B = L * A 4. 返回 W = B * B^T
std::vector<std::vector<double>>std::mt19937/******************************************************
* File: wishart.h
******************************************************/
#ifndef WISHART_H
#define WISHART_H
#include <vector>
/* 从 Wishart 分布生成样本矩阵 */
class WishartSampler {
public:
WishartSampler(int dim, int df);
std::vector<std::vector<double>> sample(const std::vector<std::vector<double>>& sigma);
private:
int p; // 维度
int n; // 自由度
std::vector<std::vector<double>> choleskyDecomposition(
const std::vector<std::vector<double>>& A);
std::vector<std::vector<double>> multiply(
const std::vector<std::vector<double>>& A,
const std::vector<std::vector<double>>& B);
std::vector<std::vector<double>> transpose(const std::vector<std::vector<double>>& A);
};
#endif
/******************************************************
* File: wishart.cpp
******************************************************/
#include "wishart.h"
#include <cmath>
#include <random>
#include <stdexcept>
/* 构造函数 */
WishartSampler::WishartSampler(int dim, int df) : p(dim), n(df) {
if (n < p) throw std::invalid_argument("df must be >= dimension");
}
/* Cholesky 分解（下三角） */
std::vector<std::vector<double>> WishartSampler::choleskyDecomposition(
const std::vector<std::vector<double>>& A) {
std::vector<std::vector<double>> L(p, std::vector<double>(p, 0.0));
for (int i = 0; i < p; ++i) {
for (int j = 0; j <= i; ++j) {
double sum = 0.0;
for (int k = 0; k < j; ++k)
sum += L[i][k] * L[j][k];
if (i == j)
L[i][j] = std::sqrt(A[i][i] - sum);
else
L[i][j] = (A[i][j] - sum) / L[j][j];
}
}
return L;
}
/* 矩阵乘法 */
std::vector<std::vector<double>> WishartSampler::multiply(
const std::vector<std::vector<double>>& A,
const std::vector<std::vector<double>>& B) {
std::vector<std::vector<double>> C(p, std::vector<double>(p, 0.0));
for (int i = 0; i < p; ++i)
for (int j = 0; j < p; ++j)
for (int k = 0; k < p; ++k)
C[i][j] += A[i][k] * B[k][j];
return C;
}
/* 转置 */
std::vector<std::vector<double>> WishartSampler::transpose(
const std::vector<std::vector<double>>& A) {
std::vector<std::vector<double>> T(p, std::vector<double>(p));
for (int i = 0; i < p; ++i)
for (int j = 0; j < p; ++j)
T[j][i] = A[i][j];
return T;
}
/* Wishart 采样 */
std::vector<std::vector<double>> WishartSampler::sample(
const std::vector<std::vector<double>>& sigma) {
std::mt19937 gen(std::random_device{}());
std::normal_distribution<> normal(0.0, 1.0);
/* Step 1: Cholesky 分解 */
auto L = choleskyDecomposition(sigma);
/* Step 2: 构造 Bartlett 矩阵 A */
std::vector<std::vector<double>> A(p, std::vector<double>(p, 0.0));
for (int i = 0; i < p; ++i) {
double chi2 = 0.0;
for (int k = 0; k < n - i; ++k) {
double z = normal(gen);
chi2 += z * z;
}
A[i][i] = std::sqrt(chi2);
for (int j = 0; j < i; ++j)
A[i][j] = normal(gen);
}
/* Step 3: B = L * A */
auto B = multiply(L, A);
/* Step 4: W = B * B^T */
auto BT = transpose(B);
return multiply(B, BT);
}
/******************************************************
* File: main.cpp
******************************************************/
#include <iostream>
#include "wishart.h"
int main() {
int p = 3;
int df = 6;
std::vector<std::vector<double>> sigma = {
{1.0, 0.5, 0.2},
{0.5, 1.0, 0.3},
{0.2, 0.3, 1.0}
};
WishartSampler sampler(p, df);
auto W = sampler.sample(sigma);
std::cout << "Wishart sample matrix:\n";
for (const auto& row : W) {
for (double v : row)
std::cout << v << " ";
std::cout << "\n";
}
return 0;
}
samplecholeskyDecomposition通过本项目，我们：
该实现：
A：这是 Wishart 分布存在的必要条件。
A：是的，Bartlett 分解理论上保证正定。
A：对 Wishart 矩阵取逆即可。
微信公众号「极客日志」，在微信中扫描左侧二维码关注。展示文案：极客日志 zeeklog
使用加密算法（如AES、TripleDES、Rabbit或RC4）加密和解密文本明文。 在线工具，加密/解密文本在线工具，online
生成新的随机RSA私钥和公钥pem证书。 在线工具，RSA密钥对生成器在线工具，online
基于 Mermaid.js 实时预览流程图、时序图等图表，支持源码编辑与即时渲染。 在线工具，Mermaid 预览与可视化编辑在线工具，online
将字符串编码和解码为其 Base64 格式表示形式即可。 在线工具，Base64 字符串编码/解码在线工具，online
将字符串、文件或图像转换为其 Base64 表示形式。 在线工具，Base64 文件转换器在线工具，online
将 Markdown（GFM）转为 HTML 片段，浏览器内 marked 解析；与 HTML 转 Markdown 互为补充。 在线工具，Markdown 转 HTML在线工具，online