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

Question

背景：

通常，R 会为众所周知的分布提供分位数。在这些分位数中，较低的 2.5% 到较高的 97.5% 覆盖了这些分布下 95% 的区域。

问题：

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

注意：显然，最高密度是“模式”（下图中的虚线）。所以我想，一个人必须从“模式”转向尾巴。

这是我的 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

你试过这个包吗：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))

并且

您甚至可以将它与多变量分布一起用于视觉 2D 高密度区域：

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 2

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 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]