Julia：分位数和置信区间答案

【问题标题】：Julia: quantiles and confidence intervalsJulia：分位数和置信区间
【发布时间】：2017-11-03 17:52:22
【问题描述】：

Julia 是否具有计算分布中包含p% 的密度点的功能？

类似于answer 中提到的scipy.stats norm.ppf 函数

示例：2 边 95% 置信区间：

> norm.ppf(1-(1-0.95)/2)
1.96

> norm.ppf(1-(1+0.95)/2)
-1.96

【问题讨论】：

P.S.我在StatsBase 包中发现了一个zscore 函数。
是的，我会建议 :-)。您应该将其发布为您自己问题的答案，除非您对此不满意？
@MichaelK.Borregaard，我不明白如何使用它。大多数情况下我得到MethodError。如果这对您来说很容易，请发布答案。谢谢！
我明白你的意思了，我被标题弄糊涂了。 zscore 通过减去向量的平均值并除以标准差来转换向量的所有值。它在统计中非常有用。我看到你在概率截止之后 - 这确实是quantile，正如下面的 Dan Getz 所回答的那样。
抱歉措辞不佳，迈克尔。随意编辑问题以澄清它！ :-)

标签： julia

【解决方案1】：

只是为了给答案添加另一个相关的增强，特别是对于贝叶斯后验的用户，我们可以定义medianinterval如下：

medianinterval(d,p = 0.95) = quantile(d,1-(1+p)/2),quantile(d,(1+p)/2)

并且拥有：

julia> medianinterval(Normal())
(-1.9599639845400576, 1.9599639845400576)

但有时更有效（即更短）的间隔将围绕分布的模式。为了解决这个问题，我们可以定义：

function modeinterval(d,p=0.95)
    mcdf = cdf(d,mode(d))
    endpoints = mcdf < p/2 ? (0,p) : mcdf > 1-p/2 ? (1-p,1) : (mcdf-p/2,mcdf+p/2)
    return map(x->quantile(d,x), endpoints)
end

对于正态分布，它并不重要，因为众数也是中位数，但对于其他分布，例如 Beta，我们可以：

julia> modeinterval(Beta(2,8),0.2)
(0.09639068616673087, 0.15355172436770012)

julia> medianinterval(Beta(2,8),0.2)
(0.1498495815725847, 0.21227857915644155)

julia> 0.15355172436770012 - 0.09639068616673087
0.05716103820096925

julia> 0.21227857915644155 - 0.1498495815725847
0.06242899758385684

众数区间以较短的长度覆盖分布的相同部分。相关讨论见Credible interval。

【讨论】：

非常整洁！自我注意：最后两行代码是间隔的“宽度”。

【解决方案2】：

Distributions 包中的 quantile 函数可能是您正在寻找的 (95% CI)。对于您拥有的正态分布：

julia> using Distributions

julia> quantile(Normal(0.0, 1.0),1-(1+0.95)/2)
-1.9599639845400576

julia> quantile(Normal(0.0, 1.0),1-(1-0.95)/2)
1.9599639845400576

同样的函数quantile可以用于其他发行版。

【讨论】：

确实很可能 - 我误解了这个问题:-)
我最终使用了这样的东西：norm_ppf2(; p = .95) = quantile(Normal(0.0, 1.0), 1-(1+p)/2), quantile(Normal(0.0, 1.0), 1-(1-p)/2)
请注意，quantile 可以用于许多发行版。所以更通用的方法是定义medianinterval(d,p=0.95) = quantile(d,1-(1+p)/2),quantile(d,(1+p)/2)。这样你也有medianinterval(Exponential()) ≈ (0.0253,3.6889)。
...和Normal() == Normal(0.0,1.0)，所以根据上面的定义，medianinterval(Normal()) 将给出标准 (-1.96,1.96)
@DanGetz 非对称分布的好例子，谢谢！