diff --git a/README.md b/README.md index 6090215..80aac1a 100644 --- a/README.md +++ b/README.md @@ -1,13 +1,13 @@ # bayes.js - MCMC and Bayes in the browser -bayes.js A small toy javascript MCMC framework that can be used fit Bayesian models in the browser. I call it "toy" because I would use it for fun, but not in production... +bayes.js is small toy javascript MCMC framework that can be used fit Bayesian models in the browser. I call it "toy" because I would use it for fun, but not in production... The two major files are: -* __mcmc.js__ Implements a MCMC framework which can be used to fit Bayesian model with both discrete and continuous parameters. Currently the only algorithm that is implemented is a version of the *adaptive Metropolis within Gibbs* (AMWG) algorithm presented by [Roberts and Rosenthal (2009) ](http://probability.ca/jeff/ftpdir/adaptex.pdf). Loading this file in the browser creates the global object `mcmc`. -* __distributions.js__ A collection of log density functions that can be used to construct Bayesian models. Follows the naming scheme `ld.*` (for example, `ld.norm` and `ld.pois`) and uses the same parameters as the `d*` density functions in R. Loading this file in the browser creates the global object `ld`. +* __mcmc.js__ - Implements a MCMC framework which can be used to fit Bayesian model with both discrete and continuous parameters. Currently the only algorithm that is implemented is a version of the *adaptive Metropolis within Gibbs* (AMWG) algorithm presented by [Roberts and Rosenthal (2009) ](http://probability.ca/jeff/ftpdir/adaptex.pdf). Loading this file in the browser creates the global object `mcmc`. +* __distributions.js__ - A collection of log density functions that can be used to construct Bayesian models. Follows the naming scheme `ld.*` (for example, `ld.norm` and `ld.pois`) and uses the same parameters as the `d*` density functions in R. Loading this file in the browser creates the global object `ld`. -In addition to this the whole thing is wrapped within an [Rstudio](https://www.rstudio.com/) project as I've use R and JAGS to write some tests. +In addition to this the whole thing is wrapped in an [Rstudio](https://www.rstudio.com/) project as I've use R and JAGS to write some tests. Minimal example of fitting a Normal distribution -------------------- @@ -23,7 +23,7 @@ var params = { mu: {type: "real"}, sigma: {type: "real", lower: 0} }; -var norm_post = function(state, data) { +var log_post = function(state, data) { var log_post = 0; // Priors log_post += ld.norm(state.mu, 0, 100); @@ -36,7 +36,7 @@ var norm_post = function(state, data) { }; // Initializing the sampler -var sampler = new mcmc.AmwgSampler(params, norm_post, data); +var sampler = new mcmc.AmwgSampler(params, log_post, data); // Burning some samples to the MCMC gods and sampling 5000 draws. sampler.burn(1000) var samples = sampler.sample(5000) @@ -45,15 +45,15 @@ var samples = sampler.sample(5000) And here is a plot of the resulting sample made in javascript using the [plotly.js](https://plot.ly/javascript/) library: -![Normal model posterior](media/normal_model_plotly.png?raw=true) +[![Normal model posterior](media/normal_model_plotly.png?raw=true)](http://codepen.io/rasmusab/pen/LpaKep?editors=001) How to implement a custom Bayesian model ----------------------------------------- -In __mcmc.js__ the currently sole sampler is `AmwgSampler`: An adaptive, general purpose MCMC sample that can be used to fit a wide range of models *in theory*. In practice `AmwgSampler` can work well as long as there are not too many parameters. Here it's hard to specify what is "too many" but generally <10 parameters should be fine but >100 might be too many, but it depends. As AmwgSampler is a Gibbs sampler it also becomes crippled when the parameters are correlated. +In __mcmc.js__ the currently sole sampler is `AmwgSampler`: An adaptive, general purpose MCMC sampler that can be used to fit a wide range of models *in theory*. In practice `AmwgSampler` can work well as long as there are not too many parameters. Here it's hard to specify what is "too many" but generally <10 parameters should be fine but >100 might be too many, but it depends. As `AmwgSampler` is a Gibbs sampler it also becomes crippled when the parameters are correlated. -Here is how to create an instance of `AmwgSampler` and to produce a number of samples from a model defined by `params` and `log_post`: +Here is how to create an instance of `AmwgSampler` and produce a number of samples from a model defined by `params` and `log_post`: ```JavaScript var sampler = new mcmc.AmwgSampler(params, log_post, data); @@ -86,8 +86,8 @@ There are five parameter properties: * `type` - a string that defines the type of parameter. Can be `"real"`, `"int"` or `"binary"`. Note that this has nothing to do with javascript types, it's just to inform the sampler how to handle the parameter. * `dim` - an array giving the dimensions of the parameter. For example `dim: [1]` would define a scalar (a single number), `dim: [3]` a vector of length 3, and `dim: [2,2]` a 2 by 2 matrix. -* `lower` - A number giving the lower support for the parameter. -* `upper` - A number giving the upper support for the parameter. +* `lower` - A number giving the lower support of the parameter. +* `upper` - A number giving the upper support of the parameter. * `init` - A number or an array giving the initial value of the parameter when starting the sampler. Not all parameter properties need to be filled in and when left out will be replaced by defaults, for example, the following will result in the same parameter definition as above: @@ -96,7 +96,7 @@ Not all parameter properties need to be filled in and when left out will be repl params = {mu: {}, sigma: {lower:0}} ``` -If no `type` is given it defaults to `"real"` and Otherwise the default parameter properties depend on the `type`. The following... +If no `type` is given it defaults to `"real"` and otherwise the default parameter properties depend on the `type`. The following... ```JavaScript params = { @@ -141,7 +141,7 @@ state = { bin_mat: [[1, 1, 0], [0, 0, 1], [1, 1, 1]]; ``` -In order to calculate the log posterior it's useful to have access to functions that return the log density of a bunch of probability distributions. This can be found in the **distributions.js** file which when imported created global object `ld` that contains a number of log density functions, for example `ld.norm`, `ld.beta`, and `ld.pois`. +In order to calculate the log posterior it's useful to have access to functions that return the log density of a bunch of probability distributions. This can be found in the **distributions.js** file which when imported creates a global object `ld` that contains a number of log density functions, for example `ld.norm`, `ld.beta`, and `ld.pois`. The design pattern for crafting a log posterior function goes something like this (here exemplified by a standard beta Bernoulli model) : @@ -151,7 +151,7 @@ var log_post = function(state, data) { var log_post = 0; // Add the log densities of the priors to log_post. - // If Priors are not specified that is the same as assuming flat uniform priors. + // If Priors are not specified it is the same as assuming flat uniform priors. log_post += ld.beta(state.theta, 2, 2); // Loop through the data and add each resulting log density to log_post. @@ -173,18 +173,18 @@ In this section I've used the variable names `state` and `data` but you are free ### Defining the data (`data`) -This is easy, whatever you pass in as the `data` argument when creating AmwgSampler will be passed on to the log posterior function (`log_post`). So it's up to you what kind of data structure you want to work with in `log_post`. +This is easy, whatever you pass in as the `data` argument when instantiating an `AmwgSampler` will be passed on to the log posterior function (`log_post`). So it's up to you what kind of data structure you want to work with in `log_post`. ### Defining the options (`options`) There are a bunch of options that can be given to `AmwgSampler` as an options object, for example: -* `thin` - an integer specifying number of steps between each saved sample. Defaults to 1.0. -* `monitor` - an array of strings specifying what parameters to monitor and return sampled from. Defaults to `null` which means that all parameters will be monitored. +* `thin` - an integer specifying the number of steps between each saved sample. Defaults to 1.0. +* `monitor` - an array of strings specifying what parameters to monitor and return samples from. Defaults to `null` which means that all parameters will be monitored. -It is also possible to control how the sampler adapts its proposal, where the two most important options are: +It is also possible to control how the sampler adapts it's proposal, where the most important options are: -* `prop_log_scale` - The initial log(SD) of the one-dimensional proposal (using the same init for each parameter). Defaults to 0.0 which implies a SD of exp(0) == 1.0. +* `prop_log_scale` - The initial log(SD) of the one-dimensional proposal (using the same initial log(SD) for each parameter). Defaults to 0.0 which implies a SD of exp(0) == 1.0. * `target_accept_rate` - The acceptance rate of the proposals the sampler tries to achieve. Defaults to 0.44. * `batch_size` - The number of steps between each update of the `prop_log_scale`. Defaults to 50. * `max_adaptation` - The maximum amount `prop_log_scale` is changed each batch update. Defaults to 0.33. @@ -219,14 +219,14 @@ Here are some interactive examples implemented as Codepens that should run in th * [Simple linear regression](http://codepen.io/rasmusab/pen/wMzyGE?editors=001) * [Multivariate logistic regression](http://codepen.io/rasmusab/pen/eJdVBm?editors=001) * [the "Pump" example from the BUGS project](http://codepen.io/rasmusab/pen/jWMYxK?editors=001) -* [An hierarchical model with varying slope and intercept from Gelman and Hill ](http://codepen.io/rasmusab/pen/VewLoM?editors=001) (This is stretching the limits of what the simple `AmwgSampler` can do...) +* [An hierarchical model with varying slope and intercept from Gelman and Hill (2006) ](http://codepen.io/rasmusab/pen/VewLoM?editors=001) (This is stretching the limits of what the simple `AmwgSampler` can do...) -These demos rely on the [plotly library](https://plot.ly/javascript/) and I haven't tested them extensively on different platforms/browsers. You should be able to change the data and model definition on the fly (but if you change some stuff, like adding multidimensional variables, the plotting might stop working). +These demos rely on the [plotly.js library](https://plot.ly/javascript/) and I haven't tested them extensively on different platforms/browsers. You should be able to change the data and model definition on the fly (but if you change some stuff, like adding multidimensional variables, the plotting might stop working). Available distributions in distributions.js -------------------------------------------- -The file distributions.js implements log densities for most of the common probability distributions. The log density functions use the naming scheme of R when possible but are not vectorized in any way. Some o the implemented distributions are: +The file distributions.js implements log densities for most of the common probability distributions. The log density functions use the naming scheme of R when possible but are not vectorized in any way. Some of the implemented distributions are: * Bernoulli distribution: `ld.bern(x, prob)` * Binomial distribution: `ld.binom(x, size, prob)` @@ -238,12 +238,26 @@ The file distributions.js implements log densities for most of the common probab For the full list of distributions just check the source of **distributions.js**. +FAQ +------------------ + +* When is a javascript MCMC sampler useful? + - Well, for starters, it's *not* particularly useful if you want to do serious Bayesian data analysis. Then you should use a serious tool like [JAGS](http://mcmc-jags.sourceforge.net/) or [STAN](http://mc-stan.org/). It could, however, be useful if you would want to put a demo of a Bayesian model online, but don't want to / can't run the computations on a server. It could also be useful as a part of a javascript application making use of Bayesian computation at some point. +* How good is the sampler that bayes.js uses? + - bayes.js implements the *adaptive Metropolis within Gibbs* described by [Roberts and Rosenthal (2009) ](http://probability.ca/jeff/ftpdir/adaptex.pdf) which is a good algorithm in that (1) it's adaptive and works out-of-the-box without you having to set a lot of tuning parameters, (2) it can handle both continuous and discrete parameters, (3) it is easy to implement. The downside with the sampler is that (1) it only works well with a small number of parameters, (2) it's a Gibbs sampler so it's going to struggle with correlated parameters. +* What is "a small number of parameters"? + - Depends heavily on context but <10 is probably "small" here. But depending on the model you're trying to run you might get away with 100+ parameters. +* How fast is it? + - Also super context dependent. On simple models it's pretty fast, for example, fitting a standard Normal model on 1000 datapoints producing a sample of 20,000 draws takes ~0.5 s. on my computer. Also, when I've been playing around with different browsers I've seen order-of-magnitude changes in performance when changing seemingly arbitrary things. For example, inlining the definition of the Normal density in the function calculating the log posterior rather than using `ld.norm` defined in **distributions.js** resulted in 10x slower sampling on Firefox 37. +* Are there any alternatives if I want to do Bayes in the browser? + - There's a probibalistic programming language called [webppl](http://webppl.org/) that's implemented in javascript and that looks great, but I don't have any experience using it. + Some notes about the implementation of bayes.js ---------------------------------------- -The sampler as defined in **mcmc.js** is somewhat over-engineered as it includes a little "framework"" of classes for implementing new sampling algorithms. The two main type of classes are `Stepper`s and `Sampler`s. A `Stepper` is a class responsible for moving one or more parameters around in parameter space. A `Sampler` is a more "user facing" class that handles the setup of a collection of `Stepper`s that together takes a full step in the model's parameter space. Here `AmwgSampler` is so far the only implemented `Sampler`. For more details just check the source of **mcmcm.js** which is not completely absent of comments. +The sampler as defined in **mcmc.js** is somewhat over-engineered as it includes a little "framework" of classes for implementing new sampling algorithms. The two main type of classes are `Stepper`s and `Sampler`s. A `Stepper` is a class responsible for moving one or more parameters around in parameter space. A `Sampler` is a more "user facing" class that handles the setup of a collection of `Stepper`s that together takes a full step in the model's parameter space. Here `AmwgSampler` is so far the only implemented `Sampler`. For more details just check the source of **mcmcm.js** which is not completely absent of comments. -bayes.js is wrapped inside an [Rstudio](https://www.rstudio.com/) and that is because I've used R to compare and test the performance of **mcmc.js** and **distributions.js**. I've also implemented a number of unit tests, however many of those tests checks that the samples of bayes.js have the right distribution but because the samples are random draws these tests actually fail from time to time. Having tests that fail even if noting is wrong is not ideal, but it's better than no tests at all... +bayes.js is wrapped inside an [Rstudio](https://www.rstudio.com/) project and that's because I've used R to compare and test the performance of **mcmc.js** and **distributions.js**. I've also implemented a number of unit tests, however, many of those tests check that the samples of bayes.js have the right distribution but because the samples are random draws these tests actually fail from time to time. Having tests that fail even if noting is wrong is not ideal, but it's better than no tests at all I figured... References --------------------