Я пытаюсь вычислить 95% вероятный интервал следующего апостериорного распределения. Я не смог найти функцию в R для нее, но правильный ли подход ниже?
x <- seq(0.4,12,0.4)
px <- c(0,0, 0, 0, 0, 0, 0.0002, 0.0037, 0.018, 0.06, 0.22 ,0.43, 0.64,0.7579, 0.7870, 0.72, 0.555, 0.37, 0.24, 0.11, 0.07, 0.02, 0.009, 0.005, 0.0001, 0,0.0002, 0, 0, 0)
plot(x,px, type="l")
mm <- sum(x*px)/sum(px)
var <- (sum((x)^2*px)/sum(px)) - (mm^2)
cat("95% credible interval: ", round(mm -1.96*sqrt(var),3), "-", round(mm + 1.96*sqrt(var),3),"\n")
Ответы:
Как отметил Генри , вы предполагаете нормальное распределение, и это совершенно нормально, если ваши данные соответствуют нормальному распределению, но будут неправильными, если вы не можете принять нормальное распределение для них. Ниже я опишу два разных подхода, которые вы могли бы использовать для неизвестного распределения, учитывая только точки данных
x
и сопутствующие оценки плотностиpx
.Первое, что нужно учитывать, это то, что именно вы хотите обобщить, используя ваши интервалы. Например, вас могут интересовать интервалы, полученные с использованием квантилей, но вас также может интересовать область с наивысшей плотностью (см. Здесь или здесь ) вашего распределения. Хотя в простых случаях, таких как симметричные унимодальные распределения, это не должно иметь большого значения (если таковое имеется), это будет иметь значение для более «сложных» распределений. Как правило, квантили дают интервал, содержащий вероятностную массу, сконцентрированную вокруг медианы (середина от вашего распределения), в то время как область с наивысшей плотностью является областью вокруг мод100α% распределения. Это будет более понятным, если вы сравните два графика на рисунке ниже - квантили «обрезают» распределение по вертикали, в то время как область с самой высокой плотностью «обрезает» его по горизонтали.
Следующее, что нужно рассмотреть, - это как справиться с тем фактом, что у вас есть неполная информация о распределении (при условии, что мы говорим о непрерывном распределении, у вас есть только куча точек, а не функция). Что вы можете сделать с этим, так это взять значения «как есть» или использовать некоторую интерполяцию или сглаживание для получения «промежуточных» значений.
Один из подходов состоит в том, чтобы использовать линейную интерполяцию (см.
?approxfun
В R) или, альтернативно, что-то более гладкое, как сплайны (см.?splinefun
В R). Если вы выбираете такой подход, вы должны помнить, что алгоритмы интерполяции не знают предметной области ваших данных и могут возвращать неверные результаты, такие как значения ниже нуля и т. Д.Второй подход, который вы могли бы рассмотреть, - это использовать плотность ядра / распределение смеси, чтобы приблизить ваше распределение, используя имеющиеся у вас данные. Здесь самое сложное - выбрать оптимальную полосу пропускания.
Далее вы найдете интересующие вас интервалы. Вы можете продолжить численно или путем моделирования.
1a) Выборка для получения квантильных интервалов
1b) Отбор проб для получения области наибольшей плотности
2a) Найти квантили численно
2b) Найти область наибольшей плотности численно
Как вы можете видеть на графиках ниже, в случае унимодального, симметричного распределения оба метода возвращают один и тот же интервал.
Конечно, вы также можете попытаться найти интервал вокруг некоторого центрального значения, такого как и использовать некоторую оптимизацию для поиска подходящего , но два описанных выше подхода, кажется, используются чаще и являются более интуитивными.100α% Pr(X∈μ±ζ)≥α ζ
источник