# Parametric bootstrap for uncertainty of parameter

Problem: Write a parametric bootstrap algorithm to compute the uncertainty in $$\hat{\tau}_{MM}$$ using 500 bootstrap samples.

Now, $$\hat{\tau}$$ is an estimate of the parameter in my PDF of a Rayleigh distribution: $$f(x)=\tau x \exp\left(\frac{-\tau x^2}{2}\right)$$ calculated using the method of moments.

From what I know, parametric bootstrap is hypothesis testing, but I am not sure how to approach this problem. Things I've done include:

• Calculating $$\hat{\tau}_{MM}$$ using method of moments:

$$\hat{\tau}_{MM}=\frac{\pi}{2 \bar{x}^2}$$

• Creating a vector of 500 samples, using Inverse Transform method in R:

 set.seed(3)
# Number of random data points drawn
n <- 500
# Pseudo-random number generator algorithm used to produce points from U(0,1)
u <- runif(n)
x_INVERSE <- (10*log(1/(1-u)))^(1/2)


I could bootstrap this sample 500 times to give me 500 bootstrap samples, but I don't know what to do with it.

From what I know, parametric bootstrapping is hypothesis testing, but I am not sure how to approach this problem.

Not really. It’s typically used to provide estimates of various statistics (such as the mean, or in your case $\hat{\tau}_{MM}$), by sampling with replacement. In nonparametric bootstrapping we would resample from existing data samples (interesting concept, right?). In your case of parametric bootstrapping, we sample from a distribution.

You seem to already have got the “parametric” part of “parametric bootstrap” thanks to your sampler for $x$:

u <- runif(n)                            # Pseudo-random number generator algorithm used to produce points from U(0,1)
x_INVERSE<- (10*log(1/(1-u)))^(1/2)      # Your parametric sampling


And you are right to say that you would repeat this procedure for a certain number of iterations (typically at least a thousand). Within each iteration, you would calculate the value of $\hat{\tau}_{MM}^*$ (where the $*$ superscript means it is bootstrapped) based on those samples, and collect them for later:

set.seed(3)
bootstrapped_bar_xs <- numeric(0)
bootstrapped_hat_tau_mms <- numeric(0)
n <- 500                                 # Number of random data points drawn
it <- 1000                               # Number of iterations

for (i in 1:it) {
u <- runif(n)                            # Pseudo-random number generator algorithm used to produce points from U(0,1)
x_INVERSE<- (10*log(1/(1-u)))^(1/2)      # Your parametric sampling

bar_x <- mean(x_INVERSE)
hat_tau_mm <- pi / (2 * bar_x^2)

bootstrapped_bar_xs <- c(bootstrapped_bar_xs, bar_x)
bootstrapped_hat_tau_mms <- c(bootstrapped_hat_tau_mms, hat_tau_mm)
}


(Note that the R code I’ve given you is rather inefficient because it uses looping. Vectorizing the code would make it more efficient, at the cost of making it a tiny bit less readable for this pedagogical purpose.)

Let's take stock of what we have so far. We have $1000$ sets of resamples. Here's the first hundred, with each line representing a smooth histogram over 500 samples:

From those, we have calculated $1000$ estimates of the sample mean, each corresponding to one of those $1000$ sets (of $500$ samples per set):

And for each of those $1000$ bootstrapped means $\bar{x}^*$, we have estimated your statistic of interest $\hat{\tau}_{MM}^* = \frac{\pi}{2 {\bar{x}^*}^2}$:

Now comes the payoff.

One measure of uncertainty for the parameter $\hat{\tau}_{MM}$ is its standard error. And the estimated standard error is simply the standard deviation of the boostrapped values, $\hat{\text{s.e.}}(\hat{\tau}_{MM}^*) = 0.0093$

se_hat_tau_mm <- sd(bootstrapped_hat_tau_mms)
print(paste("Standard error of the estimator, estimated by bootstrapping, is", se_hat_tau_mm))


Where you go from here depends on what exactly your problem is - what you already have and what you are looking for.

For example, if you don't have your own estimate for $\hat{\tau}_{MM}$ you could calculate a bootstrapped estimate $\bar{\hat{\tau}_{MM}^*} = 0.2005$:

bar_hat_tau_mm <- mean(bootstrapped_hat_tau_mms)
print(paste("Bootstrapped method of moments estimator is", bar_hat_tau_mm))


Then you could assume that $\hat{\text{s.e.}}(\hat{\tau}_{MM}^*)$ is normally-distributed, so an approximate 95% confidence interval is situated at $\bar{(\hat{\tau}_{MM}^*)} \pm 2 \times \hat{\text{s.e.}}(\hat{\tau}_{MM}^*)$ - i.e. from $0.2005 - 2 \times 0.0093$ to $0.2005 + 2 \times 0.0093$:

Similarly, you could get other confidence intervals (90%, 99%, etc) by scaling the s.e. up or down according to the quantiles of a Normal distribution (or a t-distribution if you prefer).

Another measure of uncertainty would be to build an empirical confidence interval. For this, you would need to have your own value for the estimator, say $\hat{\tau}_{MM} = 0.2$ exactly. Then you could calculate the differences between the bootstrapped estimates and your own estimate, $\delta^* = \hat{\tau}_{MM}^* - \hat{\tau}_{MM}$:

my_hat_tau_mm <- 0.2
differences <- bootstrapped_hat_tau_mms - my_hat_tau_mm


which would look like this:

To calculate a confidence interval, you would now pull out quantiles from the empirical distribution of $\delta^*$ and combine them with your estimate. Taking the example of a 95% confidence interval around $\hat{\tau}_{MM}^*$, you would use the 2.5%th and 97.5%th quantiles of $\delta^*$. The answer would be slightly different:

quantiles <- quantile(differences, c(0.025,0.975))
print(paste("2.5% and 97.5% quantiles of differences are", paste(quantiles, collapse=", ")))
interval_bounds <- my_hat_tau_mm + quantiles
print(paste("95% confidence interval of my estimate is", paste(interval_bounds, collapse=" to ")))


The previous approach resulted in a symmetric confidence interval, while this one results in a (more accurate) asymmetric one:

References: