如何解决计算并绘制 Distributions.jl 中分布的中心可信度和最高后验密度区间
我想 (i) 计算和 (ii) 绘制 distributions.jl 库中分布的中心可信区间和最高后验密度区间。 理想情况下,可以编写自己的函数来计算 CI 和 HPD,然后使用 Plots.jl 绘制它们。但是,我发现实现非常棘手(免责声明:我是 Julia 的新手)。 关于库/要点/存储库的任何建议可以让计算和绘制它们更容易?
上下文
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
预期的最终结果总结在下图或 p.第 33 个 BDA3。
目前找到的资源:
-
gist:虽然使用 PyPlot
解决方法
感谢您更新问题;它带来了新的视角。
要点是正确的;只有它使用早期版本的 Julia。
因此,linspace
应替换为 LinRange
。使用 using PyPlot
代替 using Plots
。
我会将绘图部分更改为以下内容:
plot(cred_x,pdf(B,cred_x),fill=(0,0.9,:orange))
plot!(x,x),title="pdf with 90% region highlighted")
乍一看,CI 的计算似乎是正确的。 (就像来自 Closed Limelike Curves 的答案或来自 [there][1] 问题的答案)。对于 HDP,我同意 Closed Limelike Curves。只有我想补充一点,您可以根据要点代码构建您的 HDP 函数。我还会有一个已知分布的后验版本(如您的参考文档第 33 页,图 2.2),因为您不需要采样。另一个带有像 Closed Limelike Curves 所示的采样。
,OP 编辑了问题,所以我给出了一个新答案。
对于中心可信区间,答案很简单:取每个点的分位数:
lowerBound = quantile(Normal(0,1),.025)
upperBound = quantile(Normal(0,.975)
这将为您提供一个区间,其中 x
的概率低于下限 0.025,与上限类似,加起来为 0.05。
HPD 更难计算。此外,它们往往不太常见,因为它们有一些奇怪的特性,这些特性不是由中央可信区间共享的。最简单的方法可能是使用蒙特卡罗算法。使用 randomSample = rand(Normal(0,2^12)
从正态分布中抽取 2^12
个样本。 (或者不管你想要多少样本,更多的样本会得到更准确的结果,受随机机会的影响更小。)然后,对于每个随机点,使用 pdf.(randomSample)
评估该随机点的概率密度。然后,选择概率密度最高的 95% 的点;包括最高密度区间中的所有这些点,以及它们之间的任何点(我假设您正在处理像正态分布一样的单模分布)。
对于正态分布有更好的方法可以做到这一点,但它们更难概括。
,您正在寻找 ArviZ.jl,以及 Turing.jl 的 MCMCChains。 MCMCChains 将为您提供非常基本的绘图功能,例如从每个链估计的 PDF 的图。 ArviZ.jl(Python ArviZ 包的包装器)添加了更多绘图。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。