Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Paper: Parameter Estimation Using the Python Package pymcmcstat #466

merged 39 commits into from Jul 3, 2019
Changes from 1 commit
Show all changes
39 commits
Select commit Hold shift + click to select a range
added initial rst and bib files
prmiles May 13, 2019
added doi numbers for code and articles
prmiles May 13, 2019
made changes - successfully compiled
prmiles May 13, 2019
added additional references
prmiles May 16, 2019
added methodology section
prmiles May 16, 2019
updated line formatting
prmiles May 16, 2019
first complete rough draft
prmiles May 20, 2019
added images
prmiles May 20, 2019
final edits for first draft submission
prmiles May 21, 2019
updated reference to pymcmcstat 1.7.1
prmiles May 21, 2019
updated intro
prmiles Jun 21, 2019
reordered methodology
prmiles Jun 21, 2019
updated abstract
prmiles Jun 21, 2019
fixed typo
prmiles Jun 21, 2019
added section
prmiles Jun 21, 2019
updated joss ref
prmiles Jun 21, 2019
added description of key pieces
prmiles Jun 21, 2019
updated key components list
prmiles Jun 22, 2019
explain MCMC relation to Bayes
prmiles Jun 22, 2019
added ref for Metropolis
prmiles Jun 22, 2019
updated initial example
prmiles Jun 22, 2019
added basic visualization example
prmiles Jun 22, 2019
add pairwise correlation
prmiles Jun 23, 2019
add pc image
prmiles Jun 23, 2019
added uq intervals
prmiles Jun 23, 2019
fixed captions
prmiles Jun 23, 2019
add Metropolis section
prmiles Jun 23, 2019
removed monodomain section
prmiles Jun 23, 2019
minor update to conclusion
prmiles Jun 23, 2019
minor upd to several sections
prmiles Jun 23, 2019
fixed bug and upd case studies
prmiles Jun 23, 2019
upd visco section
prmiles Jun 24, 2019
upd rsl section
prmiles Jun 24, 2019
updated conclusions
prmiles Jun 24, 2019
add ref for mcmcplot
prmiles Jun 24, 2019
added URLs via footnotes
prmiles Jun 26, 2019
moved caption to above table
prmiles Jun 26, 2019
fixed typo in visc. section
prmiles Jun 26, 2019
minor edits
prmiles Jun 26, 2019
File filter...
Filter file types
Jump to…
Jump to file or symbol
Failed to load files and symbols.


Just for now


added initial rst and bib files

  • Loading branch information...
prmiles committed May 13, 2019
commit eeef0738e2cd83ca54574b7e9bc8db4c0ee6d827
@@ -0,0 +1,38 @@
title={An Adaptive Metropolis Algorithm},
author={Haario, Heikki and Saksman, Eero and Tamminen, Johanna and others},
publisher={Bernoulli Society for Mathematical Statistics and Probability}

title={\mbox{DRAM}: Efficient Adaptive MCMC},
author={Haario, Heikki and Laine, Marko and Mira, Antonietta and Saksman, Eero},
journal={Statistics and computing},

author = {Paul Miles},
title = {prmiles/pymcmcstat: v1.6.0},
month = aug,
year = 2018,
doi = {10.5281/zenodo.1407136},
url = {}

title={Uncertainty Quantification: Theory, Implementation, and Applications},
author={Smith, Ralph C},
@@ -0,0 +1,280 @@
:author: Paul R. Miles
:institution: Department of Mathematics
:institution: North Carolina State University

:author: Ralph C. Smith
:institution: Department of Mathematics
:institution: North Carolina State University

:bibliography: mybib

Parameter Estimation Using the Python Package pymcmcstat

.. class:: abstract

The Python package pymcmcstat \cite{pymcmcstat2018v1.6.0} provides a robust
platform for a variety of engineering inverse problems. Bayesian
statistical analysis is a powerful tool for parameter estimation,
and many algorithms exist for numerical approaches that utilize
Markov Chain Monte Carlo (MCMC) methods \cite{smith2014uncertainty}.

In pymcmcstat, the user is provided with a suite of Metropolis based
algorithms, with the primary approach being Delayed Rejection Adaptive
Metropolis (DRAM) \cite{haario2006dram, haario2001adaptive}. A simple
procedure of adding data, defining model parameters and settings, and
setting up simulation options provides the user with a wide variety of
computational tools for considering inverse problem. This approach to
inverse problems utilizes data to provide insight into model limitations
and provide accurate estimation of the underlying model and observation

As many Python packages currently exist for performing MCMC simulations,
we had several goals in developing this code. To our knowledge, no
current package contains the $n$-stage delayed rejection algorithm,
so pymcmcstat was intended to fill this gap. Furthermore, many
researchers in our community have extensive experience using the MATLAB
toolbox mcmcstat. Our implementation provides a similar user environment,
while exploiting Python structures. We hope to decrease dependence on
MATLAB in academic communities by advertising comparable tools in Python.

This package has been applied to a wide variety of engineering problems,
including radiation source localization as well as constitutive model
development of smart material systems. This is not an exhaustive listing
of scientific problems that could be analyzed using pymcmcstat, and more
details regarding the program methodology can be found via the project

Localization of special nuclear material in urban environments poses a
very important task with many challenges. Accurate representation of
radiation transport in a three-dimensional domain that includes various
forms of construction materials presents many computational challenges.
For a representative domain in Ann Arbor, Michigan we can construct
surrogate models using machine learning algorithms based on Monte Carlo
N-Particle (MCNP) simulations. The surrogate models provide a
computationally efficient approach for subsequent inverse model
calibration, where we consider the source location ($x, y, z$) as our
model parameters. We will demonstrate the viability of using pymcmcstat
for localization problems of this nature.

Many smart material systems depend on robust constitutive relations for
applications in robotics, flow control, and energy harvesting. To fully
characterize the material or system behavior, uncertainty in the model must
be accurately represented. By using experimental data in conjunction with
pymcmcstat, we can estimate the model parameter distributions and visualize
how that uncertainty propagates through the system. We will consider
specific examples in viscoelastic modeling of dielectric elastomers and
also continuum approximations of ferroelectric monodomain crystal

.. class:: keywords

Markov Chain Monte Carlo (MCMC), Delayed Rejection Adaptive Metropolis (DRAM),
Parameter Estimation, Bayesian Inference


Twelve hundred years ago |---| in a galaxy just across the hill...

Lorem ipsum dolor sit amet, consectetur adipiscing elit. Vestibulum sapien
tortor, bibendum et pretium molestie, dapibus ac ante. Nam odio orci, interdum
sit amet placerat non, molestie sed dui. Pellentesque eu quam ac mauris
tristique sodales. Fusce sodales laoreet nulla, id pellentesque risus convallis
eget. Nam id ante gravida justo eleifend semper vel ut nisi. Phasellus
adipiscing risus quis dui facilisis fermentum. Duis quis sodales neque. Aliquam
ut tellus dolor. Etiam ac elit nec risus lobortis tempus id nec erat. Morbi eu
purus enim. Integer et velit vitae arcu interdum aliquet at eget purus. Integer
quis nisi neque. Morbi ac odio et leo dignissim sodales. Pellentesque nec nibh
nulla. Donec faucibus purus leo. Nullam vel lorem eget enim blandit ultrices.
Ut urna lacus, scelerisque nec pellentesque quis, laoreet eu magna. Quisque ac
justo vitae odio tincidunt tempus at vitae tortor.

Of course, no paper would be complete without some source code. Without
highlighting, it would look like this::

def sum(a, b):
"""Sum two numbers."""

This comment has been minimized.

Copy link

nudro Jun 11, 2019

Although I appreciate the transparency via the maths you provide, I'm concerned that this will be misconstrued as "technical jargon" unless the goal and utility for why you lay out the functions are relevant - particularly to emphasize the goal of how approximating these functions directly are intractable (as you mentioned the quadtature)... and therefore are better suited for mcmc - enter pymcmcstat and metropolis algorithms.

This comment has been minimized.

Copy link

prmiles Jun 24, 2019

Author Contributor

Thank you for pointing this out. I have broken the Methodology section into two subsections: 1) "Bayesian Framework" and 2) "Basic Example". I believe these two sections compliment each other, and help give context to the technical jargon. I also added a separate section to discuss Metropolis algorithms, including a basic description of how they work.

return a + b

With code-highlighting:

.. code-block:: python
def sum(a, b):
"""Sum two numbers."""
return a + b
Maybe also in another language, and with line numbers:

.. code-block:: c

int main() {
for (int i = 0; i < 10; i++) {
/* do something */
return 0;
Or a snippet from the above code, starting at the correct line number:

.. code-block:: c
:linenostart: 2

for (int i = 0; i < 10; i++) {
/* do something */
Inline code looks like this: :code:`chunk of code`.

Important Part

It is well known that Spice grows on the planet Dune. Test
some maths, for example :math:`e^{\pi i} + 3 \delta`. Or maybe an
equation on a separate line:

.. math::
g(x) = \int_0^\infty f(x) dx
or on multiple, aligned lines:

.. math::
:type: eqnarray
g(x) &=& \int_0^\infty f(x) dx \\
&=& \ldots
The area of a circle and volume of a sphere are given as

.. math::
:label: circarea
A(r) = \pi r^2.

This comment has been minimized.

Copy link

Mahdisadjadi Jun 12, 2019

The axes don't have labels/units (should be similar to Fig 1).

This comment has been minimized.

Copy link

prmiles Jun 24, 2019

Author Contributor

In my updated draft I have removed this figure. Thank you for pointing out the mistake.

.. math::
:label: spherevol
V(r) = \frac{4}{3} \pi r^3
We can then refer back to Equation (:ref:`circarea`) or
(:ref:`spherevol`) later.

Mauris purus enim, volutpat non dapibus et, gravida sit amet sapien. In at
consectetur lacus. Praesent orci nulla, blandit eu egestas nec, facilisis vel
lacus. Fusce non ante vitae justo faucibus facilisis. Nam venenatis lacinia
turpis. Donec eu ultrices mauris. Ut pulvinar viverra rhoncus. Vivamus
adipiscing faucibus ligula, in porta orci vehicula in. Suspendisse quis augue
arcu, sit amet accumsan diam. Vestibulum lacinia luctus dui. Aliquam odio arcu,
faucibus non laoreet ac, condimentum eu quam. Quisque et nunc non diam
consequat iaculis ut quis leo. Integer suscipit accumsan ligula. Sed nec eros a
orci aliquam dictum sed ac felis. Suspendisse sit amet dui ut ligula iaculis
sollicitudin vel id velit. Pellentesque hendrerit sapien ac ante facilisis
lacinia. Nunc sit amet sem sem. In tellus metus, elementum vitae tincidunt ac,
volutpat sit amet mauris. Maecenas [#]_ diam turpis, placerat [#]_ at adipiscing ac,
pulvinar id metus.

.. [#] On the one hand, a footnote.
.. [#] On the other hand, another footnote.
.. figure:: figure1.png

This is the caption.:code:`chunk of code` inside of it. :label:`egfig`

.. figure:: figure1.png
:align: center
:figclass: w

This is a wide figure, specified by adding "w" to the figclass. It is also
center aligned, by setting the align keyword (can be left, right or center).
This caption also has :code:`chunk of code`.

.. figure:: figure1.png
:scale: 20%
:figclass: bht

This is the caption on a smaller figure that will be placed by default at the
bottom of the page, and failing that it will be placed inline or at the top.
Note that for now, scale is relative to a completely arbitrary original
reference size which might be the original size of your image - you probably
have to play with it. :label:`egfig2`

As you can see in Figures :ref:`egfig` and :ref:`egfig2`, this is how you reference auto-numbered

.. table:: This is the caption for the materials table. :label:`mtable`

| Material | Units |
| Stone | 3 |
| Water | 12 |
| Cement | :math:`\alpha` |

We show the different quantities of materials required in Table

.. The statement below shows how to adjust the width of a table.
.. raw:: latex


.. table:: This is the caption for the wide table.
:class: w

| This | is | a | very | very | wide | table |

Unfortunately, restructuredtext can be picky about tables, so if it simply
won't work try raw LaTeX:

.. raw:: latex

\multirow{2}{*}{Projection} & \multicolumn{3}{c|}{Area in square miles}\tabularnewline
& Large Horizontal Area & Large Vertical Area & Smaller Square Area\tabularnewline
Albers Equal Area & 7,498.7 & 10,847.3 & 35.8\tabularnewline
Web Mercator & 13,410.0 & 18,271.4 & 63.0\tabularnewline
Difference & 5,911.3 & 7,424.1 & 27.2\tabularnewline
Percent Difference & 44\% & 41\% & 43\%\tabularnewline
\caption{Area Comparisons \DUrole{label}{quanitities-table}}

Perhaps we want to end off with a quote by Lao Tse [#]_:

*Muddy water, let stand, becomes clear.*

.. [#] :math:`\mathrm{e^{-i\pi}}`
ProTip! Use n and p to navigate between commits in a pull request.
You can’t perform that action at this time.