saddle.distn.RdApproximate an entire distribution using saddlepoint methods. This
function can calculate simple and conditional saddlepoint distribution
approximations for a univariate quantity of interest. For the simple
saddlepoint the quantity of interest is a linear combination of
W where W is a vector of random variables. For the
conditional saddlepoint we require the distribution of one linear
combination given the values of any number of other linear
combinations. The distribution of W must be one of multinomial,
Poisson or binary. The primary use of this function is to calculate
quantiles of bootstrap distributions using saddlepoint approximations.
Such quantiles are required by the function control to
approximate the distribution of the linear approximation to a
statistic.
saddle.distn(A, u = NULL, alpha = NULL, wdist = "m", type = "simp", npts = 20, t = NULL, t0 = NULL, init = rep(0.1, d), mu = rep(0.5, n), LR = FALSE, strata = NULL, ...)
| A | This is a matrix of known coefficients or a function which returns
such a matrix. If a function then its first argument must be the
point |
|---|---|
| u | If |
| alpha | The alpha levels for the quantiles of the distribution which should be returned. By default the 0.1, 0.5, 1, 2.5, 5, 10, 20, 50, 80, 90, 95, 97.5, 99, 99.5 and 99.9 percentiles are calculated. |
| wdist | The distribution of W. Possible values are |
| type | The type of saddlepoint to be used. Possible values are
|
| npts | The number of points at which the saddlepoint approximation should be calculated and then used to fit the spline. |
| t | A vector of points at which the saddlepoint approximations are
calculated. These points should extend beyond the extreme quantiles
required but still be in the possible range of the bootstrap
distribution. The observed value of the statistic should not be
included in |
| t0 | If |
| init | When |
| mu | The vector of parameter values for the distribution. The default is that the components of W are identically distributed. |
| LR | A logical flag. When |
| strata | A vector giving the strata when the rows of A relate to stratified
data. This is used only when |
| ... | When |
The returned value is an object of class "saddle.distn". See the help
file for saddle.distn.object for a description of such
an object.
The range at which the saddlepoint is used is such that the cdf
approximation at the endpoints is more extreme than required by the
extreme values of alpha. The lower endpoint is found by
evaluating the saddlepoint at the points t0[1]-2*t0[2],
t0[1]-4*t0[2], t0[1]-8*t0[2] etc. until a point is
found with a cdf approximation less than min(alpha)/10, then a
bisection method is used to find the endpoint which has cdf
approximation in the range (min(alpha)/1000,
min(alpha)/10). Then a number of, equally spaced, points are
chosen between the lower endpoint and t0[1] until a total of
npts/2 approximations have been made. The remaining
npts/2 points are chosen to the right of t0[1] in a
similar manner. Any points which are very close to the centre of the
distribution are then omitted as the cdf approximations are not
reliable at the centre. A smoothing spline is then fitted to the
probit of the saddlepoint distribution function approximations at the
remaining points and the required quantiles are predicted from the
spline.
Sometimes the function will terminate with the message
"Unable to find range". There are two main reasons why this may
occur. One is that the distribution is too discrete and/or the
required quantiles too extreme, this can cause the function to be
unable to find a point within the allowable range which is beyond the
extreme quantiles. Another possibility is that the value of
t0[2] is too small and so too many steps are required to find
the range. The first problem cannot be solved except by asking for
less extreme quantiles, although for very discrete distributions the
approximations may not be very good. In the second case using a
larger value of t0[2] will usually solve the problem.
Booth, J.G. and Butler, R.W. (1990) Randomization distributions and saddlepoint approximations in generalized linear models. Biometrika, 77, 787--796.
Canty, A.J. and Davison, A.C. (1997) Implementation of saddlepoint approximations to resampling distributions. Computing Science and Statistics; Proceedings of the 28th Symposium on the Interface 248--253.
Davison, A.C. and Hinkley, D.V. (1997) Bootstrap Methods and their Application. Cambridge University Press.
Jensen, J.L. (1995) Saddlepoint Approximations. Oxford University Press.
# The bootstrap distribution of the mean of the air-conditioning # failure data: fails to find value on R (and probably on S too) air.t0 <- c(mean(aircondit$hours), sqrt(var(aircondit$hours)/12)) if (FALSE) saddle.distn(A = aircondit$hours/12, t0 = air.t0) # alternatively using the conditional poisson saddle.distn(A = cbind(aircondit$hours/12, 1), u = 12, wdist = "p", type = "cond", t0 = air.t0)#> Warning: non-integer x = 11.344718#> Warning: non-integer x = 0.655282#> Warning: non-integer x = 11.344718#> Warning: non-integer x = 0.655282#> Warning: non-integer x = 11.832240#> Warning: non-integer x = 0.167760#> Warning: non-integer x = 11.832240#> Warning: non-integer x = 0.167760#> Warning: non-integer x = 7.444538#> Warning: non-integer x = 4.555462#> Warning: non-integer x = 7.444538#> Warning: non-integer x = 4.555462#> Warning: non-integer x = 5.494449#> Warning: non-integer x = 6.505551#> Warning: non-integer x = 5.494449#> Warning: non-integer x = 6.505551#> Warning: non-integer x = 1.594269#> Warning: non-integer x = 10.405731#> Warning: non-integer x = 1.594269#> Warning: non-integer x = 10.405731#> Warning: non-integer x = 3.544359#> Warning: non-integer x = 8.455641#> Warning: non-integer x = 3.544359#> Warning: non-integer x = 8.455641#> Warning: non-integer x = 4.519404#> Warning: non-integer x = 7.480596#> Warning: non-integer x = 4.519404#> Warning: non-integer x = 7.480596#> Warning: non-integer x = 11.561394#> Warning: non-integer x = 0.438606#> Warning: non-integer x = 11.561394#> Warning: non-integer x = 0.438606#> Warning: non-integer x = 11.290549#> Warning: non-integer x = 0.709451#> Warning: non-integer x = 11.290549#> Warning: non-integer x = 0.709451#> Warning: non-integer x = 11.019703#> Warning: non-integer x = 0.980297#> Warning: non-integer x = 11.019703#> Warning: non-integer x = 0.980297#> Warning: non-integer x = 10.748857#> Warning: non-integer x = 1.251143#> Warning: non-integer x = 10.748857#> Warning: non-integer x = 1.251143#> Warning: non-integer x = 10.478011#> Warning: non-integer x = 1.521989#> Warning: non-integer x = 10.478011#> Warning: non-integer x = 1.521989#> Warning: non-integer x = 10.207165#> Warning: non-integer x = 1.792835#> Warning: non-integer x = 10.207165#> Warning: non-integer x = 1.792835#> Warning: non-integer x = 9.936320#> Warning: non-integer x = 2.063680#> Warning: non-integer x = 9.936320#> Warning: non-integer x = 2.063680#> Warning: non-integer x = 9.665474#> Warning: non-integer x = 2.334526#> Warning: non-integer x = 9.665474#> Warning: non-integer x = 2.334526#> Warning: non-integer x = 8.785225#> Warning: non-integer x = 3.214775#> Warning: non-integer x = 8.785225#> Warning: non-integer x = 3.214775#> Warning: non-integer x = 8.175822#> Warning: non-integer x = 3.824178#> Warning: non-integer x = 8.175822#> Warning: non-integer x = 3.824178#> Warning: non-integer x = 7.566419#> Warning: non-integer x = 4.433581#> Warning: non-integer x = 7.566419#> Warning: non-integer x = 4.433581#> Warning: non-integer x = 6.957016#> Warning: non-integer x = 5.042984#> Warning: non-integer x = 6.957016#> Warning: non-integer x = 5.042984#> Warning: non-integer x = 6.347613#> Warning: non-integer x = 5.652387#> Warning: non-integer x = 6.347613#> Warning: non-integer x = 5.652387#> Warning: non-integer x = 5.738210#> Warning: non-integer x = 6.261790#> Warning: non-integer x = 5.738210#> Warning: non-integer x = 6.261790#> Warning: non-integer x = 5.128807#> Warning: non-integer x = 6.871193#> Warning: non-integer x = 5.128807#> Warning: non-integer x = 6.871193#> #> Saddlepoint Distribution Approximations #> #> Call : #> saddle.distn(A = cbind(aircondit$hours/12, 1), u = 12, wdist = "p", #> type = "cond", t0 = air.t0) #> #> Quantiles of the Distribution #> #> 0.1% 27.4 #> 0.5% 35.4 #> 1.0% 39.7 #> 2.5% 46.7 #> 5.0% 53.5 #> 10.0% 62.5 #> 20.0% 75.3 #> 50.0% 104.5 #> 80.0% 139.0 #> 90.0% 158.8 #> 95.0% 175.9 #> 97.5% 191.2 #> 99.0% 209.6 #> 99.5% 222.4 #> 99.9% 249.5 #> #> Smoothing spline used 20 points in the range 9.8 to 304.7.# Distribution of the ratio of a sample of size 10 from the bigcity # data, taken from Example 9.16 of Davison and Hinkley (1997). ratio <- function(d, w) sum(d$x *w)/sum(d$u * w) city.v <- var.linear(empinf(data = city, statistic = ratio)) bigcity.t0 <- c(mean(bigcity$x)/mean(bigcity$u), sqrt(city.v)) Afn <- function(t, data) cbind(data$x - t*data$u, 1) ufn <- function(t, data) c(0,10) saddle.distn(A = Afn, u = ufn, wdist = "b", type = "cond", t0 = bigcity.t0, data = bigcity)#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> Warning: non-integer counts in a binomial glm!#> #> Saddlepoint Distribution Approximations #> #> Call : #> saddle.distn(A = Afn, u = ufn, wdist = "b", type = "cond", t0 = bigcity.t0, #> data = bigcity) #> #> Quantiles of the Distribution #> #> 0.1% 1.070 #> 0.5% 1.092 #> 1.0% 1.104 #> 2.5% 1.122 #> 5.0% 1.139 #> 10.0% 1.158 #> 20.0% 1.184 #> 50.0% 1.237 #> 80.0% 1.304 #> 90.0% 1.348 #> 95.0% 1.392 #> 97.5% 1.436 #> 99.0% 1.494 #> 99.5% 1.537 #> 99.9% 1.636 #> #> Smoothing spline used 20 points in the range 1.014 to 1.96.# From Example 9.16 of Davison and Hinkley (1997) again, we find the # conditional distribution of the ratio given the sum of city$u. Afn <- function(t, data) cbind(data$x-t*data$u, data$u, 1) ufn <- function(t, data) c(0, sum(data$u), 10) city.t0 <- c(mean(city$x)/mean(city$u), sqrt(city.v)) saddle.distn(A = Afn, u = ufn, wdist = "p", type = "cond", t0 = city.t0, data = city)#> Warning: non-integer x = 0.866400#> Warning: non-integer x = 8.511350#> Warning: non-integer x = 0.622251#> Warning: non-integer x = 0.866400#> Warning: non-integer x = 8.511350#> Warning: non-integer x = 0.622251#> Warning: non-integer x = 3.210844#> Warning: non-integer x = 3.107208#> Warning: non-integer x = 3.681949#> Warning: non-integer x = 3.210844#> Warning: non-integer x = 3.107208#> Warning: non-integer x = 3.681949#> Warning: non-integer x = 2.038622#> Warning: non-integer x = 5.809279#> Warning: non-integer x = 2.152100#> Warning: non-integer x = 2.038622#> Warning: non-integer x = 5.809279#> Warning: non-integer x = 2.152100#> Warning: non-integer x = 1.452511#> Warning: non-integer x = 7.160314#> Warning: non-integer x = 1.387175#> Warning: non-integer x = 1.452511#> Warning: non-integer x = 7.160314#> Warning: non-integer x = 1.387175#> Warning: non-integer x = 1.159455#> Warning: non-integer x = 7.835832#> Warning: non-integer x = 1.004713#> Warning: non-integer x = 1.159455#> Warning: non-integer x = 7.835832#> Warning: non-integer x = 1.004713#> Warning: non-integer x = 3.155629#> Warning: non-integer x = 0.115412#> Warning: non-integer x = 6.728960#> Warning: non-integer x = 3.155629#> Warning: non-integer x = 0.115412#> Warning: non-integer x = 6.728960#> Warning: non-integer x = 0.205022#> Warning: non-integer x = 2.133273#> Warning: non-integer x = 7.661705#> Warning: non-integer x = 0.205022#> Warning: non-integer x = 2.133273#> Warning: non-integer x = 7.661705#> Warning: non-integer x = 1.722845#> Warning: non-integer x = 1.033105#> Warning: non-integer x = 7.244049#> Warning: non-integer x = 1.722845#> Warning: non-integer x = 1.033105#> Warning: non-integer x = 7.244049#> Warning: non-integer x = 1.787431#> Warning: non-integer x = 6.388294#> Warning: non-integer x = 1.824275#> Warning: non-integer x = 1.787431#> Warning: non-integer x = 6.388294#> Warning: non-integer x = 1.824275#> Warning: non-integer x = 2.415407#> Warning: non-integer x = 4.940756#> Warning: non-integer x = 2.643837#> Warning: non-integer x = 2.415407#> Warning: non-integer x = 4.940756#> Warning: non-integer x = 2.643837#> Warning: non-integer x = 3.043383#> Warning: non-integer x = 3.493218#> Warning: non-integer x = 3.463399#> Warning: non-integer x = 3.043383#> Warning: non-integer x = 3.493218#> Warning: non-integer x = 3.463399#> Warning: non-integer x = 3.671359#> Warning: non-integer x = 2.045680#> Warning: non-integer x = 4.282961#> Warning: non-integer x = 3.671359#> Warning: non-integer x = 2.045680#> Warning: non-integer x = 4.282961#> Warning: non-integer x = 4.299336#> Warning: non-integer x = 0.598142#> Warning: non-integer x = 5.102523#> Warning: non-integer x = 4.299336#> Warning: non-integer x = 0.598142#> Warning: non-integer x = 5.102523#> Warning: non-integer x = 3.433935#> Warning: non-integer x = 4.409277#> Warning: non-integer x = 2.156788#> Warning: non-integer x = 3.433935#> Warning: non-integer x = 4.409277#> Warning: non-integer x = 2.156788#> Warning: non-integer x = 3.360158#> Warning: non-integer x = 3.271008#> Warning: non-integer x = 3.368834#> Warning: non-integer x = 3.360158#> Warning: non-integer x = 3.271008#> Warning: non-integer x = 3.368834#> Warning: non-integer x = 3.319252#> Warning: non-integer x = 2.639889#> Warning: non-integer x = 4.040859#> Warning: non-integer x = 3.319252#> Warning: non-integer x = 2.639889#> Warning: non-integer x = 4.040859#> Warning: non-integer x = 3.278346#> Warning: non-integer x = 2.008770#> Warning: non-integer x = 4.712884#> Warning: non-integer x = 3.278346#> Warning: non-integer x = 2.008770#> Warning: non-integer x = 4.712884#> Warning: non-integer x = 3.237440#> Warning: non-integer x = 1.377650#> Warning: non-integer x = 5.384909#> Warning: non-integer x = 3.237440#> Warning: non-integer x = 1.377650#> Warning: non-integer x = 5.384909#> Warning: non-integer x = 3.196534#> Warning: non-integer x = 0.746531#> Warning: non-integer x = 6.056935#> Warning: non-integer x = 3.196534#> Warning: non-integer x = 0.746531#> Warning: non-integer x = 6.056935#> Warning: non-integer x = 3.155629#> Warning: non-integer x = 0.115412#> Warning: non-integer x = 6.728960#> Warning: non-integer x = 3.155629#> Warning: non-integer x = 0.115412#> Warning: non-integer x = 6.728960#> Warning: non-integer x = 2.734728#> Warning: non-integer x = 0.299661#> Warning: non-integer x = 6.965612#> Warning: non-integer x = 2.734728#> Warning: non-integer x = 0.299661#> Warning: non-integer x = 6.965612#> Warning: non-integer x = 2.228786#> Warning: non-integer x = 0.666383#> Warning: non-integer x = 7.104831#> Warning: non-integer x = 2.228786#> Warning: non-integer x = 0.666383#> Warning: non-integer x = 7.104831#> #> Saddlepoint Distribution Approximations #> #> Call : #> saddle.distn(A = Afn, u = ufn, wdist = "p", type = "cond", t0 = city.t0, #> data = city) #> #> Quantiles of the Distribution #> #> 0.1% 1.216 #> 0.5% 1.236 #> 1.0% 1.248 #> 2.5% 1.272 #> 5.0% 1.301 #> 10.0% 1.340 #> 20.0% 1.393 #> 50.0% 1.502 #> 80.0% 1.618 #> 90.0% 1.680 #> 95.0% 1.732 #> 97.5% 1.777 #> 99.0% 1.830 #> 99.5% 1.866 #> 99.9% 1.938 #> #> Smoothing spline used 20 points in the range 1.182 to 2.061.