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=4,
n.adapt=100, thin=1, trace=1)
plot(mcmc_run, xlim=c(0, 20), parameters="mean")
plot(mcmc_run, xlim=c(0, 10), parameters="sd")
plot(x=mcmc_run, what="MarkovChain", ylim=c(0, 15), parameters="mean",
col.posterior = colorRampPalette(c("blue", "grey"), alpha = 0.001)
(mcmc_run$parametersMCMC$n.chains))
plot(mcmc_run, what="MarkovChain", ylim=c(0, 10), parameters="sd",
col.posterior = colorRampPalette(c("blue", "grey"), alpha = 0.001)
(mcmc_run$parametersMCMC$n.chains))
plot(mcmc_run, what="LnL",
col.posterior = colorRampPalette(c("blue", "grey"), alpha = 0.001)
(mcmc_run$parametersMCMC$n.chains))
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=3.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=10))
plot(mcmc_run, parameters="p", xlim=c(-3,3),
breaks=seq(from=-3, to =3, by=0.05), legend=c(x=1, y= 3.10))
plot(mcmc_run, parameters="p", what="MarkovChain")
plot(mcmc_run, parameters="p", transform=invlogit, what="MarkovChain",
ylim=c(0, 1))
plot(mcmc_run, what="lnl", col.posterior = "black")
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=2.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=6))
}
Run the code above in your browser using DataLab