DISCRETE AND CONTINUOUS
DYNAMICAL SYSTEMS SERIES B

doi:10.3934/dcdsb.2017133

ADAPTIVE METHODS FOR STOCHASTIC DIFFERENTIAL
EQUATIONS VIA NATURAL EMBEDDINGS ANDREJECTION SAMPLING WITH MEMORY

Christopher Rackauckas
Department of Mathematics, Center for Complex Biological Systems
University of California, Irvine
CA 92697, USA

Qing Nie∗
Department of Mathematics, Center for Complex Biological Systems
University of California, Irvine
CA 92697, USA(Communicated by Yuan Lou)
Abstract. Adaptive time-stepping with high-order embedded Runge-Kutta
pairs and rejection sampling provides efficient approaches for solving differential equations. While many such methods exist for solving deterministic
systems, little progress has been made for stochastic variants. One challenge
in developing adaptive methods for stochastic differential equations (SDEs) isthe construction of embedded schemes with direct error estimates. We present
a new class of embedded stochastic Runge-Kutta (SRK) methods with strong
order 1.5 which have a natural embedding of strong order 1.0 methods. This
allows for the derivation of an error estimate which requires no additional
function evaluations. Next we derive a general method to reject the time
steps without losing information about the future Brownian path termed Rejection Sampling with Memory (RSwM). This method utilizes a stack data
structure to do rejection sampling, costing only a few floating point calculations. We show numerically that the methods generate statistically-correct
and tolerance-controlled solutions. Lastly, we show that this form of adaptivity
can be applied to systems of equations, and demonstrate that it solves a stiff
biological model 12.28x faster than common fixed timestep algorithms. Our
approach only requires the solution to a bridging problem and thus lends itself
to natural generalizations beyond SDEs.

1. Introduction. Explicit methods for solving Ordinary Differential Equations(ODEs) with adaptive time-stepping algorithms, such as the Dormand-Prince andCash-Karp algorithms, have become efficient and popular solvers for numerical simulations of ODEs due to their ease of use and accuracy [8, 24, 6]. These algorithms
consist of two major components: an embedded higher-order and lower-order pairof temporal integrators to allow for estimating the local error (i.e. the error indicator), and an adaptive time-stepping algorithm for choosing the next stepsizeand performing the update [21, 1]. Together, the algorithm is able to effectively
2010 Mathematics Subject Classification. Primary: 65C30, 60H10; Secondary: 68P05.
Key words and phrases. Stochastic differential equations, adaptive methods, rejection sampling, embedded algorithms, stochastic Runge-Kutta, strong approximation.
∗ Corresponding author: xxxx.

1

2

CHRISTOPHER RACKAUCKAS AND QING NIE

solve most equations by automatically adapting to the solutions by adjusting thetimestep.
To derive adaptive algorithms for stochastic ordinary differential equations(SODEs) of the form
dXt = f (t, Xt )dt + g(t, Xt )dWt ,(1)

where f is the drift coefficient and g is the diffusion coefficient, a natural choice
is to use embedded Stochastic Runge-Kutta (SRK) methods. Kloeden and Platen
developed a set of stochastic Taylor expansions which allowed for the derivation ofhigher order (greater than order 1.0) SRK methods [15, 13]. One interesting featureof this class of algorithms is that, since no analytical solution exists for iterated
stochastic integrals, they are usually replaced by sufficiently exact approximations
due to Wiktorsson [25, 23]. While strong order 2.0 SRK methods have been studied[3, 15], existing works mainly focus on strong order 1.5 Stochastic Runge-Kutta
methods, such as a series of methods for Stratanovich SDEs [4, 2] and for Ito SDEs
[14].
One recent advance in SRK methods is Rößler’s methods for Ito SDEs [22]. This
class of methods was shown to be efficient by reducing the number of function
evaluations required per step. Their format via extended Butcher tableaus allows
them to be easily amendable. However, these methods have not been used foradaptive algorithms. Here, we will show that for certain Rößler-SRK methods,
there exists a naturally embedded order 1.0 method such that no extra function
evaluations are needed in order to evaluate its difference from the order 1.5 method,presenting an efficient way to both perform the steps and approximate the error.
With the embedded SRK methods, the next challenge is to deal with the choiceof time steps. For ODEs, one usually performs a trial step in order to calculate anerror estimate and, using this estimate, either choose to step or reject the step andtry again with a smaller timestep. The problem with such an approach for SODEs
is that when performing the trial step, one takes a random number to simulate
the future Brownian path, and if one rejects the timestep, this will change thesample properties of the Brownian path [5]. This is particularly an issue because
larger variations in the Brownian path are correlated with more numerical error andthus such an algorithm would have a strong tendency to reduce the variance of thenoise process. Lastly, it’s not clear how to effectively sub-sample a solution along a
Brownian path after one decides it’s too coarse.
Most of the current approaches for adaptive time steps attempt to perform error
calculations before stepping in order circumvent rejection sampling. One approach[9] utilizes an order 1.0 algorithm with a Brownian path stored at the h, h2 , h4 timepoints, etc. stored as a tree. The algorithm solves for an error equation using astochastic Taylor series and then utilizes this new equation to solve a linearization
backwards in time in order to produce an error estimate for a future timestep without having to do rejection sampling. This method requires using halving/doublingof stepsizes and, since it’s not able to account for which timesteps the Brownian
path has the largest future jumps, has to be conservative in order to achieve its chosen error estimate. It also requires the user to provide many derivatives of the drift
and diffusion functions f and g, or alternatively perform a numerical estimation ofthese derivatives at each timestep. Lastly, in their paper they proved that any variable stepsize implementation must use a numerical method which has strong order

NATURAL EMBEDDINGS AND REJECTION SAMPLING WITH MEMORY

3of at least 1, meaning that any variable stepsize method must go beyond Euler-type
methods.
Another algorithm may be limited to low-order methods due to its requirementof having a small “fixed number of observations”, and requires the calculation ofderivatives of both f and g [19]. Lamba [16] presented an order 1.0 algorithm whichused properties of the Brownian bridge restricted to a grid of multiples of 13 . It
included the novel property of having separate error estimators for deterministic
and noise induced error for finer control of the solution. There is one adaptive
stepping algorithm which does not have restrictions on the stepsize and utilizes
embedded strong order 1.0/1.5 SRK methods for Stratanovich SDEs [5] . Theirembedded method requires two additional function evaluations in order to derive
an error estimate, and requires solving for a dense matrix and performing a matrix
multiplication each time an interval is subdivided. Their results show that generaltime-stepping tends to be vastly more efficient than grid-based methods from before.
The main results of the paper are the development of strong order 1.0/1.5 embedded Runge-Kutta pairs for error estimation and a new adaptive time-stepping
algorithm termed Rejection Sampling with Memory (RSwM). This extension to therejection sampling algorithm is diagrammed in Figure 1 and gives an overview of
how the various developments come together to form one complete algorithm.

Update t:=t+h

Use the future information
to determine W(t+qh)

Update h:=q*h
Accept if q≥1

Compute a step ofsize h using two
different integration
methods (e.g. Eq 2)

Estimate the error
E using the difference
between the methods(Sections 2 and 3)

Use E tocalculateq (Eq 21)

Update h:=q*h

Save the information forthe change of W overthe interval (t+qh,t+h)

Use the Brownian Bridge
to determine a value
for W(t+qh) (Eq 22)

Reject if q<1

Figure 1. Outline of the adaptive SODE algorithm based on Rejection Sampling with Memory. This depicts the general schema
for the Rejection Sampling with Memory algorithm (RSwM). The
green arrow depicts the path taken when the step is rejected,
whereas the red arrow depicts the path that is taken when the stepis accepted. The blue text denotes steps which are specific to thestochastic systems in order to ensure correct sampling properties
which are developed in section 4.
In section 2 we construct the method ESRK1 which is a strong order 1.0/1.5
embedded pair of algorithms which allows for high-order solving and error estimation. In section 3 we prove the existence of a natural order 1.0/1.5 embedded pairof Rößler SRK algorithms that allows for the efficient calculation of an error estimate for high-order SRK methods with no additional function evaluations. Using
this error estimate, we construct the rejection sampling methods for SDEs whichallows for efficient general adaptive time-stepping with almost no extra floating

4

CHRISTOPHER RACKAUCKAS AND QING NIE

point operations. In section 4 we first elucidate a simple solution to general adaptive time-stepping, discuss how one could attempt to improve the efficiency of thealgorithms, and then develop more efficient adaptive stepping algorithms.
Lastly, we perform numerical experiments to show the effectiveness of the algorithms. In section 5, we show that our methods are able to reproduce the proper
sample statistics for the stochastic differential equations and the user chosen local
error parameter  is able to control the global error. Then in section 6, we demonstrate the ability for the adaptive parameters to control for the coarseness of thesolution to nonlinear systems of equations, and show that when applied to a large
stiff system derived from biology, our method is able to stability solve a Monte
Carlo experiment of 10000 simulations without diverging 12.28x faster than a few
popular fixed timestep methods. We conclude by discussing how these algorithms
can be generalized to other problems like Stratanovich SDEs, jump diffusions, andSPDEs.
2. Construction of an order 1.0/1.5 embedded pair. Using a colored root
tree analysis, Rößler developed a set of strong order 1.5 stochastic Runge-Kutta
schemes [22]. This multi-step method resulted in less computational steps than the
KPS schemes, and the number of steps grows much slower as the Ito dimension
increases than in the KPS schemes [15]. The structure of Rößler methods are
relatively simple, making the method both easy to implement and fast in runtime.
For simplicity, we describe the equations for a single Ito dimension, though theresults generalize to the case of multiple Ito dimensions.
The Runge-Kutta methods are strong order 1.5 if they satisfy a set of orderconditions in A, taking the following form:

Un+1

=

Un +

s
X



(0)
(0)
∆t +
αi f tn + ci ∆t, Hi

(2)

i=1

 
s 

X
(1)
(3) I(1,0)(4) I(1,1,1)(1)(1)(2) I(1,1)
,
+ βi
+ βi
g tn + ci ∆t, Hi
βi I(1) + βi √
∆t
∆t
∆t
i=1

with stages(0)Hi= Un +

s
X



(0)
(0)
(0)
Aij f tn + cj ∆t, Hj
∆t

(3)

j=1

+
(1)

Hi

s
X



(0)
(1)
(1) I(1,0)
Bij g tn + cj ∆t, Hj
,
∆t
j=1

= Un +

s
X



(1)
(0)
(0)
Aij f tn + cj ∆t, Hj
∆t

(4)

j=1

+

s
X

√

(1)
(1)
(1)
Bij g tn + cj ∆t, Hj
∆t.

j=1

The iterated stochastic integrals I(1,1) , I(1,0) , and I(1,1,1) are evaluated via the approximations due to Wiktorsson as noted in C. These methods are referred to as the
SRI algorithms and are indicated by the tuple of 44 coefficients A(i) , B (i) , β (i) , α .
Similarly, for additive noise the Runge-Kutta methods are strong order 1.5 if they

NATURAL EMBEDDINGS AND REJECTION SAMPLING WITH MEMORY

satisfy the order conditions contained B, taking the form
s


X
(0)
(0)
Un+1 = Un +αi f tn + ci ∆t, Hi
∆t +

5

(5)

i=1


s 
X
(1)
(2) I(1,0)(1)βi I(1) + βig(tn + ci ∆t),∆ti=1with stages(0)Hi= Un +

s
X

s



I
X
(1,0)
(0)
(0)
(0)
(0)
(1)
. (6)
Aij f tn + cj ∆t, Hj∆t +Bij g tn + cj ∆t
∆t
j=1
j=1

Given the efficiency of order 1.5 SRK methods for stochastic problems, we willdevelop Runge-Kutta-Fehlburg-like embedded algorithms [1], which step along anorder 1.5 path using error estimates from an order 1.0 path.
We start with the order 1.5 method SRIW1 [22] denoted by coefficients A(i) , B (i) ,T
β (i) , α . Notice that the order 1.0 conditions from A only enforce β (3) B (1) e = 0
T
(3)
(4)
(3)
(4)
and β (4) B (1) e = 0 on β and β . Thus let β̃ = β̃ = 0. This gives an order
1.0 method A(i) , B (i) , β̃ (i) , α .
Let X̄ be the numerical solution of SRIW1 on the general SODE Equation 1, and
X̃ be the numerical solution from the algorithm with coefficients A(i) , B (i) , β̃ (i) , α .
By Equation 2, the difference between the numerical solutions on the next step is
 
s 

X
(4) I(1,1,1)(1)(1)(3) I(1,0)+ βig tn + ci ∆t, Hi.
(7)X̄ − X̃ =
βi
∆t
∆t
i=1
Note that we could keep β (3) not equal to zero and still have an order 1.0 method,but zeroing β (3) gives an estimate for the size of ∆Z (the second Brownian pathused to estimate I(1,0) as noted in C) which may need to be controlled as well.
However, this error estimate does not take into account stiffness
P in the drift
coefficient f . Notice
that
the
only
conditiononthecoefficientsαisαi = 1. Thus
1 1
let α̃ = 2 , 2 , 0, 0 . This gives another order 1.0 Method (A0 , B0 , βi , α̃). If we let X̂
be the numerical solution to Equation 1 using the coefficients (A0 , B0 , βi , α̃), thenby Equation 2 notice ∆t 

∆t 
(0)
(0)
(0)
(0)
+
.
(8)
X̄ − X̂ = − f tn + c1 ∆t, H1
f tn + c2 ∆t, H2
6
6
Notice that for this algorithm the difference is only due to the stiffness of the drift
coefficient. Using both the β and the α changes together gives a pair of strong order
1.0 / order 1.5 methods which is robust to errors in both coefficients. We designate
this algorithm Embedded SRK 1 (ESRK1) as the set of coefficients defined in Table 1
for which steps are taken according to Equation 2, Equation 3, Equation 4.
Adding the two previous error equations Equation 7 and Equation 8, gives E,
the error estimation via the difference between the order 1.5 and order 1.0 method,as ∆t 

∆t 
(0)
(0)
(0)
(0)
E = − f tn + c1 ∆t, H1+
f tn + c2 ∆t, H2
(9)
6
6
 
s 

X
(3) I(1,0)(4) I(1,1,1)
(1)
(1)
+
βi
+ βi
g tn + ci ∆t, Hi
.
∆t
∆t
i=1

6

CHRISTOPHER RACKAUCKAS AND QING NIE

c(0)

A(0)

B (0)

c(1)

A(1)

B (1)

αT

β (1)

T

β (2)

β (3)

T

β (4)

T

β̃ (4)β̃ (3)

α̃T

T

T

T

(a) Legend Table
0
3
4

3
4

3
2

0

0

0

0

0

0

0

0

0

0

0

0

0
1
4

1
4

1

1

0

1
4

0

0

1
4

1
3

2
3

0

1
2

1
2

1
2

0

0

0

-1

0

-5

3

1
2

-1

4
3

2
3

0

-1

4
3

− 13

0

2

− 43

− 23

0

−2

5
3

− 23

1

0

0

0

0

0

0

0

0

(b) Coefficients Table

Table 1. ESRK1. Table (a) shows the legend for how the numbers in in Table (b) correspond to the coefficient arrays/matricesc(i) , A(i) , B (i) , α, β (i) , α̃, and β˜(i) . For example, these tables showthat αT = ( 13 , 23 , 0, 0). Note that the matrices A(i) and B (i) arelower triangular since the method is explicit.

Thus ESRK1 is a pair of order 1.0/order 1.5 embedded algorithms. Notice thatevery term in the error calculation Equation 9 is also part of the calculation forthe order 1.5 timestep from Equation 2 and thus no additional function evaluations

NATURAL EMBEDDINGS AND REJECTION SAMPLING WITH MEMORY

7

are required for this error estimate. Thus ESRK1 is an extension to the order 1.5
algorithm which allows for an efficient error estimation.
3. Error estimation via natural embeddings. The construction from section 2
is part of a more general construction which we call natural embeddings in Rößler 
SRK methods. Take a Rößler SRK method determined by coefficients A(i) , B (i) , β (i) , α .
For this method, we can generate an embedded order 1.0 method via a parameterδ by changing α to α̃ where we subtract δ from two coefficients and add δ to twocoefficients. Thus by Equation 2, this gives the deterministic error estimate
 



(0)
(0)
(0)
(0)
ED = ∆t f tn + ck1 ∆t, Hk1 + f tn + ck2 ∆t, Hk2
(10)
 



(0)
(0)
(0)
(0)
− f tn + ck3 ∆t, Hk3 + f tn + ck4 ∆t, Hk4,
If we then change β (3) and β (4) to zero on the indices in the index set I2 , then thisgives the noise error estimate X  (3) I(1,0)(4) I(1,1,1)
(1)
(1)
+ βi
g tn + ci ∆t, Hi
.
EN =
βi
∆t
∆t

(11)

i∈I2

By the order conditions of A, these coefficient changes together give an order 1.0
method whose difference from the order 1.5 method A(i) , B (i) , β (i) , α is a boundederror estimate which is calculated directly from the coefficients of the order 1.5
method. This result is summarized in the following theorem:
Theorem 3.1 (Natural Embedding). For any order 1.5
coeffi SRK method with


(i)
(i)
(i)
(i)
(i)
cients A , B , βi , α , there exists a δ ∈ R such that A , B , β̃ , α̃ defines anatural order 1.0 embedded SRK method, where β̃ (3) and β̃ (4) are equivalent to β (3)and β (4) respectively except the indices in I2 are zeroed, and α̃ differs from α by δat each coefficient in I1 by a sign difference of (−1)σ(i) (σ(i) = 1 for at least one iand σ(i) = 2 for at least one i), and |I1 | even. The
 one-step difference between the

(i)
(i)
(i)
(i)
(i)
(i)
methods A , B , β , α and A , B , β̃ , α̃ is bounded by
E = δED + EN ,

(12)

where
ED

=

∆t

X
i∈I1

EN

=




X 
(0)
(0)
(0)
(0)
(−1)σ(i) f tn + ci ∆t, Hi≤ ∆t
f tn + ci ∆t, Hi
,(13)
i∈I1

 

X  (3) I(1,0)(4) I(1,1,1)
(1)
(1)
βi
+ βi
g tn + ci ∆t, Hi
.
∆t
∆t
i∈I

(14)

2

Notice that from this theorem if we take the coefficients from the SRIW1 method
and let δ = 16 then we arrive at the ESRK1 method from section 2. Also, note that
this theorem be extended to the case where |I1 | = 3 by letting α̃ differ from α by δδin one coefficient
P and 2 in the two others. However, in order to satisfy the Order
0.5 constraint
αi = 1, we cannot simply modify a single component.
This implies that there is no need to specifically derive embedded pairs of methods. Instead, for any order 1.5 Rößler SRK method that one constructs, there exists
a naturally embedded order 1.0 method such that E is an error estimate, and thus
the error estimator can be utilized without ever explicitly constructing the order 1.0
method. In addition, the calculation of E only uses values found in the order 1.5

8

CHRISTOPHER RACKAUCKAS AND QING NIE

timestep from Equation 2. As a result, the error estimate is naturally obtained without any additional function evaluations. Lastly, the equation Equation 12 also hasa adjustable parameter δ, which allows the user to scale the relative contributionsof the deterministic and noise terms to the overall error estimate.
We note that all of these equations generalize to systems by considering the variables and functions as vectors and vector functions respectively, with the absolutevalue generalizing to an appropriate norm. This result also generalizes to multiple
Ito dimensions, i.e., SODEs of the form
dXt = f (Xt , t)dt +dX

gk (Xt , t)dWt ,

(15)

k

in the case where the noise function g is diagonal by a similar construction on theRößler SRID, giving the error estimates
ED

X




X 
(0)
(0)
(0)
(0)
,(16)
≤ ∆t
f tn + ci ∆t, Hi(−1)σ(i) f tn + ci ∆t, Hi

=

∆t

=

 
d X

X
(1)
(4) I(1,1,1)(1)(3) I(1,0).
+ βi
gk tn + ci ∆t, Hi
βi
∆t
∆t
i∈I

i∈I1

EN

k

i∈I1

(17)

2

Note that in the case of SRA-type methods (additive noise), a similar theorem also
holds. By Equation 5 on the multiple Ito dimension SDE Equation 15 we utilize
the same construction and receive a natural embedding with the error estimation
Equation 12 containing the same deterministic error estimation Equation 13, witha new noise-error estimation:d X
sX(2) I(1,0)
(1)
βi
EN =
gk (tn + ci ∆t) .
(18)
∆t
i=1
k

4. Three algorithms for rejection sampling with memory.
4.1. Rejection sampling for stochastic differential equations. As in deterministic rejection sampling algorithms [1, 21, 11], definev
u n  u 1 X Ei
e=t
,
(19)
n i
sci
where
(i)

sci = abs + Xt rel ,(20)
is the absolute/relative mixed error with  is the user-chosen error absolute andrelative tolerances respectively, Ei is the error estimation of the ith element computed with the methods from section 3, and n is the size of the system. Using thismeasure of the system error with respect to tolerance, it is common in deterministic
systems to define 1/(O+1)
1
,
(21)
q=
γe
where γ is a penalty factor (in deterministic methods it is often taken as γ = 2)
and O is the minimum order of the numerical methods used in the error estimators.
For solving deterministic equations, the rejection sampling method for adaptive
time-stepping can be summarized as:
1. If q < 1, reject the initial choice of h and repeat the calculation with qh

NATURAL EMBEDDINGS AND REJECTION SAMPLING WITH MEMORY

9

2. If q ≥ 1, then accept the computed step and change h to min (hmax , qh) where
hmax is chosen by the user.
For solving SODEs, this approach needs to be modified since rejection of future
steps involves throwing away information about the future Wiener process whichbiases the sample statistics for the path. Our solution is to change rejection samplingso that all information about the future path is kept, and adapt the sampling of
the Brownian path to accommodate for this information. This is summarized asfollows:
1. If q < 1, reject the initial choice of h, store the values for the future steps ofthe Wiener process ∆Wt , and repeat the calculation with qh.
2. If q ≥ 1, then accept the computed step and change h to min (hmax , qh) where
hmax is chosen by the user. Use future information to determine the next step.
The algorithm for adapting the random number sampling to account for future
information is as follows. As noted in C, due to the term I(1,0) in the SRK algorithm,
we must consider two Brownian paths: Wt , the path from the SODE, and Zt , anindependent Brownian path for the approximation. First, propose a step with ∆W P
and ∆Z P for a timestep h. If these are rejected, change the step size to qh. Thus
we need to sample a new value at W (tn + qh) using the known values of W (tn )and W (tn + h). By the properties of the Brownian Bridge [20, 17], if W (0) = 0,
W (h) = L, and q ∈ (0, 1), thenW (qh) ∼ N (qL, (1 − q)qh) .

(22)

Thus propose a step of size qh and take the random numbers ∆W = W (qh) and∆Z = Z(qh) from the distribution Equation 22.
Notice that by the Markov property of Brownian motion, the immediate future
and current timestep values fully characterize the distribution of the Brownian
Bridge, and using the Brownian Bridge we can interpolate the Brownian process towhatever time-point we deem necessary. Thus if the every time we take a sampleof the future Brownian process (i.e. take a random number) we save the random
number (even if rejecting the step), we can recover the correct distribution forinterpolations of the Brownian process using Equation 22, giving all of the necessary
information for computing a step of size qh.
4.2. Computational issues. First we deduce the proper form to save the future
information upon rejection and the necessary details to save. Assume that we areat a time t and the Brownian process changed by an amount ∆W P in the timestep
of size h. If there was no future information, this would be the value given by arandom number generator with the distribution N (0, h). When a step is rejected, a
new stepsize is chosen as qh where q < 1. Knowing the Brownian process’ value at t
and t + h, we can determine a value for W (t + qh) using Equation 22. Labeling ∆W
as the change in W over our new step of size qh, we have future knowledge thatthe Brownian path changes over the interval (t + qh, t + h) by ∆W = ∆W P − ∆W(note that for the same reasons, we have that the secondary Brownian motion Z
used as part of the high-order integrator must also satisfy ∆Z = ∆Z P − ∆Z). To
store this future information, we save the 3-tuple (∆W , ∆Z, (1 − q)h). The firsttwo values denote how much the Brownian motions changed and the last denotes
the length of the time interval over which the changes occurred. From these values
we have the information necessary for to use Equation 22 in order to interpolate
the change in Brownian motion from t + qh to any value in [t + qh, t + h]. Since we10

CHRISTOPHER RACKAUCKAS AND QING NIE

know the value of the Brownian motions at t + qh, we can perform rejections using
this algorithm recursively to see that we can re-reject to any value in [t, t + qh] andstill accurately re-create the Brownian motions at any point in the full [t, t + h].
Therefore saving this 3-tuple at each rejection suffices for the full reconstructionof the Brownian processes. Since this algorithm requires very few floating point
operations, this algorithm has a good potential for efficiency.
At this point a naive implementation runs into many practical computational
issues. The first implementation which comes to mind is to store this 3-tuple offuture information into an array. The storage must be a mutable multi-valued data
structure because it is possible for consecutive rejections to occur, meaning that we
may have to store many pieces of future information. However, when using general
stepsizes, it is possible for a step to go past some values of future information, andif this step is then rejected and the future information is added to the end of thearray, the array is no longer necessarily in time order. The naive approach to solve
this problem is to expand to saving a 4-tuple which also includes the timepoint,
and at every timestep either search the array or re-sort the array given new inputs.
While the search can be made O(log n) by using bisection (where n is the numberof future information points), insertion into an arbitrary point in an array still
constitutes O(n) operations since the other elements in the array must be moved.
Since this is the cost per insertion and many insertions can required before accepting
a timestep, it is essential to the efficiency of the algorithm to make this process more
performant.
Knowing that we will need to insert arbitrary values into the middle of the array,
one may then propose using a linked list data structure. However, if the searchis done on a linked list data structure, then the bisection still constitutes O(n)operations due to the traversal structure. To alleviate these problems, we propose
an algorithm using stacks. Stacks are computationally efficient data structures
which allow for O(1) insertion and retrieval, but with the restriction that only thevalue at the top of the stack is readily available. By achieving this complexity,
we minimize the computational cost of rejections, allowing for one to not have to
be overly conservative with the choice of stepsize. In the following subsection wewill show how to achieve general timestepping with O(1) information insertions andretrievals while only requiring the necessary floating point operations and ordering
the data in a manner which is efficient for caching.
4.3. The RSwM algorithms. Due to the LIFO (last in, first out) structure of
the stack data structure, if the algorithm iteratively adds the future portions tothe stack, then the tuple at the top of the stack always denotes the values for thetime interval of the immediate future. Thus, after a successful step, if the stack is
non-empty, we use the values from the top of the stack to propose our next step.
If the stack is ever empty, then this indicates that the algorithm currently has no
future information. In this case, take new random numbers and propose a new stepin the same manner as the non-adaptive algorithm. By using the stack in such amanner, we are able to perform all of the memory operations in O(1) with the only
additional floating point operations due to adaptation being the few operations to
calculate the Li , resulting in a computationally efficient algorithm.
A psudocode version of the algorithm is given as Algorithm 1. While intuitive,
RSwM1 is not the most efficient since it does not always use the estimated “beststepsize” qh. If the stack is non-empty, then the step size is taken to be the firstcomponent from the tuple on the top of the stack (L1 ) which may be significantly

NATURAL EMBEDDINGS AND REJECTION SAMPLING WITH MEMORY

11

Algorithm 1 RSwM1
1: Set the values , hmax , T
2: Set t = 0, W = 0, Z = 0, X = X0
3: Take an initial h, ∆Z, ∆W ∼ N (0, h)
4: while t < T do
5:
Attempt a step with h, ∆W , ∆Z to calculate Xtemp according to (2)Calculate E according to (9)Update q using (21)if (q < 1) then
. % Reject the Step
9:
Take ∆W̃ ∼ N (q∆W, (1 − q)qh) and ∆Z̃ ∼ N (q∆Z, (1 − q)qh)10:
Calculate ∆W = ∆W − ∆ W̃ and ∆Z = ∆Z − ∆Z̃
11:Push (1 − q) h, ∆W , ∆Z into stack S
12:Update h := qh
13:Update ∆W := ∆W̃ , ∆Z := ∆Z̃14:else. % Accept the Step
15:
Update t := t + h, W := W + ∆W , Z := Z + ∆Z, X = X + Xtemp
16:
if (S is empty) then17:
Update c := min(hmax , qh), h := min(c, T − t)18:
Take ∆W, ∆Z ∼ N (0, h)
19:
else
20:
Pop the top of S as L
21:
Update h := L1 , ∆W := L2 , ∆Z := L3
22:
end if23:
end if
24: end while
6:
7:
8:

less than qh. This constrains the algorithm to step to any timepoint t which was
the value of a previously rejected timestep. Thus, RSwM1 is not an entirely general
stepping algorithm.
To extend the previous algorithm, we wish to change the stepping algorithm
when the stack is not empty to allow for unraveling the stack multiple times per
step. To do so, let the algorithm take enough items off the stack so that the total
length is close to the proposed timestep qh, and adjust for the difference.
A psudocode version of the algorithm is given as Algorithm 4 in Appendix D.
The resulting algorithm, RSwM2, is more complex than RSwM1. However, theadvantage of this algorithm is that if we propose a step of size hnew = qh, we will
always step with a size hnew regardless of the number of items on the stack. In theexample simulations in section 5, this algorithm is shown to be more efficient.
However, RSwM2 is not correct. The problem comes from “re-rejections”. If
the algorithm pops a more than one item off the stack and then rejects the step,these elements are lost. Thus when it tries to re-step with the proposed qh, it does
not properly recreate our Brownian path since we will have lost some of the future
information. In section 5 we see that this case is rare enough that the results areusually statistically correct, though we wish to extend this algorithm to a truly
correct algorithm.
To ensure correctness we need to save the popped elements until an accepted
step passes the point which they encode. The idea is to use a deque where future
information is added to the back and re-popped information is added to the front.

12

CHRISTOPHER RACKAUCKAS AND QING NIE

However, whenever a step is accepted, we will wish to drop all of the elements addedfrom the front. This means that the most natural algorithm would be similar to thecommon decomposition of the deque into two stacks. Instead of utilizing a dequedirectly, we will use one stack S1 for the future information and another S2 for the
re-popped information. Whenever a value is popped from the future information
stack S1 , it will be placed into S2 . If the step is rejected, those values will be popped
back over to S1 , and if the step is accepted they will be discarded.
The algorithm RSwM3 is explained in more detail by example. Assume we rejecta timestep which uses the first n elements from the stack, and that the new proposed
timestep qh uses what would have only been the first k. The fix would be to add
the last n − k − 1 tuples of future information back to the stack (in reverse timeorder), solve the Brownian bridge problem with the k and k + 1st items straddlingthe proposed step qh, and add the future information from qh to the k + 1st items
to the stack. The algorithm should then only “discard” items after they have been
passed by a successful step.
If we saved the n tuples as they popped from the stack S1 into a different stackS2 , then we can think of S2 as the stack which contains all future information that iscurrently being used in a proposed timestep. Since this will add the elements to S2
in chronological order, the elements of stack S2 will pop out in reverse time order.
Thus after rejecting the timestep, the n−k elements we wish to save will be the firstn − k elements of the stack S2 . Therefore, after a rejection, the algorithm should
pop off the first n − k − 1 elements of S2 , re-add them to S1 , and pop one moreelement to solve a Brownian Bridge problem. Then after any successful timestep,
the algorithm can simply discard all of the information in S2 .
A psudocode version of the algorithm is given as Algorithm 2. The resulting
algorithm, RSwM3, is an extension to RSwM2 which keeps the property of always
using the estimated “best timestep” qh but is now able to ensure correctness.
5. Numerical verification of correctness and efficiency. In order to evaluate
the efficiency of the adaptive methods, we implemented ESRK1 with the three
rejection sampling with memory algorithms RSwM1, RSwM2, and RSwM3. We
chose the three example equations Equation 31, Equation 33, and Equation 35
defined in E.
These three examples cover a variety of behaviors seen in SODEs. Example 1 isa linear stochastic differential equation. Since α > 0, the solution tends to infinity
and thus presents a case where the algorithms must decrease the timesteps in orderto ensure the accuracy. Example 2 has a periodic solution and will be used to showthat the error of our method decreases even in cases where the process is stable.
Example 3 will be used to test the algorithm’s ability to handle additive noise. For
all of the tests, the relative tolerance was set to 0 and the absolute tolerance wasadjusted in order to test various algorithm behaviors.
5.1. Correctness. We first checked for the correctness of the algorithm. In order to study the correctness of each implementation, we ran RSwM1, RSwM2,
and RSwM3 200 times on Equation 31, and tested the distribution of W2 . For√ 2 ∼ N (0, 1) should be a standard norBrownian motion, Wt ∼ N (0, t) and thus W
2
mal random variable. In order to test that our rejection sampling algorithms did
not change the distributions for the underlying Brownian processes, we used theKolmogorov-Smirnov Test. The 1-Sample Kolmogorov-Smirnov test is a standard

NATURAL EMBEDDINGS AND REJECTION SAMPLING WITH MEMORY

13

Algorithm 2 RSwM3
1: Set the values , hmax , T
2: Set t = 0, W = 0, Z = 0, X = X0
3: Take an initial h, ∆Z, ∆W ∼ N (0, h)
4: while t < T do
5:
Attempt a step with h, ∆W , ∆Z to calculate Xtemp according to (2)6:Calculate E according to (9)7:Update q using (21)8:
if (q < 1) then
. % Reject the Step9:
Set hs = 0, ∆W = 0, ∆Z = 0
10:while S2 is not empty do
11:
Pop the top of S2 as L
12:
if hs + L1 < (1 − q)h then13:
Update hs := hs + L1
14:
Update ∆Wtmp := ∆Wtmp + L2 , ∆Ztmp := ∆Ztmp + L3
15:
else
16:
Push L onto S2 and break17:
end if18:
end while19:
Set hK = h − hs , K2 = ∆W − ∆Wtmp , K3 = ∆Z − ∆Ztmp
20:
Set qK = hqh
K

21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:

Take ∆W̃ ∼ N (qK K2 , (1 − qK )qK L1 )Take ∆Z̃ ∼ N (qK K3 , (1 − qK )qK L1 )
Pop ((1 − qK )hK , K2 − ∆W̃ , K3 − ∆Z̃) onto S1
Pop (qK L1 , ∆W̃ , ∆Z̃) onto S2Update ∆W = ∆W̃ , ∆Z = ∆Z̃, h = qh
else
. % Accept the Step
Update t := t + h, W := W + ∆W , Z := Z + ∆Z, X = X + Xtemp
Empty S2
Update c := min(hmax , qh), h := min(c, T − t)
Set hs = 0, ∆W = 0, ∆Z = 0
while S is not empty do
Pop the top of S1 as L
if (hs + L1 < h) then
. % Temporary not far enoughUpdate hs := hs + L1 , ∆W := ∆W + L2 , ∆Z := ∆Z + L3
Push a copy of L onto S2
else. % Final part of step from stacks
Set qtmp = h−hL1Let ∆W̃ ∼ N (qtmp L2 , (1 − qtmp )qtmp L1 )Let ∆Z̃ ∼ N (qtmp L3 , (1 − qtmp )qtmp L1 )Push ((1 − qtmp )L1 , L2 − ∆W̃ , L3 − ∆Z̃) onto S1
Update ∆W := ∆W + ∆W̃ , ∆Z := ∆Z + ∆Z̃, hs := hs + qtmp L1
end ifend while. % Update for last portion to step. Note zero if final part is from stack44:
if (h − hs is not zero) then45:
Let ηW , ηZ ∼ N (0, h − hs )46:
Update ∆W = ∆W + ηW , ∆Z = ∆Z + ηZ
47:
Push (h − hs , ηW , ηZ ) onto S248:
end if49:
end if
50: end while
38:
39:
40:
41:
42:
43:

14

CHRISTOPHER RACKAUCKAS AND QING NIE

non-parametric test of distributions which is known for being very sensitive and easily rejecting the null-hypothesis of standard normality [10]. A plot of the p-values
for the Kolmogorov-Smirnov test on 20 trials tolerances of 10−1 , 10−3 , 10−5 is shownin Figure 2. We note that due to multiple sampling issues, one would expect 1/20
tests to fail at a p-value of .05. Indeed, for the RSwM1 and RSwM3 algorithms our
results have around the expected number of failures. We see that RSwM2 has 6
failures, which shows that it also generally passes the tests even though there is no
guarantee for its correctness. However, it is clear that the issue magnifies at lower
tolerances. These results indicate that the normalized Brownian process’s distribution was standard normal and thus were not significantly altered by the RSwM
sampling algorithms. Therefore RSwM algorithms generate statistically correct results.

Figure 2. Adaptive algorithm Kolmogorov-Smirnov Tests. Equation 31 was solved from t = 0 to t = 2. Scatter plots of the p-values
from Kolmogorov-Smirnov tests against the normal distribution.
At the end of each run, a Kolmogorov-Smirnov test was performedon the values at end of the Brownian path for 200 simulations. The
x-axis is the absolute tolerance (with relative tolerance set to zero)and the y-axis is the p-value of the Kolmogorov Smirnov tests.

As an additional verification, we included a Quantile-Quantile Plot (QQplot)√ 2 against 10,000 standard normal random variables generated via MATLAB’sof W
2randn command. The simulation results are shown in Figure 3. The estimated√ 2 as generated by RSwM 1 through 3 line up on the x = y axis,
quantiles for W2√ 2 follows a standard normal distribution. To test for edge efindicating that W
2fects, we ran the same calculation with higher/lower numbers of paths. This showsthat as the number of paths increases, the number of quantiles which align also
increases. These results show that although we cannot guarantee the correctnessof the RSwM2, its sample statistics cannot be distinguished from the correct algorithms. This indicates that the edge case that the algorithm omits is a rareoccurrence.

NATURAL EMBEDDINGS AND REJECTION SAMPLING WITH MEMORY

15

Figure 3. Adaptive algorithm correctness checks. QQplots of the
distribution of the Brownian path at the end time T = 2 over 10,000
paths. The x-axis is the quantiles of the standard normal distribution while√ the y-axis is the estimated quantiles for the distributionof W2 / 2. Each row is for a example equation for Examples 1-3
respectively, and each column is for the algorithm RSwM1-3 respectively.  = 10−4 . The red dashed line represents x = y, meaning
the quantiles of a 10,000 standard normal random variables equate
with the quantiles of the sample.
The blue circles represent the√quantile estimates for W (2)/ 2 which should be distributed as astandard normal.

5.2. Accuracy and efficiency. Next we tested the accuracy and efficiency of thethree algorithms RSwM1, RSwM2, RSwM3 on the three example SDEs Equation 31,
Equation 33, and Equation 35. The results, as summarized in Figure 4 and Figure 5,
show that in almost every case, the improved algorithms RSwM2 and RSwM3 tend

16

CHRISTOPHER RACKAUCKAS AND QING NIE

to have lower runtimes and better convergence of the error at the final timepoint.
Notice that the RSwM2 and RSwM3 algorithms tend to have error convergencesmuch closer to log-linear, whereas RSwM1 tends to flatten out at lower tolerances.
Also notice how in all cases, the RSwM2 and RSwM3 algorithms have comparable
runtimes which are much smaller than those of RSwM1. This suggests that, while
RSwM1 is vastly simpler to implement, the more complex RSwM2 and RSwM3 are
more efficient.

Figure 4. Adaptive Algorithm Error Comparison on Equation 31,
Equation 33, and Equation 35. Comparison of the rejection sampling with memory algorithms on examples 1-3. Along the x-axis
is  which is the user chosen local error tolerance. The y-axis is
the average l2 error along the path. These values are taken as themean for 100 runs. Both axes are log-scaled.

5.3. Summary of numerical correctness and efficiency. The results of thenumerical experiments can be summarized as follows. From subsection 5.1 we seethat by simulating Equation 31 Equation 33, and Equation 35 by ESRK1 with thetime-stepping algorithms RSwM1, RSwM2, and RSwM3, the resulting Brownain
paths have the appropriate sample statistics. This indicates that in all of thesecases the algorithms produce the correct results. Notably, this shows that RSwM2,
whose correctness we could not guarantee from its derivation, generates results
which are can be statistically correct at low tolerances. Given that the KolmogorovSmirnov test is a stringent test for normality, this indicates that the edge case notappropriately dealt with in RSwM2 is a rare enough occurrence that the algorithm

NATURAL EMBEDDINGS AND REJECTION SAMPLING WITH MEMORY

17

Figure 5. Adaptive Algorithm Timing Comparison. Comparisonof the rejection sampling with memory algorithms. Along the xaxis is  which is the user chosen local error tolerance. In the y-axistime is plotted (in seconds). The values are the elapsed time for a1000 Both axes are log-scaled.

can likely be statistically indistinguishable from being correct in non demanding
circumstances.
From subsection 5.2 we see that the RSwM2 and RSwM3 algorithms are themost efficient on the test equations Equation 31 and Equation 33 in terms of runtime and accuracy. Together, the results show that the three algorithms are correct
and that there exists a tradeoff between them. The simplest algorithm to implement is RSwM1, but on some classes of problems it can perform slower than theother algorithms. RSwM2 is by a slight margin the fastest algorithm but is only
psudo-correct. RSwM3 is an extension of RSwM2 which, at the cost of adding
correctness, adds complexity in the implementation and runs with only slightly
decreased efficiency.
6. Numerical experiments with systems of SDEs. Next we wished to show
the applicability of the adaptive algorithms for solving large systems of nonlinear
equations. These examples better match user cases, though no analytical solution
exists for the error analysis. Instead, we wish to show how the error tolerance parameters allow for specifying the granularity of the numerical approximation, anddemonstrate that this algorithm can be used to solve problem which would be computationally intractable without the adaptivity.

18

CHRISTOPHER RACKAUCKAS AND QING NIE

Figure 6. Solution to the Lorenz system Equation 39 on t ∈ [0, 10]with additive noise and parameters α = 10, ρ = 28, σ = 3,and β = 8/3 at varying tolerances. The system was solved using ESRK+RSwM3 with the relative tolerance fixed at zero andvarying absolute tolerances.
6.1. Control of numerical granularity. To test the ability for the error tolerance
parameters to control the granularity of the solution, the ESRK1 algorithm with
RSwM3 was used to solve the well-known Lorenz system with additive noise. The
parameters were chosen so the deterministic solution is chaotic and generates thewell-known Lorenz attractor. The relative tolerance was fixed at zero and equation
was solved at a variety of absolute tolerances. The solution was plotted as Figure 6.
One can see that as the tolerance is decreased, the solution becomes more exact.
Thus the RSwM tolerance parameters allow the users to adjust the trade-off between
the exactness of the numerical approximation and computational time on nonlinear
systems of equations.
6.2. Large nonlinear system on cell differentiation. Next we demonstrated
the ability of the ESRK+RSwM algorithms to allow for solving large systems ofnonlinear equations. As a test, we chose a model of stochastic cell differentiationin the epithelial-to-mesenchymal transitions of somatic cells. The system of 19
equations is given in Appendix F. This system serves as an ideal test case foradaptive algorithms since the stochastic switches give large intervals of relative
stability with small intervals of extreme stiffness. Thus the ability to be able toautomatically adjust the timestep depending on whether a stochastic transitionis occurring is fundamental for making the system computationally tractable. In
Figure 7 we show the solution on a time interval from 0 to 500 for two indicators of
the cell states, Ecad and Vim, and the corresponding stepsizes the ESRK+RSwM
algorithm used to solve the equations throughout the simulation. This figure showsNATURAL EMBEDDINGS AND REJECTION SAMPLING WITH MEMORY

19

that the ESRK+RSwM algorithm is able to reproduce the scientific results, andthat the adaptive algorithm is able to adjust the timestep between 10−4 and 10−11
to ensure accuracy and stability.
Unlike the deterministic case, whether a given solution is numerically stable for agiven stepsize is dependent on the random Brownian trajectory and thus one cannot
as easily classify the methods as stable for a given stepsize on a given problem.
Instead, we define a method as sufficiently stable at a stepsize h if the method has
no unstable trajectories in a sufficiently large ensemble. The reason for this choice
is because it gives a measure for how fast a Monte Carlo simulation can correctly
be computed. The condition that no trajectories can fail is because if the noisiest(and thus most unstable) trajectories are always discarded from the ensemble, theresulting Monte Carlo solution would be biased. This is a major concern for manymodels in practice. For example, the chosen epithelial-to-mesenchymal transition
model has the most numerical instability during the stochastic switching events, butthese events are the quality of interest and thus must be accurate in the resulting
simulations.
Therefore, to test how efficiently the algorithms could produce statistically correct results on a stiff system, we solved cell model from Appendix F 10,000 times
via the Euler-Maruyama, Runge-Kutta Milstein (RK-Mil) [15], and Rößler SRI
methods using increments of h = 2−i and calculated the number of runs which
diverged (failed) and the elapsed time for the 10,000 runs (note that the algorithm
exited and denoted a run as diverged when encountering a NaN). These results arecompared to the adaptive algorithm in Table 2. The adaptive algorithm solved
10,000 simulations with no failures in 187 seconds. The comparable value for fixed
timestep methods, the shortest time for which there were no failures, was the EulerMaruyama algorithm with h = 2−20 which took 2286 seconds. This shows that theadaptive algorithm performed 12.28x as fast as the fastest fixed time-stepping algorithm, indicating that the advantages of adaptive time-stepping far out-weighted
the overheads of the adaptive algorithm. Note that this testing method also far
underestimates the advantages of the adaptive method: from Figure 7 we see thatthe most stiff portions are during the middle of stochastic transition events which
would not happen during 1 second simulations. However, testing multiple solutionsof the non-adaptive algorithms at such a low h in order to converge on the fullt ∈ [0, 500] was not computationally feasible. Indeed, one large advantage of theadaptive method is that the tests to determine the appropriate h are computationally expensive themselves. Also note that the tolerance for the adaptivity had to be
set low in order to ensure stability. The advantage of the adaptive method could bemade even larger by using a numerical method with a larger stability region. The
application of the RSwM algorithms to more stable solvers for stiff systems andhigher weak order solvers is a subject for future research.
7. Implementation issues. For deploying these algorithms, there are a few extra
implementation issues that should be addressed. First of all, the choice of q is due
to calculations for the optimal stepsize for deterministic equations [7] and thus does
not directly apply to stochastic systems. We instead use the deterministic q as a
guideline. We heuristically modify it by noting that due to the extra overhead of
the stack implementations, it may be wise to over-accelerate the stepsize changes
and thus choose a higher power. In practice, we found20

CHRISTOPHER RACKAUCKAS AND QING NIE

Figure 7. Stochastic cell differentiation model solutions. (A)Timeseries of the concentration of [Ecad]. The solution is plotted once every 100 accepted steps due to memory limitations. (B)Timeseries of the concentration of [V im]. The solution is plotted
once every 100 accepted steps due to memory limitations. (C) Accepted h over time. The h values were taken every 100 accepted
steps due to memory limitations. (D) Elapsed time of the EulerMaruyama and ESRK+RSwM3 algorithms on the stochastic cell
model. Each algorithm was used to solve the model on t ∈ [0, 1]
10,000 times. The elapsed time for the fixed timestep methods forgiven h’s are shown as the filled lines, while the dashed line is theelapsed time for the adaptive method. The red circles denote theminimal h and times for which the method did not show numericalinstability in the ensemble.

2
1
q=
,
(23)
γe
to be an efficient choice. Future research should investigate this issue in more detail.
In addition, it is wise in practice to ensure that q takes acceptable values [11].
We found that the following tends to work well in practice:q = min(qmax, max(q, qmin)),(24)

where qmax and qmin are parameters set by the user. While for deterministic
systems it’s suggested that qmax is between 4 and 10 [11], we note that our algorithm has a higher overhead than in the deterministic case, and thus further
restricting q to reduce the number of rejections had the potential to increase theefficiency to of the algorithm. This was tested in Table 3 where we found a suitable

NATURAL EMBEDDINGS AND REJECTION SAMPLING WITH MEMORY

Euler-Maruyama
Fails (/10,000) Time (s)

∆t
2−16
2−17
2−18
2−19
2−20
2−21

137
39
3
1
0
0

Runge-Kutta Milstein
Fails (/10,000) Time (s)133.35269.09580.141138.412286.35
4562.20

131
26
6
1
0
0

211.92428.28861.011727.913439.906891.35

21

Rößler SRI
Fails (/10,000) Time (s)
78
17
0
0
0
0

609.271244.06
2491.374932.709827.16
19564.16

Table 2. Fixed timestep method fails and runtimes. The fixed
timestep algorithms and ESRK+RSwM3 algorithms were used tosolve the stochastic cell model on t ∈ [0, 1] 10,000 times. Failures
were detected by checking if the solution contained any NaN values.
During a run, if any NaNs were detected, the solver would instantly
end the simulations and declare a failure. The runtime for theadaptive algorithm (with no failures) was 186.81 seconds.

qmax
1 + 2−5
1 + 2−4
1 + 2−3
1 + 2−2
1 + 2−1
1 + 201 + 21
1 + 22

Example 1
Time (s)
Error
37.00
34.73
49.14
39.33
38.22
82.76
68.16
48.23

2.57e-8
2.82e-8
3.14e-8
3.59e-8
3.82e-8
4.41e-89.63e-8
1.01e-7

Example 2Time (s)Error

Example 3
Time (s)
Error

60.87
32.40
132.3333.90
159.94
34.41
33.98
33.97

67.71
66.68
65.94
66.33
68.16
568.2287.50
69.78

2.27e-7
3.10e-7
8.85e-7
1.73e-6
2.58e-6
3.58e-6
6.06e-6
9.74e-6

3.42e-93.43e-93.44e-93.44e-93.44e-93.44e-93.22e-93.44e-9

Cell Model
Time (s)229.83
196.36186.81205.57249.77337.99418.78571.59

Table 3. qmax determination tests. Equation 31, Equation 33,
and Equation 35 were solved using the ESRK+RSwM3 algorithmwith a relative tolerance of 0 and absolute tolerance of 2−14 . The
elapsed time to solve a Monte Carlo simulation of 100,000 simulations to T = 1 was saved and the mean error at T = 1
was calculated. The final column shows timing results for using
ESRK+RSwM3 on the stochastic cell model from F solved withthe same tolerance settings as in subsection 6.2 to solve a Monte
Carlo simulation of 10,000 simulations.

value qmax = 1.125 which we used as the default throughout this paper (except in
subsection 5.2 where qmax = 10 was used in order to better evaluate the ability
to accurately unravel the stack). Further research should look into the usage of
PI-controlled stepsize choices to further reduce the number of rejections.
Also, one should note that elements entering the stack should be sanitized inorder to eliminate issues due to round-off error. In some cases the mathematicalalgorithm may specify very small intervals to go on the stack. This is particularly
the case when q,qtmp , or qK are close to 1 when the step size is very small. When
this item from the stack is eventually used to calculate qtmp , there is the possibilityof overflow which can lead to NaNs and stall the algorithm. To avoid this issue, it22

CHRISTOPHER RACKAUCKAS AND QING NIE

is wise to discard any interval smaller than some size. We found 10−14 to be a good
cutoff.
Lastly, it is common for many adaptive implementations to give a heuristic determination of an initial stepsize. Using the simple approximationthat with 99%√√
probability a N (0, σ) random variable is in the interval (−3 σ, 3 σ), we adapt
the scheme from [11] to take in account this a priori estimate, adjust due to thesymmetry of the normal distribution, assume the next error power is 12 larger, andarrive at the following 3.
Algorithm 3 Initial h Determination
1: Let d0 = kX0 k
2: Calculate f0 = f (X0 , t) and σ0 = 3g(X0 , t)
3: Let d1 = kmax(|f0 + σ0 |, |f0 − σ0 |)k
4: if d0 < 10−5 or d1 < 10−5 then5:
Let h0 = 10−6
6: else7:
Let h0 = 0.01(d0 /d1 )8: end if
9: Calculate an Euler step: X1 = X0 + h0 f0
10: Calculate new estimates: f1 = f (X1 , t) and σ0 = 3g(X1 , t)
11: Determine σ1M = max(|σ0 + σ1 |), |σ0 − σ1 |)12: Let d2 = kmax(|f1 − f0 + σ1M |, |f1 − f0 − σ1M |)k/h0
13: if max(d1 , d2 ) < 10−15 then14:
Let h1 = max(10−6 , 10−3 h0 )15: else
16:
Let h1 = 10−(2+log10 (max(d1 ,d2 ))/(order+0.5)
17: end if
18: Let h = min(100h0 , h1 )Note that k · k is the norm as defined in Equation 19 with Equation 20. In
practice we found this to be a conservative estimate which would “converge” tosome approximate stepsize in < 10 steps, making it useful in practice.
8. Conclusion and discussion. In this paper we have developed a class of adaptive time-stepping methods for solving stochastic differential equations. The approach is based on a natural creation of high Strong-order embedded algorithms(with implicit error calculations) and generalized rejection sampling techniques.
These algorithms are highly efficient: the error estimators require no extra evaluations of the drift and diffusion terms, while the rejection sampling with memoryalgorithms are O(1) for rejecting and interpolating the Brownian motion, with
RSwM2 and RSwM3 allowing for an unconstrained choice of timestep. Throughnumerical testing, we demonstrated the algorithm’s ability to adapt the error to
the user-chosen parameter , the prescribed local error threshold. We also haveshown the efficiency of the algorithms and their ability to scale to efficiently solve
large stiff systems of equations.
Our results show that rejection sampling with memory can be used to speed up
the solution of large systems of equations, particularly those coming from biologicalmodels. Many of these models are created to show interesting behavior which isdependent on the stochasticity, such as stochastic switching, but can be difficult toNATURAL EMBEDDINGS AND REJECTION SAMPLING WITH MEMORY

23

simulate due to the stiffness associated with these properties. Using such an adaptive algorithm allows the simulation to focus its computational time to accurately
capture the rare but important events. Further research into the application ofthese methods with stiff solvers could further improve the effectiveness. Also, theadaptive methods reduce the practical amount of time to conduct such experiments
since one major time-consuming activity can be finding a small enough stepsize sothat large Monte Carlo simulations can run without any path diverging and biasing the result. Further research should look into using rejection sampling withmemory combined with higher order weak convergence methods to improve thecomputational efficiency when investigating moment properties from Monte Carloexperiments.
One of the unique features of the adaptive methods is that the embedded algorithm does not require an explicit construction of the order 1.0 method. It also does
not have any extra requirements on the order 1.5 algorithm. This allows one to use
any order 1.5 Rößler SRK algorithm, and thus as new coefficients for these algorithms are found to have “better” properties, such as improved stability or lower
highest-order truncation error, our construction shows the existence of embedded
versions for adaptive algorithms by default.
There are many implementation details that can be examined in order to more
efficiency utilize these algorithms. For example, in the algorithm RSwM3, one could
simulate discarding the second stack S2 by explicitly controlling the pointer on anarray and thus not have to deallocate memory every accepted timestep. Further
research can investigate such implementation details for even more computationallyefficient versions. Note that the calculation of the tuple for the stack requires only 3
floating point calculations and the generation of the new random numbers requires
only 3x2 floating point calculations. Thus in total the algorithm requires only
≈ 10 floating point operations per rejection. Therefore the adaptive time-stepping
algorithms require a minimal number of floating point operations, making them
have a theoretically small impact on the computational effort.
The time-stepping method can be easily extended to other types of equations.
Notice that as stated this algorithm controls the strong or pathwise error. There
are many cases where one may be interested in the weak or average error of a Monte
Carlo simulation. Our algorithm can be used to create a weak error estimate as well.
To do so, instead of solving N independent size n systems, combine the systems
into an equivalent system of size nN . Notice that by the definition of e, the error
estimate obtained by this algorithm will be an average error over all simulations,and thus give a weak estimate. Using this weak estimate along with the RSwM
algorithms on high weak order methods could be a path of future research.
At its core, the time-stepping algorithms are methods to sub-sample a continuoustime stochastic process as needed and effectively reproduce the sample properties
at every point. Thus for any continuous-time martingale which defines a type ofdifferential equation, this method can be applied by simply knowing the sample
statistics of a bridge-type problem. Thus no part of the derivation required thatthe underlying process was an Ito process. Therefore, if one uses an embedded
pair of methods for Stratonovich SODEs, the adaptive time-stepping algorithms
still apply. Also, this includes problems like stochastic partial differential equations(SPDEs). For example, for SPDEs with space-time white noise, if one imposes a
space-discretization then one arrives at a system of SODEs for which the Brownian
bridge statistics are the same as any standard system of SODEs and thus our24

CHRISTOPHER RACKAUCKAS AND QING NIE

methods apply. If one needs to adapt the spatial discretization, then this also
follows the sample properties of the Brownian bridge and thus one can upsample
the process along space until the desired error is met. For a Hilbert space valued
Brownian motion like a Q-Wiener process, a Brownian bridge-like problem will need
to be solved but the same general algorithm holds.
Notice that since the stepping algorithms only require the values on the stack,the RSwM form of adaptivity allows one to solve an equation while only storing
the current values and the Brownian stack values. This is important for solving
large systems of SPDEs which become memory-bound if one attempts to store
the solution of the Brownain path at every timepoint. During our study of theimplementation we saw the stack reaching maximum sizes of 20, and thus only useda modest amount of memory. For problems which are more memory bound, one
can effectively decrease the maximum stack size by changing around some of theparameters, for example lowering the q exponent or decreasing qmax.
This method also provides an effective way to develop adaptive algorithms forvariable-rate Markovian switching and jumps [18]. These problems are usually difficult to numerically simulate because one must effectively estimate the times forswitching or else incur a large penalty due to the discontinuity. One can normallyestimate whether a jump has occurred in an interval to a certain degree of accuracy, and using our adaptive algorithm one can hone in on the time at which thejump occurred, effectively solve the continuous problem to that time, and apply thejump. Therefore, the embedded time-stepping method along with the new rejection
sampling algorithm developed in this work provides a general framework which one
could build on for adaptive algorithms in solving many other stochastic problems.
Lastly, we wish to note that implementations of these algorithms are being released as part of a package DifferentialEquations.jl. DifferentialEquations.jl is a Julia library for solving ODEs, SDEs, DDEs, DAEs, and certain classes of PDEs andSPDEs, using efficient solvers in an easy-to-use scripting language. This package,developed by the authors, is freely available and includes additional functionalities such as parallelizing Monte Carlo experiments using the discussed methods on
HPCs.
Acknowledgments. This work was partially supported by NIH grants P50GM
76516 and R01GM107264 and NSF grants DMS1562176 and DMS1161621. This
material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1321846, the National Academiesof Science, Engineering, and Medicine via the Ford Foundation, and the National
Institutes of Health Award T32 EB009418. Its contents are solely the responsibilityof the authors and do not necessarily represent the official views of the NIH.
The authors would like to thank Tom Breloff for his responsiveness and development of the Plots.jl Julia package which was used extensively in this publication.
We would also like to thank Catherine Ta and Tian Hong for sharing their code forthe stochastic cell model. We thank Gerrit Ansmann and the reviewers for theircomments. Lastly, I would like to thank the Julia community for their continuing
help and support.
Appendix.
A. Order conditions for Rößler-SRI methods. The coefficients A(i) ,
B (i) , β (i) , α must satisfy the following order conditions to achieve order .5:

NATURAL EMBEDDINGS AND REJECTION SAMPLING WITH MEMORY
T

1. αT e = 1T
2. β (1) e = 1

25T

3. β (2) e = 0
T
4. β (3) e = 0

5. β (4) e = 0

additionally, for order 1:TT

1. β (1) B (1) e = 0T
2. β (2) B (1) e = 1

3. β (3) B (1) e = 0
T
4. β (4) B (1) e = 0

and lastly for order 1.5:
1. αT A(0) e = 21
2. αT B (0) e = 1
2
3. αT B (0) e = 23
T
4. β (1) A(1) e = 1
T
5. β (2) A(1) e = 0
T
6. β (3) A(1) e = −1
T
7. β (4) A(1) e = 0
2T
8. β (1) B (1) e = 116.

9. β (2)
10. β

T

(3)T
(4)T

11. β
T
12. β (1)T
13. β (2)T
14. β (3)T
15. β (4)B (1) e

2

=0

(1)

2

= −1

B

(1)

e

2

B e =2
B (1) B (1) e = 0B (1) B (1) e = 0B (1) B (1) e = 0B (1) B (1) e = 1

1 (1)T  (1)  (0)  1 (3)T  (1)  (0) βAB e + βA
B e =023where f, g ∈ C 1,2 (I × Rd , Rd ), c(i) = A(i) e, e = (1, 1, 1, 1)T [22].

B. Order conditions
for Rößler-SRA methods. The coefficientsA(i) , B (i) , β (i) , α must satisfy the conditions for order 1:
1. αT e = 1

TT

3. β (2) e = 0

2. β (1) e = 1

and the additional conditions for order 1.5:
1. αT B (0) e = 1
2. αT A(0) e = 12

3. αT B (0) e

2

= 23

T

4. β (1) c(1) = 1
T
5. β (2) c(1) = −1

where c(0) = A(0) e with f ∈ C 1,3 (I × Rd , Rd ) and g ∈ C 1 (I, Rd ) [22]. From these
conditions he proposed the following strong order 1.5 scheme found in Table 4 known
as SRA1.

26

CHRISTOPHER RACKAUCKAS AND QING NIE

0

c(0)

A(0)

B (0)

T

(1)T

α

β

β

3
4

(2)T

3
4
1
3

(a) Legend Table

1
2
2
3

1

0

-1

1

(b) Coefficients Table

Table 4. SRA1. Table (a) shows the legend for how the numbers in in Table (b) correspond to the coefficient arrays/matricesc(i) , A(i) , B (i) , α, and β (i) . Note that the matrices A(i) and B (i) arelower triangular since the method is explicit.

C. Wiktorsson iterated stochastic integral approximations. In the Ito expansions for the derivation of SRK methods of order greater than 1.0, iterated
stochastic integrals appear. We denote these as:

Z tn+1
I(1)

=

dWs ,

(25)

tn
Z tn+1 Z s

I(1,1)=

dWu dWs ,tn
tn
Z tn+1 Z

I(1,1,...,1k )

Z s

...

=
tn

(26)

dWu1 . . . dWs .

(27)

tn

Note that for our purposes we are using a single Brownian path, though thisextends to multiple Ito dimensions. For a timestep h = tn+1 −tn , the approximation
due to Wiktorsson [25, 22] is as follows:I(1,1)=

I(1,1,1)=I(1,0)=

1 2(I − h),2 (1)1 3(I(1) − 3hI(1) ),
6 

1
1
h I(1) + √ ζ .
2
3

(28)
(29)
(30)

I(1) is one timestep of the Brownian path Wt , which is referred to as ∆W ∼N (0, h). ζ ∼ N (0, h) is a standard normal random variable which is independentof ∆W . If we collect all of the ζ from [0, T ], then its cumulative sum is itself anapproximation to a Brownian path which we denote Zt with discrete steps ∆Z = ζ.
Zt is independent of Wt and is the second Brownian path considered in the RSwM
algorithms.

NATURAL EMBEDDINGS AND REJECTION SAMPLING WITH MEMORY

27

D. RSwM2 algorithm specification. The algorithm for RSwM2 is specified as
Algorithm 4.

Algorithm 4 RSwM2
1: Set the values , hmax , T
2: Set t = 0, W = 0, Z = 0, X = X0
3: Take an initial h, ∆Z, ∆W ∼ N (0, h)
4: while t < T do
5:
Attempt a step with h, ∆W , ∆Z to calculate Xtemp according to (2)6:Calculate E according to (9)Update q using (21)if (q < 1) then
. % Reject the Step
9:
Take ∆W̃ ∼ N (q∆W, (1 − q)qh) and ∆Z̃ ∼ N (q∆Z, (1 − q)qh)10:
Calculate ∆W = ∆W − ∆ W̃ and ∆Z = ∆Z − ∆Z̃
11:Push (1 − q) h, ∆W , ∆Z into stack S
12:Update h := qh
13:Update ∆W := ∆W̃ , ∆Z := ∆Z̃14:else. % Accept the Step
15:
Update t := t + h, W := W + ∆W , Z := Z + ∆Z, X = X + Xtemp
16:
Update c := min(hmax , qh), h := min(c, T − tn )17:
Set hs = 0, ∆W = 0, ∆Z = 018:while S is not empty do
19:
Pop the top of S as L
20:
if (hs + L1 < h) then
. % Temporary not far enough21:
Update hs := hs + L1 , ∆W := ∆W + L2 , ∆Z := ∆Z + L3
22:else
. % Final part of step from stacks
23:
Set qtmp = h−hL124:
Let ∆W̃ ∼ N (qtmp L2 , (1 − qtmp )qtmp L1 )25:
Let ∆Z̃ ∼ N (qtmp L3 , (1 − qtmp )qtmp L1 )26:
Push ((1 − qtmp )L1 , L2 − ∆W̃ , L3 − ∆Z̃) onto S
27:
Update ∆W := ∆W + ∆W̃ , ∆Z := ∆Z̃ , hs := hs + qtmp L1
28:
end if29:
end while. % Update for last portion to step. Note zero if final part is from stack30:
if (h − hs is not zero) then31:
Let ηW , ηZ ∼ N (0, h − hs )32:
Update ∆W = ∆W + ηW , ∆Z = ∆Z + ηZ
33:
end if34:
end if
35: end while
7:
8:

28

CHRISTOPHER RACKAUCKAS AND QING NIE

E. Example equations. The three example equations are:Example 1.
dXt = αXt dt + βXt dWt , X0 =

1
,
2

(31)

1
1
and β = 20
. Actual Solution:
where α = 10

Xt = X0 e

2

β− α2



t+αWt

.

(32)

Example 2.
dXt = −110

2

sin (Xt ) cos3 (Xt ) dt +Actual Solution:
Xt = arctan11
cos2 (Xt ) dWt , X0 = ,
10
2


1
Wt + tan (X0 ) .
10

(33)

(34)

Example 3.



1
1
β
αβ
−
Xt dt + √dWt , X0 = ,dXt = √
2
(1
+
t)
2
1+t
1+t
1
1
where α = 10 and β = 20 . Actual Solution:

(35)

1
β
X0 + √(t + αWt ) .
1+t
1+t

(36)

dXt

= α(Yt − Xt ) + σdWt ,

(37)

dYt

= Xt (ρ − Zt ) − Yt + σdWt ,

(38)

dZt

= Xt Yt − βZt + σdWt ,(39)

Xt = √
Example 4.

where X0 = Y0 = Z0 = 0, α = 10, ρ = 28, σ = 3, and β = 8/3.
F. Stochastic cell differentiation model. The stochastic cell differentiation
model is given by the following system of SDEs which correspond to a chemical
reaction network modeled via mass-action kinetics with Hill functions for the feedbacks. This model was introduced in [12] to model stochastic transitions between
Epithelial and Mesenchymal states.
A=(([T GF ] + [T GF 0]) /J0snail )n0snail + ([OV OL2] /J1snail )n1snail

d [snail1]tdt=

k0snail + ksnaild [SN AIL]dt
d [miR34]
dt
d [SR]
dt
d [zeb]dt−(([T GF ] + [T GF 0]) /J0snail )n0snail(1 + A) (1 + [SN AIL] /J2snail )kdsnail ([snail1] − [SR]) − kdSR [SR]=kSN AIL ([snail1] − [SR]) − kdSN AIL [SN AIL]=

kO34 +k34−

1 + ([SN AIL] /J134 )n134 + ([ZEB] /J234 )n234
kd34 ([miR34] − [SR]) − (1 − λSR ) kdSR [SR]=T k (KSR ([snail1] − [SR]) ([miR34] − [SR]) − [SR])=k0zeb + kzeb([SN AIL] /J1zeb )n1zeb
1 + ([SN AIL] /J1zeb )n1zeb + ([OV OL2] /J2zeb )n2zeb

NATURAL EMBEDDINGS AND REJECTION SAMPLING WITH MEMORY−[zeb] −

kdzeb

5
X

!
C5i [ZR]−=d [miR200]dt=

k0200 +−

kd200

[zeb] −kZEB

5X

kdZRi C5i [ZRi ]

i=1

i=1

d [ZEB]dt

5X!
C5i [ZRi ]− kdZEB [ZEB]

i=1

k200

1 + ([SN AIL] /J1200 )n1200 + ([ZEB] /J2200 )n2200!
5
X
i[miR200] −iC5 [ZRi ] − [T R]
i=1

−

5
X

(1 − λi ) kdZRi C5i i [ZRi ] − (1 − λT R ) kdT R [T R]

i=1d [ZR1 ]
dt

=

Tk

K1

[miR200] −5X!
iC5i [ZRi ] − [T R]

i=1[zeb] −

5X!
C5i [ZRi ]!
− [ZR1 ]i=1d [ZR2 ]dt=

d [ZR3 ]
dt

=

Tk

K2

[miR200] −5X!
iC5i [ZRi ] − [T R]!
[ZR1 ] − [ZR2 ]

i=1

Tk

K3

[miR200] −5X!
iC5i [ZRi ] − [T R]!
[ZR1 ] − [ZR3 ]i=1d [ZR4 ]
dt

=

Tk

K4

[miR200] −5X!
iC5i [ZRi ] − [T R]!
[ZR1 ] − [ZR4 ]i=1d [ZR5 ]
dt

=

Tk

K5

[miR200] −5X!
iC5i [ZRi ] − [T R]!
[ZR1 ] − [ZR5 ]i=1d [tgf ]dt
d [T GF ]dt
d [T R]dt=

ktgf − kdtgf ([tgf ] − [T R]) − kdT R [T R]=

k0T GF + kT GF ([tgf ] − [T R]) − kdT GF [T GF ]!
!
5
X
T k KT R [miR200] −iC5i [ZRi ] − [T R] ([tgf ] − [T R]) − [T R]=

i=1

d [Ecad]
dt

=

k0E +

B

=

kV 1

d [V im]dt
d [OV OL2]
dt

=
=

kE1
1 + ([SN AIL] /J1E )n1E([SN AIL] /J1V )n1V
n1V

+

kE2
1 + ([ZEB] /J2E )n2E

+ kV 2− kdE [Ecad]([ZEB] /J2V )n2V

1 + ([SN AIL] /J1V )
1 + ([ZEB] /J2V )n2V
B
k0V +− kdV [V im](1 + [OV OL2] /J3V )1
k00 + k0
− kdO [OV OL2]
1 + ([ZEB] /J0 )n0

where
5
X

iC5i [ZRi ]=5 [ZR1 ] + 20 [ZR2 ] + +30 [ZR3 ] + 20 [ZR4 ] + 5 [ZR5 ] ,=5 [ZR1 ] + 10 [ZR2 ] + 10 [ZR3 ] + 5 [ZR4 ] + [ZR5 ] .

i=1
5
X

C5i [ZRi ]i=1

The parameter values are given in Table 5.

29

30

CHRISTOPHER RACKAUCKAS AND QING NIE

Parameter

Value

Parameter

Value

Parameter

Value

Parameter

Value

J1200
J2200
J134
J234
JO
J0snailJ1snailJ2snailk0snail
n1200
n2200
n134
n234
nO
n0snail
kO
kzeb
kdZR1
kdZR5
kdSR

3
0.2
0.15
0.35
0.9
0.6
0.5
1.8
0.0005
3
2
2
2
2
2
1.2
0.06
0.5
0.5
0.9

J1E
J2E
J1V
J2V
J3V
J1zeb
J2zeb
K1
k0zeb
n1snail
n1E
n2E
n1V
n2V
n2zeb
k200
kT GF
kdZR2
kdO
kdE

0.1
0.3
0.4
0.4
2
3.5
0.9
1
0.003
2
2
2
2
2
6
0.02
1.5
0.5
1.0
0.05

K2
K3
K4
K5
KT R
KSR
T GF 0
Tk
λ1
λ2
λ3
λ4
λ5
λSR
λT R
k34
kSN AIL
kdZR3
kd200
kdV

1
1
1
1
20
100
0
1000
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.01
16
0.5
0.035
0.05

k0O
kO200
kO34
kdsnail
kdtgf
kdzeb
kdT GF
kdZEB
k0T GF
k0E
k0V
kE1
kE2
kV 1
kV 2
ktgf
kZEB
kdZR4
kd34

0.35
0.0002
0.001
0.09
0.1
0.1
0.9
1.66
1.1
5
5
15
5
2
5
0.05
16
0.5
0.035

Table 5. Table of Parameter Values for the Stochastic Cell Model.

REFERENCES[1] R. L. Burden and J. D. Faires, Numerical analysis, 9th edn, Brooks/Cole, Cengage Learning,
Boston, MA, 2011.
[2] K. Burrage and P. Burrage, General order conditions for stochastic runge-kutta methods forboth commuting and non-commuting stochastic ordinary differential equation systems, Appl.
Numer. Math., 28 (1998), 161–177.
[3] K. Burrage and P. M. Burrage, High strong order explicit runge-kutta methods for stochastic
ordinary differential equations, Applied Numerical Mathematics, 22 (1996), 81–101.
[4] P. M. Burrage, Runge-kutta methods for stochastic differential equations, Thesis, The University of Queensland Brisbane, 1999.
[5] P. M. Burrage and K. Burrage, A variable stepsize implementation for stochastic differential
equations, SIAM Journal on Scientific Computing, 24 (2002), 848–864.
[6] J. R. Cash and A. H. Karp, A variable order runge-kutta method for initial value problems
with rapidly varying right-hand sides, ACM Transactions on Mathematical Software, 16(1990), 201–222.
[7] F. Ceschino, Modification de la longueur du pas dans l’integration numerique parles methodesa pas lies, Chiffres, 4 (1961), 101–106.
[8] J. R. Dormand and P. J. Prince, A family of embedded runge kutta formulae, Journal ofComputational and Applied Mathematics, 6 (1980), 19–26.
[9] J. G. Gaines and T. J. Lyons, Variable step size control in the numerical solution of stochastic
differential equations, SIAM Journal on Applied Mathematics, 57 (1997), 1455–1484.
[10] A. Ghasemi and S. Zahediasl, Normality tests for statistical analysis: A guide for nonstatisticians, Int J Endocrinol Metab, 10 (2012), 486–489.
[11] E. Hairer, S. P. N orsett and G. Wanner, Solving Ordinary Differential Equations I, 2nd rev.
edn, Springer series in computational mathematics. Springer-Verlag, Berlin, New York, 1993.
[12] T. Hong, K. Watanabe, C. H. Ta, A. Villarreal-Ponce, Q. Nie and X. Dai, An ovol2-zeb1
mutual inhibitory circuit governs bidirectional and multi-step transition between epithelial
and mesenchymal states, PLoS Comput Biol, 11 (2015), e1004569.
[13] S. M. Iacus, Simulation and Inference for Stochastic Differential Equations: With R Examples(Springer Series in Statistics), Springer Publishing Company, Incorporated, 2008.
[14] J. Kaneko, Explicit order 1.5 runge kutta scheme for solutions of ito stochastic differential
equations, S/urikaisekikenky/usho K/oky/uroku, 932 (1995), 46–60.

NATURAL EMBEDDINGS AND REJECTION SAMPLING WITH MEMORY

31[15] P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations,
Springer Berlin Heidelberg, 1992.
[16] H. Lamba, An adaptive timestepping algorithm for stochastic differential equations, Journalof Computational and Applied Mathematics, 161 (2003), 417–430.
[17] T. Liggett, Continuous Time Markov Processes: An Introduction, American MathematicalSociety, 2010.
[18] X. Mao and C. Yuan, Stochastic Differential Equations with Markovian Switching, Imperial
College Press, London, 2006.
[19] N. Hofmann, T. Muller-Gronbach and K. Ritter, Optimal approximation of stochastic differential equations by adaptive step-size control, Mathematics of Computation, 69 (2000),
1017–1034.
[20] B. Oksendal, Stochastic Differential Equations: An Introduction with Applications, Springer,
2003.
[21] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes 3rd
Edition: The Art of Scientific Computing, Cambridge University Press, 2007.
[22] A. Rößler, Runge kutta methods for the strong approximation of solutions of stochastic differential equations, SIAM Journal on Numerical Analysis, 48 (2010), 922–952.
[23] T. Ryden and M. Wiktorsson, On the simulation of iterated ito integrals, Stochastic Processes
and their Applications 91 (2001), 151–168.
[24] L. F. Shampine, Some practical runge-kutta formulas, Mathematics of Computation, 46(1986), 135–150.
[25] M. Wiktorsson, Joint characteristic function and simultaneous simulation of iterated ito integrals for multiple independent brownian motions, The Annals of Applied Probability, 11(2001), 470–487.

Received August 2016; revised December 2016.
E-mail address: contact@chrisrackauckas.com
E-mail address: qnie@math.uci.edu