确定 R 中分布的高密度区域

Question

背景：

通常，R 给出 well-known 分布的分位数。在这些分位数中，较低的 2.5% 到较高的 97.5% 覆盖了这些分布下 95% 的区域。

问题：

假设我有一个 F 分布 (df1 = 10, df2 = 90)。在 R 中，我如何确定此分布下的 95% 区域，以便这 95% 仅覆盖高密度区域，而不是 R 通常给出的 95%（参见 我的 R 代码下面)?

注：显然，密度最高的是"mode"（下图中的虚线）。所以我想，必须从 "mode" 向尾巴移动。

这是我的 R 代码：

curve(df(x, 10, 90), 0, 3, ylab = 'Density', xlab = 'F value', lwd = 3)

Mode = ( (10 - 2) / 10 ) * ( 90 / (90 + 2) )

abline(v = Mode, lty = 2)

CI = qf( c(.025, .975), 10, 90)

arrows(CI[1], .05, CI[2], .05, code = 3, angle = 90, length = 1.4, col= 'red' )

points(Mode, .05, pch = 21, bg = 'green', cex = 3)

Answer 1

DBDA2E 的第 25.2 节给出了完整的 R 代码，用于确定以三种方式指定的分布的最高密度区间：累积密度函数、网格近似或样本。对于累积密度函数，该函数称为 HDIofICDF()。它位于本书网站（上面的链接）的实用程序脚本 DBDA2E-utilities.R 中。这是代码：

HDIofICDF = function( ICDFname , credMass=0.95 , tol=1e-8 , ... ) {
  # Arguments:
  # ICDFname is R’s name for the inverse cumulative density function
  # of the distribution.
  # credMass is the desired mass of the HDI region.
  # tol is passed to R’s optimize function.
  # Return value:
  # Highest density interval (HDI) limits in a vector.
  # Example of use: For determining HDI of a beta(30,12) distribution, type
  # > HDIofICDF( qbeta , shape1 = 30 , shape2 = 12 )
  # Notice that the parameters of the ICDFname must be explicitly named;
  # e.g., HDIofICDF( qbeta , 30 , 12 ) does not work.
  # Adapted and corrected from Greg Snow’s TeachingDemos package.
  incredMass = 1.0 - credMass
  intervalWidth = function( lowTailPr , ICDFname , credMass , ... ) {
    ICDFname( credMass + lowTailPr , ... ) - ICDFname( lowTailPr , ... )
  }
  optInfo = optimize( intervalWidth , c( 0 , incredMass ) , ICDFname=ICDFname ,
    credMass=credMass , tol=tol , ... )
  HDIlowTailPr = optInfo$minimum
  return( c( ICDFname( HDIlowTailPr , ... ) ,
    ICDFname( credMass + HDIlowTailPr , ... ) ) )
}

Answer 2

你试过这个包吗：https://github.com/robjhyndman/hdrcde？

按照你的例子：

library(hdrcde)
hdr.den(rf(1000,10,90),prob=95)

您可以使用各种高密度区域，它适用于多模式 pdf。

hdr.den(c(rf(1000,10,90),rnorm(1000,4,1)),prob=c(50,75,95))

并且

你甚至可以将它与多变量分布一起使用到视觉二维高密度区域：

hdrs=c(50,75,95)
x=c(rf(1000,10,90),rnorm(1000,4,1))
y=c(rf(1000,5,50),rnorm(1000,7,1) )
par(mfrow=c(1,3))
hdr.den(x,prob=hdrs,xlab="x")
hdr.den(y,prob=hdrs,xlab="y")
hdr.boxplot.2d(x,y,prob=hdrs,shadecol="red",xlab="x",ylab="y")

Answer 3

使用`stat.extend`包中的`HDR.f`函数

stat.extend package 为 R 中的所有基本发行版及其扩展包中的一些发行版提供了 HDR 函数。它使用基于分布的分位数函数的方法，并自动调整分布的形状（单峰、双峰等）。以下是如何使用该函数计算您感兴趣的 HDR。

#Load library
library(stat.extend)

#Compute HDR for an F distribution
HDR.f(cover.prob = 0.9, df1 = 10, df2 = 20)

        Highest Density Region (HDR) 
 
90.00% HDR for F distribution with 10 numerator degrees-of-freedom and 
20 denominator degrees-of-freedom 
Computed using nlm optimisation with 9 iterations (code = 3) 

[0.220947190373167, 1.99228812929142]

确定 R 中分布的高密度区域

Determining High Density Region for a distribution in R

statistics

r

distribution

bayesian

confidence-interval

使用`stat.extend`包中的`HDR.f`函数

确定 R 中分布的高密度区域

Determining High Density Region for a distribution in R

statistics

r

distribution

bayesian

confidence-interval

使用stat.extend包中的HDR.f函数

使用`stat.extend`包中的`HDR.f`函数