# PAWL-package: PARALLEL ADAPTIVE WANG-LANDAU In PAWL: Implementation of the PAWL algorithm

## 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``` ``` demo(discreteexample) demo(gaussianexample) demo(mixture2kexample) ```

### Example output

```Loading required package: mvtnorm

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

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

> #######
> #######
> #mhparameters <- tuningparameters(nchains = 10, niterations = 1000, storeall = TRUE)
> #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")
>
>
> ######
> ######
> N <- 10

> T <- 2000

> # here we disable the adaptive proposal to highlight the automatic binning
> # mechanism
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...