3

我想 (i) 计算和 (ii) 绘制 Distributions.jl 库中分布的中心可信区间和最高后验密度区间。理想情况下,可以编写自己的函数来计算 CI 和 HPD,然后使用 Plots.jl 绘制它们。但是,我发现实现非常棘手(免责声明:我是 Julia 的新手)。有什么关于库/gists/repo 的建议可以让计算和绘制它们更容易吗?

语境

using Plots, StatsPlots, LaTeXStrings
using Distributions

dist = Beta(10, 10)
plot(dist)  # thanks to StatsPlots it nicely plots the distribution

# missing piece 1: compute CI and HPD
# missing piece 2: plot CI and HPD

预期的最终结果总结在下图或第页。BDA3的第 33 条。

在此处输入图像描述

目前找到的资源:

4

3 回答 3

3

感谢您更新问题;它带来了一个新的视角。

要点是正确的;只有它使用早期版本的 Julia。因此linspace应替换为LinRange. 而不是using PyPlot使用using Plots. 我会将绘图部分更改为以下内容:

plot(cred_x, pdf(B, cred_x), fill=(0, 0.9, :orange))
plot!(x,pdf(B,x), title="pdf with 90% region highlighted")

乍一看,CI 的计算似乎是正确的。(就像 Closed Limelike Curves 的答案或问题 [there][1] 的答案)。对于 HDP,我同意 Closed Limelike Curves。只有我要补充一点,您可以根据 gist 代码构建 HDP 函数。我也有一个已知分布的后验版本(如您的参考文档第 33 页,图 2.2 中),因为您不需要采样。另一个像封闭的石灰样曲线所示的采样。

于 2021-06-03T20:39:23.723 回答
1

OP编辑了这个问题,所以我给出了一个新的答案。

对于中央可信区间,答案很简单:取每个点的分位数:

lowerBound = quantile(Normal(0, 1), .025)
upperBound = quantile(Normal(0, 1), .975)

这将为您提供一个区间,其中位于下限以下 0.025 的概率x,同样对于上限,加起来为 0.05。

HPD 更难计算。此外,它们往往不太常见,因为它们具有一些中央可信区间不共享的奇怪属性。最简单的方法可能是使用蒙特卡洛算法。用于从正态分布randomSample = rand(Normal(0, 1), 2^12)中抽取样本。2^12(或者,无论您想要多少样本,更多的样本会提供更准确的结果,这些结果受随机机会的影响较小。)然后,对于每个随机点,使用 评估该随机点的概率密度pdf.(randomSample)。然后,选择概率密度最高的 95% 的点;包括最高密度区间中的所有这些点,以及它们之间的任何点(我假设您正在处理像正态一样的单模分布)。

对于正态分布,有更好的方法可以做到这一点,但它们更难概括。

于 2021-06-05T18:59:23.713 回答
0

您正在寻找 ArviZ.jl 以及 Turing.jl 的 MCMCChains。MCMCChains 将为您提供非常基本的绘图功能,例如从每个链估计的 PDF 绘图。ArviZ.jl(Python ArviZ 包的包装器)添加了更多绘图。

于 2021-06-04T01:20:05.540 回答