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:
