Last chance! 50% off unlimited learning
Sale ends in
Plot the results within a mcmcComposite object.
If scale.prior is TRUE, another scale is shown at right.
legend can take these values:
FALSE, TRUE, topleft, topright, bottomleft, bottomright, c(x=, y=)
# S3 method for mcmcComposite
plot(
x,
...,
chain = 1,
parameters = 1,
transform = NULL,
scale.prior = TRUE,
legend = "topright",
ylab = "Posterior density",
las = 1,
show.prior = TRUE,
col.prior = "red",
lty.prior = 1,
lwd.prior = 1,
col.posterior = "white",
lty.posterior = 1,
lwd.posterior = 1,
ylab.prior = "Prior density"
)None
A mcmcComposite object
Graphical parameters to be sent to hist()
The chain to use
Name of parameters or "all"
Function to be used to transform the variable
If TRUE, the prior is scaled at the same size as posterior
If FALSE, the legend is not shown; see description
y-label for posterior
las parameter (orientation of y-axis graduation)
whould the prior be shown?
Color for prior curve
Type of line for prior curve
Width of line for prior curve
Color for posterior histogram
Type of line for posterior histogram
Width of line for posterior histogram
y-label for prior
Marc Girondot marc.girondot@gmail.com
plot.mcmcComposite plots the result of a MCMC search
Other mcmcComposite functions:
MHalgoGen(),
as.mcmc.mcmcComposite(),
as.parameters(),
as.quantiles(),
merge.mcmcComposite(),
plot.PriorsmcmcComposite(),
setPriors(),
summary.mcmcComposite()
if (FALSE) {
library(HelpersMG)
require(coda)
x <- 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(1, 1),
Min=c(-3, 0), Max=c(100, 10), Init=c(10, 2), stringsAsFactors = FALSE,
row.names=c('mean', 'sd'))
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=x,
adaptive = TRUE,
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=1)
plot(mcmc_run, xlim=c(0, 20))
plot(mcmc_run, xlim=c(0, 10), parameters="sd")
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=50000, parameters=parameters_mcmc, data=x,
adaptive = TRUE,
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, data=x,
adaptive = TRUE,
likelihood=dnormx, n.chains=1, n.adapt=0, thin=1, trace=1)
#########################
## Example with transform
#########################
x.1<-rnorm(6000, 2.4, 0.6)
x.2<-rlnorm(10000, 1.3,0.1)
X<-c(x.1, x.2)
hist(X,100,freq=FALSE, ylim=c(0,1.5))
Lnormlnorm <- function(par, val) {
p <- invlogit(par["p"])
return(-sum(log(p*dnorm(val, par["m1"], abs(par["s1"]), log = FALSE)+
(1-p)*dlnorm(val, par["m2"], abs(par["s2"]), log = FALSE))))
}
# Mean 1
m1=2.3; s1=0.5
# Mean 2
m2=1.3; s2=0.1
# proportion of category 1 - logit transform
p=0
par<-c(m1=m1, s1=s1, m2=m2, s2=s2, p=p)
result2<-optim(par, Lnormlnorm, method="BFGS", val=X,
hessian=FALSE, control=list(trace=1))
lines(seq(from=0, to=5, length=100),
dnorm(seq(from=0, to=5, length=100),
result2$par["m1"], abs(result2$par["s1"])), col="red")
lines(seq(from=0, to=5, length=100),
dlnorm(seq(from=0, to=5, length=100),
result2$par["m2"], abs(result2$par["s2"])), col="green")
p <- invlogit(result2$par["p"])
paste("Proportion of Gaussian data", p)
lines(seq(from=0, to=5, length=100),
p*dnorm(seq(from=0, to=5, length=100),
result2$par["m1"], result2$par["s1"])+
(1-p)*dlnorm(seq(from=0, to=5, length=100),
result2$par["m2"], result2$par["s2"]), col="blue")
parameters_mcmc <- data.frame(Density=c('dunif', 'dunif', 'dunif', 'dunif', 'dunif'),
Prior1=c(0, 0.001, 0, 0.001, -3),
Prior2=c(10, 10, 10, 10, 3),
SDProp=c(1, 1, 1, 1, 1),
Min=c(0, 0.001, 0, 0.001, -3),
Max=c(10, 10, 10, 10, 3),
Init=result2$par, stringsAsFactors = FALSE,
row.names=c('m1', 's1', 'm2', 's2', 'p'))
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, val=X,
parameters_name = "par",
adaptive = TRUE,
likelihood=Lnormlnorm, n.chains=1,
n.adapt=100, thin=1, trace=100)
plot(mcmc_run, parameters="m1", breaks=seq(from=0, to =10, by=0.1),
legend=c(x=6, y=0.10))
plot(mcmc_run, parameters="p", transform=invlogit, xlim=c(0,1),
breaks=seq(from=0, to=1, by=0.01), legend=c(x=0.6, y=0.10))
plot(mcmc_run, parameters="p", xlim=c(-3,3),
breaks=seq(from=-3, to =3, by=0.05), legend=c(x=1, y= 0.10))
parameters_mcmc <- data.frame(Density=c('dunif', 'dunif', 'dunif', 'dunif', 'dnorm'),
Prior1=c(0, 0.001, 0, 0.001, 0.5),
Prior2=c(10, 10, 10, 10, 1),
SDProp=c(1, 1, 1, 1, 1),
Min=c(0, 0.001, 0, 0.001, -3),
Max=c(10, 10, 10, 10, 3),
Init=result2$par, stringsAsFactors = FALSE,
row.names=c('m1', 's1', 'm2', 's2', 'p'))
mcmc_run_pnorm <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, val=X,
parameters_name = "par",
adaptive = TRUE,
likelihood=Lnormlnorm, n.chains=1,
n.adapt=100, thin=1, trace=100)
plot(mcmc_run_pnorm, parameters="m1", breaks=seq(from=0, to =10, by=0.1),
legend=c(x=6, y=0.10))
plot(mcmc_run_pnorm, parameters="p", transform=invlogit, xlim=c(0,1),
breaks=seq(from=0, to=1, by=0.01), legend=c(x=0.6, y=0.10))
plot(x=mcmc_run_pnorm, parameters="p", xlim=c(-3,3),
breaks=seq(from=-3, to =3, by=0.05), legend=c(x=1, y= 0.10))
# Note that it is more logic to use beta distribution for p as a
# proportion. However p value must be checked to be used in optim
# The use of logit transform can be a problem because it can stuck
# the p value to 1 or 0 during fit.
Lnormlnorm <- function(par, val) {
p <- par["p"]
return(-sum(log(p*dnorm(val, par["m1"], abs(par["s1"]), log = FALSE)+
(1-p)*dlnorm(val, par["m2"], abs(par["s2"]), log = FALSE))))
}
# Example of beta distribution
# Mean is alpha/(alpha+beta)
# Variance is (alpha*beta)/((alpha+beta)^2*(alpha+beta+1))
alpha = 5
beta = 9
plot(x = seq(0.0001, 1, by = .0001),
y = dbeta(seq(0.0001, 1, by = .0001), alpha, beta),
type = "l", ylab="Density", xlab="p", bty="n")
points(x=alpha/(alpha+beta), y=0, pch=4)
segments(x0=alpha/(alpha+beta)-sqrt((alpha*beta)/((alpha+beta)^2*(alpha+beta+1))),
x1=alpha/(alpha+beta)+sqrt((alpha*beta)/((alpha+beta)^2*(alpha+beta+1))),
y0=0, y1=0)
# Use of optim with L-BFGS-B to limit p between 0 and 1 and s > 0
# Mean 1
m1=2.3; s1=0.5
# Mean 2
m2=1.3; s2=0.1
# proportion of category 1 - logit transform
p=0.5
par <- c(m1=m1, s1=s1, m2=m2, s2=s2, p=p)
result2 <- optim(par, Lnormlnorm, method="L-BFGS-B", val=X,
lower = c(-Inf, 0, -Inf, 0, 0),
upper = c(Inf, Inf, Inf, Inf, 1),
hessian=FALSE, control=list(trace=1))
parameters_mcmc <- data.frame(Density=c('dunif', 'dunif', 'dunif', 'dunif', 'dbeta'),
Prior1=c(0, 0.001, 0, 0.001, 5),
Prior2=c(10, 10, 10, 10, 9),
SDProp=c(1, 1, 1, 1, 1),
Min=c(0, 0.001, 0, 0.001, 0),
Max=c(10, 10, 10, 10, 1),
Init=c('m1' = 2.4,
's1' = 0.6,
'm2' = 1.3,
's2' = 0.1,
'p' = 0.5), stringsAsFactors = FALSE,
row.names=c('m1', 's1', 'm2', 's2', 'p'))
mcmc_run_pbeta <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, val=X,
parameters_name = "par",
adaptive = TRUE,
likelihood=Lnormlnorm, n.chains=1,
n.adapt=100, thin=1, trace=100)
plot(mcmc_run_pbeta, parameters="m1", breaks=seq(from=0, to =10, by=0.1),
legend=c(x=6, y=0.10))
plot(mcmc_run_pbeta, parameters="p", xlim=c(0,1),
breaks=seq(from=0, to=1, by=0.01), legend=c(x=0.6, y=2))
}
Run the code above in your browser using DataLab