fit.sampling: fit.sampling
In khaors/pumpingtest: R Package for the analysis and evaluation of pumping tests

Description Usage Arguments Value Author(s) References See Also Examples

Function to estimate the aquifer parameters from a pumping test using Markov Chain Monte Carlo simulation.

1
2
3

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: amcmc: Adaptative Markov Chain Monte Carlo proposed by Brooks and Rosenthal(2009). This method requires the specification of the variances of the proposal distributions ( see proposal.sigma below). amcmc1: Adaptative Markov Chain Monte Carlo proposed by Bai(2009). This method does not require the specification of any additional parameter. twalk: transverse-walk method proposed by Christen and Fox (2010). This methods does not require any additional parameter.
`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)

Oscar Garcia-Cabrejo khaors@gmail.com

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)

## End(Not run)