Function to estimate the aquifer parameters from a pumping test using Markov Chain Monte Carlo simulation.
fit.sampling(ptest, model, method = "amcmc", prior.pdf, prior.parameters, iterations = 10000, burnIn = 0.8 * iterations, seed = 12345, proposal.sigma = NULL, iter_update_par = 100, cov.corr = FALSE)
| ptest | A pumping_test object. |
|---|---|
| model | A character string specifying the model used in the parameter estimation. |
| method | A character string specifying the sampling method used in the parameter estimation. Currently the following methods are supported:
|
| prior.pdf | A character vector with the names of the probability density functions used for each hydraulic parameter. Currently only 'unif' and 'norm' are implemented. |
| prior.parameters | A matrix with the parameters of the probability density functions of the hydraulic parameters used in the drawdown calculations. For the uniform distribution, this parameters are the minimun and maximum, whereas for the normal distribution the mean and standard deviation must be specified. |
| iterations | An integer specifying the number of iterations required for sampling. |
| burnIn | An integer specifying the length of the burn-in period of the Markov Chain. |
| seed | A random seed. |
| proposal.sigma | A numeric vector with the standard deviation of the normal distributions used as proposal distribution. These standard deviations are only used during the first stage of the sampling where the proposals come from a set of independent normal distributions. |
| iter_update_par | An integer specifying the number of iterations to update the covariance matrix of the proposal distribution. |
| cov.corr | A logical flag indicating if the covariance matrix needs to be corrected for positive definiteness. |
A list with the hydraulic parameters of the model (includes transmissivity, storage coefficient and radius of influence) and the fitted parameters (includes a and t0)
Bai. An Adaptive Directional Metropolis-within-Gibbs algorithm. Technical Report, Department of Statistics, University of Toronto, 2009. http://www.utstat.toronto.edu/wordpress/WSFiles/technicalreports/0903.pdf
Roberts, G. O. & Rosenthal, J. S. Examples of adaptive MCMC Journal of Computational and Graphical Statistics, 2009, 18, 349-367.
Christen, J. A. & Fox, C. A general purpose sampling algorithm for continuous distributions (the t-walk). Bayesian Analysis, 2010, 5, 2, 263-281.
Other base functions: additional.parameters<-,
confint.pumping_test,
confint_bootstrap,
confint_jackniffe,
confint_wald, estimated<-,
evaluate, fit.optimization,
fit.parameters<-, fit,
hydraulic.parameter.names<-,
hydraulic.parameters<-,
model.parameters, model<-,
plot.pumping_test,
plot_model_diagnostic,
plot_sample_influence,
plot_uncert,
print.pumping_test,
pumping_test, simulate,
summary.pumping_test
# NOT RUN { # Bayesian estimation of hydraulic parameters for a confined aquifer data("theis") # Create a pumping_test object ptest.theis <- pumping_test("Well1", Q = 0.01388, r = 250, t = theis$t, s = theis$s) # Define prior distributions of hydraulic parameters prior.pdf <- rep("unif", 3) prior.parameters <- matrix(0, 3, 2) # # Transmissivity PDF par prior.parameters[1, ] <- c(1e-5, 1e-2) # Storage coefficient PDF par prior.parameters[2, ] <- c(1e-7, 1e-3) # Std deviation of residuals PDF par prior.parameters[3, ] <- c(1e-05, 0.1) proposal.sigma <- c(0.1, 5, 0.001) # Parameter estimation using AMCMC ptest.theis.mcmc <- fit.sampling(ptest = ptest.theis, model = "theis", method = "amcmc", prior.pdf = prior.pdf, prior.parameters = prior.parameters, iterations = 30000, burnIn = 25000, seed = 12355, proposal.sigma = proposal.sigma, iter_update_par = 300, cov.corr = TRUE) # Assign parameters hydraulic.parameters(ptest.theis) <- ptest.theis.mcmc$hydraulic.parameters[25000:30000,] hydraulic.parameter.names(ptest.theis) <- ptest.theis.mcmc$hydraulic.parameters.names estimated(ptest.theis) <- TRUE # Plot Results cex <- 1.0 plot(ptest.theis, type = "uncertainty", cex = cex, cex.axis = cex, cex.lab = cex, cex.main = cex) # Bayesian estimation using twalk ptest.theis1 <- ptest.theis ptest.theis.mcmc1 <- fit.sampling(ptest = ptest.theis1, model = "theis", method = "twalk", prior.pdf = prior.pdf, prior.parameters = prior.parameters, iterations = 30000) # Assign parameters hydraulic.parameters(ptest.theis1) <- ptest.theis.mcmc1$hydraulic.parameters[25000:30000,] hydraulic.parameter.names(ptest.theis1) <- ptest.theis.mcmc1$hydraulic.parameters.names estimated(ptest.theis1) <- TRUE # Plot Results cex <- 1.0 plot(ptest.theis1, type = "uncertainty", cex = cex, cex.axis = cex, cex.lab = cex, cex.main = cex) # Bayesian estimation of hydraulic parameters for Phreatic aquifer data("boulton") ptest.boulton <- pumping_test("Well1", Q = 0.03, r = 20, t = boulton$t, s = boulton$s) ptest.boulton.fit <- fit(ptest.boulton, 'boulton') prior.pdf <- rep("unif", 5) prior.parameters <- matrix(0, 5, 2) # Transmissivity prior.parameters[1, ] <- c(1e-5, 1e-1) # Storage coefficient prior.parameters[2, ] <- c(1e-6, 1e-1) # Omegad (drainage porosity) prior.parameters[3,] <- c(1e-4, 1e-1) # phi prior.parameters[4, ] <- c(1e-05, 0.1) # Sigma prior.parameters[5, ] <- c(1e-05, 0.1) proposal.sigma <- c(0.1, 5, 5, 0.001, 0.001) ptest.boulton.mcmc <- fit.sampling(ptest = ptest.boulton, model = 'boulton', method = 'amcmc', prior.pdf = prior.pdf, prior.parameters = prior.parameters, iterations = 10000, burnIn = 9000, seed = 12345, proposal.sigma = proposal.sigma, iter_update_par = 200, cov.corr = TRUE) # Assign results hydraulic.parameters(ptest.boulton) <- ptest.boulton.mcmc$hydraulic.parameters[9000:10000,] hydraulic.parameter.names(ptest.boulton) <- ptest.boulton.mcmc$hydraulic.parameters.names estimated(ptest.boulton) <- TRUE # Plot Results cex <- 1.0 plot(ptest.boulton, type = "uncertainty", cex = cex, cex.axis = cex, cex.lab = cex, cex.main = cex) # }