diff --git a/_quarto.yml b/_quarto.yml index 4cef2c0..7631883 100644 --- a/_quarto.yml +++ b/_quarto.yml @@ -25,6 +25,8 @@ book: - montyDSL.qmd - composability.qmd - samplers.qmd + - runners.qmd + - advanced_monty.qmd - part: "Odin & Monty" chapters: - inference.qmd diff --git a/advanced_monty.qmd b/advanced_monty.qmd new file mode 100644 index 0000000..5dbb364 --- /dev/null +++ b/advanced_monty.qmd @@ -0,0 +1,9 @@ +# Advanced monty + +## Observers + +## Working with other packages + +Exporting chains to work with other packages + +## Nested models \ No newline at end of file diff --git a/composability.qmd b/composability.qmd index 5b4977f..196c563 100644 --- a/composability.qmd +++ b/composability.qmd @@ -1 +1,24 @@ -# Composability \ No newline at end of file +# Composability + +```{r} +library(monty) +``` + +```{r} +# A simple example; a model that contains something of interest, +# and a simple prior from monty_dsl +likelihood <- monty_example("banana") +prior <- monty_dsl({ + alpha ~ Normal(0, 1) + beta ~ Normal(0, 10) +}) +posterior <- likelihood + prior +posterior + +# # The same thing, more explicitly: +# monty_model_combine(likelihood, prior) +# +# # Control properties of the combined model: +# monty_model_combine(likelihood, prior, +# monty_model_properties(has_gradient = FALSE)) +``` \ No newline at end of file diff --git a/model.qmd b/model.qmd index 38795a1..6f2c160 100644 --- a/model.qmd +++ b/model.qmd @@ -4,4 +4,42 @@ We often refers to using models. ## SimCity vs Statistical models -## Normalised vs unormalised densities \ No newline at end of file +## Normalised vs unormalised densities + +In Bayesian inference, we often encounter probability densities that can either be normalised or unnormalised, and understanding the distinction between these is crucial for implementing Monte Carlo methods effectively. Normalisation is particularly important when interpreting probabilities and integrating over parameter spaces, a concept with parallels in physics, where distributions are sometimes normalised to ensure they represent meaningful physical quantities. + +### Unnormalised Density + +An **unnormalised density** represents a probability density function that is proportional to the target distribution but has not been scaled to integrate to 1 over the entire parameter space. In Bayesian inference, we are often interested in the posterior distribution of a parameter $\theta$ given data $y$. The posterior density is given by: + +$$p(\theta | y) = \frac{p(y|\theta) p(\theta)}{p(y)}$$ + +where: + +- $p(y | \theta))$ is the likelihood of the data given the parameter, +- $p(\theta)$ is the prior distribution, and +- $p(y) = \int p(y | \theta) p(\theta) \, d\theta \)$ is the marginal likelihood, or the evidence, which acts as a normalising constant. + +In many cases, computing $p(y)$ is challenging or even intractable, so we work with the unnormalised posterior density \( p(y \mid \theta) p(\theta) \). The unnormalised density is sufficient for methods like **Markov Chain Monte Carlo (MCMC)**, where the absolute probability values are less relevant than their relative proportions across different values of \( \theta \). + +### Normalised Density + +A **normalised density** ensures that the total probability over all possible parameter values is 1. Normalisation is achieved by dividing by the integral over the entire density, as shown above with \( p(y) \) in the denominator. This step is essential when absolute probability values are required or when we need to make direct probabilistic interpretations. + +In Bayesian inference, although we often rely on unnormalised densities, there are cases—such as when comparing models using Bayes factors—where we require the fully normalised posterior. + +### Link to Physics: Partition Function and Probability Density + +The concept of normalisation has parallels in statistical mechanics. In physics, the **partition function** \( Z \) normalises the probability density over all microstates \( x \) of a system, enabling calculations of thermodynamic properties. Here, the probability density for a state \( x \) with energy \( E(x) \) at temperature \( T \) is proportional to \( e^{-E(x)/kT} \), where \( k \) is Boltzmann's constant. The partition function \( Z = \sum_{x} e^{-E(x)/kT} \) or its integral form over a continuous state space normalises this density, ensuring probabilities are meaningful when summing over all possible states. + +In Monte Carlo algorithms, this parallel becomes especially useful. For example, in the **Metropolis-Hastings algorithm**, we only need the ratio of probabilities between two states. This mirrors physical systems where absolute energies are less important than energy differences between states, similar to unnormalised densities in Bayesian inference. + +### Monte Carlo Methods with Unnormalised Densities + +Monte Carlo methods like **Metropolis-Hastings** and **Importance Sampling** often operate with unnormalised densities, since only the relative probabilities affect the algorithm's behaviour. This allows us to bypass the need to calculate the normalising constant, which can be computationally expensive or intractable in high-dimensional spaces. + +In Importance Sampling, for example, we estimate expectations with respect to the target distribution by drawing samples from a simpler proposal distribution and weighting them according to the unnormalised target density. + + + +## Properties of monty models \ No newline at end of file diff --git a/monty.qmd b/monty.qmd index fdf6fc1..218d9e2 100644 --- a/monty.qmd +++ b/monty.qmd @@ -1,16 +1,13 @@ # Getting started with monty -* the philosophie -* three options - -`monty` is a R package design to work with statistical models in a *modular* way. It uses Monte Carlo methods to facilitate drawing from these statistical models and is currently built to create complex models through four different routes: +`monty` is a R package designed to work with statistical models in a modular way. It uses Monte Carlo methods to sample from these statistical models. It is currently built to create complex models through four different routes: * simple R code * a dedicated DSL * odin models * composing two monty models. -The simple example in this chapter is built from a simple R function and is the most basic way of defining a monty model. We will cover more complex examples in the model section. Section ** deals with the monty DSL, while section ** explains how to obtain new monty models through composition of two monty models. monty models obtains from odin models are dealt in the odin & monty chapter. +The simple example in this chapter comes from a simple R function and is the most basic way of defining a monty model. We will cover more complex examples in the model section. Section ** deals with the monty DSL, while section ** explains how to obtain new monty models through the composition of two monty models. monty models obtained from odin models are dealt in the odin & monty chapter. ## A simple example @@ -24,16 +21,43 @@ For this we start by defining ```{r} fn <- function(l, p, m1, m2) { - p*dnorm(l,mean = m1) + (1-p)*dnorm(l, mean = m2) + log(p*dnorm(l,mean = m1) + (1-p)*dnorm(l, mean = m2)) } ``` ```{r} -l_distribution <- monty_model_function(fn, fixed = list(p = 0.75, m1 = 5, m2=7)) +l_distribution <- monty_model_function(fn, fixed = list(p = 0.75, m1 = 3, m2=7)) +``` + +We have just created a `monty` model. + +```{r} +l <- seq(from = 0, to = 10, by = 0.1) +plot(l, + exp(Vectorize(l_distribution$density)(l)), + type = "l", + ylab = "density") ``` + ## Sampling from our example distribution + +```{r} +sampler <- monty_sampler_random_walk(matrix(1)) +samples <- monty_sample(l_distribution, sampler, 2000, initial = 3, + n_chains = 4) +``` + +```{r} +matplot(samples$density, type = "l", lty = 1, + xlab = "log posterior density", ylab = "sample", col = "#00000055") +``` + +```{r} +hist(samples$pars["l",,], breaks=100) +``` + ## Further work Before showing how to fit odin models to data, we'll focus for a bit on monty itself. If you are anxious to get on and fit the model from @sec-data-filter, you might skip ahead to @sec-inference, where we resume this. diff --git a/runners.qmd b/runners.qmd new file mode 100644 index 0000000..d25e818 --- /dev/null +++ b/runners.qmd @@ -0,0 +1,3 @@ +# Runners + +One of the focus of `monty` is tailoring implementation to optimise usage of available hardware while ensuring reproducible results accross platforms. \ No newline at end of file diff --git a/samplers.qmd b/samplers.qmd index 71c6c4f..318fe05 100644 --- a/samplers.qmd +++ b/samplers.qmd @@ -1,4 +1,4 @@ -# Samplers +# Samplers and inference This vignette will describe the samplers available in monty.