# Exploring Code Style Transfer with Neural Networks

Karl Munson*, Anish Savla*, Chih-Kai Ting*, Serenity Wade*
University of California, Santa Cruz
{kmunson, asavla, cting3<EMAIL_ADDRESS>Kiran Kate  Kavitha Srinivas
IBM Research

###### Abstract

Style is a significant component of natural language text, reflecting a change
in the tone of text while keeping the underlying information the same. Even
though programming languages have strict syntax rules, they also have style.
Code can be written with the same functionality but using different language
features. However, programming style is difficult to quantify, and thus as
part of this work, we define style attributes, specifically for Python. To
build a definition of style, we utilized hierarchical clustering to capture a
style definition without needing to specify transformations. In addition to
defining style, we explore the capability of a pre-trained code language model
to capture information about code style. To do this, we fine-tuned pre-trained
code-language models and evaluated their performance in code style transfer
tasks.

##  1 Introduction

###  1.1 Motivation

Recent years have seen a sharp increase in research on code language modeling,
with work ranging from pre-training large language models over code (Ahmad et
al. (2021), Wang et al. (2021), Feng et al. (2020)), as well as related work
examining the application of other natural language processing tasks applied
to code, such as code machine translation or code generation. Szafraniec et
al. (2022), Chen et al. (2021). However, accompanying this rise has been a
lack of work studying how well these models understand the human notion of
code style.

Code style is an integral element of code, so many organizations regulate code
style through code reviews to improve readability and maintainability. Code
style is, broadly, the appearance of code to a viewer and is thus - as much
the same as writing style - something very individual to its writer. Due to
the strict syntactic grammars governing how code is structured, many aspects
of code style can be quantified automatically. Yun Peng (2021)
programmatically studies the usage of Python language features across publicly
available repositories.

We assert that the content of code (its functionality) and the style of code
(its aesthetics) are, like the content and style of text, aspects whose latent
representations can be extricated and individually manipulated through neural
style transfer.

In this work, we defined an embedding space whose dimensions are extricable
metrics measuring different aspects of code style. We then explored how to
analyze code style through clustering in this space and how well we can
classify code into these clusters. From there, we generated data and fine-
tuned, pre-trained code language models to perform code style transfer along
one aspect of code style.

###  1.2 Contribution

The main contributions of this work are fourfold. First, we applied clustering
algorithms over code scripts embedded in a code style embedding space. We did
this to find code clusters with distinct styles and evaluate the quality and
predictability of the discovered clusters. Second, from publicly available
code repositories, we generated a suite of parallel corpora; each corpus
represents an atomic style transformation, such as casing or comments. Third,
we used these parallel corpora to fine-tune pre-trained code language models
on sequence-to-sequence individual style transfer tasks. These models can
perform precise stylistic alterations one at a time, such as modifying casing
or adding comments. We also trained a joint model that can apply multiple
style transfers at once. Finally, whereas the existing workaround code
language models focus on the syntactic integrity of processed code, our work
focuses on something novel: code style.

##  2 Related Work / Background

###  2.1 Code Style

Due to the subjective nature of code style, the first background area of
interest to address is what constitutes code style. This paper focuses
primarily on style aspects of Python as a result of the large number of
stylistic choices available in Python, large number of publicly available code
samples, and familiarity with the language.

Many of the python language features that we focused on were selected based on
the analysis done in An Empirical Study for Common Language Features Used in
Python Projects Yun Peng (2021). This paper tries to study the use and impact
of Python language features in real-world Python projects. The authors analyze
Python language features and automatically identify the use of 22 kinds of
common Python language features in 6 categories in Python source code. This
study was conducted over 35 popular Python projects from eight application
domains covering 4.3 million lines of code; to investigate the usage of these
language features in the project. Inheritance, decorator use, keyword
arguments, for loops, and nested classes are listed as the top 5 used language
features.

###  2.2 Code Language Models

Recent years have seen the natural language processing literature dominated by
the pre-train and fine-tune paradigm. In line with this, we will primarily use
two pre-trained language models for code for our experiments, PLBART and
CodeT5. Pre-trained model usage will eliminate the need for extensive
resources for training a code language model from scratch.

CodeT5 is an encoder-decoder model pre-trained on masked span prediction,
identifier tagging, and masked identifier prediction for bimodal, programming,
and natural language, dual generation Wang et al. (2021). It uses
CodeSearchNet data which consists of unimodal (programming languages only) and
bimodal (programming language and natural language) data on six programming
languages. Developed by Salesforce, it is available on HuggingFace in
different model sizes.

PLBART is another encoder-decoder model we explored. The model is pre-trained
in Java, Python, and natural language on reconstruction and denoising Ahmad et
al. (2021). Like CodeT5, it is publicly available on HuggingFace in different
model sizes, and has checkpoints fine-tuned on various code-related tasks.

Even though language models for code are used for many tasks, such as bug
detection, code translation, and code generation; to the best of our
knowledge, our work is the first one to address the task of code style
transfer.

###  2.3 Text Style Transfer

Text style transfer is the task of reformatting a natural language text from
one style to a text with the same meaning in a different style, while
maintaining the textual content of the source.(See figure 1) There are a
number of ways to approach this task, and a significant amount of work has
been done to use neural models for text style transfer (Riley et al. (2021),
Tikhonov et al. (2019), Reif et al. (2021), Jin et al. (2020)). This existing
work in text style transfer in natural language serves as a blueprint on how
to apply the same style transfer to code.

Figure 1:  Example of text style transfer. Riley et al. (2021)

Most approaches to this task utilize either an encoder-decoder framework or a
generative adversarial network with different variations to those two basic
frameworks, such as the one pictured below (2).

Figure 2:  Example of GAN style transfer architecture. Tikhonov et al. (2019)

##  3 Requirements / Proposal

###  3.1 Defining Code Styles

The core underlying topic of this research is code style, which is, most
broadly, the appearance of the code. Code style is vital to the readability of
code. Many organizations value code styles which adhere to a centralized style
guide and conduct code reviews to ensure that the style of written code
matches a certain quality threshold. Though code style is subjective and very
complex, we develop a framework for analyzing specific aspects of code style
through quantifiable metrics and use this framework to conduct experiments
with both clustering and sequence-to-sequence model training.

Listed below are definitions of the aspects of code style we decided to focus
on. These categories describe the main attributes of code style we chose to
focus on for parallel corpora generation. We also used additional attributes
during clustering experiments, mostly related to the above-described traits.

####  3.1.1 Casing

For this metric, we consider the various types of casing found within code:
snake case, lower and upper camel case, and lowercase. Also, we consider how
these casing types appear in different varieties of identifier names, such as
variable, method, or class names.

####  3.1.2 Comment

Comments are vital elements of code style as they greatly aid the readability
of the code. We consider many attributes of code comments to be style-
defining. Comment word length, comment density, and the ratio of commented
lines in a script all act as approximate measures of how much a script has
been annotated with comments.

####  3.1.3 Docstrings

Docstrings, like comments, are written by the programmer to aid the
readability of the code and are integral to code style. They serve as headers
to inform readers of the intended function, input and output arguments, and
types of a method, class, or module. They often contain other information as
well, relating to copyright or authorship.

####  3.1.4 List Comprehensions

We additionally focus on Python language-specific features, such as
comprehensions. List comprehensions are one such language feature, used to
iterate and operate on the items within a data structure. A programmer’s
decision between using a for loop or comprehension to do this is a Python-
specific stylistic choice.

####  3.1.5 Usage of Classes

Being an object-oriented language, Python allows programmers to define
classes. How one chooses to organize classes, design the methods within
classes, or use inheritance, are all additional hallmarks of code style.

###  3.2 Code Style Transfer

Each code sample has a style as defined in the previous section. We want to be
able to map a target style on to a sample of code which does not have that
style. This is the same task as text style transfer but applied to code
snippets.

##  4 Design

###  4.1 Clustering

Code style is difficult to define, so we selected several features (See Table
1) to embed Python code into a multidimensional space. Then we used
unsupervised clustering algorithms to identify the types of code that are most
similar to each other. From these clusters, we define spatial buckets that
represent different and distinct code styles.

To cluster style features, we utilized the HDBSCANMalzer and Baum (2020)
clustering model to group code samples together based on similar counts of
certain style features. HDBSCAN is a hierarchical clustering algorithm that we
approximated by a random forest classifier to perform the prediction of
cluster labels on new data.

Code Style Features for Clustering
---
Feature Type | Features
Casing | Snake case variable usage ratio
Snake case function name usage ratio
Snake case class usage ratio
Upper camel case variable usage ratio
Upper camel case function name usage ratio
Upper camel case class usage ratio
Lower camel case variable usage ratio
Lower camel case function name usage ratio
Lower camel case class usage ratio
Documentation | Docstring Density
Comment Density
Function/Class | Average function decorators usage count
Average class decorators usage count
Average class inheritance usage count
Python Features | Average list comprehension usage count
Average generators usage count
Average lambda function usage count
Table 1: The selected code style features which were used for clustering
experiments.

###  4.2 Style Transfer GAN

We adopt a similar style transfer architecture to Tikhonov et al., the state-
of-the-art text style transfer on the Yelp small review sentiment transfer
task. The model will be slightly different from the original paper, so the
pre-trained encoder-decoder is uninterrupted by inserting style between the
components. The input code is concatenated with the style code and reduced
using a linear layer so that it fits into the pre-trained encoder. The pre-
trained encoder decoder will take this compressed input and produce a styled
code sample that goes into the discriminator. We used the discriminator loss
to modify the style representation concatenated to the input.

##### Substitute Pre-trained Model

We replace the encoder from Tikhonov et al. with the encoder of our pre-
trained model CodeT5 or PLBART. We use the pre-trained encoder to generate a
latent representation from the input code snippet. For the model generator, we
will use the decoder of CodeT5 or PLBART. The generator will output decoded
code. This code passes to the discriminator, a 3-layer multi-layer perceptron;
again following from Tikhonov et al.. The MLP will learn to represent the
style of the generated code, and a style loss passes back to the encoder.

###  4.3 Style Classification with Pre-Trained Models

The GAN-trained model did not generate intelligible Python code let alone
correctly styled code. This begs the question, do pre-trained code language
models even have a notion of code styles?

As a result, we built a classification model from the pre-trained language
models we considered to see if the pre-trained models could capture any style
information from the HBDSCAN clusters. The inputs for this model are code
samples, and the outputs for the model are the labels generated from the
clustering experiments. We implemented two versions of the classifier, one for
PLBART and one for CodeT5. The PLBART model is a pre-trained PLBART model with
a sequence classification head. The CodeT5 model uses a CodeT5 model and fine-
tunes the generation model to produce the labels via generation instead of
through a classification head.

###  4.4 Seq2Seq Style Transfer

The style classification results showed pre-trained code language models have
some ability to capture style elements. As a result, we constructed a
sequence-to-sequence model to do style transformations. Instead of
transferring from one style to another based on complex features, we focused
on individual style transfers like going from a for loop to an equivalent list
comprehension. To do this, we constructed parallel corpora for five transfer
tasks, generating comments, generating docstrings, adding casing, transforming
code to a class structure, and converting for loops to list comprehensions
when relevant. The details for parallel corpora generation are in the data
section.

The sequence-to-sequence model is a pre-trained CodeT5 small instance with 60
million parameters fine-tuned on the parallel corpora. For each of the
individual transformations, we fine-tuned a sequence-to-sequence model to see
if it was capable of performing specified single style transforms. In
addition, we also fine-tuned a combined sequence-to-sequence model on all of
the parallel corpora. The resulting multi-task model can do all the
transformations with only one instance of CodeT5 instead of five. To fine-tune
the model on multiple tasks, the dataloader randomly sample batches from each
of the five parallel corpora. The dataloader samples proportionally to the
size of the parallel corpora meaning that larger corpora are sampled more
based on their size.

The combined sequence-to-sequence model also had a modification to its
training data. For each code sample, a natural language prompt, separated from
the code by a special token, was used to indicate to the model which
transformation the model should perform (figure:4). This allows the desired
transformation to be specified during the inference of the model and for
prompting of multiple transformations to be performed on a single code sample.

##  5 Data

###  5.1 Py150k

We used Py150k Raychev et al. (2016), which consists of parsed Abstract Syntax
Trees from Python 2.7 scripts collected from a set of Github repositories with
permissive and non-viral licenses. The dataset creators removed duplicate
files and project forks, keeping only programs that parse and have at most
30,000 nodes in the AST. The set encompassed code from 5958 unique Github
users and 8422 Github repositories. The dataset is split into two parts:
100,000 files for training and 50,000 for evaluation.

###  5.2 BigQuery

We also used Google BigQuery for additional sources of Python code. The set
contains 1.3 million python files drawn from Github with the watch count
available. The data is split into two parts, samples from the CuBERTKanade et
al. (2020) dataset and new code samples which were added to BigQuery after the
CuBERT paper.

The files from the CuBERT dataset are already deduplicated to avoid
contamination of dataset with py150, and duplicate scripts Kanade et al.
(2020). For non-deduplicated files, we did the deduplication via their paths
in the initial query to remove direct copies of code samples. In order to
ensure deduplication of the remaining files, we calcluated an MD5 checksum on
the contents of each file, and removed duplicates from the dataset. The query
used to collect the big query data was based on the one used in the CuBERT
paper and will be made available in addition to other resources used in this
research.

###  5.3 Parallel Corpora

For the sequence-to-sequence transformation, we generated parallel corpora for
each individual style transformation for training. For all transformations
other than the comment transformation, the input and labels for code samples
had no comments as the AST module used for transformations strips comments; so
for the transformed data to match the original, the original code also needs
comments removed. These parallel corpora transformations are easy to produce
programmatically, but are difficult to undo. The model needs to identify when
to make transformations and correctly produce them which is not a trivial
solution.

####  5.3.1 Casing

For casing, we selected identifiers and method names using the AST module,
lower-cased all of the data and stripped out underscores. This data with
uncased code is the input, and the label is the code sample with the original
casing. The model thus learns how to format identifier names. There are
700,000 examples in this dataset.

####  5.3.2 Docstring

Code samples were split using the AST module isolating methods with
corresponding docstrings. The input to the sequence-to-sequence model has the
docstrings stripped; the label is the original code with the docstring. There
are 2,500,000 examples in this dataset.

####  5.3.3 Comments

We took code samples and parsed them, then unparsed the code using the AST
module. This parse strips the comments and returns the same code with only
comments removed. The input is the uncommented code; the label is the original
commented code. There are 700,000 examples in this dataset.

####  5.3.4 Classes

We took code samples with classes and removed the class definitions and all
references to self in methods. The input for this model is the code without
class-related syntax, and the gold label is the original code. There are
400,000 examples in this dataset.

####  5.3.5 List Comprehensions

Code samples with list comprehensions in them are selected via the AST. To
make the versions without list comprehensions, each list comprehension is
converted to an equivalent function. The input for this model is the code with
list comprehensions converted to for loops, and the gold label is the original
code. This dataset has 36,000 examples.

##  6 Experiments

###  6.1 Code Style Clustering

We collected metrics(table 1) for the scripts based on the style features
specified in the design section and then used HDBScan Malzer and Baum (2020)
to find and evaluate clusters of code. HBDSCAN has three parameters that we
tuned to reduce the number of clusters while keeping clusters distinct.

In addition to measuring cluster distinctness, we also measured the
predictability of the cluster results using several non-neural classifiers.
The experiment tested whether the features have enough signals to build a
high-quality classification model since the code style could be sophisticated
and difficult to characterize. We tested SVM, decision trees, random forest,
naive Bayes, and logistic regression models’ ability to predict cluster labels
given the constructed style space for a script.

###  6.2 Style Transfer GAN

To test the viability of this model, we implemented it as described in the
design section with two different configurations of pre-trained models. In one
case, we used CodeT5; in the other, we used PLBART as the pre-trained
component. In both cases, the GAN model could not output valid python code,
let alone perform the style transfer task, despite tuning parameters and
troubleshooting the model. As a result, we aborted the experiments with this
architecture.

###  6.3 Style Cluster Classification

The failure of the Style Transfer GAN could indicate that neural network
architectures could not capture the styles defined with our clusters. As a
result of this possibility, we tested if CodeT5 or PLBART were capable of
capturing code style at all. In order to test this, we tested different
classifiers built on pre-trained code language models. We tested the CodeT5
small, 60 million parameters, and base, 120 million parameters. For PLBART, we
tested with the “uclanlp/plbart-multi_task-python” checkpoint on HuggingFace
as well as “uclanlp/plbart-base.” Both PLBART models have the exact parameter
count but with different fine-tuning. For the hyperparameters of both models,
we used a learning rate of 1e-4, weight decay of 0.01, and trained for four
epochs. We chose these parameters as similar pre-trained classifiers used
these parameters. Four epochs of training were finished as validation loss was
not improving much. The rest of the hyperparameters used are default
HuggingFace parameters for each model.

###  6.4 Seq2Seq Style Transfer

Style cluster classification showed the ability of pre-trained models to
capture some elements of style. CodeT5 had the best performance for the
classification task, so we next trained a set of sequence-to-sequence models
to perform each style transfer. All these experiments were based on CodeT5
small because the base version had marginally better performance. However,
more parameters make it less efficient to train multiple sequence-to-sequence
models. We built a baseline with CodeT5 architecture but trained it from
scratch to evaluate how much CodeT5 can learn with versus without pre-training
code.

####  6.4.1 Individual Style Transfer

We fine-tuned CodeT5 for each parallel corpora we generated and trained all
models with the same hyperparameters. Each model was fine-tuned for four
epochs for docstring, comments, casing, and class corpora. List comprehensions
corpora had only 36,000 examples, and the corresponding task was thus fine-
tuned for six epochs. Docstring was fine-tuned on a 400,000 example subset of
the 2.5 million examples to reduce training time and make dataset size
comparable to other tasks. For the baseline of this model, we used a
randomized weights version of CodeT5 to train on the task. This model would
have same architecture but none of the information about code learned in pre-
training. By doing this we could show what parts of code style are captured in
pre-training.

####  6.4.2 Multiple Style Transfer

The combined model was trained with a subset of each parallel corpora. For
docstring transfer, 250,000 examples were used; for casing 200,000 examples
were used; for comments, 350,000 examples were used; for class, 350,000
examples were used; and for list comprehensions, 36,000 examples were used.
Each data sample had a natural language prompt attached to the front of the
input to indicate the transformation to be performed. The model samples
batches proportionally to the size of training data, and each batch only
contained examples from one parallel corpus. The model was trained for 12
epochs over this dataset with the same hyperparameters as the individual
models. The model was tested every 4 epochs for performance on each of the
individual tasks using the evaluation set.

In addition to evaluating individual tasks, the model was evaluated using the
evaluation set on doing multiple transformations at once. The model was given
a natural language prompt(figure 4 and figure 5) that specified that multiple
transformations should be done. Examples from the evaluation set where
multiple transformations could be done were selected. Several combinations of
features were performed, with at most five features transferred at once. All
of these multiple transformations are zero-shot tasks, with the idea being to
test the model’s ability to use what it learned in the individual tasks
without creating different combinations of parallel corpora for fine-tuning
multiple transformations.

##  7 Evaluation

We report our evaluation on multiple components: clustering to define code
style and style transfer using sequence to sequence generation.

###  7.1 Code Style Clustering

For clustering, we evaluated the internal and external validity of the
clusters.

####  7.1.1 Metrics

##### Davies-Boules Index

For the internal metric, the Davies-Boules Index(DBI)(Davies and Bouldin,
1979) was selected; DBI measures the density of a cluster by calculating
similarity across its containing data points. The score is considered better
for the smaller value of DBI.

##### Purity

For the external metric, we used Purity(Manning, 2008) to measure how well the
clusters mapped to authorship of code scripts, with the assumption that a
given author or a given organization will have the same style. Purity
evaluates if a cluster can be recognized as a single class.

####  7.1.2 Qualitative Analysis

##### Hyperparameter Tuning

Evaluation and Setup of our clusterer |
---|---
Model | Hyperparameters | Evaluation
| Min Cluster Size | Min Samples | Cluster Epsilon | Purity | DBI
Baseline | - | - | - | 0.246 | -
Our clusterer | 400 | 100 | 0.01 | 0.404 | 1.06
Table 2: The performance and setup for the final selected clusterer.

For the results of the tuning experiments of HDBScan, we tuned three
hyperparameters, which are the minimum of each cluster size(Min Cluster Size),
minimum samples for calculating the distance to the centroid(Min Samples), and
Cluster Epsilon for the vague of the boundary between each cluster. The best
setup gave the best performance on both Purity and DBI. For Purity, the table
2 shows that our clustering model performed a significant improvement from the
random sampling baseline, randomly assigning data points to clusters and
calculating the Purity with author labels.

Figure 3: t-SNE(Van der Maaten and Hinton, 2008) visualization of the final
clusters

Given the selected hyperparameters, we generated the final clusters from the
entire dataset. Figure 3 shows the t-SNE visualization of the data metrics
with colors representing the cluster of the data point. The light green points
indicate the data detected as outliers, around 741K examples. Other than
those, around 234K examples were clustered out of 975K.

##### Pair Authorship Analysis

Cluster Authorship Pair Comparison
---
Authorship 1 | Authorship 2 | Purity | DBI | Purity-Improve | DBI-Improve
anhstudios | plotly | 0.945 | 0.845 | 0.232 | 579.528
anhstudios | AppScale | 0.764 | 0.838 | 0.036 | 75.124
anhstudios | dropbox | 0.884 | 1.125 | 0.029 | 88.904
anhstudios | mozilla | 0.859 | 1.003 | 0.020 | 17.735
anhstudios | google | 0.767 | 1.082 | 0.015 | 32.572
anhstudios | openstack | 0.635 | 0.881 | 0.015 | 94.766
Table 3: The comparison on authorship pair for the quality of the clusters.

Furthermore, we compared the pairwise authorship across all clusters to
evaluate the most distinctive coding style when assuming each author
represents a style. The score will only consider the data points that belong
to the pair of authors(authorship $1$ and authorship $2$). The improved score
indicates how much it is improved from the random sampling baseline, randomly
picking data points and assigning them to the clusters. The baseline score can
be high simply because of the high authorship distribution. So the improvement
score is considered for targeting the most valuable comparison by the HDBScan
results. The table 3 shows the top comparing evaluation results. It suggests
the anhstudios and plotly is having the most distinctive code style from
others. Figure 6 shows an example of how their coding styles are different
from each other.

##### Traditional Machine Learning models

Performance on Traditional ML Classifier
---
Classifier | Accuracy | Precision | Recall | F1
SVM | 0.777 | 0.796 | 0.777 | 0.784
Logistic Regression | 0.736 | 0.798 | 0.736 | 0.752
Decision Tree | 0.980 | 0.980 | 0.980 | 0.980
Random Forest | 0.989 | 0.989 | 0.989 | 0.989
Naive Bayes | 0.755 | 0.571 | 0.755 | 0.650
Table 4: This is the evaluation for testing how well the quality of clusters
is by fitting the cluster labels and data metrics with various traditional
machine learning classifiers.

We also evaluated the clusters using a traditional classification algorithm
with F-1 score, precision, recall, and accuracy. This evaluation tested the
quality of clusters based on whether simple classifiers can map data points
based on the same features we used for clustering into proposed clusters.
Table 4 shows the performance score for measuring the predictability of the
cluster for five different classifiers. The Naive Bayes was not performing
well, randomly guessing the data points to the clusters. The logistic
regression and SVM classifiers did relatively better but failed to recognize
many cases. However, the tree-based classifiers performed exceptionally well,
especially for the random forest classifier, which reached almost 99% score
for each metric. The reason depends on how the HDBScan splits the data points
using hierarchical tree parsing. The results indicate that our clustering
results perform good predictability on simple classifiers, but how will
neural-based models do even without directly looking at the data features on
which the clusters are based?

##### Pre-trained models Cluster Classification

Performance for Pre-trained Code Language Model
---
Model | Accuracy | Precision | Recall | F1
PLBart-Base | 0.8 | 0.8 | 0.8 | 0.8
PLBart-Multi-task-Python | 0.79 | 0.79 | 0.79 | 0.79
CodeT5-Small | 0.919 | 0.92 | 0.92 | 0.92
CodeT5-Base | 0.92 | 0.92 | 0.92 | 0.92
Table 5: The classification report for the pre-trained language model fine-
tuned on cluster classification.

To ensure that pre-trained models could identify aspects of code style, we
used F-1, precision, recall, and accuracy to evaluate if the pre-trained code
language models, PLBart and CodeT5, could classify the clusters based purely
on an input of source code. The table 5 shows the classification performance
of different checkpoints on these two architectures. CodeT5 performed a lot
better results than PLBART on classifying the style clusters. The small
checkpoint and the base checkpoint only had a very marginal difference in
their scores, which means the compact version of CodeT5 was enough to capture
the code style we defined with clusters. The better scores could contribute to
the pre-training methods of CodeT5, which target token-level prediction,
especially for identifiers. Identifiers defined by users contain much more
style information about the code than other keyword tokens. It also has the
pre-training task of bi-modal dual generation, which will understand the
alignment between the code and the documentation. Overall, CodeT5 can capture
the code style according to our proposed clusters even without seeing style
features. However, will CodeT5 be able to perform style transfer?

###  7.2 Code style transfer using Seq2Seq Generation

We used the following metrics to evaluate sequence to sequence generation.

####  7.2.1 Metrics

##### CodeBLEU

The sequence-to-sequence generation task is evaluated with CodeBLEU(Ren et
al., 2020), which is considered better than BLEU for code because it also
compares code syntax and semantics. It leverages keywords of programming
languages, the abstract syntax tree structure, and the semantics of programs.
The final CodeBLEU score is an average value over BLEU score (N-gram),
weighted BLEU score (Weighted N-gram), the syntax score (Syntax Match), and
the semantic score (Dataflow Match).

##### BLEU

For natural language tasks such as docstring transfer and comment transfer,
evaluating them with the standard BLEU score(Papineni et al., 2002) is more
appropriate.

##### DiffBLEU

Since the CodeBLEU is too optimistic when evaluating the transfer task, we
proposed DiffBLEU, which is BLEU only on the difference between reference code
$X$ and prediction code $\hat{Y}$ and the difference between input code $X$
and label code $Y$. CodeBLEU sometimes shows a high score which simply
contributes to the successfully reconstructed code from the input other than
the part that should be transferred. Especially for list comprehension
transfer, the portion that needs to be modified by the model is significantly
small. The CodeBLEU can still show high performance even without attempting
any transfer. DiffBLEU can evaluate only the difference in the code that
should be performed in the transfer task and applies to all transfer tasks.
For natural language tasks such as docstring transfer and comment transfer,
evaluating them with the standard BLEU score(Papineni et al., 2002) is more
appropriate.

$\textit{DiffBLEU}=\textit{BLEU}(\textit{Diff}(X,\hat{Y}),\textit{Diff}(X,Y))$
(1)

The equation 1 shows the concept of the DiffBLEU calculation. The method of
extracting difference, Diff, is using the “SequenceMatcher” function from the
built-in library “difflib”. However, the DiffBLEU is not precise in capturing
the difference when multiple positions and a large portion of change happen.
DiffBLEU can only serve as a proxy for evaluating relative performance on the
transfer part and will need a more precise metric, such as a parsed-tree-based
difference score.

##### Parsability

We have also evaluated the validity of the generated code with the parsability
metric. Parsability is the accuracy of whether the code is still parsable from
the AST module in Python3. Any error occurred during parsing will be
considered as failing to parse.

####  7.2.2 Individual Style Transfer

Overall Performance for Seq2Seq Generation on Individual Style Transfer
---
Task/Model | CodeBLEU | BLEU-NL | DiffBLEU | Parsability
Comment Transfer
Baseline | 0.031 | 0 | 2.20E-82 | 0.000
Individual Fine-tuned | 0.792 | 0.14 | 0.283 | 0.91
Combined Fine-tuned | 0.779 | 0.049 | 0.23 | 0.896
Class Transfer
Baseline | 0.028 | - | 1.52E-234 | 0.000
Individual Fine-tuned | 0.848 | - | 0.45 | 0.98
Combined Fine-tuned | 0.807 | - | 0.37 | 0.918
Docstring Transfer
Baseline | 0.027 | 0 | 3.43E-234 | 0.000
Individual Fine-tuned | 0.656 | 0 | 0.008 | 0.847
Combined Fine-tuned | 0.715 | 0 | 0.002 | 0.843
Casing Transfer
Baseline | 0.026 | - | 2.53E-235 | 0.001
Individual Fine-tuned | 0.967 | - | 0.728 | 0.947
Combined Fine-tuned | 0.894 | - | 0.382 | 0.925
List Comprehension Transfer
Baseline | 0.021 | - | 0.000 | 0.000
Individual Fine-tuned | 0.982 | - | 0.862 | 0.869
Combined Fine-tuned | 0.924 | - | 0.499 | 0.793
Table 6: The overall performance table on individual feature transfer.

According to the table 6, the baseline was reasonably performing scores close
to zero since it is hard to transfer style without any pre-trained
understanding of the programming language. The scores show that almost all
individual fine-tuned models were performing the best to maintain generated
code quality and perform transfer tasks.(See figure 7 for example.) The
combined fine-tuned model was trained with multi-task learning, giving noises
affecting the individual transfer results.

However, for the docstring transfer, the combined model got better performance
on CodeBLEU. One of the analyses is the combined model learned from a greater
variety of training data than just the parallel corpus on docstring transfer.
So its ability of code reconstruction can perform on a more extensive
distribution of code compared to the individual model.

##### Manual Evaluation on Docstring Transfer

Manual Evaluation on Docstring Transfer
---
| Sensible | Meaning
Task/Model | Avg. | Agreement | Avg. | Agreement
Docstring Transfer
Individual Fine-tuned | 0.848 | 0.848 | 1.283 | 0.435
Table 7: The manual evaluation on sensible and meaning of docstrings. The
agreement score represent the percentage of how two annotators agree to each
other.

The manual evaluation was necessary for some natural language generation tasks
such as docstring/comments transfer. Automatic scores for the natural language
can mislead the true performance by giving a low score without considering the
meaning of the generated text and the Sensibility of the model. So we
evaluated the docstring transfer with manual evaluation on two metrics:
Sensibility and Meaning.

Meaning metric was a $1$ to $3$ scale depending on how well the meaning of the
predicted docstrings matched the reference docstrings. The table 9 shows the
description of each scale. If a docstring did not have a corresponding
docstring in the ground truth, it was considered as not having a similar
meaning.

Sensibility was measured on a $0$ or $1$ basis to see if the docstrings
produced for each code snippet are grammatical and not logical nonsense.

The manual evaluation(table 7) indicates that the model can generate
docstrings by sensing the contents of the code. However, the produced meanings
are too simple and do not match well with the targets.

####  7.2.3 Multiple Style Transfer

Overall Performance for Seq2Seq Generation on Multiple Style Transfer
---
Task/Model | CodeBLEU | BLEU-Comment | BLEU-Docstring | DiffBLEU | Parsability
Comment + Docstring
Baseline | 0.035 | 0 | 0 | 4.53E-82 | 0
Combined Fine-tuned | 0.560 | 0.003 | 0.005 | 0.003 | 0.839
Casing + Class
Baseline | 0.026 | - | - | 1.51E-234 | 0
Combined Fine-tuned | 0.669 | - | - | 0.097 | 0.894
List Comp + Casing + Class
Baseline | 0.019 | - | - | 0 | 0
Combined Fine-tuned | 0.667 | - | - | 0.065 | 0.827
List Comp + Casing + Class + Docstring + Comment
Baseline | 0.025 | 0 | 0 | 1.84E-157 | 0
Combined Fine-tuned | 0.418 | 0.022 | 0.002 | 0.002 | 0.929
Table 8: The overall performance table on multiple feature transfer.

According to the table 8, the combined model performed significantly better
than the baseline. However, the scores are performing in a lot worse range
compared to the individual transfer. Even though the model can still generate
the grammatical output, many predictions are the same as the input, not
performing any transfer. Some only performed a subset of the transfer instead
of the full combined transfer. However, there were also some successful
cases(see figure 8) in which a full combined transfer was attempted. The
observation indicates that the CodeT5-small revealed the capability of
multiple code style transfers at once.

##  8 Conclusion

In this work, we performed a novel task - code style transfer with individual
and multiple styles. We defined code styles through clustering and found that
the style clusters reached decent quality and predictability. We explored pre-
trained code language models, CodeT5 and PLBART, and found CodeT5 has the
solid capability for capturing code style. We generated parallel corpora for
our proposed tasks and fine-tuned CodeT5. The results strongly suggest that
such a model can successfully learn the notion of code style and perform
multiple style transfers at once.

##### Future Work

For further experiments, GAN training with the combined fine-tuned model as
the generator will be an interesting test. Moreover, we can explore more
powerful pre-trained language models such as
BLOOM111https://huggingface.co/bigscience/bloom and Codex(Chen et al., 2021)
with prompting to examine whether and how they recognize the code styles.

## References

  * Ahmad et al. (2021) Wasi Uddin Ahmad, Saikat Chakraborty, Baishakhi Ray, and Kai-Wei Chang. 2021.  Unified pre-training for program understanding and generation.  _CoRR_ , abs/2103.06333.
  * Chen et al. (2021) Mark Chen, Jerry Tworek, Heewoo Jun, Qiming Yuan, Henrique Ponde de Oliveira Pinto, Jared Kaplan, Harrison Edwards, Yuri Burda, Nicholas Joseph, Greg Brockman, Alex Ray, Raul Puri, Gretchen Krueger, Michael Petrov, Heidy Khlaaf, Girish Sastry, Pamela Mishkin, Brooke Chan, Scott Gray, Nick Ryder, Mikhail Pavlov, Alethea Power, Lukasz Kaiser, Mohammad Bavarian, Clemens Winter, Philippe Tillet, Felipe Petroski Such, Dave Cummings, Matthias Plappert, Fotios Chantzis, Elizabeth Barnes, Ariel Herbert-Voss, William Hebgen Guss, Alex Nichol, Alex Paino, Nikolas Tezak, Jie Tang, Igor Babuschkin, Suchir Balaji, Shantanu Jain, William Saunders, Christopher Hesse, Andrew N. Carr, Jan Leike, Joshua Achiam, Vedant Misra, Evan Morikawa, Alec Radford, Matthew Knight, Miles Brundage, Mira Murati, Katie Mayer, Peter Welinder, Bob McGrew, Dario Amodei, Sam McCandlish, Ilya Sutskever, and Wojciech Zaremba. 2021.  Evaluating large language models trained on code.  _CoRR_ , abs/2107.03374.
  * Davies and Bouldin (1979) David L Davies and Donald W Bouldin. 1979.  A cluster separation measure.  _IEEE transactions on pattern analysis and machine intelligence_ , (2):224–227.
  * Feng et al. (2020) Zhangyin Feng, Daya Guo, Duyu Tang, Nan Duan, Xiaocheng Feng, Ming Gong, Linjun Shou, Bing Qin, Ting Liu, Daxin Jiang, and Ming Zhou. 2020.  Codebert: A pre-trained model for programming and natural languages.
  * Jin et al. (2020) Di Jin, Zhijing Jin, Zhiting Hu, Olga Vechtomova, and Rada Mihalcea. 2020.  Deep learning for text style transfer: A survey.
  * Kanade et al. (2020) Aditya Kanade, Petros Maniatis, Gogul Balakrishnan, and Kensen Shi. 2020.  Learning and evaluating contextual embedding of source code.  In _International Conference on Machine Learning_ , pages 5110–5121. PMLR.
  * Van der Maaten and Hinton (2008) Laurens Van der Maaten and Geoffrey Hinton. 2008.  Visualizing data using t-sne.  _Journal of machine learning research_ , 9(11).
  * Malzer and Baum (2020) Claudia Malzer and Marcus Baum. 2020.  A hybrid approach to hierarchical density-based cluster selection.  In _2020 IEEE International Conference on Multisensor Fusion and Integration for Intelligent Systems (MFI)_ , pages 223–228. IEEE.
  * Manning (2008) Christopher D Manning. 2008.  _Introduction to information retrieval_.  Syngress Publishing,.
  * Papineni et al. (2002) Kishore Papineni, Salim Roukos, Todd Ward, and Wei-Jing Zhu. 2002.  Bleu: a method for automatic evaluation of machine translation.  In _Proceedings of the 40th Annual Meeting of the Association for Computational Linguistics_ , pages 311–318, Philadelphia, Pennsylvania, USA. Association for Computational Linguistics.
  * Raychev et al. (2016) Veselin Raychev, Pavol Bielik, and Martin Vechev. 2016.  Probabilistic model for code with decision trees.  _SIGPLAN Not._ , 51(10):731–747.
  * Reif et al. (2021) Emily Reif, Daphne Ippolito, Ann Yuan, Andy Coenen, Chris Callison-Burch, and Jason Wei. 2021.  A recipe for arbitrary text style transfer with large language models.
  * Ren et al. (2020) Shuo Ren, Daya Guo, Shuai Lu, Long Zhou, Shujie Liu, Duyu Tang, Neel Sundaresan, Ming Zhou, Ambrosio Blanco, and Shuai Ma. 2020.  Codebleu: a method for automatic evaluation of code synthesis.  _arXiv preprint arXiv:2009.10297_.
  * Riley et al. (2021) Parker Riley, Noah Constant, Mandy Guo, Girish Kumar, David Uthus, and Zarana Parekh. 2021.  Text{settr}: Label-free text style extraction and tunable targeted restyling.
  * Szafraniec et al. (2022) Marc Szafraniec, Baptiste Roziere, Hugh Leather, Francois Charton, Patrick Labatut, and Gabriel Synnaeve. 2022.  Code translation with compiler representations.
  * Tikhonov et al. (2019) Alexey Tikhonov, Viacheslav Shibaev, Aleksander Nagaev, Aigul Nugmanova, and Ivan P. Yamshchikov. 2019.  Style transfer for texts: Retrain, report errors, compare with rewrites.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_. Association for Computational Linguistics.
  * Wang et al. (2021) Yue Wang, Weishi Wang, Shafiq Joty, and Steven C. H. Hoi. 2021.  Codet5: Identifier-aware unified pre-trained encoder-decoder models for code understanding and generation.
  * Yun Peng (2021) Mingzhe Hu Yun Peng, Yu Zhang. 2021.  An empirical study for common language features used in python projects.

Figure 4: Example individual transfer prompt Figure 5: Example multi-transfer
prompt Figure 6: Example of a pair of scripts with distinctive code style: the
code from plotly, which is used as the component for visualization tools,
tends to be object-oriented designed; on the other hand, the code from
anhstudio, which is used for game server, is having the style of single
execution script. They are having difference on class, casing, and docstring
usages. Meaning Scale Table
---
Scale | Description
1 | No predicted docstrings have any similar meaning to label docstrings.
2 | Some predicted docstrings have similar meaning to labels.
3 | More than half predicted docstrings have similar meaning to labels.
Table 9: Scale Table for Meaning Metric Figure 7: List comprehension transfer
examples: the red text is the original text and the green text is the changed
part. Figure 8: List comprehension + Casing + Class transfer examples: the red
text is the original text and the green text is the changed part.


# Single chain differential evolution Monte-Carlo for self-tuning Bayesian
inference

Willem Bonnaffé1

1\. Big Data Institute, University of Oxford, Old Road Campus, Oxford OX3 7LF

Abstract

1\. Bayesian inference is difficult because it often requires time consuming
tuning of samplers. Differential evolution Monte-Carlo (DEMC) is a self-tuning
multi-chain sampling approach which requires minimal input from the operator
as samples are obtained by taking the difference of the current position of
multiple randomly selected chains. However, this can also make DEMC more
computationally intensive than single chain samplers. 2\. We provide a single-
chain adaptation of the DEMC algorithm by taking samples according to the
difference in previous states of the chain, rather than the current state of
multiple chains. This minimises computational costs by requiring only one
posterior evaluation per step, while retaining the self-adaptive property of
DEMC. We test the algorithm by sampling a bivariate normal distribution and by
estimating the posterior distribution of parameters of an ODE model fitted to
an artificial prey-predator time series. In both cases we compare the quality
of DEMC generated samples to those obtained by a standard adaptive Markov
chain Monte-Carlo sampler (AMC). 3\. In both case studies, DEMC is as accurate
as AMC in estimating posterior distributions, while being an order of
magnitude faster due to simpler computations. DEMC also provides a higher
effective samples size than AMC, and lower initial samples autocorrelations.
4\. Its low computational cost and self-adaptive property make single chain
DEMC particularly suitable for fitting models that are costly to evaluate,
such as ODE models. The simplicity of the algorithm also makes it easy to
implement in base R, hence offering a simple alternative to STAN.

Keywords: MCMC; DEMC; DREAM; Differential evolution; Bayesian inference; ODE;
Fitting; Self-adaptive;

Emails<EMAIL_ADDRESS>
##  1 Introduction

Bayesian inference differs from frequentist inference in that it provides
parameter distributions instead of point estimates, which allow for a thorough
and intuitive quantification of uncertainty. This also makes Bayesian models
more computationally intensive to fit, as they require numerical sampling
schemes to estimate the posterior density distribution of model parameters.
However, modern computers have greatly reduced the computational costs of
fitting Bayesian models. As a result, Bayesian inference has gained popularity
in ecology over the past decade.

Gradient-free Markov chain Monte-Carlo (MCMC) is the main family of methods
for estimating parameter distributions. These methods do not rely on the
calculation of the derivative of the posterior distribution with respect to
model parameters. The chain moves in the parameter space according to a jump
kernel, also known as proposal distribution ([10]). In most cases, the
proposal distribution is defined as a multivariate Gaussian distribution,
which handles covariations between parameters ([10]). The new parameter vector
is accepted with a probability equal to the Metropolis-Hastings (MH) ratio,
namely the ratio of the log posterior density of the proposed and current
parameter vectors ([10]).

Gradient-free MCMC implementations have evolved considerably over the past two
decades. Adaptive Monte-Carlo (AMC) introduced adaptation of the jump kernel
to better match the shape of the target distribution ([9]). In the case of a
Gaussian jump kernel, this amounts to updating the variance-covariance matrix
with previous samples. Later implementations introduced delayed rejections,
where rejected samples are temporarily added to the chains, as a way to
encourage further explorations of the parameter space ([8]). Several
approaches were also developed to approximate the target distribution, either
with PCA ([11]), exponential family functions ([18]), linear, quadratic, or
Gaussian processes ([6]), in order to move more efficiently in the parameter
space and reduce computational costs. More recently, the no-U-turn sampler
(NUTS), has gained popularity among ecologists, due to its capacity to handle
high dimensional and hierarchical models ([13]). Yet, these algorithms are not
equally well-suited to solve different ecological problems ([14]), which often
requires a trial and error search for a suitable implementation.

A further limitation of these methods is that the quality of the sampling is
contingent on an adequate tuning of the hyperparameters of the jump kernel,
which is a difficult and time consuming task. Ter Braak et al. developed a
self-tuning approach, differential evolution Monte-Carlo (DEMC, [4]), also
known as DREAM (differential evolution adaptive Monte-Carlo, [5]), by
combining differential evolution (DE), a powerful genetic algorithm for global
optimisation, and MCMC sampling. Multiple chains explore the parameter space
by choosing new samples via a linear combination of the current state of the
other chains ([4]). To improve chain mixing, Ter Braak et al. further allowed
new samples to be taken according to previous samples, in addition to the
current state of the chains ([5]). The great strength of the approach is that
very little input from the operator is required to perform the inference, as
the algorithm itself tunes the jump kernel ([4]).

Though DEMC is computationally competitive compared to MCMC, the need for
multiple chains still increases computational cost. We hence propose a single
chain adaptation of the DEMC algorithm. Using a single chain means that we
only need to evaluate the posterior once per step, which makes the approach
computationally inexpensive. The shape of the posterior is learned by DE from
previous states of the chain, with no gradient nor tuning required from the
operator. This makes the approach especially suitable for sampling costly
mathematically intractable posterior distributions, such as that of ODE model
parameters.

In spite of successful implementations in other fields (e.g. [19, 20, 17]),
DEMC remains seldom used by ecologists (but see [12]). STAN and JAGS have
increasingly been used to perform Bayesian inference of ecological models
([13]). Yet, they require a specific language, which makes it hard to
customise models and interface with R ([13]). DEMC is a simple and robust
sampling algorithm which can be implemented in only a few lines of code in
base R, hence easing Bayesian inference. Therefore, we believe it useful to
reiterate the main strengths of DEMC in a context that may appeal to
ecologists.

In this paper, we first introduce a single chain adaptation of the DEMC
algorithm. We provide a self-standing and ready to use implementation of the
algorithm in R, Rcpp, and python. In a first case study, we assess how well
the sampler estimates the density of a bivariate normal distribution, for
which true samples are known. In a second case study, we provide a typical
example of population dynamics analysis, by fitting an ODE model to prey-
predator oscillations generated by the Lotka-Volterra system. In both cases
studies, we benchmark the approach against a standard AMC sampler. We show
that single chain DEMC is faster, and has a larger effective sample size
(ESS), than AMC. Our work establishes that single chain DEMC is a simple,
fast, and robust sampling algorithm, which provides ecologists who wish to
perform inference in base R with a simple, self-standing, alternative to STAN,
or JAGS.

##  2 Method

Population DEMC

Differential evolution Monte-Carlo (DEMC) is an ensemble sampling scheme
([4]). The ensemble consists of a population of chains that interact with each
other to explore the target distribution. In DEMC, the steps are performed by
differential evolution (DE), where a given chain uses the position of other
randomly selected chains to move more efficiently than a random walker. More
specifically, the step of a chain is proportional to the difference between
the position of two other randomly selected chains in the parameter space,
plus an error term. The error term encourages the chain to jump to unexplored
areas of the parameter space. DEMC hence does not require the specification
and tuning of a proposal distribution, as the shape of the target is
implicitly approximated by DE the jump kernel. This relieves the operator from
the need to tune the algorithm. However, DEMC requires at least 3 chains, and
preferably more ([5]), which can make it more computationally intensive than
simple adaptive MCMC.

Single chain DEMC

We propose here a single chain adaptation of the population DEMC algorithm.
The motivation behind this is that a single chain algorithm requires only a
single evaluation of the posterior per step, which is useful for
computationally intensive problems. In addition, this makes the algorithm
extremely simple, and hence easy to implement in any programming language. The
main difference between the population algorithm and the single chain
algorithm is that instead of using other chains to perform a jump, we use
previous states of the chain:

$\theta_{p}=\theta^{(k)}+\gamma\left(\theta^{(u)}-\theta^{(v)}\right)+\varepsilon,$
(1)

where $\theta_{p}$ is the proposed state, $\theta^{(k)}$ the current state,
$\theta^{(u)}$ and $\theta^{(v)}$ are the $u^{th}$ and $v^{th}$ past states of
the chain, and $\gamma$ and $\varepsilon$ are tuning parameters. As with the
population algorithm, $\gamma=2.38/\sqrt{2d}$ and
$\varepsilon\sim\mathcal{U}(-\delta,\delta)$, where $\delta=0.001$ should
yield optimal results ([4]).

Algorithm

The following algorithm can be used to sample the target distribution
$\pi(.)$:

  * •

Initialisation $k=0$

    * –

Set $\theta^{(0)}$ such that $\pi(\theta^{(0)})>0$

    * –

Set $\gamma=2.38/\sqrt{2d}$ and $\delta=0.001$

  * •

For each iteration $0<k\leq k_{max}$:

    * –

Draw $\theta^{(u)}$ and $\theta^{(v)}$ with $u\sim\mathcal{U}(0,k-1)$ and
$v\sim\mathcal{U}(0,k-1)$

    * –

Draw $\varepsilon\sim\mathcal{U}(-\delta,\delta)$

    * –

Compute
$\theta_{p}=\theta^{(k)}+\gamma\left(\theta^{(u)}-\theta^{(v)}\right)+\varepsilon$

    * –

Draw $q\sim\mathcal{U}(0,1)$

    * –

If $\pi(\theta_{p})/\pi(\theta^{(k)})>q$ then set $\theta^{(k+1)}=\theta_{p}$

    * –

Else set $\theta^{(k+1)}=\theta^{(k)}$

Detailed balance

It is necessary to show that single chain DEMC retains the main properties of
Monte-Carlo (MC) methods, namely that detailed balance and ergodicity are
preserved. Given the similarity of the population and single chain DEMC
algorithms, we can use the proofs used by Ter Braak in his original paper
([4]).

Detailed balance can be proven by showing that the transition kernel is
symmetrical, which implies that each step is reversible. The probability
density of the proposed value conditional on the current and past states of
the chain is

$q(\theta_{p}|\theta^{(k)},\theta^{(u)},\theta^{(v)})=f_{U}(\theta_{p}|\theta^{(k)}+\gamma(\theta^{(u)}-\theta^{(v)}),\delta),$
(2)

where $q(\theta_{p}|.)$ is the transition kernel to $\theta_{p}$ given
$\theta^{(k)}$, $\theta^{(u)}$, $\theta^{(v)}$, and $f_{U}$ the probability
density function of a uniform distribution centred on
$\theta^{(k)}+\gamma(\theta^{(u)}-\theta^{(v)})$ with spread $\delta$.

The total probability of jumping from the current to the proposed value is
obtained by summing the transition kernel over all $k-1$ past states of the
chain

$q(\theta_{p}|\theta^{(k)})=\sum_{u}^{k-1}\sum_{v\neq
u}^{k-1}f_{U}(\theta_{p}|\theta^{(k)}+\gamma(\theta^{(u)}-\theta^{(v)}),\delta).$
(3)

Similarly the total probability of the reverse jump is

$\displaystyle q(\theta^{(k)}|\theta_{p})$
$\displaystyle=\sum_{u}^{k-1}\sum_{v\neq
u}^{k-1}f_{U}(\theta^{(k)}|\theta_{p}-\gamma(\theta^{(u)}-\theta^{(v)}),\delta).$
(4)

By noting that setting $u=v$ and $v=u$ yields a mathematically identical
transition kernel we show that
$q(\theta_{p}|\theta^{(k)})=q(\theta^{(k)}|\theta_{p})$ which satisfies
detailed balance.

Ergodicity

As with any adaptive schemes, ergodicity is more difficult to establish given
that the proposal distribution changes with each iteration ([15]). Ergodicity
can be demonstrated by showing that the chain is not transient, aperiodic,
irreducible ([4]). Nontransiency is guaranteed if the algorithm converges to a
stationary distribution, which requires a stationary proposal distribution.
For adaptive schemes this condition is achieved by progressively reducing the
adaptation ([9]). In practice this is done by sampling from an increasing part
of the history of the chain (e.g. by taking samples from iteration $k/2$ to
$k$, $k$ being the current iteration, [9]). Aperiodicity requires that the
chain can return to its current state in a finite amount of steps, while
irreducibility requires that any state can be reached within a finite amount
of steps. Both conditions are guaranteed here due to the error term
$\varepsilon$, which generates a random walk behaviour, thus theoretically
allowing for the exploration of the entire parameter space ([4]).

##  3 Case studies

In order to assess the reliability and efficiency of our single chain DEMC
algorithm we perform two tests. The first test consists in estimating the
density of a bivariate correlated normal distribution. The second test is to
recover the probability density of the parameters of a differential equation
system. In both tests, the algorithm is compared to a standard adaptive
Metropolis-Hastings algorithm (AMC).

Bivariate normal distribution

We test the capacity of our algorithm to recover the probability density of a
highly correlated bivariate normal distribution. The target distribution is
given by

$\pi(\theta)=(2\pi)^{-1}det(\Sigma)^{-\frac{1}{2}}\exp\left(-\frac{1}{2}\theta^{T}\Sigma^{-1}\theta\right),$
(5)

where $\pi(\theta)$ refers to the posterior probability density of the
parameter vector $\theta=[\theta_{1},\theta_{2}]^{T}$, $\Sigma$ is the
covariance matrix, the diagonal elements of which are
$\sigma_{1}=\sigma_{2}=1$ and the correlation coefficient is set to
$\rho=0.99$.

We run three AMC chains of 10000 iterations with a Gaussian transition kernel.
The scale parameter is set to $\alpha=2.38/\sqrt{4}$ and the covariance matrix
of the transition kernel is initially set to a two dimensional identity matrix
$\Sigma_{p}=I_{2}$. Initially, 250 samples are taken without adaptation. A
further 10000 samples are taken with adaptation of the covariance matrix at
every step using an increasing part of the chain (i.e. samples from $k/2$ to
$k$, $k$ being the current iteration). We also run three DEMC chains of 10000
iterations. The scale parameter is set to $\gamma=2.38/\sqrt{4}$ and the noise
parameter is set to $\delta=0.001$. Both DEMC and AMC chains are initiated at
random starting points
$\theta_{i}^{(0)}\overset{iid}{\sim}\mathcal{N}(0,1),~{}~{}i=\\{1,2\\}$.

Differential equation system

We test the capacity of single chain DEMC to recover the probability density
of parameters of a differential equation system. Here we use the Lotka-Voltera
differential equation system, which consists in a set of two coupled
differential equations

$\left\\{\begin{aligned} \frac{dX}{dt}&=\alpha X-\beta XY\\\
\frac{dY}{dt}&=\gamma XY-\delta Y\\\ \end{aligned}\right.$ (6)

where the dynamics of the system $dZ/dt=[dX/dt,dY/dt]^{T}$ are governed by the
state vector is $Z=[X,Y]^{T}$, $Z\in\mathbb{R}_{+}^{2}$ and the parameter
vector is $\theta=[\alpha,\beta,\gamma,\delta]^{T}$,
$\theta\in\mathbb{R}_{+}^{4}$.

We use this system to simulate an artificial time series that we fit the model
to. To do this we set
$\theta_{true}=[\alpha=1,\beta=0.1,\gamma=0.1,\delta=1]^{T}$ and integrate the
system on $t\in[0,10]$. Then we simply sample the predictions of the model
every time step $\Delta t=1$.

We define a simple Bayesian model to fit the differential equation system

$\pi(\theta|Z)\propto\prod_{k}p(Z_{k}|\bar{Z}(k|\theta))p(\theta),$ (7)

where $Z_{k}$ is the observed state and
$\bar{Z}(t|\theta)=Z_{0}+\int_{0}^{k}d\bar{Z}$ the expected state of the
system at the observed time step $k=\\{0,\ldots,10\\}$. The likelihood is
defined as a gamma distribution around the expected state of the system
$Z_{k}\overset{iid}{\sim}\Gamma(\bar{Z}(k|\theta),1),~{}\forall~{}k\in{0,\ldots,10}$.
The prior distribution of the parameter vector is defined as a multivariate
uniform distribution $\theta\sim\mathcal{U}_{4}(0,2)$.

We run three AMC chains of 10000 iterations. The scale parameter is set to
$\alpha=2.38/\sqrt{8}$ the covariance matrix of the transition kernel is
initiated as a multivariate identity matrix $\Sigma_{p}=I_{4}$. Initially, 250
samples are taken without adaptation. A further 10000 samples are taken with
adaptation of the covariance matrix at every step using an increasing part of
the chain (i.e. samples from $k/2$ to $k$, $k$ being the current iteration).
We also run three DEMC chains of 10000 iterations. As before the scale
parameter is set to $\gamma=2.38/\sqrt{8}$ and the noise parameter is set to
$\delta=0.001$ All chains are initiated at points randomly sampled from the
prior distribution.

Analysis and comparison of algorithms

We compare the AMC and DEMC algorithms on three levels. First, we observe the
trace of the chains to assess the rate of convergence of the algorithms.
Second, we compute the autocorrelation in the chain samples up to the
$50^{th}$ order and effective sample size ([10]). The first 1000 iterations of
the chains are removed to eliminate the transient states. Finally, we compare
the estimated quantiles of the stationary probability distributions obtained
from the two algorithms. In this case the chains are further thinned down to
1000 independent samples. Results are presented below.

##  4 Results

Bivariate normal distribution

Single chain DEMC and AMC converge to virtually identical posterior
distributions, which approximate well the shape of the bivariate normal
distribution (Figure 1). In doing so, our single chain DEMC algorithm is about
an order of magnitude faster than AMC (Table 1). The autocorrelation in the
samples of DEMC drops faster than that of AMC which remains slightly above 0
for longer (Figure 2). Samples estimated by single chain DEMC match up with
both those generated by AMC and by the true distribution (Figure 3).

Table 1:  Summary results of sampling the bivariate normal distribution with AMC and DEMC. The Dim column corresponds to the dimension of the bivariate normal, the MaP, mean and sd are the maximum a posteriori, mean and standard deviation estimated from the samples, $\hat{r}$ is the convergence statistic (as in [10]), and ESS denotes the effective sample size. The time column corresponds to the time in seconds required to run 3 chains of 10000 iterations.  Algorithm | Dim. | MaP | mean | sd | $\hat{r}$ | Time | ESS
---|---|---|---|---|---|---|---
AMC | $\theta_{1}$ | -0.0042 | -0.0105 | 1.0001 | 1.0024 | 77.163 | 186.4
| $\theta_{2}$ | -0.0087 | -0.0097 | 1 | 1.0022 | - | 183.56
DE-MC | $\theta_{1}$ | -0.0328 | -0.0013 | 0.9893 | 1.0006 | 4.007 | 223.96
| $\theta_{2}$ | -0.0332 | -0.006 | 0.9928 | 1.0005 | - | 222.47
Figure 1:  Estimated posterior distributions of the parameters of the
bivariate normal model. Comparison of the estimated posterior distributions of
parameters $\theta_{1}$ and $\theta_{2}$ of the bivariate normal model,
obtained by AMC (left) and DEMC (right). dLogPost corresponds to the log-
density of the bivariate normal distribution. V1 and V2 correspond to the
first and second parameter, $\theta_{1}$ and $\theta_{2}$, respectively. The
shade of grey, from dark to light, is proportional to the posterior density.
The histograms correspond to the marginal probability density of the
parameters. Final samples are obtained by burning and thinning the chains to
retain 1000 samples.

Figure 2:  Sample autocorrelations of the bivariate normal model. Comparison
of autocorrelations of samples generated by AMC (red) and DEMC (blue). Each
row corresponds to one of three chains. Each column corresponds to each
parameter of the bivariate Gaussian. Graph a., c., and e., correspond to the
three chains for parameter $\theta_{1}$, while graph b., d., and f. correspond
to parameter $\theta_{2}$. Autocorrelations are calculated on the chains
before thinning.

Figure 3:  Samples of the bivariate normal model. Comparison of true samples
taken from the bivariate normal distribution and samples obtained by AMC (a.
and c.) and DEMC (b. and d.) for parameter $\theta_{1}$ (a. and b.) and
$\theta_{2}$ (c. and d.). Different colours correspond to different chains.

ODE system

In the ODE case study, the chains of single chain DEMC also converge faster to
the target distribution than AMC (Table 2). The posterior distributions
estimated by AMC and DEMC are virtually identical (Figure 4). The
autocorrelation in the samples of single chain DEMC drops faster than AMC
(Figure 5). Finally, parameter samples obtained by AMC and DEMC are virtually
identical (Figure 6).

Table 2:  Summary results of fitting the ODE model with the AMC and DEMC. The Dim column corresponds to the parameters of the ODE model, the MaP, mean and sd are the maximum a posteriori, mean and standard deviation estimated from the samples, $\hat{r}$ is the convergence statistic (as in [10]), and ESS denotes the effective sample size. The time column corresponds to the time in seconds required to run 3 chains of 10000 iterations.  Algorithm | Dim. | MaP | mean | sd | $\hat{r}$ | Time | ESS
---|---|---|---|---|---|---|---
AMC | $\alpha$ | 1.0029 | 0.999 | 0.0353 | 1.002 | 389.223 | 151.31
| $\beta$ | 0.1007 | 0.1014 | 0.0105 | 1.0006 | - | 154.55
| $\gamma$ | 0.0984 | 0.1014 | 0.0134 | 1.0022 | - | 147.44
| $\delta$ | 0.9688 | 0.975 | 0.0891 | 1.0032 | - | 144.15
DE-MC | $\alpha$ | 1.0078 | 0.9972 | 0.036 | 1.0007 | 115.251 | 182.21
| $\beta$ | 0.1017 | 0.1011 | 0.0109 | 1.0006 | - | 172.34
| $\gamma$ | 0.097 | 0.1019 | 0.0136 | 1.0009 | - | 178.13
| $\delta$ | 0.9608 | 0.9743 | 0.0864 | 1 | - | 183.54
Figure 4:  Estimated posterior distributions of the parameters of the ODE
model. Comparison of the estimated posterior distributions of parameters
$\alpha$, $\beta$, $\gamma$, and $\delta$ of the Lotka-Volterra ODE system,
obtained by AMC (left) and DEMC (right). dLogPost corresponds to the log-
density of the bivariate normal distribution. V1, V2, V3, and V4 correspond to
the ODE parameters, $\alpha$, $\beta$, $\gamma$, and $\delta$, respectively.
The shade of grey, from dark to light, is proportional to the posterior
density. The histograms correspond to the marginal probability density of the
parameters. Final samples are obtained by burning and thinning the chains to
retain 1000 samples. The fit of the ODE models to the artificial time series
is shown in graph a. (AMC) and b. (DEMC). Solid lines, either red for the
prey, X, or blue for the predator, Y, correspond to the best fit, and the
shaded area is the 90% confidence interval on predictions.

Figure 5:  Sample autocorrelations of the ODE model. Comparison of
autocorrelations of samples generated by AMC (red) and DEMC (blue). Each row
corresponds to one of three chains. Each column corresponds to each parameter
of the bivariate Gaussian. Graph a., c., and e., correspond to the three
chains for parameter $\theta_{1}$, while graph b., d., and f. correspond to
parameter $\theta_{2}$. Autocorrelations are calculated on the chains before
thinning.

Figure 6: Samples of the ODE model. Comparison of samples obtained by AMC (y
axis) and DEMC (x axis) for each of the four parameters of the ODE system,
$\alpha$ (a.), $\beta$ (b.), $\gamma$ (c.), $\delta$ (d.). Different colours
correspond to different chains.

##  5 Discussion

Bayesian inference is challenging because of the need for sampling the
posterior distributions of model parameters. Differential evolution Monte-
Carlo (DEMC) is a self-adaptive sampling method, which relieves the operator
from the need of tuning the sampler. The method however requires an ensemble
of chains, which can make more computationally intensive to use than single
chain alternatives. Here we provide a single chain adaptation of DEMC, making
it computationally inexpensive and simple to implement. We show that our
method is accurate by sampling a bivariate normal distribution where true
samples are known. We also show that the method is accurate in fitting ODE
models to prey-predator time series data. In both case studies we compare
results of single chain DEMC to a standard adaptive Monte-Carlo (AMC). We find
that our method is quantitatively identical to AMC while being much faster.
Single chain DEMC hence provides a simple, accurate, and robust alternative to
performing Bayesian inference in STAN and JAGs.

Comparison of single chain DEMC to AMC, DREAM, and HMC

We find that parameter estimation by single chain DEMC and AMC are identical,
for both the bivariate normal model and ODE model. However, single chain DEMC
is much faster than AMC. This is due to the fact that tuning the jump kernel
of AMC involves computing a covariance matrix of previous samples ([9]), which
is computationally expensive. In contrast, DEMC only requires taking the
difference between two previous samples, which is inexpensive ([4]). In
addition, sample autocorrelations of DEMC are lower than AMC, which translates
into a higher effective sample size for DEMC. This may be due to the fact that
the jump kernel of DEMC has no predetermined functional form, hence making it
better at approximating different target shapes than AMC.

The original implementation of DEMC (also known as DREAM) involves multiple
chains ([4]). The new location of the chain in the parameter space is
determined by taking the difference of the current location of the other
chains, chosen at random. Using multiple chains offers several advantages over
a single chain implementation. First, it reduces the risk of getting trapped
in a local maximum. Ter Braak also noted that chain mixing increased with the
number of elements in the ensemble ([5]). However, computational cost also
increases with the number of chains. So that, even if single chain DEMC is
potentially less powerful than multiple chain DEMC, it makes up for it by
taking more samples per unit of time.

We choose not to compare HMC and DEMC because we think they serve different
purposes. HMC is only suitable for targets that are not too costly to evaluate
and for which gradients are available ([18]). DEMC is more suitable for
computationally expensive and mathematically intractable targets. Regardless,
even in cases where HMC may be more suitable, DEMC would be useful for fast
prototyping the model and running pilot chains, due to its simplicity and the
fact that it does not require gradients.

Suitability of single chain DEMC for parameter estimation in ODE models

Fitting ODE models to time series data is a powerful technique to identify the
drivers of population dynamics ([7, 16, 1, 3, 2]). For instance, ODE
parameters can describe the birth, death rate of individuals, as well as the
effect of competitors and predators in population growth. Estimating ODE
parameters is challenging because the predictions of ODE models are obtained
from numerical simulation (e.g. [3, 2]). This makes it computationally
expensive to evaluate the posterior distribution, and prevents the calculation
of gradients. Single chain DEMC is highly suitable for fitting these models.
This is because only one evaluation of the posterior is required per step,
hence reducing computational costs as much as possible. The self tuning
property of single DEMC also alleviates the need for running pilot chains and
tuning the algorithm, thus saving the operator time. We demonstrate this in
our second case study, where we fit an ODE model to a prey-predator time
series generated by a Lotka-Volterra model. Single chain DEMC provides
accurate parameter estimates compared to ground truth, while being much faster
than AMC. Single chain DEMC hence takes the hassle away from the fitting,
allowing the operator to focus on designing and comparing ODE models.

Limits and prospects

In our single chain DEMC algorithm, a new sample is taken by computing the
difference between two previous samples selected at random. This, in essence,
is a linear interpolation/extrapolation of the previous states of the chain.
While we show that this performs well on correlated/nonlinear distributions, a
nonlinear combination of previous states may perform better. To do this, we
could use nonlinear combinations of the previous states, by using polynomial
functions, or even artificial neural networks (ANN). One difficulty with this
is that is requires extra tuning parameters to control the nonlinearity of the
jump kernel. However, those could be treated as ordinary model parameters and
learned by differential evolution as the chain is running. We view neural
differential evolution Monte-Carlo as a particularly promising prospect given
the powerful approximation capacity of ANNs.

Finally, we only considered a simple implementation of single chain DEMC in R,
Rcpp, and python. We did not exploit the potential for CPU and GPU
parallelization of the algorithm. This could result in order of magnitude
faster inference.

Conclusion

In conclusion, we demonstrate that our single chain adaptation of the DEMC
algorithm is easy to use and implement, computationally efficient, and
improves sampling efficiency. The differential evolution step alleviates the
need for trial and error tuning of the sampler, allowing the operator to focus
on model design and comparison. The simplicity of the algorithm makes it easy
to implement in base R in a few lines of code, hence providing an easy
alternative to more complicated Bayesian inference tools (e.g. STAN, JAGS).
The single chain nature of the algorithm makes it computationally inexpensive
compared to multiple chain DEMC. Overall, our algorithm is particularly
suitable for Bayesian inference of computationally expensive and complicated
models, such as ODE models, that require an algorithm that minimises model
evaluations.

Statement of authorship

Willem Bonnaffé designed the method, performed the analysis, wrote the
manuscript;

Acknowledgments

We thank warmly Tim Coulson and the Ecological and Evolutionary Lab for
comments and feedback on early versions of the work. This work was supported
by the Oxford-Oxitec scholarship and the NERC DTP.

Data accessibility

All data and code is available at https://github.com/WillemBonnaffe/DEMC.

## References

  * [1] Matthew P. Adams et al.  “Informing management decisions for ecological networks, using dynamic models calibrated to noisy time-series data”  In _Ecology Letters_ 23 Blackwell Publishing Ltd, 2020, pp. 607–619  DOI: 10.1111/ele.13465
  * [2] Willem Bonnaffé and Tim Coulson  “Fast fitting of neural ordinary differential equations by Bayesian neural gradient matching to infer ecological interactions from time series data”  In _arXiv_ , 2022, pp. 1–57  URL: http://arxiv.org/abs/2209.06184
  * [3] Willem Bonnaffé, Ben C. Sheldon and Tim Coulson  “Neural ordinary differential equations for ecological and evolutionary time series analysis”  In _Methods in Ecology and Evolution_ 2, 2021, pp. 1–46  DOI: 10.1111/2041-210x.13606
  * [4] Cajo J. F. Ter Braak  “A Markov Chain Monte Carlo version of the genetic algorithm Differential Evolution: easy Bayesian computing for real parameter spaces”  In _Statistics and Computing_ 16, 2006, pp. 239–249  DOI: 10.1007/s11222-006-8769-1
  * [5] Cajo J.F. Ter Braak and Jasper A. Vrugt  “Differential Evolution Markov Chain with snooker updater and fewer chains”  In _Statistics and Computing_ 18, 2008, pp. 435–446  DOI: 10.1007/s11222-008-9104-9
  * [6] Patrick R. Conrad, Youssef M. Marzouk, Natesh S. Pillai and Aaron Smith  “Accelerating Asymptotically Exact MCMC for Computationally Intensive Models via Local Approximations”  In _Journal of the American Statistical Association_ 111 American Statistical Association, 2016, pp. 1591–1607  DOI: 10.1080/01621459.2015.1096787
  * [7] V. Demyanov, S. N. Wood and T. J. Kedwards  “Improving ecological impact assessment by statistical data synthesis using process-based models”  In _Journal of the Royal Statistical Society. Series C: Applied Statistics_ 55, 2006, pp. 41–62  DOI: 10.1111/j.1467-9876.2005.00527.x
  * [8] Heikki Haario, Marko Laine, Antonietta Mira and Eero Saksman  “DRAM: Efficient adaptive MCMC”  In _Statistics and Computing_ 16, 2006, pp. 339–354  DOI: 10.1007/s11222-006-9438-0
  * [9] Heikki Haario, Eero Saksman and Johanna Tamminen  “An Adaptive Metropolis Algorithm”  In _Bernoulli_ 7, 2001, pp. 223  DOI: 10.2307/3318737
  * [10] N Thompson Hobbs and Mevin B Hooten  “A Statistical Primer for Ecologists”  Princeton University Press, 2015  DOI: doi:10.1515/9781400866557
  * [11] David A. Kennedy, Vanja Dukic and Greg Dwyer  “Combining principal component analysis with parameter line-searches to improve the efficacy of Metropolis–Hastings MCMC”  In _Environmental and Ecological Statistics_ 22 Kluwer Academic Publishers, 2015, pp. 247–274  DOI: 10.1007/s10651-014-0297-0
  * [12] Dan Lu, Daniel Ricciuto, Anthony Walker, Cosmin Safta and William Munger  “Bayesian calibration of terrestrial ecosystem models: A study of advanced Markov chain Monte Carlo methods”  In _Biogeosciences_ 14 Copernicus GmbH, 2017, pp. 4295–4314  DOI: 10.5194/bg-14-4295-2017
  * [13] Cole C. Monnahan and Kasper Kristensen  “No-U-turn sampling for fast Bayesian inference in ADMB and TMB: Introducing the adnuts and tmbstan R packages”  In _PLoS ONE_ 13 Public Library of Science, 2018  DOI: 10.1371/journal.pone.0197954
  * [14] Lauren C. Ponisio, Perry Valpine, Nicholas Michaud and Daniel Turek  “One size does not fit all: Customizing MCMC methods for hierarchical models using NIMBLE”  In _Ecology and Evolution_ 10 John WileySons Ltd, 2020, pp. 2385–2416  DOI: 10.1002/ece3.6053
  * [15] Gareth Roberts and Jeffrey Rosenthal  “Coupling and Ergodicity of Adaptive Markov Chain Monte Carlo Algorithms”  In _Journal of Applied Probability_ 44, 2012, pp. 458–475
  * [16] Benjamin Rosenbaum, Michael Raatz, Guntram Weithoff, Gregor F. Fussmann and Ursula Gaedke  “Estimating parameters from multiple time series of population dynamics using bayesian inference”  In _Frontiers in Ecology and Evolution_ 6, 2019, pp. 1–14  DOI: 10.3389/fevo.2018.00234
  * [17] Erin M. Shockley, Jasper A. Vrugt and Carlos F. Lopez  “PyDREAM: High-dimensional parameter inference for biological models in python”  In _Bioinformatics_ 34 Oxford University Press, 2018, pp. 695–697  DOI: 10.1093/bioinformatics/btx626
  * [18] Heiko Strathmann, Dino Sejdinovic, Samuel Livingstone, Zoltan Szabo and Arthur Gretton  “Gradient-free Hamiltonian Monte Carlo with Efficient Kernel Exponential Families”  In _arXiv_ , 2015, pp. 1–9  URL: https://github.com/karlnapf/kernel_hmc
  * [19] Brandon M. Turner and Per B. Sederberg  “Approximate Bayesian computation with differential evolution”  In _Journal of Mathematical Psychology_ 56, 2012, pp. 375–385  DOI: 10.1016/j.jmp.2012.06.004
  * [20] Brandon M. Turner, Per B. Sederberg, Scott D. Brown and Mark Steyvers  “A Method for efficiently sampling from distributions with correlated dimensions”  In _Psychological Methods_ 18, 2013, pp. 368–384  DOI: 10.1037/a0032222

##  6 Supplementary

##  Appendix A Traces of AMC chains fitted to the multivariate Gaussian

Figure S1: Trace of full AMC chains sampling the bivariate normal model.
Figure S2:  Trace of downsampled AMC chains sampling the bivariate normal
distribution. Downsampling is performed by burning the first 2000 iterations
(out of 10000) and then retaining 1000 evenly spaced samples of the burnt
chain.

##  Appendix B Traces of DEMC chains fitted to the multivariate Gaussian

Figure S3: Trace of full DEMC chains sampling the bivariate normal model.
Figure S4:  Trace of downsampled DEMC chains sampling the bivariate normal
distribution. Downsampling is performed by burning the first 2000 iterations
(out of 10000) and then retaining 1000 evenly spaced samples of the burnt
chain.

##  Appendix C Traces of AMC chains fitted to the ODE system

Figure S5: Trace of full AMC chains sampling the ODE model. Figure S6:
Continued  Figure S7:  Trace of downsampled AMC chains sampling the ODE model.
Downsampling is performed by burning the first 2000 iterations (out of 10000)
and then retaining 1000 evenly spaced samples of the burnt chain.  Figure S8:
Continued

##  Appendix D Traces of DEMC chains fitted to the ODE system

Figure S9: Trace of full DEMC chains sampling the ODE model. Figure S10:
Continued Figure S11:  Trace of downsampled DEMC chains sampling the ODE
model. Downsampling is performed by burning the first 2000 iterations (out of
10000) and then retaining 1000 evenly spaced samples of the burnt chain.
Figure S12:  Continued