🧄Python分析绘图蒙特卡洛和马尔可夫链蒙特卡洛样本
Python | 蒙特卡洛 | 马尔可夫链蒙特卡洛 | 条件概率 | 概率分布 | 核密度估计(KDE) | 边缘化 | 密度 | 样本 | 高斯分布 | 相关长度 | NumPy | scipy
本文概述了使用核密度估计 (KDE) 计算边缘化一维和二维密度的方法，分析和绘图蒙特卡洛和马尔可夫链蒙特卡洛样本。
蒙特卡洛
蒙特卡罗是一种计算技术，它基于为问题构建随机过程并通过从具有预定概率分布的随机数字序列中进行 N 倍抽样来进行数值实验。
$X$ - 随机变量
$$ \hat{\mathrm{X}}=\frac{1}{\mathrm{~N}} \sum_{\mathrm{i}=1}^{\mathrm{N}} \mathrm{x}_{\mathrm{i}} $$
$\hat{\mathrm{X}}$ - $X$ 的估计或样本均值
$\bar{X}$ - $X$ 的期望值或真实平均值
如果可以给问题一个概率解释，那么可以使用随机数对其进行建模。
算法组成部分
概率分布函数 (pdf) - 物理（或数学）系统必须由一组 pdf 描述。
随机数生成器 - 必须有一个均匀分布在单位间隔上的随机数源。
抽样规则 - 从指定的 pdf 抽样的处方，假设在单位间隔上随机数的可用性。
计分（或计分）- 结果必须累积成总体计分或感兴趣数量的分数。
误差估计 - 必须确定作为试验次数和其他数量函数的统计误差（方差）的估计。
方差减少技术 - 减少估计解中的方差以减少蒙特卡罗模拟计算时间的方法。
并行化和矢量化 - 高效使用先进的计算机架构。
马尔可夫链蒙特卡洛
马尔可夫链是一个数学过程，它经历了从一种状态到另一种状态的转换。 马尔可夫过程的关键特性是它是随机的，并且过程中的每一步都是“无记忆的”； 换句话说，未来状态仅取决于过程的当前状态，而不取决于过去。
这些步骤的连续是马尔可夫链。
这个随机过程可以用条件概率来描述：
$$ \operatorname{Pr}\left(X_{n+1}=x \mid X_{1}=x_{1}, X_{2}=x_{2}, \ldots, X_{n}=x_{n}\right)=\operatorname{Pr}\left(X_{n+1}=x \mid X_{n}=x_{n}\right) $$
$X_i$ 的可能值形成可数集 S，即链的状态空间。 马尔可夫链经常用有向图表示（与我们通常的有向无环图相反），其中边用从一种状态 (S) 到另一种状态的概率标记。
马尔可夫链蒙特卡罗 (MCMC) 模拟允许参数估计，例如均值、方差、预期值和贝叶斯模型的后验分布的探索。 为了评估“后验”的属性，从该分布中抽取许多具有代表性的随机值。
Metropolis 算法是一种常用的 MCMC 过程。 该算法产生所谓的“随机游走”，其中以小步长重复采样分布； 独立于之前的移动，因此是无记忆的。 然后使用此过程，例如，描述分布或计算期望值。
核密度估计 (KDE)
使用最大似然交叉验证的内核构建和带宽优化。
核函数的第一个属性是它必须是对称的。 这意味着 u 和 –u 的核函数值相同，如下图所示。 这可以在数学上表示为 K (-u) = K ( u)。 核函数的对称特性使函数的最大值 (max(K(u)) 位于曲线的中间。
函数曲线下的面积必须等于 1。
在数学上，这个性质表示为
$$ \int_{-\infty}^{+\infty} K(u) d u=1 $$
使用高斯密度函数作为核函数是因为高斯密度曲线下的面积是1，而且也是对称的。
核函数的值，即密度，不能为负，K(u) ≥ 0 对于所有 -∞ < u < ∞。
Python分析和绘制
# Get some random samples for demonstration:
# make random covariance, then independent samples from Gaussian
ndim = 4
nsamp = 10000
random_state = np.random.default_rng(10) # seed random generator
A = random_state.random((ndim,ndim))
cov = np.dot(A, A.T)
samps = random_state.multivariate_normal([0]*ndim, cov, size=nsamp)
A = random_state.random((ndim,ndim))
cov = np.dot(A, A.T)
samps2 = random_state.multivariate_normal([0]*ndim, cov, size=nsamp)# Get objects for the samples, specifying same parameter
# names and labels; if not specified weights are assumed to all be unity
names = ["x%s"%i for i in range(ndim)]
labels = ["x_%s"%i for i in range(ndim)]
samples = Samples(samples=samps,names = names, labels = labels)
samples2 = Samples(samples=samps2,names = names, labels = labels, label='Second set')# Triangle plot
g = plots.get_subplot_plotter()
g.triangle_plot([samples, samples2], filled=True)Last updated