RandomFromHessianOrMCMC: Random numbers based on Hessian matrix or MCMC

In HelpersMG: Tools for Environmental Analyses, Ecotoxicology and Various R Functions

Random numbers based on Hessian matrix or MCMC

Description

If it is very long, use silent parameter to check if something goes wrong.
If replicates is NULL or is 0, or if method is NULL, parameters are just copied into data.frame.

Usage

RandomFromHessianOrMCMC(
  se = NULL,
  Hessian = NULL,
  mcmc = NULL,
  chain = 1,
  regularThin = TRUE,
  MinMax = NULL,
  fitted.parameters = NULL,
  fixed.parameters = NULL,
  method = NULL,
  probs = c(0.025, 0.5, 0.975),
  replicates = 10000,
  fn = NULL,
  silent = FALSE,
  ParTofn = "par",
  ...
)

Arguments

se	A named vector with SE of parameters
Hessian	An Hessian matrix
mcmc	A result from MHalgogen()
chain	MCMC chain to be used
regularThin	If TRUE, use regular thin for MCMC
MinMax	A data.frame with at least two columns: Min and Max and rownames being the variable names
fitted.parameters	The fitted parameters
fixed.parameters	The fixed parameters
method	Can be NULL, "SE", "Hessian", "MCMC", or "PseudoHessianFromMCMC"
probs	Probability for quantiles
replicates	Number of replicates to generate the randoms
fn	The function to apply to each replicate
silent	Should the function display some information
ParTofn	Name of the parameter to send random values to fn
...	Parameters send to fn function

Details

RandomFromHessianOrMCMC returns random numbers based on Hessian matrix or MCMC

Value

Returns a list with three data.frames named random, fn, and quantiles

Author(s)

Marc Girondot [e-mail]

Examples

## Not run: 
library(HelpersMG)
val <- rnorm(100, mean=20, sd=5)+(1:100)/10
# Return -ln L of values in val in Gaussian distribution with mean and sd in par
fitnorm <- function(par, data) {
  -sum(dnorm(data, par["mean"], abs(par["sd"]), log = TRUE))
}
# Initial values for search
p<-c(mean=20, sd=5)
# fit the model
result <- optim(par=p, fn=fitnorm, data=val, method="BFGS", hessian=TRUE)
# Using Hessian
df <- RandomFromHessianOrMCMC(Hessian=result$hessian, 
                              fitted.parameters=result$par, 
                              method="Hessian")$random
hist(df[, 1], main="mean")
hist(df[, 2], main="sd")
plot(df[, 1], df[, 2], xlab="mean", ylab="sd", las=1, bty="n")

# Using MCMC
parameters_mcmc <- data.frame(Density=c('dnorm', 'dlnorm'), 
Prior1=c(10, 0.5), Prior2=c(2, 0.5), SDProp=c(0.35, 0.2), 
Min=c(-3, 0), Max=c(100, 10), Init=c(10, 2), stringsAsFactors = FALSE, 
row.names=c('mean', 'sd'))
# Use of trace and traceML parameters
# trace=1 : Only one likelihood is printed
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=val, 
parameters_name = "par", 
likelihood=fitnorm, n.chains=1, n.adapt=100, thin=1, trace=1)
df <- RandomFromHessianOrMCMC(mcmc=mcmc_run, fitted.parameters=NULL, 
                              method="MCMC")$random
hist(df[, 1], main="mean")
hist(df[, 2], main="sd")
plot(df[, 1], df[, 2], xlab="mean", ylab="sd", las=1, bty="n")

# Using a function fn
fitnorm <- function(par, data, x) { 
  y=par["a"]*(x)+par["b"]
  -sum(dnorm(data, y, abs(par["sd"]), log = TRUE))
}
p<-c(a=0.1, b=20, sd=5)
# fit the model
x <- 1:100
result <- optim(par=p, fn=fitnorm, data=val, x=x, method="BFGS", hessian=TRUE)
# Using Hessian
df <- RandomFromHessianOrMCMC(Hessian=result$hessian, fitted.parameters=result$par, 
                              method="Hessian", 
                              fn=function(par) (par["a"]*(x)+par["b"]))
plot(1:100, val)
lines(1:100, df$quantiles["50%", ])
lines(1:100, df$quantiles["2.5%", ], lty=2)
lines(1:100, df$quantiles["97.5%", ], lty=2)

## End(Not run)

HelpersMG documentation built on Sept. 12, 2024, 9:35 a.m.