HelpersMG (version 4.2)

RandomFromHessianOrMCMC: Random numbers based on Hessian matrix or MCMC

Description

A data.frame with one column for each parameter

Usage

RandomFromHessianOrMCMC(
  Hessian = NULL,
  mcmc = NULL,
  chain = 1,
  fitted.parameters = NULL,
  fixed.parameters = NULL,
  probs = c(0.025, 0.5, 0.975),
  replicates = 10000,
  fn = NULL,
  silent = FALSE,
  ...
)

Arguments

Hessian

An Hessian matrix

mcmc

A result from MHalgogen()

chain

MCMC chain to be used

fitted.parameters

The fitted parameters

fixed.parameters

The fixed parameters

probs

Probability for quantiles

replicates

Number of replicates to generate

fn

The function to apply to each replicate

silent

Should the function display some information

...

Parameters send to fn function

Value

Returns a list with two data.frames named df_random and df_fn

Details

RandomFromHessianOrMCMC returns random numbers based on Hessian matrix or MCMC

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)$df_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)$df_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
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, 
                              fn=function(par) (par["a"]*(x)+par["b"]))
plot(1:100, val)
lines(1:100, df$quantile["50%", ])
lines(1:100, df$quantile["2.5%", ], lty=2)
lines(1:100, df$quantile["97.5%", ], lty=2)
# }