Bestimmung Region hoher Dichte für eine Verteilung in R

Question

Hintergrund:

Normalerweise R Quantile für bekannte Verteilungen gibt. Von diesen Quantilen decken die unteren 2,5% bis zu den oberen 97,5% 95% der Fläche unter diesen Verteilungen ab.

Frage:

Angenommen, ich habe eine F-Verteilung (df1 = 10, DF2 = 90). In R, wie kann ich die 95% der Fläche unter dieser Verteilung so bestimmen, dass diese 95% nur den Bereich mit hoher Dichte abdecken, nicht die 95%, die R normalerweise gibt (siehe mein R-Code unter)?

Hinweis: Die höchste Dichte ist eindeutig der "Modus" (gestrichelte Linie im Bild unten). Ich denke, man muss vom "Modus" in Richtung der Schwänze gehen.

Hier ist mein R-Code:

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

Abschnitt 25.2 von DBDA2E vollständigen R-Code gibt drei Möglichkeiten zur Bestimmung der höchsten Dichte Intervalle für Verteilungen angegeben: als kumulative Dichtefunktionen, wie Gitter Annäherungen oder als Proben. Für eine kumulative Dichtefunktion heißt die Funktion HDIofICDF(). Es ist in der Dienstprogramme Skript, DBDA2E-utilities.R auf der Website des Buches (oben verlinkt). Hier ist der Code:

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 , ...))) 
}