PAWL-package: PARALLEL ADAPTIVE WANG-LANDAU
In PAWL: Implementation of the PAWL algorithm
Description
Details
Author(s)
Examples
Description
The package implements the Parallel Adaptive Wang-Landau algorithm on various examples.
The provided demos allow to reproduce the figures of the article.
Details
Package: PAWL
Type: Package
Version: 1.0
Date: 2011-08-11
License: GPL (>= 2)
LazyLoad: yes
Depends: mvtnorm
Suggests: ggplot2
The main function is pawl. It takes algorithmic parameters in arguments (see the help of the pawl function),
as well a target distribution. Look at the demos to learn how to specify a target distribution.
Author(s)
Luke Bornn <bornn@stat.harvard.edu>,
Pierre E. Jacob <pierre.jacob.work@gmail.com>
Examples
1
2
3
Example output
Loading required package: mvtnorm
Loading required package: foreach
Loading required package: reshape
Loading required package: ggplot2
demo(discreteexample)
---- ~~~~~~~~~~~~~~~
> #rm(list = ls())
> try(detach(package:PAWL, unload = TRUE), silent = TRUE)
> library(PAWL)
> ### toy example: the state space is made of three states
> # target distribution:
> targetpdf <- c(0.5, 0.4, 0.1)
> # target log probability:
> parameters <- list(logpi = log(targetpdf))
> logdensity <- function(x, parameters){
+ parameters$logpi[x]
+ }
> # we have to specify a proposal mechanism
> # for the Metropolis-Hastings kernel
> # since the (default) continuous
> # gaussian random walk cannot be used
> transitionmatrix <- t(matrix(c(0.5, 0.5, 0.0,
+ 0.3, 0.5, 0.2,
+ 0.0, 0.6, 0.4), ncol = 3))
> proposalparam <- list(transitionmatrix = transitionmatrix, card = 3)
> # function that generates proposal:
> rproposal <- function(states, proposalparam){
+ for (index in 1:length(states)){
+ states[index] <- sample(x = 1:proposalparam$card,
+ size = 1, prob = proposalparam$transitionmatrix[states[index],])
+ }
+ return(list(states = states))
+ }
> # function to compute the density of the proposal kernel
> # (necessary to compute the acceptance rate)
> dproposal <- function(states, ys, proposalparam){
+ for (index in 1:(length(states))){
+ states[index] <- log(transitionmatrix[states[index], ys[index]])
+ }
+ return(states)
+ }
> proposalinstance <- proposal(rproposal = rproposal,
+ dproposal = dproposal,
+ proposalparam = proposalparam)
adaptiveproposal unspecified: set to FALSE
sigma_init unspecified: set 1
> # function to draw starting points for the MCMC algorithms:
> rinit <- function(size) return(rep(1, size))
> # define the target
> discretetarget <- target(name = "discrete toy example", dimension = 1, type = "discrete",
+ rinit = rinit, logdensity = logdensity, parameters = parameters)
dinit has to be specified if you want to use SMC
> # specify Metropolis-Hastings tuning parameters:
> mhparameters <- tuningparameters(nchains = 1, niterations = 10000, storeall = TRUE)
> # Rprof(tmp <- tempfile())
> amhresults <- adaptiveMH(discretetarget, mhparameters, proposalinstance)
[1] "Launching Adaptive Metropolis-Hastings algorithm with parameters:"
Object of class tuningparameters.
*number of parallel chains: 1
*number of iterations: 10000
*compute mean: FALSE
*compute mean (burnin): 0
*save every nth iteration: 1
saving every 1 iterations
hence saving a vector of size 10001 x 1 x 1 = 10001
Iteration 200 / 10000
Iteration 400 / 10000
Iteration 600 / 10000
Iteration 800 / 10000
Iteration 1000 / 10000
Iteration 1200 / 10000
Iteration 1400 / 10000
Iteration 1600 / 10000
Iteration 1800 / 10000
Iteration 2000 / 10000
Iteration 2200 / 10000
Iteration 2400 / 10000
Iteration 2600 / 10000
Iteration 2800 / 10000
Iteration 3000 / 10000
Iteration 3200 / 10000
Iteration 3400 / 10000
Iteration 3600 / 10000
Iteration 3800 / 10000
Iteration 4000 / 10000
Iteration 4200 / 10000
Iteration 4400 / 10000
Iteration 4600 / 10000
Iteration 4800 / 10000
Iteration 5000 / 10000
Iteration 5200 / 10000
Iteration 5400 / 10000
Iteration 5600 / 10000
Iteration 5800 / 10000
Iteration 6000 / 10000
Iteration 6200 / 10000
Iteration 6400 / 10000
Iteration 6600 / 10000
Iteration 6800 / 10000
Iteration 7000 / 10000
Iteration 7200 / 10000
Iteration 7400 / 10000
Iteration 7600 / 10000
Iteration 7800 / 10000
Iteration 8000 / 10000
Iteration 8200 / 10000
Iteration 8400 / 10000
Iteration 8600 / 10000
Iteration 8800 / 10000
Iteration 9000 / 10000
Iteration 9200 / 10000
Iteration 9400 / 10000
Iteration 9600 / 10000
Iteration 9800 / 10000
Iteration 10000 / 10000
> # Rprof()
> # print(summaryRprof(tmp))
> # unlink(tmp)
> chains <- ConvertResults(amhresults)
[1] "converting chains before plotting..."
[1] "...done"
> cat("AMH: target probabilities:", targetpdf, "\n")
AMH: target probabilities: 0.5 0.4 0.1
> amhcount <- tabulate(chains$X1, nbins = 3)
> cat("AMH: obtained frequencies:", amhcount / sum(amhcount), "\n")
AMH: obtained frequencies: 0.4884512 0.4087591 0.1027897
> # we bin such that states 1 and 2 are in bin 1, and state 3 is in bin 2
> getPos <- function(points, logdensity) 2 - (points <= 2)
> # we further specify some parameters, like the bins,
> # the desired frequency in each bin...
> positionbinning <- binning(position = getPos,
+ name = "position",
+ bins = c(1, 2),
+ desiredfreq = c(0.8, 0.2),
+ useLearningRate = FALSE)
fhthreshold unspecified: set to 0.5
splitThreshold unspecified: set to 0.1
minSimEffort unspecified: set 200
learningrate unspecified: set to x -> x^(-1)
useFH unspecified: set to FALSE
autobinning unspecified: set to FALSE
> pawlresults <- pawl(discretetarget, binning = positionbinning, AP = mhparameters, proposalinstance)
[1] "Launching Particle Wang-Landau algorithm ..."
saving every 1 iterations
hence saving a vector of size 10001 x 1 x 1 = 10001
iteration 200 / 10000
iteration 400 / 10000
iteration 600 / 10000
iteration 800 / 10000
iteration 1000 / 10000
iteration 1200 / 10000
iteration 1400 / 10000
iteration 1600 / 10000
iteration 1800 / 10000
iteration 2000 / 10000
iteration 2200 / 10000
iteration 2400 / 10000
iteration 2600 / 10000
iteration 2800 / 10000
iteration 3000 / 10000
iteration 3200 / 10000
iteration 3400 / 10000
iteration 3600 / 10000
iteration 3800 / 10000
iteration 4000 / 10000
iteration 4200 / 10000
iteration 4400 / 10000
iteration 4600 / 10000
iteration 4800 / 10000
iteration 5000 / 10000
iteration 5200 / 10000
iteration 5400 / 10000
iteration 5600 / 10000
iteration 5800 / 10000
iteration 6000 / 10000
iteration 6200 / 10000
iteration 6400 / 10000
iteration 6600 / 10000
iteration 6800 / 10000
iteration 7000 / 10000
iteration 7200 / 10000
iteration 7400 / 10000
iteration 7600 / 10000
iteration 7800 / 10000
iteration 8000 / 10000
iteration 8200 / 10000
iteration 8400 / 10000
iteration 8600 / 10000
iteration 8800 / 10000
iteration 9000 / 10000
iteration 9200 / 10000
iteration 9400 / 10000
iteration 9600 / 10000
iteration 9800 / 10000
iteration 10000 / 10000
> pawlchains <- ConvertResults(pawlresults)
[1] "converting chains before plotting..."
[1] "...done"
> cat("PAWL: desired frequencies:", positionbinning@desiredfreq, "\n")
PAWL: desired frequencies: 0.8 0.2
> pawlcount <- tabulate(getPos(pawlchains$X1, pawlchains$logdens), nbins = 2)
> cat("PAWL: obtained frequencies:", pawlcount / sum(pawlcount), "\n")
PAWL: obtained frequencies: 0.80012 0.19988
> # show the trace plot of log theta:
> #PlotLogTheta(pawlresults)
>
> counts <- tabulate(pawlchains$X1, nbins = 3)
> counts <- counts[1:2]
> counts <- counts / sum(counts)
> cat("PAWL: obtained proportions inside bin 1:", counts, "\n")
PAWL: obtained proportions inside bin 1: 0.5644839 0.4355161
> cat("PAWL: compared to:", targetpdf[1:2] / sum(targetpdf[1:2]), "\n")
PAWL: compared to: 0.5555556 0.4444444
demo(gaussianexample)
---- ~~~~~~~~~~~~~~~
> # remove all objects
> #graphics.off()
> #rm(list = ls())
> # try to detach the package if it was already loaded
> try(detach(package:PAWL, unload = TRUE), silent = TRUE)
> # load the package
> library(PAWL)
> # starting points for MCMC algorithms
> rinit <- function(size) rnorm(size)
> # target log density function: a gaussian distribution N(mean = 2, sd = 3)
> parameters <- list(mean = 2, sd = 3)
> logdensity <- function(x, parameters) dnorm(x, parameters$mean, parameters$sd, log = TRUE)
> # creating the target object
> gaussiantarget <- target(name = "gaussian", dimension = 1,
+ rinit = rinit, logdensity = logdensity,
+ parameters = parameters)
dinit has to be specified if you want to use SMC
> # setting a seed for the RNG
> set.seed(17)
> #######
> ## Adaptive Metropolis-Hastings
> #######
> #mhparameters <- tuningparameters(nchains = 10, niterations = 1000, storeall = TRUE)
> #amhresults <- adaptiveMH(gaussiantarget, mhparameters)
> #range(amhresults$alllogtarget)
> ## check that it's working
> #PlotHist(results = amhresults, component = 1)
> #curve(dnorm(x, mean = gaussiantarget@parameters$mean,
> # sd = gaussiantarget@parameters$sd), add = TRUE, lwd = 2, col = "red")
>
>
> ######
> # Parallel Adaptive Wang-Landau
> ######
> N <- 10
> T <- 2000
> # here we disable the adaptive proposal to highlight the automatic binning
> # mechanism
> proposal <- createAdaptiveRandomWalkProposal(N, gaussiantarget@dimension,
+ adaptiveproposal = FALSE)
adaptationrate unspecified: set to x -> x^(-0.6)
sigma_init unspecified: set 1
> pawlparameters <- tuningparameters(nchains = N, niterations = T, storeall = TRUE)
> print(pawlparameters)
Object of class tuningparameters.
*number of parallel chains: 10
*number of iterations: 2000
*compute mean: FALSE
*compute mean (burnin): 0
*save every nth iteration: 1
> #########
> # PAWL where we bin along the state space
> getPos <- function(points, logdensity) points
> ## we further specify some parameters, like the bins,
> ## the desired frequency in each bin
> ncuts <- 2
> positionbinning <- binning(position = getPos,
+ name = "position",
+ binrange = c(-20, -3),
+ ncuts = ncuts,
+ useLearningRate = TRUE,
+ autobinning = TRUE)
desiredfreq unspecified: set to 1 / nbins for each bin
fhthreshold unspecified: set to 0.5
splitThreshold unspecified: set to 0.1
minSimEffort unspecified: set 200
learningrate unspecified: set to x -> x^(-1)
useFH unspecified: set to TRUE
> print(positionbinning)
Object of class binning.
*name: position
*autobinning: TRUE
*bins: -Inf -20 -3
*Flat Histogram threshold: 0.5
*split threshold: 0.1
*use Flat Histogram criterion: TRUE
*use learning rate: TRUE
*min number of iterations between FH: 200
> pawlresults <- pawl(gaussiantarget, binning = positionbinning, AP = pawlparameters,
+ proposal = proposal)
[1] "Launching Particle Wang-Landau algorithm ..."
saving every 1 iterations
hence saving a vector of size 2001 x 10 x 1 = 20010
iteration 100 / 2000
/** diagnosis at iteration 100
* current number of bins = 3
desired freq - proportions:
1 0.05841584 -1.058416
* skewed bins: 2
* bins with enough points to be split: 2 3
* bins to split: 2
*/
iteration 200 / 2000
/** diagnosis at iteration 200
* current number of bins = 4
desired freq - proportions:
1 1 -2.516 -0.242
* skewed bins: 3
* bins with enough points to be split: 3 4
* bins to split: 3
*/
iteration 300 / 2000
/** diagnosis at iteration 300
* current number of bins = 5
desired freq - proportions:
1 1 -0.308 -2.48 -0.803
* skewed bins: 3
* bins with enough points to be split: 3 4 5
* bins to split: 3
*/
iteration 400 / 2000
/** diagnosis at iteration 400
* current number of bins = 6
desired freq - proportions:
1 -0.374 -1.616 -2.768 -0.428 -0.158
* skewed bins: 2 5
* bins with enough points to be split: 2 3 4 5 6
* bins to split: 2
*/
iteration 500 / 2000
/** diagnosis at iteration 500
* current number of bins = 7
desired freq - proportions:
1 1 -4.46 -0.128 0.4 -0.884 0.052
* skewed bins: 3 5 6
* bins with enough points to be split: 3 6 7
* bins to split: 3
*/
iteration 600 / 2000
/** diagnosis at iteration 600
* current number of bins = 8
desired freq - proportions:
1 -4.088 -0.728 -1.496 1 1 0.076 0.031
* skewed bins: 2 3 7
* bins with enough points to be split: 2 4 8
* bins to split: 2
*/
iteration 700 / 2000
/** diagnosis at iteration 700
* current number of bins = 9
desired freq - proportions:
0.007 -0.824 -2.84 -0.536 -0.32 0.832 0.928 0.64 0.178
* skewed bins: 2 4 5 7
* bins with enough points to be split: 1 3 9
*/
iteration 800 / 2000
/** diagnosis at iteration 800
* current number of bins = 9
desired freq - proportions:
-0.371 -0.188 -1.172 -0.008 -0.296 0.916 0.904 0.556 0.2125
* skewed bins: 2 4 5
* bins with enough points to be split: 1 3 5 9
* bins to split: 5
*/
iteration 900 / 2000
/** diagnosis at iteration 900
* current number of bins = 10
desired freq - proportions:
-1.1 1 -0.104 0.28 0.952 0.088 0.952 0.952 1 0.4
* skewed bins: 3 5 6 8
* bins with enough points to be split: 1 10
*/
iteration 1000 / 2000
/** diagnosis at iteration 1000
* current number of bins = 10
desired freq - proportions:
-1.232 0.976 0.352 0.616 0.976 0.544 0.976 0.976 1 0.4
* skewed bins: 2 5 6 8
* bins with enough points to be split: 1 10
*/
iteration 1100 / 2000
/** diagnosis at iteration 1100
* current number of bins = 10
desired freq - proportions:
-1.288 0.984 0.568 0.744 0.984 0.696 0.984 0.984 1 0.4
* skewed bins: 2 5 6 8
* bins with enough points to be split: 1 10
*/
iteration 1200 / 2000
/** diagnosis at iteration 1200
* current number of bins = 10
desired freq - proportions:
-1.21175 0.256 0.61 0.772 0.988 0.772 0.988 0.988 1 0.4
* skewed bins: 5 6 8
* bins with enough points to be split: 1 2 10
*/
iteration 1300 / 2000
/** diagnosis at iteration 1300
* current number of bins = 10
desired freq - proportions:
-1.0136 -0.1136 0.1312 0.2224 0.7024 0.6736 0.9904 0.9904 1 0.4
* skewed bins: 3 6 8
* bins with enough points to be split: 1 2 3 4 10
*/
iteration 1400 / 2000
/** diagnosis at iteration 1400
* current number of bins = 10
desired freq - proportions:
-0.844 -0.188 -0.04 -0.012 0.224 0.128 0.896 0.988 0.936 0.3825
* skewed bins: 8
* bins with enough points to be split: 1 2 3 4 6 10
*/
iteration 1500 / 2000
/** diagnosis at iteration 1500
* current number of bins = 10
desired freq - proportions:
-0.7318571 -0.2171429 -0.05942857 0.07428571 0.28 0.1977143 0.8354286 0.9691429 0.9331429 0.2684286
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 10
*/
iteration 1600 / 2000
/** diagnosis at iteration 1600
* current number of bins = 10
desired freq - proportions:
-0.5975 -0.089 0.061 0.082 0.256 0.238 0.793 0.889 0.73 0.167125
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 9 10
*/
iteration 1700 / 2000
/** diagnosis at iteration 1700
* current number of bins = 10
desired freq - proportions:
-0.4866667 0.032 0.1653333 0.184 0.3386667 0.3226667 0.736 0.7493333 0.528 0.08
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 9 10
*/
iteration 1800 / 2000
/** diagnosis at iteration 1800
* current number of bins = 10
desired freq - proportions:
-0.398 0.124 0.1696 0.1696 0.3424 0.1792 0.5656 0.5944 0.4588 0.0478
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
Flat histogram criterion met at iteration 1850 !
iteration 1900 / 2000
iteration 2000 / 2000
> # histogram of the binned coordinate
> PlotHistBin(pawlresults, positionbinning)
[1] "converting chains before plotting..."
[1] "...done"
> #########
> ## PAWL where we bin along the energy (= - log density)
> getPos <- function(points, logdensity) - logdensity
> # we can get the range using pre exploratory mcmc
> preexpresults <- preexplorationAMH(gaussiantarget, N, 1000)
Pre exploration to get log energy range
[1] "Launching Adaptive Metropolis-Hastings algorithm with parameters:"
Object of class tuningparameters.
*number of parallel chains: 10
*number of iterations: 1000
*compute mean: FALSE
*compute mean (burnin): 0
*save every nth iteration: -1
adaptationrate unspecified: set to x -> x^(-0.6)
sigma_init unspecified: set 1
Iteration 100 / 1000
Iteration 200 / 1000
Iteration 300 / 1000
Iteration 400 / 1000
Iteration 500 / 1000
Iteration 600 / 1000
Iteration 700 / 1000
Iteration 800 / 1000
Iteration 900 / 1000
Iteration 1000 / 1000
> binrange <- preexpresults$SuggestedRange
> # or specify it ourselves
> #binrange <- c(2.1, 20)
> # we further specify some parameters, like the bins,
> # the desired frequency in each bin...
> ncuts <- 2
> energybinning <- binning(position = getPos,
+ name = "energy",
+ binrange = binrange,
+ ncuts = ncuts,
+ useLearningRate = TRUE,
+ autobinning = TRUE,
+ alongenergy = TRUE)
desiredfreq unspecified: set to 1 / nbins for each bin
fhthreshold unspecified: set to 0.5
splitThreshold unspecified: set to 0.1
minSimEffort unspecified: set 200
learningrate unspecified: set to x -> x^(-1)
useFH unspecified: set to TRUE
> print(energybinning)
Object of class binning.
*name: energy
*autobinning: TRUE
*bins: -Inf 2.027298 4.683023
*Flat Histogram threshold: 0.5
*split threshold: 0.1
*use Flat Histogram criterion: TRUE
*use learning rate: TRUE
*min number of iterations between FH: 200
> # launching the algorithm...
> pawlresults <- pawl(gaussiantarget, binning = energybinning, AP = pawlparameters,
+ proposal = proposal)
[1] "Launching Particle Wang-Landau algorithm ..."
saving every 1 iterations
hence saving a vector of size 2001 x 10 x 1 = 20010
right end bin reached: disactivating bin diagnosis
iteration 100 / 2000
iteration 200 / 2000
Flat histogram criterion met at iteration 200 !
iteration 300 / 2000
iteration 400 / 2000
Flat histogram criterion met at iteration 400 !
iteration 500 / 2000
iteration 600 / 2000
Flat histogram criterion met at iteration 600 !
iteration 700 / 2000
iteration 800 / 2000
Flat histogram criterion met at iteration 800 !
iteration 900 / 2000
iteration 1000 / 2000
Flat histogram criterion met at iteration 1000 !
iteration 1100 / 2000
iteration 1200 / 2000
Flat histogram criterion met at iteration 1200 !
iteration 1300 / 2000
iteration 1400 / 2000
Flat histogram criterion met at iteration 1400 !
iteration 1500 / 2000
iteration 1600 / 2000
Flat histogram criterion met at iteration 1600 !
iteration 1700 / 2000
iteration 1800 / 2000
Flat histogram criterion met at iteration 1800 !
iteration 1900 / 2000
iteration 2000 / 2000
Flat histogram criterion met at iteration 2000 !
> # histogram of the binned coordinate
> #PlotHistBin(pawlresults, energybinning)
> ##
> #########
>
> #########
> ## results
> #
> getFrequencies(pawlresults, energybinning)
/* Frequency check
*Do the obtained frequencies match the desired frequencies?
*final bins: -Inf 2.027298 4.683023
*corresponding frequencies: 0.3317841 0.332084 0.3361319
*initial bins: -Inf 2.027298 4.683023
*desired frequencies: 0.3333333 0.3333333 0.3333333
*obtained frequencies: 0.3317841 0.332084 0.3361319
[1] 0.3317841 0.3320840 0.3361319
> # Histogram of the chains
> #PlotHist(pawlresults, 1)
> # Trace plot of the log theta penalties
> #print(PlotLogTheta(pawlresults))
> # trace plot of all the variables
> # (here there is only one variable)
> #print(PlotAllVar(pawlresults))
>
>
>
>
>
>
>
demo(mixture2kexample)
---- ~~~~~~~~~~~~~~~~
> #rm(list = ls())
> try(detach(package:PAWL, unload = TRUE), silent = TRUE)
> library(PAWL)
> set.seed(17)
> mixture <- createMixtureTarget()
ncomponents unspecified: setting it to 2
mixtureparameters unspecified: putting some default parameters
$ncomponents
[1] 2
$componentweights
[1] 0.5 0.5
$componentmeans
[1] 0.0 2.5
$componentvariances
[1] 0.3 0.3
no sample provided and no sample size: setting size to 100
generating sample of size 100
> N <- 10
> T <- 5000
> betaindex <- mixture@dimension
> getBeta <- function(points, logdensity) exp(points[,betaindex])
> betabinning <- binning(position = getBeta,
+ name = "beta",
+ binrange = c(2, 16),
+ autobinning = TRUE,
+ fhthreshold = 0.5,
+ useLearningRate = TRUE)
number of cuts (ncuts) unspecified: set to 9
desiredfreq unspecified: set to 1 / nbins for each bin
splitThreshold unspecified: set to 0.1
minSimEffort unspecified: set 200
learningrate unspecified: set to x -> x^(-1)
useFH unspecified: set to TRUE
> print(betabinning)
Object of class binning.
*name: beta
*autobinning: TRUE
*bins: -Inf 2 3.75 5.5 7.25 9 10.75 12.5 14.25 16
*Flat Histogram threshold: 0.5
*split threshold: 0.1
*use Flat Histogram criterion: TRUE
*use learning rate: TRUE
*min number of iterations between FH: 200
> pawlparameters <- tuningparameters(nchains = N, niterations = T, storeall = TRUE)
> print(pawlparameters)
Object of class tuningparameters.
*number of parallel chains: 10
*number of iterations: 5000
*compute mean: FALSE
*compute mean (burnin): 0
*save every nth iteration: 1
> ### Launching the algorithm...
> pawlresults <- pawl(mixture, binning = betabinning, AP = pawlparameters)
[1] "Launching Particle Wang-Landau algorithm ..."
adaptationrate unspecified: set to x -> x^(-0.6)
sigma_init unspecified: set 1
saving every 1 iterations
hence saving a vector of size 5001 x 10 x 7 = 350070
iteration 100 / 5000
/** diagnosis at iteration 100
* current number of bins = 10
desired freq - proportions:
-3.970297 -0.3663366 -0.07920792 0.0990099 0.2871287 0.4950495 0.8514851 0.8613861 0.9207921 0.9009901
* skewed bins: 7 8
* bins with enough points to be split: 1 2 3
*/
iteration 200 / 5000
/** diagnosis at iteration 200
* current number of bins = 10
desired freq - proportions:
-3.487562 0.119403 0.1840796 0.2437811 0.2636816 0.3930348 0.4825871 0.5572139 0.6368159 0.6069652
* skewed bins: 9
* bins with enough points to be split: 1 2 3 4 5 6 7
*/
iteration 300 / 5000
/** diagnosis at iteration 300
* current number of bins = 10
desired freq - proportions:
-3.325581 -0.07641196 0.07973422 0.2724252 0.3488372 0.4252492 0.4784053 0.551495 0.6212625 0.6245847
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 400 / 5000
/** diagnosis at iteration 400
* current number of bins = 10
desired freq - proportions:
-3.244389 0.1620948 0.1795511 0.2668329 0.3017456 0.4114713 0.4164589 0.4613466 0.4962594 0.5486284
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 500 / 5000
/** diagnosis at iteration 500
* current number of bins = 10
desired freq - proportions:
-3.195609 0.1297405 0.253493 0.3313373 0.3772455 0.3752495 0.3932136 0.4291417 0.4491018 0.4570858
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 600 / 5000
/** diagnosis at iteration 600
* current number of bins = 10
desired freq - proportions:
-3.036606 -0.01830283 0.3111481 0.3327787 0.3527454 0.3610649 0.3793677 0.422629 0.437604 0.4575707
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 700 / 5000
/** diagnosis at iteration 700
* current number of bins = 10
desired freq - proportions:
-2.773181 -0.1041369 0.2097004 0.3266762 0.3238231 0.3708987 0.382311 0.403709 0.4265335 0.4336662
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 800 / 5000
/** diagnosis at iteration 800
* current number of bins = 10
desired freq - proportions:
-2.55181 -0.09113608 0.1810237 0.2509363 0.3245943 0.3245943 0.3645443 0.3895131 0.3832709 0.4244694
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 900 / 5000
/** diagnosis at iteration 900
* current number of bins = 10
desired freq - proportions:
-2.379578 -0.08102109 0.1498335 0.2519423 0.3085461 0.3229745 0.3263041 0.3529412 0.3695893 0.3784684
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1000 / 5000
/** diagnosis at iteration 1000
* current number of bins = 10
desired freq - proportions:
-2.217782 -0.07492507 0.1628372 0.2267732 0.2887113 0.2917083 0.3056943 0.3206793 0.3356643 0.3606394
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1100 / 5000
/** diagnosis at iteration 1100
* current number of bins = 10
desired freq - proportions:
-2.016349 0.02270663 0.1407811 0.1762035 0.2179837 0.2579473 0.2752044 0.2960945 0.3042688 0.3251589
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1200 / 5000
/** diagnosis at iteration 1200
* current number of bins = 10
desired freq - proportions:
-1.84846 0.08909242 0.1298918 0.1681932 0.1948376 0.2181515 0.2414654 0.253955 0.2706078 0.2822648
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1300 / 5000
/** diagnosis at iteration 1300
* current number of bins = 10
desired freq - proportions:
-1.70638 0.1145273 0.1283628 0.1583397 0.187548 0.1975404 0.2113759 0.2244427 0.2398155 0.2444274
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1400 / 5000
/** diagnosis at iteration 1400
* current number of bins = 10
desired freq - proportions:
-1.584582 0.06995004 0.1349036 0.1577445 0.1705924 0.1855817 0.1948608 0.2098501 0.2212705 0.2398287
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1500 / 5000
/** diagnosis at iteration 1500
* current number of bins = 10
desired freq - proportions:
-1.479014 0.00799467 0.07661559 0.1672219 0.1832112 0.188541 0.2018654 0.2111925 0.2225183 0.2198534
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1600 / 5000
/** diagnosis at iteration 1600
* current number of bins = 10
desired freq - proportions:
-1.386633 0.04434728 0.06183635 0.105559 0.16802 0.1773891 0.1892567 0.1986259 0.2117427 0.2298563
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1700 / 5000
/** diagnosis at iteration 1700
* current number of bins = 10
desired freq - proportions:
-1.305115 -0.01410935 0.09641387 0.1387419 0.1475603 0.1604938 0.1769547 0.191652 0.1981188 0.2092887
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1800 / 5000
/** diagnosis at iteration 1800
* current number of bins = 10
desired freq - proportions:
-1.232649 0.0294281 0.05830094 0.1071627 0.1338145 0.141588 0.1776791 0.1860078 0.197668 0.2009994
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1900 / 5000
/** diagnosis at iteration 1900
* current number of bins = 10
desired freq - proportions:
-1.167806 0.04523935 0.07627564 0.09942136 0.135718 0.1430826 0.1520252 0.162546 0.1746449 0.1788532
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2000 / 5000
/** diagnosis at iteration 2000
* current number of bins = 10
desired freq - proportions:
-1.109445 -0.006996502 0.07446277 0.109945 0.1344328 0.1529235 0.1554223 0.1589205 0.1629185 0.1674163
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2100 / 5000
/** diagnosis at iteration 2100
* current number of bins = 10
desired freq - proportions:
-1.05664 -0.05425988 0.07615421 0.104712 0.1199429 0.1413613 0.1608758 0.1651594 0.1694431 0.1732508
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2200 / 5000
/** diagnosis at iteration 2200
* current number of bins = 10
desired freq - proportions:
-1.008632 -0.0945025 0.07905498 0.1054066 0.1281236 0.1385734 0.1535666 0.1612903 0.1671967 0.1699228
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2300 / 5000
/** diagnosis at iteration 2300
* current number of bins = 10
desired freq - proportions:
-0.9647979 -0.09039548 0.0821382 0.09821817 0.1142981 0.1403738 0.1460235 0.1499348 0.1594959 0.164711
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2400 / 5000
/** diagnosis at iteration 2400
* current number of bins = 10
desired freq - proportions:
-0.9246147 -0.08663057 0.07455227 0.09829238 0.1249479 0.1316118 0.1336943 0.1474386 0.1474386 0.1532695
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2500 / 5000
/** diagnosis at iteration 2500
* current number of bins = 10
desired freq - proportions:
-0.8876449 -0.06477409 0.06517393 0.0967613 0.1107557 0.1167533 0.127549 0.1363455 0.1491403 0.14994
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2600 / 5000
/** diagnosis at iteration 2600
* current number of bins = 10
desired freq - proportions:
-0.8535179 -0.02383699 0.05690119 0.09034987 0.1022684 0.1149558 0.1199539 0.1222607 0.1318724 0.1387928
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2700 / 5000
/** diagnosis at iteration 2700
* current number of bins = 10
desired freq - proportions:
-0.8219178 0.01406886 0.06553128 0.08811551 0.09440948 0.1025546 0.1058867 0.1121807 0.1181044 0.1210663
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2800 / 5000
/** diagnosis at iteration 2800
* current number of bins = 10
desired freq - proportions:
-0.7925741 0.02677615 0.05176723 0.08782578 0.09210996 0.09603713 0.1024634 0.1074616 0.1128169 0.115316
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2900 / 5000
/** diagnosis at iteration 2900
* current number of bins = 10
desired freq - proportions:
-0.7652534 0.01378835 0.06204757 0.07583592 0.08376422 0.09410548 0.0989314 0.104102 0.1154774 0.117201
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3000 / 5000
/** diagnosis at iteration 3000
* current number of bins = 10
desired freq - proportions:
-0.7397534 0.01332889 0.05564812 0.06864379 0.0853049 0.09596801 0.1002999 0.1016328 0.1066311 0.1122959
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3100 / 5000
/** diagnosis at iteration 3100
* current number of bins = 10
desired freq - proportions:
-0.7158981 0.01289906 0.03611738 0.05449855 0.08577878 0.09126088 0.09996775 0.1057723 0.1106095 0.1189939
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3200 / 5000
/** diagnosis at iteration 3200
* current number of bins = 10
desired freq - proportions:
-0.6935333 0.02967823 0.03561387 0.0659169 0.07216495 0.08059981 0.09122149 0.100906 0.1052796 0.1121525
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3300 / 5000
/** diagnosis at iteration 3300
* current number of bins = 10
desired freq - proportions:
-0.6725235 0.03847319 0.05574068 0.06301121 0.0724023 0.07967283 0.08330809 0.08664041 0.09391094 0.09936383
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3400 / 5000
/** diagnosis at iteration 3400
* current number of bins = 10
desired freq - proportions:
-0.6527492 0.03469568 0.05116142 0.05880623 0.07115554 0.07703617 0.08056454 0.08820935 0.09085563 0.1002646
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3500 / 5000
/** diagnosis at iteration 3500
* current number of bins = 10
desired freq - proportions:
-0.6341045 -0.01770923 0.04970009 0.06798058 0.07597829 0.08169095 0.08626107 0.08968866 0.09625821 0.1042559
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3600 / 5000
/** diagnosis at iteration 3600
* current number of bins = 10
desired freq - proportions:
-0.6164954 -0.0449875 0.0441544 0.06137184 0.07470147 0.08303249 0.09136351 0.09913913 0.1027492 0.1049708
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3700 / 5000
/** diagnosis at iteration 3700
* current number of bins = 10
desired freq - proportions:
-0.5998379 -0.07079168 0.05674142 0.06727911 0.07673602 0.08430154 0.0878141 0.09159687 0.1010538 0.1051067
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3800 / 5000
/** diagnosis at iteration 3800
* current number of bins = 10
desired freq - proportions:
-0.5840568 -0.0952381 0.05945804 0.06840305 0.07945278 0.08392528 0.09181794 0.09760589 0.0973428 0.1012891
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3900 / 5000
/** diagnosis at iteration 3900
* current number of bins = 10
desired freq - proportions:
-0.5690849 -0.1133043 0.05537042 0.07331453 0.07792874 0.08561907 0.09074596 0.09510382 0.1004871 0.1038195
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 4000 / 5000
/** diagnosis at iteration 4000
* current number of bins = 10
desired freq - proportions:
-0.5548613 -0.08672832 0.0419895 0.05898525 0.07223194 0.08097976 0.08922769 0.09522619 0.1009748 0.1019745
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 4100 / 5000
/** diagnosis at iteration 4100
* current number of bins = 10
desired freq - proportions:
-0.5413314 -0.06022921 0.04437942 0.05754694 0.06827603 0.07461595 0.08217508 0.08827115 0.0919288 0.09436723
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 4200 / 5000
/** diagnosis at iteration 4200
* current number of bins = 10
desired freq - proportions:
-0.5284456 -0.03499167 0.04284694 0.05379672 0.06736491 0.07545822 0.07760057 0.080219 0.08307546 0.08307546
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 4300 / 5000
/** diagnosis at iteration 4300
* current number of bins = 10
desired freq - proportions:
-0.516159 -0.01092769 0.02999302 0.057196 0.06231109 0.06626366 0.07207626 0.07649384 0.0802139 0.08253894
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 4400 / 5000
/** diagnosis at iteration 4400
* current number of bins = 10
desired freq - proportions:
-0.5044308 0.01204272 0.03181095 0.05157919 0.0624858 0.06157691 0.066803 0.07089298 0.07498296 0.07225631
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
Flat histogram criterion met at iteration 4439 !
iteration 4500 / 5000
iteration 4600 / 5000
Flat histogram criterion met at iteration 4639 !
iteration 4700 / 5000
iteration 4800 / 5000
Flat histogram criterion met at iteration 4873 !
iteration 4900 / 5000
iteration 5000 / 5000
> chains <- ConvertResults(pawlresults)
[1] "converting chains before plotting..."
[1] "...done"
> getFrequencies(pawlresults, betabinning)
/* Frequency check
*Do the obtained frequencies match the desired frequencies?
*final bins: -Inf 2 3.75 5.5 7.25 9 10.75 12.5 14.25 16
*corresponding frequencies: 0.1443911 0.0989802 0.09906019 0.09636073 0.0939812 0.0940012 0.09374125 0.09336133 0.09312138 0.0930014
*initial bins: -Inf 2 3.75 5.5 7.25 9 10.75 12.5 14.25 16
*desired frequencies: 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
*obtained frequencies: 0.1443911 0.0989802 0.09906019 0.09636073 0.0939812 0.0940012 0.09374125 0.09336133 0.09312138 0.0930014
[1] 0.14439112 0.09898020 0.09906019 0.09636073 0.09398120 0.09400120
[7] 0.09374125 0.09336133 0.09312138 0.09300140
> ## Plot the reaction coordinate values along the binning axis,
> ## with red vertical lines denoting the endpoints of the bins
> # print(PlotHistBin(pawlresults, betabinning))
> ## Plot the log penalties associated with each bin
> # print(PlotLogTheta(pawlresults))
> ## Plot sigma versus iterations, where sigma is
> ## the vector of the standard deviations used by the MH kernel along the iterations.
> #plot(pawlresults$sigma, type = "l")
> #PlotAllVar(chains)
> #X11()
> print(PlotComp1vsComp2(chains, "X3", "X4"))
> #X11()
> #print(PlotHistBin(pawlresults, betabinning))
> #print(PlotHist(pawlresults, component = 7))
> #X11(); print(PlotHist(pawlresults, component = 3))
>
> #PlotFH(pawlresults)
> #binincrease <- data.frame(cbind(c(1, pawlresults$splitTimes, T),
> # c(pawlresults$nbins, max(pawlresults$nbins))))
> #g <- ggplot(binincrease, aes(x = X1, y = X2)) + geom_step()
> #g <- g + ylim(0, 1.5 * max(pawlresults$nbins)) + ylab("Number of bins") + xlab("iterations")
> #g <- g + theme(title = "Number of bins along the iterations")
> #print(g)
>
>
> ### Adaptive MH for comparison
> #mhparameters <- tuningparameters(nchains = N, niterations = T, storeall = TRUE)
> #### launching the algorithm...
> #amhresults <- adaptiveMH(mixture, mhparameters)
>
> #PlotAllVar(amhresults)
> #X11()
> #print(PlotComp1vsComp2(amhresults, "X3", "X4"))
> #X11(); print(PlotHist(amhresults, component = 3))
> #
PAWL documentation built on May 2, 2019, 5:58 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.
Description Details Author(s) Examples
Description
The package implements the Parallel Adaptive Wang-Landau algorithm on various examples. The provided demos allow to reproduce the figures of the article.
Details
| Package: | PAWL |
| Type: | Package |
| Version: | 1.0 |
| Date: | 2011-08-11 |
| License: | GPL (>= 2) |
| LazyLoad: | yes |
| Depends: | mvtnorm |
| Suggests: | ggplot2 |
The main function is pawl. It takes algorithmic parameters in arguments (see the help of the pawl function),
as well a target distribution. Look at the demos to learn how to specify a target distribution.
Author(s)
Luke Bornn <bornn@stat.harvard.edu>, Pierre E. Jacob <pierre.jacob.work@gmail.com>
Examples
1 2 3 |
Example output
Loading required package: mvtnorm
Loading required package: foreach
Loading required package: reshape
Loading required package: ggplot2
demo(discreteexample)
---- ~~~~~~~~~~~~~~~
> #rm(list = ls())
> try(detach(package:PAWL, unload = TRUE), silent = TRUE)
> library(PAWL)
> ### toy example: the state space is made of three states
> # target distribution:
> targetpdf <- c(0.5, 0.4, 0.1)
> # target log probability:
> parameters <- list(logpi = log(targetpdf))
> logdensity <- function(x, parameters){
+ parameters$logpi[x]
+ }
> # we have to specify a proposal mechanism
> # for the Metropolis-Hastings kernel
> # since the (default) continuous
> # gaussian random walk cannot be used
> transitionmatrix <- t(matrix(c(0.5, 0.5, 0.0,
+ 0.3, 0.5, 0.2,
+ 0.0, 0.6, 0.4), ncol = 3))
> proposalparam <- list(transitionmatrix = transitionmatrix, card = 3)
> # function that generates proposal:
> rproposal <- function(states, proposalparam){
+ for (index in 1:length(states)){
+ states[index] <- sample(x = 1:proposalparam$card,
+ size = 1, prob = proposalparam$transitionmatrix[states[index],])
+ }
+ return(list(states = states))
+ }
> # function to compute the density of the proposal kernel
> # (necessary to compute the acceptance rate)
> dproposal <- function(states, ys, proposalparam){
+ for (index in 1:(length(states))){
+ states[index] <- log(transitionmatrix[states[index], ys[index]])
+ }
+ return(states)
+ }
> proposalinstance <- proposal(rproposal = rproposal,
+ dproposal = dproposal,
+ proposalparam = proposalparam)
adaptiveproposal unspecified: set to FALSE
sigma_init unspecified: set 1
> # function to draw starting points for the MCMC algorithms:
> rinit <- function(size) return(rep(1, size))
> # define the target
> discretetarget <- target(name = "discrete toy example", dimension = 1, type = "discrete",
+ rinit = rinit, logdensity = logdensity, parameters = parameters)
dinit has to be specified if you want to use SMC
> # specify Metropolis-Hastings tuning parameters:
> mhparameters <- tuningparameters(nchains = 1, niterations = 10000, storeall = TRUE)
> # Rprof(tmp <- tempfile())
> amhresults <- adaptiveMH(discretetarget, mhparameters, proposalinstance)
[1] "Launching Adaptive Metropolis-Hastings algorithm with parameters:"
Object of class tuningparameters.
*number of parallel chains: 1
*number of iterations: 10000
*compute mean: FALSE
*compute mean (burnin): 0
*save every nth iteration: 1
saving every 1 iterations
hence saving a vector of size 10001 x 1 x 1 = 10001
Iteration 200 / 10000
Iteration 400 / 10000
Iteration 600 / 10000
Iteration 800 / 10000
Iteration 1000 / 10000
Iteration 1200 / 10000
Iteration 1400 / 10000
Iteration 1600 / 10000
Iteration 1800 / 10000
Iteration 2000 / 10000
Iteration 2200 / 10000
Iteration 2400 / 10000
Iteration 2600 / 10000
Iteration 2800 / 10000
Iteration 3000 / 10000
Iteration 3200 / 10000
Iteration 3400 / 10000
Iteration 3600 / 10000
Iteration 3800 / 10000
Iteration 4000 / 10000
Iteration 4200 / 10000
Iteration 4400 / 10000
Iteration 4600 / 10000
Iteration 4800 / 10000
Iteration 5000 / 10000
Iteration 5200 / 10000
Iteration 5400 / 10000
Iteration 5600 / 10000
Iteration 5800 / 10000
Iteration 6000 / 10000
Iteration 6200 / 10000
Iteration 6400 / 10000
Iteration 6600 / 10000
Iteration 6800 / 10000
Iteration 7000 / 10000
Iteration 7200 / 10000
Iteration 7400 / 10000
Iteration 7600 / 10000
Iteration 7800 / 10000
Iteration 8000 / 10000
Iteration 8200 / 10000
Iteration 8400 / 10000
Iteration 8600 / 10000
Iteration 8800 / 10000
Iteration 9000 / 10000
Iteration 9200 / 10000
Iteration 9400 / 10000
Iteration 9600 / 10000
Iteration 9800 / 10000
Iteration 10000 / 10000
> # Rprof()
> # print(summaryRprof(tmp))
> # unlink(tmp)
> chains <- ConvertResults(amhresults)
[1] "converting chains before plotting..."
[1] "...done"
> cat("AMH: target probabilities:", targetpdf, "\n")
AMH: target probabilities: 0.5 0.4 0.1
> amhcount <- tabulate(chains$X1, nbins = 3)
> cat("AMH: obtained frequencies:", amhcount / sum(amhcount), "\n")
AMH: obtained frequencies: 0.4884512 0.4087591 0.1027897
> # we bin such that states 1 and 2 are in bin 1, and state 3 is in bin 2
> getPos <- function(points, logdensity) 2 - (points <= 2)
> # we further specify some parameters, like the bins,
> # the desired frequency in each bin...
> positionbinning <- binning(position = getPos,
+ name = "position",
+ bins = c(1, 2),
+ desiredfreq = c(0.8, 0.2),
+ useLearningRate = FALSE)
fhthreshold unspecified: set to 0.5
splitThreshold unspecified: set to 0.1
minSimEffort unspecified: set 200
learningrate unspecified: set to x -> x^(-1)
useFH unspecified: set to FALSE
autobinning unspecified: set to FALSE
> pawlresults <- pawl(discretetarget, binning = positionbinning, AP = mhparameters, proposalinstance)
[1] "Launching Particle Wang-Landau algorithm ..."
saving every 1 iterations
hence saving a vector of size 10001 x 1 x 1 = 10001
iteration 200 / 10000
iteration 400 / 10000
iteration 600 / 10000
iteration 800 / 10000
iteration 1000 / 10000
iteration 1200 / 10000
iteration 1400 / 10000
iteration 1600 / 10000
iteration 1800 / 10000
iteration 2000 / 10000
iteration 2200 / 10000
iteration 2400 / 10000
iteration 2600 / 10000
iteration 2800 / 10000
iteration 3000 / 10000
iteration 3200 / 10000
iteration 3400 / 10000
iteration 3600 / 10000
iteration 3800 / 10000
iteration 4000 / 10000
iteration 4200 / 10000
iteration 4400 / 10000
iteration 4600 / 10000
iteration 4800 / 10000
iteration 5000 / 10000
iteration 5200 / 10000
iteration 5400 / 10000
iteration 5600 / 10000
iteration 5800 / 10000
iteration 6000 / 10000
iteration 6200 / 10000
iteration 6400 / 10000
iteration 6600 / 10000
iteration 6800 / 10000
iteration 7000 / 10000
iteration 7200 / 10000
iteration 7400 / 10000
iteration 7600 / 10000
iteration 7800 / 10000
iteration 8000 / 10000
iteration 8200 / 10000
iteration 8400 / 10000
iteration 8600 / 10000
iteration 8800 / 10000
iteration 9000 / 10000
iteration 9200 / 10000
iteration 9400 / 10000
iteration 9600 / 10000
iteration 9800 / 10000
iteration 10000 / 10000
> pawlchains <- ConvertResults(pawlresults)
[1] "converting chains before plotting..."
[1] "...done"
> cat("PAWL: desired frequencies:", positionbinning@desiredfreq, "\n")
PAWL: desired frequencies: 0.8 0.2
> pawlcount <- tabulate(getPos(pawlchains$X1, pawlchains$logdens), nbins = 2)
> cat("PAWL: obtained frequencies:", pawlcount / sum(pawlcount), "\n")
PAWL: obtained frequencies: 0.80012 0.19988
> # show the trace plot of log theta:
> #PlotLogTheta(pawlresults)
>
> counts <- tabulate(pawlchains$X1, nbins = 3)
> counts <- counts[1:2]
> counts <- counts / sum(counts)
> cat("PAWL: obtained proportions inside bin 1:", counts, "\n")
PAWL: obtained proportions inside bin 1: 0.5644839 0.4355161
> cat("PAWL: compared to:", targetpdf[1:2] / sum(targetpdf[1:2]), "\n")
PAWL: compared to: 0.5555556 0.4444444
demo(gaussianexample)
---- ~~~~~~~~~~~~~~~
> # remove all objects
> #graphics.off()
> #rm(list = ls())
> # try to detach the package if it was already loaded
> try(detach(package:PAWL, unload = TRUE), silent = TRUE)
> # load the package
> library(PAWL)
> # starting points for MCMC algorithms
> rinit <- function(size) rnorm(size)
> # target log density function: a gaussian distribution N(mean = 2, sd = 3)
> parameters <- list(mean = 2, sd = 3)
> logdensity <- function(x, parameters) dnorm(x, parameters$mean, parameters$sd, log = TRUE)
> # creating the target object
> gaussiantarget <- target(name = "gaussian", dimension = 1,
+ rinit = rinit, logdensity = logdensity,
+ parameters = parameters)
dinit has to be specified if you want to use SMC
> # setting a seed for the RNG
> set.seed(17)
> #######
> ## Adaptive Metropolis-Hastings
> #######
> #mhparameters <- tuningparameters(nchains = 10, niterations = 1000, storeall = TRUE)
> #amhresults <- adaptiveMH(gaussiantarget, mhparameters)
> #range(amhresults$alllogtarget)
> ## check that it's working
> #PlotHist(results = amhresults, component = 1)
> #curve(dnorm(x, mean = gaussiantarget@parameters$mean,
> # sd = gaussiantarget@parameters$sd), add = TRUE, lwd = 2, col = "red")
>
>
> ######
> # Parallel Adaptive Wang-Landau
> ######
> N <- 10
> T <- 2000
> # here we disable the adaptive proposal to highlight the automatic binning
> # mechanism
> proposal <- createAdaptiveRandomWalkProposal(N, gaussiantarget@dimension,
+ adaptiveproposal = FALSE)
adaptationrate unspecified: set to x -> x^(-0.6)
sigma_init unspecified: set 1
> pawlparameters <- tuningparameters(nchains = N, niterations = T, storeall = TRUE)
> print(pawlparameters)
Object of class tuningparameters.
*number of parallel chains: 10
*number of iterations: 2000
*compute mean: FALSE
*compute mean (burnin): 0
*save every nth iteration: 1
> #########
> # PAWL where we bin along the state space
> getPos <- function(points, logdensity) points
> ## we further specify some parameters, like the bins,
> ## the desired frequency in each bin
> ncuts <- 2
> positionbinning <- binning(position = getPos,
+ name = "position",
+ binrange = c(-20, -3),
+ ncuts = ncuts,
+ useLearningRate = TRUE,
+ autobinning = TRUE)
desiredfreq unspecified: set to 1 / nbins for each bin
fhthreshold unspecified: set to 0.5
splitThreshold unspecified: set to 0.1
minSimEffort unspecified: set 200
learningrate unspecified: set to x -> x^(-1)
useFH unspecified: set to TRUE
> print(positionbinning)
Object of class binning.
*name: position
*autobinning: TRUE
*bins: -Inf -20 -3
*Flat Histogram threshold: 0.5
*split threshold: 0.1
*use Flat Histogram criterion: TRUE
*use learning rate: TRUE
*min number of iterations between FH: 200
> pawlresults <- pawl(gaussiantarget, binning = positionbinning, AP = pawlparameters,
+ proposal = proposal)
[1] "Launching Particle Wang-Landau algorithm ..."
saving every 1 iterations
hence saving a vector of size 2001 x 10 x 1 = 20010
iteration 100 / 2000
/** diagnosis at iteration 100
* current number of bins = 3
desired freq - proportions:
1 0.05841584 -1.058416
* skewed bins: 2
* bins with enough points to be split: 2 3
* bins to split: 2
*/
iteration 200 / 2000
/** diagnosis at iteration 200
* current number of bins = 4
desired freq - proportions:
1 1 -2.516 -0.242
* skewed bins: 3
* bins with enough points to be split: 3 4
* bins to split: 3
*/
iteration 300 / 2000
/** diagnosis at iteration 300
* current number of bins = 5
desired freq - proportions:
1 1 -0.308 -2.48 -0.803
* skewed bins: 3
* bins with enough points to be split: 3 4 5
* bins to split: 3
*/
iteration 400 / 2000
/** diagnosis at iteration 400
* current number of bins = 6
desired freq - proportions:
1 -0.374 -1.616 -2.768 -0.428 -0.158
* skewed bins: 2 5
* bins with enough points to be split: 2 3 4 5 6
* bins to split: 2
*/
iteration 500 / 2000
/** diagnosis at iteration 500
* current number of bins = 7
desired freq - proportions:
1 1 -4.46 -0.128 0.4 -0.884 0.052
* skewed bins: 3 5 6
* bins with enough points to be split: 3 6 7
* bins to split: 3
*/
iteration 600 / 2000
/** diagnosis at iteration 600
* current number of bins = 8
desired freq - proportions:
1 -4.088 -0.728 -1.496 1 1 0.076 0.031
* skewed bins: 2 3 7
* bins with enough points to be split: 2 4 8
* bins to split: 2
*/
iteration 700 / 2000
/** diagnosis at iteration 700
* current number of bins = 9
desired freq - proportions:
0.007 -0.824 -2.84 -0.536 -0.32 0.832 0.928 0.64 0.178
* skewed bins: 2 4 5 7
* bins with enough points to be split: 1 3 9
*/
iteration 800 / 2000
/** diagnosis at iteration 800
* current number of bins = 9
desired freq - proportions:
-0.371 -0.188 -1.172 -0.008 -0.296 0.916 0.904 0.556 0.2125
* skewed bins: 2 4 5
* bins with enough points to be split: 1 3 5 9
* bins to split: 5
*/
iteration 900 / 2000
/** diagnosis at iteration 900
* current number of bins = 10
desired freq - proportions:
-1.1 1 -0.104 0.28 0.952 0.088 0.952 0.952 1 0.4
* skewed bins: 3 5 6 8
* bins with enough points to be split: 1 10
*/
iteration 1000 / 2000
/** diagnosis at iteration 1000
* current number of bins = 10
desired freq - proportions:
-1.232 0.976 0.352 0.616 0.976 0.544 0.976 0.976 1 0.4
* skewed bins: 2 5 6 8
* bins with enough points to be split: 1 10
*/
iteration 1100 / 2000
/** diagnosis at iteration 1100
* current number of bins = 10
desired freq - proportions:
-1.288 0.984 0.568 0.744 0.984 0.696 0.984 0.984 1 0.4
* skewed bins: 2 5 6 8
* bins with enough points to be split: 1 10
*/
iteration 1200 / 2000
/** diagnosis at iteration 1200
* current number of bins = 10
desired freq - proportions:
-1.21175 0.256 0.61 0.772 0.988 0.772 0.988 0.988 1 0.4
* skewed bins: 5 6 8
* bins with enough points to be split: 1 2 10
*/
iteration 1300 / 2000
/** diagnosis at iteration 1300
* current number of bins = 10
desired freq - proportions:
-1.0136 -0.1136 0.1312 0.2224 0.7024 0.6736 0.9904 0.9904 1 0.4
* skewed bins: 3 6 8
* bins with enough points to be split: 1 2 3 4 10
*/
iteration 1400 / 2000
/** diagnosis at iteration 1400
* current number of bins = 10
desired freq - proportions:
-0.844 -0.188 -0.04 -0.012 0.224 0.128 0.896 0.988 0.936 0.3825
* skewed bins: 8
* bins with enough points to be split: 1 2 3 4 6 10
*/
iteration 1500 / 2000
/** diagnosis at iteration 1500
* current number of bins = 10
desired freq - proportions:
-0.7318571 -0.2171429 -0.05942857 0.07428571 0.28 0.1977143 0.8354286 0.9691429 0.9331429 0.2684286
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 10
*/
iteration 1600 / 2000
/** diagnosis at iteration 1600
* current number of bins = 10
desired freq - proportions:
-0.5975 -0.089 0.061 0.082 0.256 0.238 0.793 0.889 0.73 0.167125
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 9 10
*/
iteration 1700 / 2000
/** diagnosis at iteration 1700
* current number of bins = 10
desired freq - proportions:
-0.4866667 0.032 0.1653333 0.184 0.3386667 0.3226667 0.736 0.7493333 0.528 0.08
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 9 10
*/
iteration 1800 / 2000
/** diagnosis at iteration 1800
* current number of bins = 10
desired freq - proportions:
-0.398 0.124 0.1696 0.1696 0.3424 0.1792 0.5656 0.5944 0.4588 0.0478
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
Flat histogram criterion met at iteration 1850 !
iteration 1900 / 2000
iteration 2000 / 2000
> # histogram of the binned coordinate
> PlotHistBin(pawlresults, positionbinning)
[1] "converting chains before plotting..."
[1] "...done"
> #########
> ## PAWL where we bin along the energy (= - log density)
> getPos <- function(points, logdensity) - logdensity
> # we can get the range using pre exploratory mcmc
> preexpresults <- preexplorationAMH(gaussiantarget, N, 1000)
Pre exploration to get log energy range
[1] "Launching Adaptive Metropolis-Hastings algorithm with parameters:"
Object of class tuningparameters.
*number of parallel chains: 10
*number of iterations: 1000
*compute mean: FALSE
*compute mean (burnin): 0
*save every nth iteration: -1
adaptationrate unspecified: set to x -> x^(-0.6)
sigma_init unspecified: set 1
Iteration 100 / 1000
Iteration 200 / 1000
Iteration 300 / 1000
Iteration 400 / 1000
Iteration 500 / 1000
Iteration 600 / 1000
Iteration 700 / 1000
Iteration 800 / 1000
Iteration 900 / 1000
Iteration 1000 / 1000
> binrange <- preexpresults$SuggestedRange
> # or specify it ourselves
> #binrange <- c(2.1, 20)
> # we further specify some parameters, like the bins,
> # the desired frequency in each bin...
> ncuts <- 2
> energybinning <- binning(position = getPos,
+ name = "energy",
+ binrange = binrange,
+ ncuts = ncuts,
+ useLearningRate = TRUE,
+ autobinning = TRUE,
+ alongenergy = TRUE)
desiredfreq unspecified: set to 1 / nbins for each bin
fhthreshold unspecified: set to 0.5
splitThreshold unspecified: set to 0.1
minSimEffort unspecified: set 200
learningrate unspecified: set to x -> x^(-1)
useFH unspecified: set to TRUE
> print(energybinning)
Object of class binning.
*name: energy
*autobinning: TRUE
*bins: -Inf 2.027298 4.683023
*Flat Histogram threshold: 0.5
*split threshold: 0.1
*use Flat Histogram criterion: TRUE
*use learning rate: TRUE
*min number of iterations between FH: 200
> # launching the algorithm...
> pawlresults <- pawl(gaussiantarget, binning = energybinning, AP = pawlparameters,
+ proposal = proposal)
[1] "Launching Particle Wang-Landau algorithm ..."
saving every 1 iterations
hence saving a vector of size 2001 x 10 x 1 = 20010
right end bin reached: disactivating bin diagnosis
iteration 100 / 2000
iteration 200 / 2000
Flat histogram criterion met at iteration 200 !
iteration 300 / 2000
iteration 400 / 2000
Flat histogram criterion met at iteration 400 !
iteration 500 / 2000
iteration 600 / 2000
Flat histogram criterion met at iteration 600 !
iteration 700 / 2000
iteration 800 / 2000
Flat histogram criterion met at iteration 800 !
iteration 900 / 2000
iteration 1000 / 2000
Flat histogram criterion met at iteration 1000 !
iteration 1100 / 2000
iteration 1200 / 2000
Flat histogram criterion met at iteration 1200 !
iteration 1300 / 2000
iteration 1400 / 2000
Flat histogram criterion met at iteration 1400 !
iteration 1500 / 2000
iteration 1600 / 2000
Flat histogram criterion met at iteration 1600 !
iteration 1700 / 2000
iteration 1800 / 2000
Flat histogram criterion met at iteration 1800 !
iteration 1900 / 2000
iteration 2000 / 2000
Flat histogram criterion met at iteration 2000 !
> # histogram of the binned coordinate
> #PlotHistBin(pawlresults, energybinning)
> ##
> #########
>
> #########
> ## results
> #
> getFrequencies(pawlresults, energybinning)
/* Frequency check
*Do the obtained frequencies match the desired frequencies?
*final bins: -Inf 2.027298 4.683023
*corresponding frequencies: 0.3317841 0.332084 0.3361319
*initial bins: -Inf 2.027298 4.683023
*desired frequencies: 0.3333333 0.3333333 0.3333333
*obtained frequencies: 0.3317841 0.332084 0.3361319
[1] 0.3317841 0.3320840 0.3361319
> # Histogram of the chains
> #PlotHist(pawlresults, 1)
> # Trace plot of the log theta penalties
> #print(PlotLogTheta(pawlresults))
> # trace plot of all the variables
> # (here there is only one variable)
> #print(PlotAllVar(pawlresults))
>
>
>
>
>
>
>
demo(mixture2kexample)
---- ~~~~~~~~~~~~~~~~
> #rm(list = ls())
> try(detach(package:PAWL, unload = TRUE), silent = TRUE)
> library(PAWL)
> set.seed(17)
> mixture <- createMixtureTarget()
ncomponents unspecified: setting it to 2
mixtureparameters unspecified: putting some default parameters
$ncomponents
[1] 2
$componentweights
[1] 0.5 0.5
$componentmeans
[1] 0.0 2.5
$componentvariances
[1] 0.3 0.3
no sample provided and no sample size: setting size to 100
generating sample of size 100
> N <- 10
> T <- 5000
> betaindex <- mixture@dimension
> getBeta <- function(points, logdensity) exp(points[,betaindex])
> betabinning <- binning(position = getBeta,
+ name = "beta",
+ binrange = c(2, 16),
+ autobinning = TRUE,
+ fhthreshold = 0.5,
+ useLearningRate = TRUE)
number of cuts (ncuts) unspecified: set to 9
desiredfreq unspecified: set to 1 / nbins for each bin
splitThreshold unspecified: set to 0.1
minSimEffort unspecified: set 200
learningrate unspecified: set to x -> x^(-1)
useFH unspecified: set to TRUE
> print(betabinning)
Object of class binning.
*name: beta
*autobinning: TRUE
*bins: -Inf 2 3.75 5.5 7.25 9 10.75 12.5 14.25 16
*Flat Histogram threshold: 0.5
*split threshold: 0.1
*use Flat Histogram criterion: TRUE
*use learning rate: TRUE
*min number of iterations between FH: 200
> pawlparameters <- tuningparameters(nchains = N, niterations = T, storeall = TRUE)
> print(pawlparameters)
Object of class tuningparameters.
*number of parallel chains: 10
*number of iterations: 5000
*compute mean: FALSE
*compute mean (burnin): 0
*save every nth iteration: 1
> ### Launching the algorithm...
> pawlresults <- pawl(mixture, binning = betabinning, AP = pawlparameters)
[1] "Launching Particle Wang-Landau algorithm ..."
adaptationrate unspecified: set to x -> x^(-0.6)
sigma_init unspecified: set 1
saving every 1 iterations
hence saving a vector of size 5001 x 10 x 7 = 350070
iteration 100 / 5000
/** diagnosis at iteration 100
* current number of bins = 10
desired freq - proportions:
-3.970297 -0.3663366 -0.07920792 0.0990099 0.2871287 0.4950495 0.8514851 0.8613861 0.9207921 0.9009901
* skewed bins: 7 8
* bins with enough points to be split: 1 2 3
*/
iteration 200 / 5000
/** diagnosis at iteration 200
* current number of bins = 10
desired freq - proportions:
-3.487562 0.119403 0.1840796 0.2437811 0.2636816 0.3930348 0.4825871 0.5572139 0.6368159 0.6069652
* skewed bins: 9
* bins with enough points to be split: 1 2 3 4 5 6 7
*/
iteration 300 / 5000
/** diagnosis at iteration 300
* current number of bins = 10
desired freq - proportions:
-3.325581 -0.07641196 0.07973422 0.2724252 0.3488372 0.4252492 0.4784053 0.551495 0.6212625 0.6245847
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 400 / 5000
/** diagnosis at iteration 400
* current number of bins = 10
desired freq - proportions:
-3.244389 0.1620948 0.1795511 0.2668329 0.3017456 0.4114713 0.4164589 0.4613466 0.4962594 0.5486284
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 500 / 5000
/** diagnosis at iteration 500
* current number of bins = 10
desired freq - proportions:
-3.195609 0.1297405 0.253493 0.3313373 0.3772455 0.3752495 0.3932136 0.4291417 0.4491018 0.4570858
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 600 / 5000
/** diagnosis at iteration 600
* current number of bins = 10
desired freq - proportions:
-3.036606 -0.01830283 0.3111481 0.3327787 0.3527454 0.3610649 0.3793677 0.422629 0.437604 0.4575707
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 700 / 5000
/** diagnosis at iteration 700
* current number of bins = 10
desired freq - proportions:
-2.773181 -0.1041369 0.2097004 0.3266762 0.3238231 0.3708987 0.382311 0.403709 0.4265335 0.4336662
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 800 / 5000
/** diagnosis at iteration 800
* current number of bins = 10
desired freq - proportions:
-2.55181 -0.09113608 0.1810237 0.2509363 0.3245943 0.3245943 0.3645443 0.3895131 0.3832709 0.4244694
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 900 / 5000
/** diagnosis at iteration 900
* current number of bins = 10
desired freq - proportions:
-2.379578 -0.08102109 0.1498335 0.2519423 0.3085461 0.3229745 0.3263041 0.3529412 0.3695893 0.3784684
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1000 / 5000
/** diagnosis at iteration 1000
* current number of bins = 10
desired freq - proportions:
-2.217782 -0.07492507 0.1628372 0.2267732 0.2887113 0.2917083 0.3056943 0.3206793 0.3356643 0.3606394
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1100 / 5000
/** diagnosis at iteration 1100
* current number of bins = 10
desired freq - proportions:
-2.016349 0.02270663 0.1407811 0.1762035 0.2179837 0.2579473 0.2752044 0.2960945 0.3042688 0.3251589
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1200 / 5000
/** diagnosis at iteration 1200
* current number of bins = 10
desired freq - proportions:
-1.84846 0.08909242 0.1298918 0.1681932 0.1948376 0.2181515 0.2414654 0.253955 0.2706078 0.2822648
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1300 / 5000
/** diagnosis at iteration 1300
* current number of bins = 10
desired freq - proportions:
-1.70638 0.1145273 0.1283628 0.1583397 0.187548 0.1975404 0.2113759 0.2244427 0.2398155 0.2444274
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1400 / 5000
/** diagnosis at iteration 1400
* current number of bins = 10
desired freq - proportions:
-1.584582 0.06995004 0.1349036 0.1577445 0.1705924 0.1855817 0.1948608 0.2098501 0.2212705 0.2398287
* skewed bins:
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1500 / 5000
/** diagnosis at iteration 1500
* current number of bins = 10
desired freq - proportions:
-1.479014 0.00799467 0.07661559 0.1672219 0.1832112 0.188541 0.2018654 0.2111925 0.2225183 0.2198534
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1600 / 5000
/** diagnosis at iteration 1600
* current number of bins = 10
desired freq - proportions:
-1.386633 0.04434728 0.06183635 0.105559 0.16802 0.1773891 0.1892567 0.1986259 0.2117427 0.2298563
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1700 / 5000
/** diagnosis at iteration 1700
* current number of bins = 10
desired freq - proportions:
-1.305115 -0.01410935 0.09641387 0.1387419 0.1475603 0.1604938 0.1769547 0.191652 0.1981188 0.2092887
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1800 / 5000
/** diagnosis at iteration 1800
* current number of bins = 10
desired freq - proportions:
-1.232649 0.0294281 0.05830094 0.1071627 0.1338145 0.141588 0.1776791 0.1860078 0.197668 0.2009994
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 1900 / 5000
/** diagnosis at iteration 1900
* current number of bins = 10
desired freq - proportions:
-1.167806 0.04523935 0.07627564 0.09942136 0.135718 0.1430826 0.1520252 0.162546 0.1746449 0.1788532
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2000 / 5000
/** diagnosis at iteration 2000
* current number of bins = 10
desired freq - proportions:
-1.109445 -0.006996502 0.07446277 0.109945 0.1344328 0.1529235 0.1554223 0.1589205 0.1629185 0.1674163
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2100 / 5000
/** diagnosis at iteration 2100
* current number of bins = 10
desired freq - proportions:
-1.05664 -0.05425988 0.07615421 0.104712 0.1199429 0.1413613 0.1608758 0.1651594 0.1694431 0.1732508
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2200 / 5000
/** diagnosis at iteration 2200
* current number of bins = 10
desired freq - proportions:
-1.008632 -0.0945025 0.07905498 0.1054066 0.1281236 0.1385734 0.1535666 0.1612903 0.1671967 0.1699228
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2300 / 5000
/** diagnosis at iteration 2300
* current number of bins = 10
desired freq - proportions:
-0.9647979 -0.09039548 0.0821382 0.09821817 0.1142981 0.1403738 0.1460235 0.1499348 0.1594959 0.164711
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2400 / 5000
/** diagnosis at iteration 2400
* current number of bins = 10
desired freq - proportions:
-0.9246147 -0.08663057 0.07455227 0.09829238 0.1249479 0.1316118 0.1336943 0.1474386 0.1474386 0.1532695
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2500 / 5000
/** diagnosis at iteration 2500
* current number of bins = 10
desired freq - proportions:
-0.8876449 -0.06477409 0.06517393 0.0967613 0.1107557 0.1167533 0.127549 0.1363455 0.1491403 0.14994
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2600 / 5000
/** diagnosis at iteration 2600
* current number of bins = 10
desired freq - proportions:
-0.8535179 -0.02383699 0.05690119 0.09034987 0.1022684 0.1149558 0.1199539 0.1222607 0.1318724 0.1387928
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2700 / 5000
/** diagnosis at iteration 2700
* current number of bins = 10
desired freq - proportions:
-0.8219178 0.01406886 0.06553128 0.08811551 0.09440948 0.1025546 0.1058867 0.1121807 0.1181044 0.1210663
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2800 / 5000
/** diagnosis at iteration 2800
* current number of bins = 10
desired freq - proportions:
-0.7925741 0.02677615 0.05176723 0.08782578 0.09210996 0.09603713 0.1024634 0.1074616 0.1128169 0.115316
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 2900 / 5000
/** diagnosis at iteration 2900
* current number of bins = 10
desired freq - proportions:
-0.7652534 0.01378835 0.06204757 0.07583592 0.08376422 0.09410548 0.0989314 0.104102 0.1154774 0.117201
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3000 / 5000
/** diagnosis at iteration 3000
* current number of bins = 10
desired freq - proportions:
-0.7397534 0.01332889 0.05564812 0.06864379 0.0853049 0.09596801 0.1002999 0.1016328 0.1066311 0.1122959
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3100 / 5000
/** diagnosis at iteration 3100
* current number of bins = 10
desired freq - proportions:
-0.7158981 0.01289906 0.03611738 0.05449855 0.08577878 0.09126088 0.09996775 0.1057723 0.1106095 0.1189939
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3200 / 5000
/** diagnosis at iteration 3200
* current number of bins = 10
desired freq - proportions:
-0.6935333 0.02967823 0.03561387 0.0659169 0.07216495 0.08059981 0.09122149 0.100906 0.1052796 0.1121525
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3300 / 5000
/** diagnosis at iteration 3300
* current number of bins = 10
desired freq - proportions:
-0.6725235 0.03847319 0.05574068 0.06301121 0.0724023 0.07967283 0.08330809 0.08664041 0.09391094 0.09936383
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3400 / 5000
/** diagnosis at iteration 3400
* current number of bins = 10
desired freq - proportions:
-0.6527492 0.03469568 0.05116142 0.05880623 0.07115554 0.07703617 0.08056454 0.08820935 0.09085563 0.1002646
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3500 / 5000
/** diagnosis at iteration 3500
* current number of bins = 10
desired freq - proportions:
-0.6341045 -0.01770923 0.04970009 0.06798058 0.07597829 0.08169095 0.08626107 0.08968866 0.09625821 0.1042559
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3600 / 5000
/** diagnosis at iteration 3600
* current number of bins = 10
desired freq - proportions:
-0.6164954 -0.0449875 0.0441544 0.06137184 0.07470147 0.08303249 0.09136351 0.09913913 0.1027492 0.1049708
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3700 / 5000
/** diagnosis at iteration 3700
* current number of bins = 10
desired freq - proportions:
-0.5998379 -0.07079168 0.05674142 0.06727911 0.07673602 0.08430154 0.0878141 0.09159687 0.1010538 0.1051067
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3800 / 5000
/** diagnosis at iteration 3800
* current number of bins = 10
desired freq - proportions:
-0.5840568 -0.0952381 0.05945804 0.06840305 0.07945278 0.08392528 0.09181794 0.09760589 0.0973428 0.1012891
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 3900 / 5000
/** diagnosis at iteration 3900
* current number of bins = 10
desired freq - proportions:
-0.5690849 -0.1133043 0.05537042 0.07331453 0.07792874 0.08561907 0.09074596 0.09510382 0.1004871 0.1038195
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 4000 / 5000
/** diagnosis at iteration 4000
* current number of bins = 10
desired freq - proportions:
-0.5548613 -0.08672832 0.0419895 0.05898525 0.07223194 0.08097976 0.08922769 0.09522619 0.1009748 0.1019745
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 4100 / 5000
/** diagnosis at iteration 4100
* current number of bins = 10
desired freq - proportions:
-0.5413314 -0.06022921 0.04437942 0.05754694 0.06827603 0.07461595 0.08217508 0.08827115 0.0919288 0.09436723
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 4200 / 5000
/** diagnosis at iteration 4200
* current number of bins = 10
desired freq - proportions:
-0.5284456 -0.03499167 0.04284694 0.05379672 0.06736491 0.07545822 0.07760057 0.080219 0.08307546 0.08307546
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 4300 / 5000
/** diagnosis at iteration 4300
* current number of bins = 10
desired freq - proportions:
-0.516159 -0.01092769 0.02999302 0.057196 0.06231109 0.06626366 0.07207626 0.07649384 0.0802139 0.08253894
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
iteration 4400 / 5000
/** diagnosis at iteration 4400
* current number of bins = 10
desired freq - proportions:
-0.5044308 0.01204272 0.03181095 0.05157919 0.0624858 0.06157691 0.066803 0.07089298 0.07498296 0.07225631
* skewed bins: 2
* bins with enough points to be split: 1 2 3 4 5 6 7 8 9 10
*/
Flat histogram criterion met at iteration 4439 !
iteration 4500 / 5000
iteration 4600 / 5000
Flat histogram criterion met at iteration 4639 !
iteration 4700 / 5000
iteration 4800 / 5000
Flat histogram criterion met at iteration 4873 !
iteration 4900 / 5000
iteration 5000 / 5000
> chains <- ConvertResults(pawlresults)
[1] "converting chains before plotting..."
[1] "...done"
> getFrequencies(pawlresults, betabinning)
/* Frequency check
*Do the obtained frequencies match the desired frequencies?
*final bins: -Inf 2 3.75 5.5 7.25 9 10.75 12.5 14.25 16
*corresponding frequencies: 0.1443911 0.0989802 0.09906019 0.09636073 0.0939812 0.0940012 0.09374125 0.09336133 0.09312138 0.0930014
*initial bins: -Inf 2 3.75 5.5 7.25 9 10.75 12.5 14.25 16
*desired frequencies: 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
*obtained frequencies: 0.1443911 0.0989802 0.09906019 0.09636073 0.0939812 0.0940012 0.09374125 0.09336133 0.09312138 0.0930014
[1] 0.14439112 0.09898020 0.09906019 0.09636073 0.09398120 0.09400120
[7] 0.09374125 0.09336133 0.09312138 0.09300140
> ## Plot the reaction coordinate values along the binning axis,
> ## with red vertical lines denoting the endpoints of the bins
> # print(PlotHistBin(pawlresults, betabinning))
> ## Plot the log penalties associated with each bin
> # print(PlotLogTheta(pawlresults))
> ## Plot sigma versus iterations, where sigma is
> ## the vector of the standard deviations used by the MH kernel along the iterations.
> #plot(pawlresults$sigma, type = "l")
> #PlotAllVar(chains)
> #X11()
> print(PlotComp1vsComp2(chains, "X3", "X4"))
> #X11()
> #print(PlotHistBin(pawlresults, betabinning))
> #print(PlotHist(pawlresults, component = 7))
> #X11(); print(PlotHist(pawlresults, component = 3))
>
> #PlotFH(pawlresults)
> #binincrease <- data.frame(cbind(c(1, pawlresults$splitTimes, T),
> # c(pawlresults$nbins, max(pawlresults$nbins))))
> #g <- ggplot(binincrease, aes(x = X1, y = X2)) + geom_step()
> #g <- g + ylim(0, 1.5 * max(pawlresults$nbins)) + ylab("Number of bins") + xlab("iterations")
> #g <- g + theme(title = "Number of bins along the iterations")
> #print(g)
>
>
> ### Adaptive MH for comparison
> #mhparameters <- tuningparameters(nchains = N, niterations = T, storeall = TRUE)
> #### launching the algorithm...
> #amhresults <- adaptiveMH(mixture, mhparameters)
>
> #PlotAllVar(amhresults)
> #X11()
> #print(PlotComp1vsComp2(amhresults, "X3", "X4"))
> #X11(); print(PlotHist(amhresults, component = 3))
> #
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.