Simulation of Inverse Gaussian Confidence Regions

2006-02-14

The following R code will simulate Wald and Log Likelihood Ratio confidence regions for the two-parameter inverse Gaussian distribution, keeping track of which regions miss the true values of the parameters. For each sample generated, the Wald and LLR regions are drawn on the same plot with the MLE indicated by a +. The true values of the parameters are indicated by a green dot which turns yellow if only the Wald region misses, blue if only the LLR region misses, and red if both regions miss.

You will need to have library(SuppDists) attached. The code for the the functions follows the examples. Note that the first and third examples do 100 simulations; you will see all 100 when they run but only the final one is shown here. The second example shows the results of 20 simulations superimposed.

> library(SuppDists)
		
> IGconfLLRWboot(n=20, mu=10,lambda=10, B=100, conflevel=95)
$"Observed Confidence"
 LLRconf Waldconf 
      95       85 

$Crosstab
       Waldmiss
LLRmiss FALSE TRUE
  FALSE    85   10
  TRUE      0    5

> IGconfLLRWboot(n=20, mu=10, lambda=10, B=20, conflevel=95, refresh=F)
$"Observed Confidence"
 LLRconf Waldconf 
      90       85 

$Crosstab
       Waldmiss
LLRmiss FALSE TRUE
  FALSE    17    1
  TRUE      0    2

> IGconfLLRWboot(n=50, mu=1, lambda=100, conflevel=90)
$"Observed Confidence"
 LLRconf Waldconf 
      91       90 

$Crosstab
       Waldmiss
LLRmiss FALSE TRUE
  FALSE    89    2
  TRUE      1    8

> IGconfLLRWboot(n=50, mu=1, lambda=100, B=1000, conflevel=90)
$"Observed Confidence"
 LLRconf Waldconf 
    88.2     88.6 

$Crosstab
       Waldmiss
LLRmiss FALSE TRUE
  FALSE   871   11
  TRUE     15  103

			

					
> IGconfLLRWboot
function (n, mu, lambda, B = 100, conflevel = 95, refresh = T, 
    ...) 
{
    mgr <- seq(max(mu - 15 * sqrt(mu^3/(n * lambda)), 0), mu + 
        15 * sqrt(mu^3/(n * lambda)), len = 101)
    lgr <- seq(max(lambda - 15 * lambda * sqrt(2/n), 0), lambda + 
        15 * lambda * sqrt(2/n), len = 101)
    LLRmiss <- NULL
    Waldmiss <- NULL
    for (BB in 1:B) {
        xx <- rinvGauss(n, mu, lambda)
        xbar <- mean(xx)
        xbar1 <- mean(1/xx)
        muhat <- xbar
        lambdahat <- 1/(xbar1 - 1/xbar)
        LLRout <- neg2LLRIG(mu, lambda, n, xbar, xbar1) > qchisq(conflevel/100, 
            2)
        LLRmiss <- c(LLRmiss, LLRout)
        Waldout <- n * ((mu - muhat)^2 * lambdahat/muhat^3 + 
            (lambda - lambdahat)^2/(2 * lambdahat^2)) > qchisq(conflevel/100, 
            2)
        Waldmiss <- c(Waldmiss, Waldout)
        contour(mgr, lgr, outer(mgr, lgr, neg2LLRIG, n = n, xbar = xbar, 
            xbar1 = xbar1), levels = qchisq(conflevel/100, 2), 
            labels = paste(conflevel, "%", sep = ""), xlab = "mu", 
            ylab = "lambda", add = (!refresh & BB > 1), ...)
        ellipsem(c(muhat, lambdahat), n * diag(c(lambdahat/muhat^3, 
            1/(2 * lambdahat^2))), qchisq(conflevel/100, 2), 
            lty = 2)
        points(muhat, lambdahat, pch = 3)
        points(mu, lambda, col = c("forestgreen", "yellow", "blue", 
            "red")[1 + 2 * LLRout + Waldout], pch = 16)
    }
    list("Observed Confidence" = c(LLRconf = 100 * (B - sum(LLRmiss))/B, 
        Waldconf = 100 * (B - sum(Waldmiss))/B), Crosstab = table(LLRmiss, 
        Waldmiss))
}

> ellipsem
function (mu, amat, c2, npoints = 100, showcentre = T, ...) 
{
    if (all(dim(amat) == c(2, 2))) {
        eamat <- eigen(amat)
        hlen <- sqrt(c2/eamat$val)
        theta <- angle(eamat$vec[1, 1], eamat$vec[2, 1])
        ellipse(hlen[1], hlen[2], theta, mu[1], mu[2], npoints = npoints, 
            ...)
        if (showcentre) 
            points(mu[1], mu[2], pch = 3)
    }
    invisible()
}

> ellipse
function (hlaxa = 1, hlaxb = 1, theta = 0, xc = 0, yc = 0, newplot = F, 
    npoints = 100, ...) 
{
    a <- seq(0, 2 * pi, length = npoints + 1)
    x <- hlaxa * cos(a)
    y <- hlaxb * sin(a)
    alpha <- angle(x, y)
    rad <- sqrt(x^2 + y^2)
    xp <- rad * cos(alpha + theta) + xc
    yp <- rad * sin(alpha + theta) + yc
    if (newplot) 
        plot(xp, yp, type = "l", ...)
    else lines(xp, yp, ...)
    invisible()
}

> angle
function (x, y) 
{
    angle2 <- function(xy) {
        x <- xy[1]
        y <- xy[2]
        if (x > 0) {
            atan(y/x)
        }
        else {
            if (x < 0 & y != 0) {
                atan(y/x) + sign(y) * pi
            }
            else {
                if (x < 0 & y == 0) {
                  pi
                }
                else {
                  if (y != 0) {
                    (sign(y) * pi)/2
                  }
                  else {
                    NA
                  }
                }
            }
        }
    }
    apply(cbind(x, y), 1, angle2)
}