Блог

Я хотел бы (i) вычислить и (ii) построить центральный вероятный интервал и самые высокие интервалы задней плотности для распределения в распределениях.библиотека jl. В идеале можно написать собственную функцию для вычисления CI и HPD, а затем использовать Plots.jl для их построения. Тем не менее, я нахожу реализацию довольно сложной (отказ от ответственности: Я новичок в Джулии). Есть какие-либо предложения по библиотекам/gist/репозиториям, которые можно проверить, чтобы упростить их вычисления и построение графиков?

Контекст

 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
 

Ожидаемый конечный результат, обобщенный на изображении ниже или на стр. 33 BDA3.

Ресурсы, найденные до сих пор:

• суть: использует PyPlot, хотя
• hdrcde Пакет R

Спасибо за обновление вопроса; это открывает новую перспективу.

Суть в некотором роде верна; только в ней используется более ранняя версия Джулии. Следовательно 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 кажется правильным. (Например, ответ из Замкнутых кривых, похожих на Лимбы, или ответ из вопроса [там][1]). Что касается HDP, я согласен с Замкнутыми кривыми, похожими на Лимбы. Только я бы добавил, что вы могли бы построить свою функцию HDP на основе кода gist. У меня также была бы версия для posterior с известным дистрибутивом (как в вашем справочном документе, стр. 33, рис. 2.2), поскольку вам не нужно пробовать. И еще один с выборкой, похожей на Замкнутые кривые, похожие на Известь.

ОП отредактировал вопрос, поэтому я даю новый ответ.

Для центральных достоверных интервалов ответ довольно прост: возьмите квантили в каждой точке:

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

Это даст вам интервал, в котором вероятность x лежать ниже нижней границы 0,025, и аналогично для верхней границы, в сумме до 0,05.

HPD сложнее рассчитать. Кроме того, они, как правило, встречаются реже, потому что у них есть некоторые странные свойства, которые не разделяются центральными надежными интервалами. Самый простой способ сделать это, вероятно, с помощью алгоритма Монте-Карло. Используется randomSample = rand(Normal(0, 1), 2^12) для извлечения 2^12 выборок из нормального распределения. (Или сколько бы образцов вы ни хотели, большее количество дает более точные результаты, на которые меньше влияет случайность.) Затем для каждой случайной точки оцените плотность вероятности в этой случайной точке с помощью pdf.(randomSample) . Затем выберите 95% точек с наибольшей плотностью вероятности; включите все эти точки в интервал наибольшей плотности и любые точки между ними (я предполагаю, что вы имеете дело с одномодовым распределением, подобным нормальному).

Есть лучшие способы сделать это для нормального распределения, но их труднее обобщить.

Вы ищете ArviZ.jl вместе с Тьюрингом.MCMCChains Джей Эл. MCMCChains предоставит вам очень простые возможности построения графиков, например, график PDF, оцененный по каждой цепочке. ArviZ.jl (оболочка вокруг пакета Python ArviZ) добавляет гораздо больше сюжетов.