MHalgoGen: Monte-Carlo Markov-chain with Metropolis-Hastings algorithm
In HelpersMG: Tools for Environmental Analyses, Ecotoxicology and Various R Functions
MHalgoGen R Documentation
Monte-Carlo Markov-chain with Metropolis-Hastings algorithm
Description
The parameters must be stored in a data.frame with named rows for each parameter with the following columns:
Density. The density function name, example dnorm, dlnorm, dunif, dbeta
Prior1. The first parameter to send to the Density function
Prior2. The second parameter to send to the Density function
SDProp. The standard error from new proposition value of this parameter
Min. The minimum value for this parameter
Max. The maximum value for this parameter
Init. The initial value for this parameter
This script has been deeply modified from a MCMC script provided by Olivier Martin (INRA, Paris-Grignon).
The likelihood function must use a parameter named parameters_name for the nammed parameters.
For adaptive mcmc, see:
Rosenthal, J. S. 2011. Optimal Proposal Distributions and Adaptive MCMC. Pages 93-112 in S. Brooks, A. Gelman,
G. Jones, and X.-L. Meng, editors. MCMC Handbook. Chapman and Hall/CRC.
Usage
MHalgoGen(
likelihood = stop("A likelihood function must be supplied"),
parameters = stop("Priors must be supplied"),
...,
parameters_name = "x",
n.iter = 10000,
n.chains = 1,
n.adapt = 100,
thin = 30,
trace = FALSE,
traceML = FALSE,
progress.bar.ini = NULL,
progress.bar = NULL,
adaptive = FALSE,
adaptive.lag = 500,
adaptive.fun = function(x) {
ifelse(x > 0.234, 1.3, 0.7)
},
intermediate = NULL,
filename = "intermediate.Rdata",
previous = NULL,
session = NULL
)
Arguments
likelihood
The function that returns -ln likelihood using data and parameters
parameters
A data.frame with priors; see description and examples
...
Parameters to be transmitted to likelihood function
parameters_name
The name of the parameters in the likelihood function, default is "x"
n.iter
Number of iterations for each chain
n.chains
Number of chains
n.adapt
Number of iteration to stabilize likelihood
thin
Interval for thinning likelihoods
trace
Or FALSE or period to show progress
traceML
TRUE or FALSE to show ML
progress.bar.ini
The command to initialize progress bar
progress.bar
The command to run the progress bar
adaptive
Should an adaptive process for SDProp be used
adaptive.lag
Lag to analyze the SDProp value in an adaptive context
adaptive.fun
Function used to change the SDProp
intermediate
Or NULL of period to save intermediate result
filename
Name of file in which intermediate results are saved
previous
The content of the file in which intermediate results are saved
session
The shiny session
Details
MHalgoGen is a function to use mcmc with Metropolis-Hastings algorithm
Value
A mcmcComposite object with all characteristics of the model and mcmc run
Author(s)
Marc Girondot marc.girondot@gmail.com
See Also
Other mcmcComposite functions:
as.mcmc.mcmcComposite(),
as.parameters(),
as.quantiles(),
merge.mcmcComposite(),
plot.PriorsmcmcComposite(),
plot.mcmcComposite(),
setPriors(),
summary.mcmcComposite()
Examples
## Not run:
library(HelpersMG)
require(coda)
val <- rnorm(30, 10, 2)
dnormx <- function(data, x) {
data <- unlist(data)
return(-sum(dnorm(data, mean=x['mean'], sd=x['sd'], log=TRUE)))
}
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,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=1)
# trace=10 : 10 likelihoods are printed
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=val,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=10)
# trace=TRUE : all likelihoods are printed
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=val,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=TRUE)
# trace=FALSE : No likelihood is printed
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=val,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=FALSE)
# traceML=TRUE : values when likelihood is better are shown
mcmc_run <- MHalgoGen(n.iter=100, parameters=parameters_mcmc, data=val,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=TRUE, traceML=TRUE)
mcmc_run <- MHalgoGen(n.iter=100, parameters=parameters_mcmc, data=val,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=FALSE, traceML=TRUE)
plot(mcmc_run, xlim=c(0, 20))
plot(mcmc_run, xlim=c(0, 10), parameters="sd")
library(graphics)
library(fields)
# show a scatter plot of the result
x <- mcmc_run$resultMCMC[[1]][, 1]
y <- mcmc_run$resultMCMC[[1]][, 2]
marpre <- par(mar=c(4, 4, 2, 6)+0.4)
smoothScatter(x, y)
# show a scale
n <- matrix(0, ncol=128, nrow=128)
xrange <- range(x)
yrange <- range(y)
for (i in 1:length(x)) {
posx <- 1+floor(127*(x[i]-xrange[1])/(xrange[2]-xrange[1]))
posy <- 1+floor(127*(y[i]-yrange[1])/(yrange[2]-yrange[1]))
n[posx, posy] <- n[posx, posy]+1
}
image.plot(legend.only=TRUE, zlim= c(0, max(n)), nlevel=128,
col=colorRampPalette(c("white", blues9))(128))
# Compare with a heatmap
x <- seq(from=8, to=12, by=0.2)
y <- seq(from=1, to=4, by=0.2)
df <- expand.grid(mean=x, sd=y)
df <- cbind(df, L=rep(0, length(nrow(df))))
for (i in 1:nrow(df)) df[i, "L"] <- -sum(dnorm(val, df[i, 1], df[i, 2], log = TRUE))
hm <- matrix(df[, "L"], nrow=length(x))
par(mar = marpre)
image.plot(x=x, y=y, z=hm, las=1)
# Diagnostic function from coda library
mcmcforcoda <- as.mcmc(mcmc_run)
#' heidel.diag(mcmcforcoda)
raftery.diag(mcmcforcoda)
autocorr.diag(mcmcforcoda)
acf(mcmcforcoda[[1]][,"mean"], lag.max=20, bty="n", las=1)
acf(mcmcforcoda[[1]][,"sd"], lag.max=20, bty="n", las=1)
batchSE(mcmcforcoda, batchSize=100)
# The batch standard error procedure is usually thought to
# be not as accurate as the time series methods used in summary
summary(mcmcforcoda)$statistics[,"Time-series SE"]
summary(mcmc_run)
as.parameters(mcmc_run)
lastp <- as.parameters(mcmc_run, index="last")
parameters_mcmc[,"Init"] <- lastp
# The n.adapt set to 1 is used to not record the first set of parameters
# then it is not duplicated (as it is also the last one for
# the object mcmc_run)
mcmc_run2 <- MHalgoGen(n.iter=1000, parameters=parameters_mcmc, x=x, data=val,
likelihood=dnormx, n.chains=1, n.adapt=1, thin=1, trace=1)
mcmc_run3 <- merge(mcmc_run, mcmc_run2)
####### no adaptation, n.adapt must be 0
parameters_mcmc[,"Init"] <- c(mean(x), sd(x))
mcmc_run3 <- MHalgoGen(n.iter=1000, parameters=parameters_mcmc, x=x, data=val,
likelihood=dnormx, n.chains=1, n.adapt=0, thin=1, trace=1)
# Here is how to use adaptive mcmc
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=val, adaptive = FALSE,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=1)
1-rejectionRate(as.mcmc(mcmc_run))
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=val, adaptive = TRUE,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=1)
1-rejectionRate(as.mcmc(mcmc_run))
# To see the dynamics :
var <- "mean"
par(mar=c(4, 4, 1, 1)+0.4)
plot(1:nrow(mcmc_run$resultMCMC[[1]]), mcmc_run$resultMCMC[[1]][, var], type="l",
xlab="Iterations", ylab=var, bty="n", las=1)
# Exemple with a progress bar
val <- rnorm(30, 10, 2)
dnormx <- function(data, x) {
data <- unlist(data)
return(-sum(dnorm(data, mean=x['mean'], sd=x['sd'], log=TRUE)))
}
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'))
# Set up the progress bar
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=val,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=FALSE,
progress.bar.ini=function(n.iter) {
assign("pb", txtProgressBar(min=0, max=n.iter, style=3),
env = parent.frame())},
progress.bar=function(iter) {setTxtProgressBar(get("pb", envir = parent.frame()), iter)})
## End(Not run)
HelpersMG documentation built on Sept. 12, 2024, 9:35 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| MHalgoGen | R Documentation |
Monte-Carlo Markov-chain with Metropolis-Hastings algorithm
Description
The parameters must be stored in a data.frame with named rows for each parameter with the following columns:
Density. The density function name, example
dnorm,dlnorm,dunif,dbetaPrior1. The first parameter to send to the
DensityfunctionPrior2. The second parameter to send to the
DensityfunctionSDProp. The standard error from new proposition value of this parameter
Min. The minimum value for this parameter
Max. The maximum value for this parameter
Init. The initial value for this parameter
This script has been deeply modified from a MCMC script provided by Olivier Martin (INRA, Paris-Grignon).
The likelihood function must use a parameter named parameters_name for the nammed parameters.
For adaptive mcmc, see:
Rosenthal, J. S. 2011. Optimal Proposal Distributions and Adaptive MCMC. Pages 93-112 in S. Brooks, A. Gelman,
G. Jones, and X.-L. Meng, editors. MCMC Handbook. Chapman and Hall/CRC.
Usage
MHalgoGen(
likelihood = stop("A likelihood function must be supplied"),
parameters = stop("Priors must be supplied"),
...,
parameters_name = "x",
n.iter = 10000,
n.chains = 1,
n.adapt = 100,
thin = 30,
trace = FALSE,
traceML = FALSE,
progress.bar.ini = NULL,
progress.bar = NULL,
adaptive = FALSE,
adaptive.lag = 500,
adaptive.fun = function(x) {
ifelse(x > 0.234, 1.3, 0.7)
},
intermediate = NULL,
filename = "intermediate.Rdata",
previous = NULL,
session = NULL
)
Arguments
likelihood |
The function that returns -ln likelihood using data and parameters |
parameters |
A data.frame with priors; see description and examples |
... |
Parameters to be transmitted to likelihood function |
parameters_name |
The name of the parameters in the likelihood function, default is "x" |
n.iter |
Number of iterations for each chain |
n.chains |
Number of chains |
n.adapt |
Number of iteration to stabilize likelihood |
thin |
Interval for thinning likelihoods |
trace |
Or FALSE or period to show progress |
traceML |
TRUE or FALSE to show ML |
progress.bar.ini |
The command to initialize progress bar |
progress.bar |
The command to run the progress bar |
adaptive |
Should an adaptive process for SDProp be used |
adaptive.lag |
Lag to analyze the SDProp value in an adaptive context |
adaptive.fun |
Function used to change the SDProp |
intermediate |
Or NULL of period to save intermediate result |
filename |
Name of file in which intermediate results are saved |
previous |
The content of the file in which intermediate results are saved |
session |
The shiny session |
Details
MHalgoGen is a function to use mcmc with Metropolis-Hastings algorithm
Value
A mcmcComposite object with all characteristics of the model and mcmc run
Author(s)
Marc Girondot marc.girondot@gmail.com
See Also
Other mcmcComposite functions:
as.mcmc.mcmcComposite(),
as.parameters(),
as.quantiles(),
merge.mcmcComposite(),
plot.PriorsmcmcComposite(),
plot.mcmcComposite(),
setPriors(),
summary.mcmcComposite()
Examples
## Not run:
library(HelpersMG)
require(coda)
val <- rnorm(30, 10, 2)
dnormx <- function(data, x) {
data <- unlist(data)
return(-sum(dnorm(data, mean=x['mean'], sd=x['sd'], log=TRUE)))
}
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,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=1)
# trace=10 : 10 likelihoods are printed
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=val,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=10)
# trace=TRUE : all likelihoods are printed
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=val,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=TRUE)
# trace=FALSE : No likelihood is printed
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=val,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=FALSE)
# traceML=TRUE : values when likelihood is better are shown
mcmc_run <- MHalgoGen(n.iter=100, parameters=parameters_mcmc, data=val,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=TRUE, traceML=TRUE)
mcmc_run <- MHalgoGen(n.iter=100, parameters=parameters_mcmc, data=val,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=FALSE, traceML=TRUE)
plot(mcmc_run, xlim=c(0, 20))
plot(mcmc_run, xlim=c(0, 10), parameters="sd")
library(graphics)
library(fields)
# show a scatter plot of the result
x <- mcmc_run$resultMCMC[[1]][, 1]
y <- mcmc_run$resultMCMC[[1]][, 2]
marpre <- par(mar=c(4, 4, 2, 6)+0.4)
smoothScatter(x, y)
# show a scale
n <- matrix(0, ncol=128, nrow=128)
xrange <- range(x)
yrange <- range(y)
for (i in 1:length(x)) {
posx <- 1+floor(127*(x[i]-xrange[1])/(xrange[2]-xrange[1]))
posy <- 1+floor(127*(y[i]-yrange[1])/(yrange[2]-yrange[1]))
n[posx, posy] <- n[posx, posy]+1
}
image.plot(legend.only=TRUE, zlim= c(0, max(n)), nlevel=128,
col=colorRampPalette(c("white", blues9))(128))
# Compare with a heatmap
x <- seq(from=8, to=12, by=0.2)
y <- seq(from=1, to=4, by=0.2)
df <- expand.grid(mean=x, sd=y)
df <- cbind(df, L=rep(0, length(nrow(df))))
for (i in 1:nrow(df)) df[i, "L"] <- -sum(dnorm(val, df[i, 1], df[i, 2], log = TRUE))
hm <- matrix(df[, "L"], nrow=length(x))
par(mar = marpre)
image.plot(x=x, y=y, z=hm, las=1)
# Diagnostic function from coda library
mcmcforcoda <- as.mcmc(mcmc_run)
#' heidel.diag(mcmcforcoda)
raftery.diag(mcmcforcoda)
autocorr.diag(mcmcforcoda)
acf(mcmcforcoda[[1]][,"mean"], lag.max=20, bty="n", las=1)
acf(mcmcforcoda[[1]][,"sd"], lag.max=20, bty="n", las=1)
batchSE(mcmcforcoda, batchSize=100)
# The batch standard error procedure is usually thought to
# be not as accurate as the time series methods used in summary
summary(mcmcforcoda)$statistics[,"Time-series SE"]
summary(mcmc_run)
as.parameters(mcmc_run)
lastp <- as.parameters(mcmc_run, index="last")
parameters_mcmc[,"Init"] <- lastp
# The n.adapt set to 1 is used to not record the first set of parameters
# then it is not duplicated (as it is also the last one for
# the object mcmc_run)
mcmc_run2 <- MHalgoGen(n.iter=1000, parameters=parameters_mcmc, x=x, data=val,
likelihood=dnormx, n.chains=1, n.adapt=1, thin=1, trace=1)
mcmc_run3 <- merge(mcmc_run, mcmc_run2)
####### no adaptation, n.adapt must be 0
parameters_mcmc[,"Init"] <- c(mean(x), sd(x))
mcmc_run3 <- MHalgoGen(n.iter=1000, parameters=parameters_mcmc, x=x, data=val,
likelihood=dnormx, n.chains=1, n.adapt=0, thin=1, trace=1)
# Here is how to use adaptive mcmc
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=val, adaptive = FALSE,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=1)
1-rejectionRate(as.mcmc(mcmc_run))
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=val, adaptive = TRUE,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=1)
1-rejectionRate(as.mcmc(mcmc_run))
# To see the dynamics :
var <- "mean"
par(mar=c(4, 4, 1, 1)+0.4)
plot(1:nrow(mcmc_run$resultMCMC[[1]]), mcmc_run$resultMCMC[[1]][, var], type="l",
xlab="Iterations", ylab=var, bty="n", las=1)
# Exemple with a progress bar
val <- rnorm(30, 10, 2)
dnormx <- function(data, x) {
data <- unlist(data)
return(-sum(dnorm(data, mean=x['mean'], sd=x['sd'], log=TRUE)))
}
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'))
# Set up the progress bar
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=val,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=FALSE,
progress.bar.ini=function(n.iter) {
assign("pb", txtProgressBar(min=0, max=n.iter, style=3),
env = parent.frame())},
progress.bar=function(iter) {setTxtProgressBar(get("pb", envir = parent.frame()), iter)})
## End(Not run)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.