In CharlotteJana/pdmppoly: Moment Approximation for Polynomial PDMPs

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-"
)

pdmppoly

The goal of the package pdmppoly is to simulate polynomial piecewise deterministic Markov processes (polynomial PDMPs) within R and to provide methods to calculate or approximate their moments.

In addition, it contains all models that are used in the (not yet published) doctoral thesis of Charlotte Jana.

Installation

You can install pmppoly from GitHub with:

# install.packages("devtools")
devtools::install_github("CharlotteJana/pdmppoly")

Polynomial PDMPs

A polynomial PDMP is a piecewise deterministic Markov process (PDMP), whose dynamics and rates are all polynomials, depending on the different variables and parameters of the PDMP.

A polynomial PDMP in pdmppoly is represented by an object of class polyPdmpModel. This class is very similar to class pdmpModel of the package pdmpsim and every method provided by pdmpsim can be applied to polyPdmpModel-objects. The polynomials appear in the new slots dynfunc and ratefunc. They are defined using the package spray and internally represented as sparse coefficient matrices.

This is a simple example modelling gene expression with positive feedback:

library(pdmppoly)

genePolyBF <- new("polyPdmpModel",
                  descr = "Gene regulation with positive feedback",
                  parms = list(b = 0.5, a0 = 1, a1 = 3, k10 = 1, k01 = 0.5), 
                  init = c(f = 1, d = 1), 
                  discStates = list(d = 0:1),
                  dynpolys = quote(list(
                    list(overall = -b*lone(1,2),
                         specific = list(a0, a1))
                  )),
                  ratepolys = quote(list(  
                    list(k01*lone(1,2), k10)
                  )),
                  jumpfunc = function(t, x, parms, jtype) {
                    c(x[1], 1 - x[2])
                  }, 
                  times = c(from = 0, to = 25, by = 0.1), 
                  solver = "lsodar")

Moment calculation

The function momApp calculates the moments of the distribution of of a polynomial PDMP, up to a given order. It solves a system of differential equations with the package deSove to obtain the results. This system is usually not closed and has to be altered by replacing moments of higher orders with fixed values. The moments to replace can be raw or central moments, depending on the argument centralize. The following closure methods are available:

• closure = "zero": Replace the moments with 0
• closure = "normal": Replace with moments of a normal distribution
• closure = "lognormal": Replace with moments of a lognormal distribution
• closure = "gamma": Replace with moments of a gamma distribution

mom <- momApp(genePolyBF, maxorder = 8,
              closure = c("zero", "normal", "zero"), 
              centralize = c(FALSE, TRUE, TRUE))
mom

plot(mom, plotorder = 1:2)

Compare with empirical moments

The results of momApp are approximated values of the moments of the distribution. To test their accuracy, it is best to simulate the model multiple times, calculate the empirical moments and compare them with the approximated moments.

simulations <- multSim(genePolyBF, seeds = 1:300)
mom <- addSimulation(mom, simulations)
plot(mom, plotorder = 1)

This simulation takes several minutes. If time consuming simulations are necessary, it is advised to define the same model as object of class pdmpModel and to use this object for simulation. Simulating with a polyPdmpModel-object is slower because of the internal representation of functions as coefficient matrices.

These are the results:

Testing if the distribution is unimodal

The results of the function momApp can also be used to test if the distribution of the PDMP is unimodal. The function is.unimodal from the package momcalc performs this test for one dimensional distributions. It is only meaningful for distributions with a known compact support. The package pdmppoly provides the function modalityTest to apply the test directly on a data.frame generated by momApp. It performs the test for every time value and every variable of the PDMP separately.

testresults <- modalityTest(mom, lower = 0, upper = 35, vars = "f")
head(testresults)

plotModalities(mom, modalities = testresults)

Be aware that both momApp and modalityTest do not depend on time-consuming simulations!

Analyse a polynomial PDMP

The package pdmppoly provides a function analyseModel that analyses a polynomial PDMP as well as possible, meaning that it performs most of the functions available on the packages pdmpsim and pdmppoly. It

• simulates the model multiple times,
• plots the simulated data with all available plot methods,
• calculates statistics such as mean, sd, median, etc. ,
• performs moment approximation with all available closure methods,
• tests if the distribution is unimodal and plots the result.

License

This package is free open source software licensed under the GNU Public License (GPL 2.0 or above). The software is provided as is and comes WITHOUT WARRANTY.

CharlotteJana/pdmppoly documentation built on Sept. 4, 2019, 4:40 p.m.