plot.mcmcComposite: Plot the result of a mcmcComposite object
In HelpersMG: Tools for Environmental Analyses, Ecotoxicology and Various R Functions

plot.mcmcComposite

R Documentation

Plot the result of a mcmcComposite object

Description

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=)

Usage

## S3 method for class '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"
)

Arguments

`x`	A mcmcComposite object
`...`	Graphical parameters to be sent to hist()
`chain`	The chain to use
`parameters`	Name of parameters or "all"
`transform`	Function to be used to transform the variable
`scale.prior`	If TRUE, the prior is scaled at the same size as posterior
`legend`	If FALSE, the legend is not shown; see description
`ylab`	y-label for posterior
`las`	las parameter (orientation of y-axis graduation)
`show.prior`	whould the prior be shown?
`col.prior`	Color for prior curve
`lty.prior`	Type of line for prior curve
`lwd.prior`	Width of line for prior curve
`col.posterior`	Color for posterior histogram
`lty.posterior`	Type of line for posterior histogram
`lwd.posterior`	Width of line for posterior histogram
`ylab.prior`	y-label for prior

Details

plot.mcmcComposite plots the result of a MCMC search

Value

None

Author(s)

Marc Girondot marc.girondot@gmail.com

Examples

## Not run: 
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))



## End(Not run)

HelpersMG documentation built on Oct. 5, 2023, 5:08 p.m.