Today

Program

• Model selection, continued
• A glimpse into:
  • Sequential Change of Measure
  • Approximate Bayesian Computation
  • Pseudo-marginal methods

Logistics

• Exercise 3 is out, due Monday April 8
• Office hour on Monday

Recall: model selection notation

• \(I\): an index over a discrete set of models.
• \(M_i\): event model \(i\) is true
• \({\mathscr{Z}}_i\) for \(i\in I\): latent space for model \(i\), i.e. \(z \in {\mathscr{Z}}_i\) when model \(i\) is selected.
• \(p_i\), \(\ell_i\), \(m_i\): prior, likelihood, and marginal likelihood densities for model \(i\), i.e.

\[ m_i = \int p_i(z) \ell_i(z) {\text{d}}z \]

Recall: key idea for Bayesian model selection

Put a prior \(p\) on \(I\), and make the uncertainty over models part of the probabilistic model.

The new joint probability density is given by:

\[ p((i, z), x) = p(i) p_i(z) \ell_i(x | z), \]

where \((i, z)\) is a member of a new latent space given by:

\[ {\mathscr{Z}}= \bigcup_{i\in I} \left( \{i\} \times {\mathscr{Z}}_i \right), \]

Recall: Bayes factor

Ratio of the marginal likelihood for two models:

\[ B_{12} = \frac{m_1(x)}{m_2(x)} \]

Values of \(B_{12}\) greater than 1.0 favor model #1 over #2. Values smaller than 1.0 favor #2 over #1.

Computation of Bayes factor / marginal likelihood

Computing \(m_i\) is hard (19 dubious ways to compute the marginal likelihood of a phylogenetic tree topology)

I will outline some key methods:

• Model saturation (done)
• Thermodynamic integration (done)
• Sequential Monte Carlo via change of measure (today)
• Reversible Jump MCMC (today)

Model saturation: spike and slab example

The idea is to build an augmented model, which can be written as a graphical model, and from which we can still approximate \(m_i(x)\).

Construction of the auxiliary latent space:

• Instead of defining the global latent space as a union of each model’s latent space, define it as a product space,
• and add to that an indicator \(\mu\) that selects which model to use to explain the data. The event \(M_1\) corresponds to \(\mu = 1\) and \(M_2\), to \(\mu = 2\).

Variation on this theme: spike and slab via saturation

• For a spike and slab with 1 parameter, a state has the form \((s,x)\) where \(s\in\{0,1\}\) determines if the variable is selected (non-zero), and \(x\) is the value of the parameter if selected.
• So the original space is \(\{(0,0)\} \cup \{(1, x) : x\in{\mathbf{R}}\}\). What are: \(p\)? \(p_i\)? \(\ell_i\)?
• Saturated space: \(\{0, 1\} \times {\mathbf{R}}\)
• Extend to \(d\) parameters with \((\{0, 1\} \times {\mathbf{R}})^d\)

Recall: Thermodynamic integration

• Here as in parallel tempering we consider a continuum of intermediate distributions from the prior \(p(z)\) to the posterior \(\pi(z)\).
• Denote the intermediate distributions \(\pi_t = \gamma_t / C_t\), \(t\in[0,1]\) given by \(\gamma_t(z) = p(z) (\ell(x|z))^t\).

We start with the following identity, which holds under the assumption that we can swap a derivative and an integral (see for example Folland, Real Analysis):

\[ \begin{aligned} \log C_1 &= \int_0^1 \int \pi_t(z) \log \ell(x|z) \\ &=: \int_0^1 g(t) \end{aligned} \]

• This can be approximated by numerical integration from the output of a Parallel Tempering run. The function \(g(t)\) at one of the grid points \(t = t_i\) is computed as a Monte Carlo average from the \(n\) samples \(z_{t,1}, z_{t,2}, \dots, z_{t,n}\) for parallel chain \(j\): \[ g(t) \approx \frac{1}{n} \sum_{j=1}^n \log \ell(x, z_{t,j}) \]
• This is one of the two methods available in Blang to automatically compute the marginal likelihood of a model

SMC: wrap-up

• The other method used in Blang for marginal likelihood estimates
• Also used to initialize the PT sampler

See files/SMC.pdf

Factoid useful soon: SMC’s estimator of the marginal likelihood is unbiased.

Pseudo-marginal methods

Setup

• Let \(\theta\) denote a parameter from an arbitrary space with prior density \(p(\theta)\).
• Let \(y\) denote data with likelihood \(p(y | \theta)\).
• Our goal is to construct an MCMC method with a stationary distribution given by \(\pi(\theta) = p(\theta | y)\).

Twist in the pseudo-marginal context (Beaumont, 2003): we assume that the likelihood \(p(y|\theta)\) is difficult or impossible to compute pointwise, but that it can be estimated in a positive unbiased fashion.

Formally: there is a function \(T(u, y, \theta)\), where \(u\) is an auxiliary variable in an arbitrary space, with density \(m(u|y,\theta)\), such that:

• \(T \ge 0\), (positive assumption)
• \(\int T(u, y, \theta) m(u|y,\theta) {\text{d}}u = p(y|\theta),\) (unbiased assumption)
• we can sample from \(m(\cdot|y,\theta)\), and evaluate \(T\) point-wise. (computability)

For example, in the context a state space model:

• \(\theta\) would be static (global) parameters,
• \(y = (y_1, \dots, y_n)\), a vector of time-indexed observations conditionally independent given:
  • \(\theta\) and
  • a latent Markov chain \(x = (x_1, \dots, x_n)\),
• \(u\), all the random variables involved in the simulation of an SMC algorithm (for the proposals and resampling steps), and
• \(T\), the SMC estimator of the marginal likelihood given static parameters
  • The fact that this estimator is positive unbiased is a standard result in the SMC literature.

Pseudo marginal methods: executive summary

• Do MCMC on \(\theta\).
• You can replace \(p(y | \theta)\) in the Metropolis-Hastings ratio by the positive unbiased estimator.
• Surprisingly, this preserves the nice theoretical guarantees of MCMC if done carefully.

Pseudo marginal methods: some details

From any positive unbiased estimator, we can now build a Metropolis-Hastings algorithms on a certain extended space. The extended target is given by: \[\begin{equation}\label{eq:extended} \tilde \pi(\theta, u) \propto p(\theta) T(u, y, \theta) m(u | y, \theta). \end{equation}\]

By Assumption (unbiased), the posterior distribution of interest is indeed a marginal of the extended target: \[ \int \tilde \pi(\theta, u) {\text{d}}u = \pi(\theta). \] Assumption (positive) and the above argument imply that the extended target is indeed a well defined density.

Now we build a Metropolis-Hastings algorithm on the extended target as follows. Let \(q(\theta'|\theta)\) denote a user-provided proposal on the parameters. We augment it to a proposal on \(u\) and \(\theta\) using \[ q(\theta', u'|\theta, u) = q(\theta'|\theta) m(u'|y, \theta'). \] The Metropolis-Hastings ratio arising from this proposal is given by:

\[ \begin{align*} \frac{\tilde \pi(\theta', u')}{\tilde \pi(\theta, u)} \frac{q(\theta, u|\theta', u')}{q(\theta', u'|\theta, u)} &= \frac{p(\theta') T(u', y, \theta') m(u' | y, \theta')}{p(\theta) T(u, y, \theta) m(u | y, \theta)} \frac{q(\theta|\theta') m(u|y, \theta)}{q(\theta'|\theta) m(u'|y, \theta')} \\ &= \frac{p(\theta') T(u', y, \theta')}{p(\theta) T(u, y, \theta)} \frac{q(\theta|\theta') }{q(\theta'|\theta)}. \end{align*} \]

From Assumption (computability), this ratio can be computed, so we have an MCMC algorithm targeting \(\tilde \pi\). By discarding the \(u\) part of the samples, we therefore get an MCMC targeting \(\pi\).

This yields the following algorithm:

• Initialize \(\theta_0\) arbitrarily.
• Sample \(u_0 \sim m(\cdot|y, \theta_0)\).
• For MCMC iteration \(\iota = 1, 2, \dots\)
  • Propose a parameter using the user-provide proposal, \(\theta^* \sim q(\cdot | \theta_{\iota-1})\).
  • Sample \(u^* \sim m(\cdot|y, \theta^*)\).
  • Sample \(v \sim \text{Uniform}(0,1)\).
  • If \[ v < \frac{p(\theta^*)}{p(\theta_{\iota - 1})} \frac{T(u^*, y, \theta^*)}{T(u_{\iota - 1}, y, \theta_{\iota - 1})} \frac{q(\theta_{\iota-1} | \theta^*)}{q(\theta^* | \theta_{\iota-1})}. \]
    • Accept: \(u_\iota = u^*\), \(\theta_\iota = \theta^*\).
  • Else
    • Reject: \(u_\iota = u_{\iota-1}\), \(\theta_\iota = \theta_{\iota-1}\).

Side note: Approximate Bayesian Computation (ABC)

• What if there is no positive unbiased estimator available for the likelihood \(p(y|\theta)\)?
• Idea:
  • approximate the likelihood by generating synthetic datasets (often easier than computing the likelihood),
  • the closer the generated data to the actual data, the higher the estimate of the likelihood
• Simplest example: two parameters, \(\theta_1\) and \(\theta_2\), discrete data \(y\)
  • Simulate \(y_{1,1}, y_{1,2}, \dots, y_{1,N} \sim p(\cdot | \theta_1)\)
  • Compute \(\hat p(y | \theta_1) = \frac{1}{N} \sum_i {{\bf 1}}[y_{1,i} = y] \to p(y|\theta_1)\)
  • Then \(y_{2,1}, y_{2,2}, \dots, y_{2,N} \sim p(\cdot | \theta_2)\)
  • Compute \(\hat p(y | \theta_2) = \frac{1}{N} \sum_i {{\bf 1}}[y_{2,i} = y] \to p(y|\theta_2)\)
• What if \(y\) is continuous/combinatorial?
  • Assume we have a distance measure on datasets \(y\), \(d(y, y')\) (often based on summary statistics)
  • Let \(\epsilon > 0\) denote a tolerance
  • Use the following estimator instead: \(\hat p(y | \theta) = \frac{1}{N} \sum_i {{\bf 1}}[d(y_{1,i}, y) < \epsilon] \to p(\{y' : d(y,y') < \epsilon\}|\theta)\)
  • Note: RHS is not equal to \(p(y|\theta)\), hence the “Approximate” in ABC
• State of the art: more than two \(\theta\)’s, automatic tuning of \(\epsilon\)
  • Sequential change of measure with \(\pi_n(\theta) \propto p(\theta) p(\{y' : d(y,y') < \epsilon_n\}|\theta)\)
  • SMC proposal: MCMC moves on \(\theta\)
  • Accept reject using MH ratio based on \(\pi_n\)
  • Use an adaptive SMC scheme to pick the sequence \(\epsilon_n\), see Del Moral, Doucet, Jasra

Reversible jump

• Stay in a “union space” of models… \[ \tilde {\mathscr{Z}}= \bigcup_{i\in I} \left( \{i\} \times \tilde {\mathscr{Z}}_i \right), \]
• but make the dimensionality of the space in the union match, dim(\(\tilde {\mathscr{Z}}_i\)) = dim(\(\tilde {\mathscr{Z}}_j\))
• How? “Pad” with auxiliary iid random variables.

Key advantage:

• We do not to instantiate all the auxiliary random variable.
• Lazy computation: only sample these auxiliary random variable when they will be needed.
• This means we can have an infinite number of auxiliary variables!
• This becomes important when \(I\) is countable infinite, e.g. for non-parametric models.

Reversible jump algorithm

• We pad a variable number of auxiliary variables in order to be able to build diffeomorphic mappings \(\Psi\) (more specifically, mappings with non-vanishing Jacobians). The mappings \(\Psi\) can be thought of as proposals.
• We many need more than one \(\Psi_j\), selected at random according to some probabilities \(\rho_{i\to j}\).

Dimensionality matching: a necessary conditions for the mapping to be diffeomorphic is that the input dimensionality of \(\Psi\) should match the output dimensionality of \(\Psi\).

Consequence: let us say that we want to “jump” from a model with \(m_1\) dimensions into one with \(m_2\) dimensions. What constraints do we have on the number \(n_1\) of auxiliary variables we add to the first model, and the number \(n_2\) we add to the second?

Notation:

• \(p(i)\) prior on model \(i\)
• \(\pi_i\) posterior given model \(i\)
• \(i,i'\) old and proposed model indices
• \(x, x'\) old and proposed model parameters
• \(u_i\): auxiliary variables before the move, input into \(\Psi_j\), with density \(g_i\)
• \(u_{i'}\): auxiliary variables after the move, output of \(\Psi_j\), with density \(g_{i'}\)

Ratio for RJMCMC:

\[ \frac{p(i')\pi_{i'}(x')}{p(i)\pi_i(x)} \frac{\rho_{i'\to i}}{\rho_{i\to i'}} \frac{g_{i'}(u_{i'})}{g_{i}(u_{i})} \left| J(x', u_2) \right| \]

Example: textbook, page 365.

Office hour materials

Set-up

• Exercise instructions
• Blang IDE/eclipse
• CLI (Command Line Interface)
• Web version
• Eclipse on Linux/PC (unsupported)

Upgrading Blang (optional)

• Only needed for eclipse/CLI
• In build.gradle locate the line that starts with compile group: 'ca.ubc.stat', name: 'blangSDK', version: ...
• Change the version to the latest available at: this page
• Go to the root of the repo, and type..
  • If using CLI: ./gradlew clean then ./gradlew installDist
  • If using eclipse: ./gradlew eclipse then in eclipse,
    • right click on project and select Refresh
    • select Project > Clean

Hierarchical model for Ariane rockets

Setting up the code

• Blang code for Ariane hierarchical model
  • For web version, create ArianeHierarchical.bl the root
  • For eclipse (i.e. Blang IDE) or CLI (Command Line Interface), in src/java/demo/ArianeHierarchical.bl
• Usual data: failure_counts.csv
• Command line arguments: copy paste from below
  • For the web version, put in configuration.txt
  • For eclipse/CLI, keep handy, will use when we run the code

--model.data data/failure_counts.csv 
--model.rocketTypes.name LV.Type
--postProcessor DefaultPostProcessor
--postProcessor.imageFormat pdf
--engine PT 
--engine.nChains 10

High-level overview of Blang

See files/blang-high-level.pdf

Running the code

• Web: click on the triangle button
• Eclipse: see bottom of this page
• CLI: use the following

./gradlew clean
./gradlew installDist
java -Xmx2g -cp build/install/blang520Assign/lib/\\* demo.ArianeHierarchical [COMMAND LINE OPTIONS]

Some quick experimentation with this model

Mini-exercise:

• Use Bayesian model selection to decide to use a flat or hierarchical model on this data.
• Compare the normalization constant estimate for Parallel Tempering and Sequential Monte Carlo
• Compare two implementations of the non-hierarchical model, one given by setting the argument --model.useHierarchical false, the other from ArianeFlat.bl available here

Deciding on flat vs hierarchical

Idea: I run the code first with --model.useHierarchical true, then with --model.useHierarchical false, and compute the Bayes factor from that.

Interpreting the standard out messages

We will go over the meaning of the output below in the office hours.

• SMC-based initialization
• SMC diagnostics:
  • annealing parameters
  • SMC effective sampling size (ESS)
• Adaptive rounds
• Parallel tempering diagnostics
  • swapSummaries [ round=8 lowest=0.574 average=0.658 ]
  • logNormalizationContantProgress [ round=8 value=-20.692 ]
  • estimatedLambda [ round=8 value=3.075 ]
  • actualTemperedRestarts [ round=8 count=34 rate=0.068 ]
  • asymptoticRoundTripBound [ round=8 count=61.347 rate=0.123 ]

Note: data source data/failure_counts.csv does not contain column failureprobabilities --- treating as missing
WARNING: There were files not up to date in the code repository (see /Users/bouchard/w/blangExample/results/all/2019-04-01-12-43-49-TokB5Hso.exec/executionInfo/code/dirty-files.txt)
Preprocess {
  17 samplers constructed with following prototypes:
    RealScalar sampled via: [RealSliceSampler]
  Initialization {
    Warning: small concentrations may cause numeric instability to Dirichlet and Beta distributions. Consider enforcing a lower bound of say 0.5 This message may also occur when slice samling outside of such constraint, you can then ignore this message. 
    Propagation [ annealParam=0.0 ess=0.900 ]
    Propagation [ annealParam=0.001 ess=0.736 ]
    Propagation [ annealParam=0.004 ess=0.563 ]
    Propagation [ annealParam=0.008 ess=0.410 ]
    Resampling [ iter=3 annealParam=0.014 logNormalization=-1.678 ]
    Propagation [ annealParam=0.014 ess=0.900 ]
    Propagation [ annealParam=0.022 ess=0.717 ]
    Propagation [ annealParam=0.034 ess=0.548 ]
    Propagation [ annealParam=0.051 ess=0.413 ]
    Resampling [ iter=7 annealParam=0.082 logNormalization=-4.323 ]
    Propagation [ annealParam=0.082 ess=0.958 ]
    Propagation [ annealParam=0.111 ess=0.793 ]
    Propagation [ annealParam=0.179 ess=0.677 ]
    Propagation [ annealParam=0.222 ess=0.463 ]
    Resampling [ iter=0 annealParam=0.305 logNormalization=-9.731 ]
    Propagation [ annealParam=0.305 ess=0.992 ]
    Propagation [ annealParam=0.333 ess=0.849 ]
    Propagation [ annealParam=0.434 ess=0.827 ]
    Propagation [ annealParam=0.444 ess=0.596 ]
    Propagation [ annealParam=0.556 ess=0.414 ]
    Resampling [ iter=0 annealParam=0.667 logNormalization=-17.196 ]
    Propagation [ annealParam=0.667 ess=0.900 ]
    Propagation [ annealParam=0.771 ess=0.889 ]
    Propagation [ annealParam=0.778 ess=0.707 ]
    Propagation [ annealParam=0.889 ess=0.548 ]
  } [ endingBlock=Initialization blockTime=338.4ms blockNErrors=1 ]
} [ endingBlock=Preprocess blockTime=587.1ms blockNErrors=1 ]
Inference {
  Round(1/9) {
    Performing 1220 moves... [ nScans=2 nChains=10 movesPerScan=61 ]
    swapSummaries [ round=0 lowest=0.0 average=0.375 ]
    logNormalizationContantProgress [ round=0 value=-142.950 ]
    estimatedLambda [ round=0 value=5.622 ]
    actualTemperedRestarts [ round=0 count=0 rate=0.0 ]
    asymptoticRoundTripBound [ round=0 count=0.151 rate=0.076 ]
  } [ endingBlock=Round(1/9) blockTime=71.62ms blockNErrors=0 ]
  Round(2/9) {
    Performing 2440 moves... [ nScans=4 nChains=10 movesPerScan=61 ]
    swapSummaries [ round=1 lowest=0.100 average=0.477 ]
    logNormalizationContantProgress [ round=1 value=-17.644 ]
    estimatedLambda [ round=1 value=4.710 ]
    actualTemperedRestarts [ round=1 count=0 rate=0.0 ]
    asymptoticRoundTripBound [ round=1 count=0.350 rate=0.088 ]
  } [ endingBlock=Round(2/9) blockTime=81.04ms blockNErrors=0 ]
...
  Round(8/9) {
    Performing 153110 moves... [ nScans=251 nChains=10 movesPerScan=61 ]
    swapSummaries [ round=7 lowest=0.595 average=0.668 ]
    logNormalizationContantProgress [ round=7 value=-20.780 ]
    estimatedLambda [ round=7 value=2.992 ]
    actualTemperedRestarts [ round=7 count=11 rate=0.044 ]
    asymptoticRoundTripBound [ round=7 count=31.436 rate=0.125 ]
  } [ endingBlock=Round(8/9) blockTime=811.1ms blockNErrors=0 ]
  Round(9/9) {
    Performing 305000 moves... [ nScans=500 nChains=10 movesPerScan=61 ]
    swapSummaries [ round=8 lowest=0.574 average=0.658 ]
    logNormalizationContantProgress [ round=8 value=-20.692 ]
    estimatedLambda [ round=8 value=3.075 ]
    actualTemperedRestarts [ round=8 count=34 rate=0.068 ]
    asymptoticRoundTripBound [ round=8 count=61.347 rate=0.123 ]
  } [ endingBlock=Round(9/9) blockTime=1.473s blockNErrors=0 ]
} [ endingBlock=Inference blockTime=3.505s blockNErrors=0 ]
Postprocess {
  Post-processing allLogDensities
  Post-processing energy
  Post-processing failureProbabilities
  Post-processing logDensity
  Post-processing m
  Post-processing s
  MC diagnostics
} [ endingBlock=Postprocess blockTime=27.55s blockNErrors=0 ]
executionMilliseconds : 31649
outputFolder : /Users/bouchard/w/blangExample/results/all/2019-04-01-12-43-49-TokB5Hso.exec

Looking at the files produced

Look at results/latest. We will go over those in the office hours. A copy of what is expected is available here.

• Main thing to look at for this exercise: monitoringPlots/logNormalizationContantProgress.pdf
• This suggest the final value, -21.017 is a reliable estimate of \(\log(m_1)\), the log probability of the data under model 1 (the hierarchical one)

Understanding the command line options

Add --help to the options. You should see the following, which we will go over:

Note: data source data/failure_counts.csv does not contain column failureprobabilities --- treating as missing
 --checkIsDAG <boolean> (default value: true)

 --engine <PosteriorInferenceEngine: SCM|PT|Forward|Exact|None|fully qualified> (default value: SCM)

 --engine.adaptFraction <double> (default value: 0.5)

 --engine.initialization <InitType: COPIES|FORWARD|SCM> (default value: SCM)

 --engine.ladder <TemperatureLadder: Geometric|EquallySpaced|Polynomial|UserSpecified|fully qualified> (default value: EquallySpaced)

 --engine.nChains <Integer> (optional)
   description: If unspecified, use the number of threads.

 --engine.nPassesPerScan <double> (default value: 3)

 --engine.nScans <int> (default value: 1_000)

 --engine.nThreads <Cores: Single|Max|Dynamic|Fixed|fully qualified> (default value: Dynamic)

 --engine.nThreads.fraction <double> (default value: 0.5)

 --engine.nThreads.ignoreUtilizedCores <boolean> (default value: true)

 --engine.nThreads.verbose <boolean> (default value: false)

 --engine.random <Random> (default value: 1)

 --engine.reversible <boolean> (default value: false)

 --engine.scmInit.maxAnnealingParameter <double> (default value: 1.0)
   description: Use higher values for likelihood maximization

 --engine.scmInit.nFinalRejuvenations <int> (default value: 5)
   description: Number of rejuvenation passes to do after the change of measure.

 --engine.scmInit.nParticles <int> (default value: 1_000)

 --engine.scmInit.nThreads <Cores: Single|Max|Dynamic|Fixed|fully qualified> (default value: Dynamic)

 --engine.scmInit.nThreads.fraction <double> (default value: 0.5)

 --engine.scmInit.nThreads.ignoreUtilizedCores <boolean> (default value: true)

 --engine.scmInit.nThreads.verbose <boolean> (default value: false)

 --engine.scmInit.random <Random> (default value: 1)
   description: Random seed used for proposals and resampling.

 --engine.scmInit.resamplingESSThreshold <double> (default value: 0.5)
   description: If the (relative) Effective Sample Size (ESS) falls below, perform a resampling round.

 --engine.scmInit.resamplingScheme <ResamplingScheme: STRATIFIED|MULTINOMIAL> (default value: STRATIFIED)

 --engine.scmInit.temperatureSchedule <TemperatureSchedule: AdaptiveTemperatureSchedule|FixedTemperatureSchedule|fully qualified> (default value: AdaptiveTemperatureSchedule)
   description: Algorithm selecting annealing parameter increments.

 --engine.scmInit.temperatureSchedule.nudgeFromZeroIfOutOfSupport <double> (default value: 1e-10)
   description: If all particles are out of support at first iteration, nudge the temperature a bit so that support constraints kick in.

 --engine.scmInit.temperatureSchedule.threshold <double> (default value: 0.9999)
   description: Annealing parameter is selected to get the (conditional) ESS decrease specified by this parameter.

 --engine.scmInit.temperatureSchedule.useConditional <boolean> (default value: true)
   description: See Zhou, Johansen and Aston (2013).

 --engine.targetAccept <Double> (optional)

 --engine.usePriorSamples <boolean> (default value: true)

 --excludeFromOutput <List: Space separated items or "file <path>" to load from newline separated file> (optional)

 --experimentConfigs.configFile <File> (optional)
   description: If set, use those arguments in provided file that do not appear in the provided arguments.

 --experimentConfigs.managedExecutionFolder <boolean> (default value: true)
   description: Automatically organize results into subdirectories of 'results/all'?

 --experimentConfigs.maxIndentationToPrint <int> (default value: inf)
   description: Use -1 to silence all HLogs output

 --experimentConfigs.recordExecutionInfo <boolean> (default value: true)
   description: Record information such as timing, main class, code version, etc for this run?

 --experimentConfigs.recordGitInfo <boolean> (default value: true)

 --experimentConfigs.saveStandardStreams <boolean> (default value: true)
   description: Save combined standard out and err into a file?

 --experimentConfigs.tabularWriter <TabularWriterFactory: CSV|Spark|fully qualified> (default value: CSV)

 --initRandom <Random> (default value: 1)

 --model <ModelBuilder: fully qualified>

 --model.data <GlobalDataSource: Path to the DataSource.>

 --model.data.reader <DataSourceReader: CSV|fully qualified> (default value: CSV)

 --model.data.reader.commentCharacter <Character> (optional)

 --model.data.reader.ignoreLeadingWhiteSpace <boolean> (default value: true)

 --model.data.reader.separator <char> (default value: ,)

 --model.data.reader.strictQuotes <boolean> (default value: false)

 --model.failureProbabilities.dataSource <DataSource: Path to the DataSource.>

 --model.failureProbabilities.dataSource.reader <DataSourceReader: CSV|fully qualified> (default value: CSV)

 --model.failureProbabilities.dataSource.reader.commentCharacter <Character> (optional)

 --model.failureProbabilities.dataSource.reader.ignoreLeadingWhiteSpace <boolean> (default value: true)

 --model.failureProbabilities.dataSource.reader.separator <char> (default value: ,)

 --model.failureProbabilities.dataSource.reader.strictQuotes <boolean> (default value: false)

 --model.failureProbabilities.name <ColumnName> (optional)
   description: Name of variable in the plate

 --model.m <RealVar: A number or NA> (optional)

 --model.numberOfFailures.dataSource <DataSource: Path to the DataSource.>

 --model.numberOfFailures.dataSource.reader <DataSourceReader: CSV|fully qualified> (default value: CSV)

 --model.numberOfFailures.dataSource.reader.commentCharacter <Character> (optional)

 --model.numberOfFailures.dataSource.reader.ignoreLeadingWhiteSpace <boolean> (default value: true)

 --model.numberOfFailures.dataSource.reader.separator <char> (default value: ,)

 --model.numberOfFailures.dataSource.reader.strictQuotes <boolean> (default value: false)

 --model.numberOfFailures.name <ColumnName> (optional)
   description: Name of variable in the plate

 --model.numberOfLaunches.dataSource <DataSource: Path to the DataSource.>

 --model.numberOfLaunches.dataSource.reader <DataSourceReader: CSV|fully qualified> (default value: CSV)

 --model.numberOfLaunches.dataSource.reader.commentCharacter <Character> (optional)

 --model.numberOfLaunches.dataSource.reader.ignoreLeadingWhiteSpace <boolean> (default value: true)

 --model.numberOfLaunches.dataSource.reader.separator <char> (default value: ,)

 --model.numberOfLaunches.dataSource.reader.strictQuotes <boolean> (default value: false)

 --model.numberOfLaunches.name <ColumnName> (optional)
   description: Name of variable in the plate

 --model.rocketTypes.dataSource <DataSource: Path to the DataSource.>

 --model.rocketTypes.dataSource.reader <DataSourceReader: CSV|fully qualified> (default value: CSV)

 --model.rocketTypes.dataSource.reader.commentCharacter <Character> (optional)

 --model.rocketTypes.dataSource.reader.ignoreLeadingWhiteSpace <boolean> (default value: true)

 --model.rocketTypes.dataSource.reader.separator <char> (default value: ,)

 --model.rocketTypes.dataSource.reader.strictQuotes <boolean> (default value: false)

 --model.rocketTypes.maxSize <Integer> (optional)

 --model.rocketTypes.name <ColumnName> (optional)

 --model.s <RealVar: A number or NA> (optional)

 --model.useHierarchical <Boolean> (optional)

 --postProcessor <PostProcessor: DefaultPostProcessor|NoPostProcessor|fully qualified> (default value: NoPostProcessor)

 --postProcessor.blangExecutionDirectory <File> (optional)
   description: When called from Blang, this will be the latest run, otherwise point to the .exec folder created by Blang

 --postProcessor.burnInFraction <double> (default value: 0.5)

 --postProcessor.essEstimator <EssEstimator: BATCH|ACT|AR> (default value: BATCH)

 --postProcessor.experimentConfigs.configFile <File> (optional)
   description: If set, use those arguments in provided file that do not appear in the provided arguments.

 --postProcessor.experimentConfigs.managedExecutionFolder <boolean> (default value: true)
   description: Automatically organize results into subdirectories of 'results/all'?

 --postProcessor.experimentConfigs.maxIndentationToPrint <int> (default value: inf)
   description: Use -1 to silence all HLogs output

 --postProcessor.experimentConfigs.recordExecutionInfo <boolean> (default value: true)
   description: Record information such as timing, main class, code version, etc for this run?

 --postProcessor.experimentConfigs.recordGitInfo <boolean> (default value: true)

 --postProcessor.experimentConfigs.saveStandardStreams <boolean> (default value: true)
   description: Save combined standard out and err into a file?

 --postProcessor.experimentConfigs.tabularWriter <TabularWriterFactory: CSV|Spark|fully qualified> (default value: CSV)

 --postProcessor.facetHeight <double> (default value: 2.0)
   description: In inches

 --postProcessor.facetWidth <double> (default value: 4.0)
   description: In inches

 --postProcessor.imageFormat <String> (default value: png)

 --postProcessor.rCmd <String> (default value: Rscript)

 --printAccessibilityGraph <boolean> (default value: false)

 --samplers.additional <SamplerSet: Fully qualified instances of blang.mcmc.Sampler>
   description: Samplers to be added.

 --samplers.excluded <SamplerSet: Fully qualified instances of blang.mcmc.Sampler>
   description: Samplers to be excluded (only useful if useAnnotation = true).

 --samplers.useAnnotation <boolean> (default value: true)
   description: If the arguments of the annotations @Samplers should be used to determine a starting set of sampler types.

 --stripped <boolean> (default value: false)
   description: Stripped means that the construction of forward simulators and annealers is skipped

 --treatNaNAsNegativeInfinity <boolean> (default value: false)

 --version <String> (optional)
   description: Version of the blang SDK to use (see https://github.com/UBC-Stat-ML/blangSDK/releases), of the form of a git tag x.y.z where x >= 2. If omitted, use the local SDK's 'master' version.

Model selection

Now I run again with --model.useHierarchical false and get \(\log(m_2) \approx -37.417\). Therefore the Bayes factor is given by:

\[ \frac{m_1}{m_2} = \exp( \log m_1 - \log m_2) = \exp((-20.692) - (-37.417)) = 18347429 \]

So the evidence for the hierarchical model vs. flat model is “decisive” (see this page for terminology, but use with a pinch of salt!)

Model selection, again

Note: the version with --model.useHierarchical false is a bit strange in the sense that the random variables \(m\) and \(s\) are still sampled over (but not used in the hierachical).

• A frequentist may worry this means that removing these parameters will increase the model score for non-hiearchical.
• Let’s check!
• Here is a straightforward flat model:

model ArianeFlat {
    
  param GlobalDataSource data
  param Plate<String> rocketTypes
  param Plated<IntVar> numberOfLaunches
  random Plated<RealVar> failureProbabilities
  random Plated<IntVar> numberOfFailures
  
  laws {
    for (Index<String> rocketType : rocketTypes.indices.filter[key.startsWith("Ariane")]) {
      failureProbabilities.get(rocketType)  ~ Beta(1.0, 1.0)
      numberOfFailures.get(rocketType)
        | RealVar failureProbability = failureProbabilities.get(rocketType),
          IntVar numberOfLaunch = numberOfLaunches.get(rocketType)
        ~ Binomial(numberOfLaunch, failureProbability)
    }
  }
}

Running this we get an estimate of \(m_2 \approx -37.589\) which is very close to our old estimate, \(-37.417\).

In fact, you can show that as the number of PT iterations and chains go to infinity, these numbers will converge to the same value.

Conclusion: this confirms we do not need to worry about parameters counting when using Bayes factors.