Applying a function to find high density area of a distribution (coding)

Question

Background:

Recently, I came across an R function called HDIofICDF (also see below) that provides two limit-values for any distribution (a unimodal curve) such that from one limit-value to the other limit-value covers 95% highly dense area of that distribution.

The code requires that user put Inverse cdf of a distribution. For R-recognized distributions, this is achieved by "qdistribution name" e.g., qf, qchisq etc., and "," and then the arguments needed for the "qdistribution name". For example for an F distribution finding the two limit values is possible via:

HDIofICDF( qf , df1 = 10 , df2 = 90 ) OR for a Chi-Square Distribution, one can use: HDIofICDF( qchisq, df = 10)

Coding Question:

Suppose I have created my own distribution and have the Inverse cdf of this distribution. My Inverse cdf (similar to a "qdistribution name" in R) is given when I use (just as an example to show what arguments my Inverse cdf uses):

invcdf.posterior(p = .025, t = 2.81, N1 = 10, N2 = 10, rscale = 1 )

Now: given my invcdf.posterior and its arguments above, how can I use the HDIofICDF function to obtain the two limit-values for my distribution?

Here is my R codes I have:

HDIofICDF = function( ICDFname , credMass=0.95 , tol=1e-8 , ... ) {

 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 , ... ) ) )
}
##########################################################################
## Example 1 of use: ##

HDIofICDF( qf , df1 = 10 , df2 = 90 ) ## This is a F distribution working OK

THIS IS HOW I TRIED TO APPLY THE HDIofICDF function TO MY OWN DISTRIBUTION:

HDIofICDF( invcdf.posterior, t = 2.81, N1 = 10 , N2 = 10, rscale = 1 ) ## Not working

The error I get is:

Error in eval(expr, envir, enclos) : argument "t" is missing, with no default Called from: eval(expr, envir, enclos)

Answer 1

Notice that your HDIofICDF function has an argument tol.

Because of partial argument matching, when you call HDIofICDF(invcdf.posterior, t = 2.81, N1 = 10 , N2 = 10, rscale = 1), the t that you specified is assumed to mean tol. (You can verify this by modifying your HDIofICDF function to print tol.) Because of this, when your invcdf.posterior is subsequently called, R complains that the t argument is missing.

Here's a toy example illustrating what is going on:

fun1 <- function(a, b) a + b
fun2 <- function(aaa, ...) fun1(...)
fun2(a = 1, b = 2)
# Error in fun1(...) : argument "a" is missing, with no default
fun2(aaa = 99, a = 1, b = 2)
# [1] 3

As suggested by the example above, to get around your problem, explicitly specify tol in your call to HDIofICDF, e.g.

HDIofICDF( invcdf.posterior, t = 2.81, N1 = 10 , N2 = 10, rscale = 1, tol = 1e-8)