HelpersMG (version 2025.12.22)

RandomFromHessianOrMCMC: 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.
If method is NULL, replicate.CI is set to 0.
If method is hessian, it will generate replicate.CI random numbers using Hessian matrix with covariance.
If method is se, it will generate replicate.CI random numbers using SE values then without covariance.
If method is mcmc, it will generate replicate.CI random numbers using random samples of the MCMC or the regularThin number if it is a number.

Usage

RandomFromHessianOrMCMC(
  se = NULL,
  Hessian = NULL,
  mcmc = NULL,
  chain = "all",
  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",
  ...
)

Value

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

Arguments

se

A named vector with SE of parameters

Hessian

An Hessian matrix

mcmc

A result from MHalgogen()

chain

MCMC chain to be used or "all"

regularThin

If TRUE, use regular thin for MCMC or use a number

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

Author

Marc Girondot

Details

RandomFromHessianOrMCMC returns random numbers based on Hessian matrix or MCMC

Examples

if (FALSE) {
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")

# Return the two first elements of the MCMC
df <- RandomFromHessianOrMCMC(mcmc=mcmc_run, fitted.parameters=NULL, 
               method="MCMC", replicates = 2, regularThin = c(1, 2))$random

# 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)
}