Backward Simulation of Stochastic Process using a Time Reverse Monte Carlo method
Shinichi Takayanagi Yukito Iba

A Stochastic Analysis of Bike Sharing Systems
     
  Shuang Tao 
 School of Operations Research and Information Engineering 
 Cornell University

 293 Rhodes Hall, Ithaca, NY 14853 
  st754@cornell.edu  
 
Jamol Pender 
 School of Operations Research and Information Engineering 
 Cornell University

 228 Rhodes Hall, Ithaca, NY 14853 
  jjp274@cornell.edu  
 
 
    December 30, 2023
==============================================================================================================================================================================================================================================================================================================================





§ INTRODUCTION


In this paper, we discuss “backward simulation,” which traces a time-reversed path from a target region A to the initial configuration (Fig. <ref>). If the outputs of the original simulation (“forward simulation”) are easily restored from those obtained by backward dynamics, we can use backward simulation as a computational tool. In particular, the time required to calculate the probability to reach A from the initial configurations can be significantly reduced when the target region A is small but the initial distribution is broad.
An example is a computation of the probability that a typhoon will hit the Tokyo area exactly under a given stochastic model (Sect. <ref>).



It is, however, difficult to design backward dynamics with the desired properties. Specifically, consider the forward dynamics of a D-dimensional stochastic process X defined by 

    X_i+1 = g(X_i) + η_i,
 
where η_i is an independent noise that obeys an arbitrary distribution and the function g: ℝ^D→ℝ^D describes noiseless forward dynamics.
Then, a naïve way to derive a time-reversed equation is 
to rearrange Eq. (<ref>) as

    X_i = g^-1(X_i+1 - η_i).
 
Here, we assume that function g is a one-to-one and onto function and
denotes the inverse function of g as g^-1. We can construct a time-reversed path iteratively, using Eq. (<ref>) and the independent realization of η_i, starting from the target region. It defines an apparently natural candidate for backward dynamics.  

Surprisingly, this naïve method does not work as expected; it does not reproduce the correct probabilities defined by the forward simulation, and the calculation of factors required to correct the bias is often computationally expensive. This becomes clear in Sect. <ref>. Furthermore, the computation of g^-1 in Eq. (<ref>) is time-consuming and reduces the efficiency of the computation.    

The aim of this research is to draw attention to these facts and propose an algorithm that partially resolves the problem. We named this algorithm the time reverse Monte Carlo method (TRMC). TRMC is based on the ideas of sequential importance sampling (SIS) <cit.> and sequential Monte Carlo (SMC) <cit.>. We discuss TRMC based on SIS in Sects. <ref> and <ref> and its improved version based on SMC in Sect. <ref>. 
There have been several studies using “path reweighting” in computational chemistry <cit.>.
These studies used reweighting for different purposes.

TRMC based on SIS is tested for a stochastic typhoon model and the Lorenz 96 model in Sect. <ref>. The improved version with SMC is also tested in Sect. <ref> for simulations with a larger number of steps. In these examples, TRMC provides unbiased estimates of the probabilities without expensive computation. In Sect. <ref>, we discuss its relation to the Bayes formula, as well as the possible improvement and limitations of TRMC. 

Time-reversed dynamics itself was discussed in several studies <cit.>, mostly from a theoretical viewpoint. On the other hand, related computational problems are  
found in data science, especially in time-series analysis using 
state-space models <cit.>. Our problem can formally be regarded as a limiting case of the “smoothing” part of these algorithms, where only one observation (“target”) is available at the end of the time series. There are, however, important differences from our problem, which are discussed in Sect. <ref>.
Studies related to the statistical inference 
on a discrete state stochastic process, such as  
genepropagation <cit.> and 
information source detection <cit.> were also reported.
These studies, however, did not consider 
dynamical systems of continuous variables. 




§ FAILURE OF NAÏVE METHOD
 

Here, we provide a detailed discussion of the naïve method and its drawbacks, 
which form the motivation for our algorithm. 
Before providing details, we formulate the problem.
Let S_T = {0 = t_0 ≤ t_1 ⋯≤ t_N = T} be a partition of the interval [0, T], and let step size Δ t= t_i+1-t_i be a constant; x_i is used to represent the value of stochastic process X at time point t_i. The transition probability density from x_i to x_i+1 defined by Eq. (<ref>) is denoted as p( x_i+1 | x_i).
We consider an estimation of the probability 

that X_N hits a small target region A in the D-dimensional space.
The probability is formally written as

    P(X_N ∈ A) 
    	=
    	∫ dx_0:N1_x_N ∈ A{∏_i=0^N-1  p( x_i+1 | x_i) }
    	p(x_0)
    	,

where 1_x ∈ A is the indicator function that takes value 1 when x ∈ A, and 0 otherwise, and 
p(x_0) is the initial distribution of the forward simulation.
Hereafter,  dx_k:l indicates dx_kdx_k+1⋯ dx_l for k ≤ l. 

A naïve method is defined as a repeated simulation 
with a uniformly distributed initial condition in 
the target region A using Eq. (<ref>).
Initially, it appears sufficient to evaluate  as .
However, there are two problems with this naïve method.
First, the exact computation of g^-1 in Eq. (<ref>) is not easy.
Computing g^-1 using numerical root-finding techniques such as the Newton-Raphson method is computationally intensive and its severity increases as the dimension increases.

Second, this computation does not reproduce the correct probability  even with the exact g^-1.
To understand this problem, we show the difference between the forward simulation and the naïve method.
Let us define 

    Y_i=X_i-η_i-1 = g(X_i-1); i ∈ [1, …, N].

Using this definition, we can rewrite Eq. (<ref>) as 

    Y_i + η_i-1 = g^-1(Y_i+1) .


Equation (<ref>) can be simplified into

    Y_i =  g^-1(Y_i+1) - η_i-1.


The probability calculated using Eq. (<ref>) corresponds to  the equation

    ∫ dy_1:N dx_N1_x_N ∈ Ap̃_f( y_N | x_N)
    	{∏_i=1^N-1p̃( y_i | y_i+1)
    	}
    	p(g^-1(y_1))
    	,

where p̃_f( y_N | x_N ) is the transition probability density from x_N to y_N defined by Eq. (<ref>) with i=N and
p̃( y_i | y_i+1) is the transition probability density from y_i+1 to y_i defined by Eq. (<ref>).
An initial condition x_N is uniformly distributed in the target region A.

We have to introduce the Jacobian of function g so that Eq. (<ref>) is consistent with Eq. (<ref>).
To show this, Eq. (<ref>) is rewritten using equations

    p(x_i|x_i-1)
    	dx_i   =
    	|det(J_g^-1(y_i+1))|
    	p̃( y_i | y_i+1)
    	dy_i, 
    
    	i    ∈ [1,⋯, N-1] 
    
    	p(x_0)
    	dx_0   =
    	|det(J_g^-1(y_1))|
    	p( g^-1(y_1) )
    	dy_1,

where |det(J_g^-1(y_i))| is the absolute value of the Jacobian of function g^-1.
As a result, the probability P(X_N ∈ A) is calculated as

	
    P(X_N ∈ A) 
    		   =
    		∫ dy_1:N
    		dx_N
    		1_x_N ∈ A
    		J(y_1, …, y_N)
    		p̃_f( y_N | x_N)
    		{∏_i=1^N-1p̃( y_i | y_i+1)
    		}
    		p( g^-1(y_1) )
    		, 
    		
    
    		J(y_1, …, y_N)
    		   =
    		{∏_i=0^N-1|det(J_g^-1(y_i+1))|
    		}.



We can obtain the correct probability using Eq. (<ref>) instead of Eq. (<ref>).
The Jacobian J_g^-1 calculation is, however, computationally expensive.

We note that the factor J(y_1, …, y_N) goes to

    exp(
    -
    ∫_0^Tdivf(x_t) dt
    )

in the limit as Δ t → 0 when we assume that g(x) = x + f(x)Δ t.
The proof of Eq. (<ref>) is given in Appendix <ref>.
This shows that we must include factor  J(y_1, …, y_N) for unbiased estimation even in the limit of infinitesimal Δ t.
We can regard the factor written in Eq. (<ref>) as the change in  the infinitesimal volume along each path (Fig. <ref>).






§ TIME REVERSE MONTE CARLO METHOD
 
To overcome these difficulties, we propose the TRMC method.
TMRC essentially involves introducing simplified backward dynamics with a weight.
This weight enables the bias of estimators to be corrected.
First, we introduce a backward transition probability q( x_i+1→ x_i) from x_i+1 to x_i.
We can choose an arbitrary probability density q, 
while the computation efficiency strongly depends on it.
Once we introduce q( x_i+1→ x_i), 
we can rewrite Eq. (<ref>) as

	
    P(X_N ∈ A) 
    		=
    		∫ dx_0:N1_x_N ∈ A/V_A{∏_i=0^N-1  
    		q( x_i+1→ x_i)
    		W_i}
    		V_A
    		p(x_0)
    		,


where 

    W_i 
    	=
    			p( x_i+1 | x_i)
    	/
    		q( x_i+1→ x_i)

is the weight required to correct the bias of estimators 
and V_A is the volume of target region A.
Suppose p(x_0) is uniformly distributed on B ⊂ℝ^D;
p(x_0)=1/V_B1_x_0 ∈ B, V_B is the volume of B.
The efficiency of our algorithm does not depend on the factor V_A when V_B is considerably large.
This is the advantage of using our algorithm.

The algorithm consists of the following steps.

TRMC Algorithm

	
	
  * Draw M samples { x_N^(1), ⋯, x_N^(M)} from the uniform distribution in V_A.
	
  * Apply the following steps for j=1, …, M, and for i = N-1, …, 0.
	
		
  * Generate a sample from x^(j)_i+1 to x^(j)_i with transition probability 
		q( x^(j)_i+1→ x^(j)_i).
		
  * Calculate weight W_i^(j) using Eq. (<ref>).
	
	
  * Evaluate the unbiased estimates of probability P(X_N ∈ A) as 
	
    P(X_N ∈ A)
    	⋍1/M∑_j=1^M 
    	W^(j)	
    	,

	where the factor
	
    W^(j)
    	=
    	{∏_i=0^N-1  
    	W_i^(j)}
    	V_A
    	p(x_0^(j))

	is attached to each simulation path.

The inputs of our algorithm are the number of Monte Carlo paths M, 
the number of time steps N, 	
the initial distribution p(x_0),
the target region A,
and the transition probability density q.
When we actually perform the simulation on our computers, we take the logarithm of these weights to prevent numerical overflow.

This algorithm provides unbiased estimates of the desired probabilities.
The idea of this scheme is a kind of SIS <cit.>.
An advantage of our method is that we do not need to calculate g^-1 or their Jacobian matrices at each i.

The remaining problem involves determining the method for choosing the transition probability q( x_i+1→ x_i). 
The basic idea is to choose the backward dynamics that generates trajectories similar to the forward dynamics defined by Eq. (<ref>).
The similarity of the trajectory is measured by W^(j) in Eq. (<ref>).



§ IMPLEMENTATION FOR STOCHASTIC DIFFERENCE EQUATION
 
To give concrete examples of the transition probability 
q( x_i+1→ x_i), 
we assume the forward dynamics to be given in the following form:

    X_i+1 = X_i + f(X_i)Δ t + ϵ_i√(Δ t).

This corresponds to the case wherein g(x) = x + f(x)Δ t in Eq. (<ref>). 
The noise ϵ_i is assumed to be i.i.d. Gaussian noise with mean zero and the variance-covariance matrix Σ = σσ^T.
This class of equations appears in a wide range of problems in many different fields such as physics <cit.>, computational chemistry <cit.>, and mathematical finance <cit.>.

In this case, as a simple choice, we can use the following backward dynamics:

    X_i = X_i+1 - f(X_i+1)Δ t + ϵ_i√(Δ t).

This approximation corresponds to 
substituting f(X_i+1) for f(X_i) in Eq. (<ref>).

With this choice, weight W_i in Eq. (<ref>) takes the form

	
    W_i   =
    		
    			p( x_i+1 | x_i)
    		/
    			q( x_i+1→ x_i)
    		
    		=
    		exp[
    			-1/2(x_i+1 - x_i -  f(x_i)Δ t )^T(ΣΔ t)^-1(x_i+1 - x_i -  f(x_i)Δ t )
    			]
    		/exp[
    			-1/2(x_i+1 - x_i -  f(x_i+1)Δ t )^T(ΣΔ t)^-1(x_i+1 - x_i -  f(x_i+1)Δ t )
    			]
    		
       =
    		exp[
    		-( f(x_i+1) - f(x_i) )^TΣ^-1(
    		(x_i+1 - x_i)
    		-
    		Δ t/2(
    		f(x_i+1) + f(x_i)
    		)
    		)
    		]
    		.


As we show in the next section, the resultant algorithm is simple yet effective compared with the forward simulation when the target region A is smaller than the support of the initial distribution p(x_0).

We note that the factor ∏_i=0^N-1 W_i goes to

    exp(
    -
    ∫_0^Tdivf(x_t) dt
    )

in the limit as Δ t → 0.
The proof of Eq. (<ref>) is given in Appendix <ref>.
Note that Eq. (<ref>) coincides with Eq. (<ref>) derived from a different assumption.







§ EXPERIMENTAL RESULTS


We present the numerical results in this section.
Forward simulations (FS) are used to check the consistency and computational efficiency of our result.

Using forward and backward dynamics, we simulate sample trajectories 
x={x_1, ⋯, x_N } generated by each model and 
compute the probability P(X_N ∈ A) from M independent simulations.

We denote a standard error of TRMC to evaluate the computational efficiency by σ_s.
We also denote the standard error of FS by σ_s^F.
Using these variables, we define a relative value of variance by 

    ρ_1
    	=
    	(
    	σ_s^F/σ_s)^2
    	.

The factor ρ_1 indicates the computational efficiency only including the effect caused by the variance of estimators for a fixed sample size.
With this definition, more complex algorithms tend to be more efficient while they require more computational time.
Then, we also define another measure of the relative computational efficiency ρ_2 as

    ρ_2
    	=
    	ρ_1τ^F/τ,

where τ is the computational time in seconds of the simulation and τ^F is the computational time of FS in seconds. 
This efficiency is defined in the sense of the actual performance considering both the computational time and the variance of the resulting estimates.



 §.§ Stochastic typhoon model


The first example is a stochastic typhoon model <cit.>, which gives an example of risk estimation by the proposed method. 
The stochastic typhoon model was designed to reproduce the statistics of typhoons in the northwestern
part of the Pacific Ocean.
This is a four-dimensional model given by

    x_i+1   = x_i + v_i, 
    
    		v_i+1   = V(x_i+1)  +
    		w ( v_i - V(x_i) )
    		+ 
    		ϵ_i
    		,
     
    		V(x_i)     = a_0 + a_1x_ϕ, i + a_2sin x_λ, i + a_3sin^2 x_λ, i,

where we use a global coordinate system defined by the geographic longitude (ϕ) and latitude (λ).
We also define the two-dimensional position x = (x_ϕ, x_λ), speed v = (v_ϕ, v_λ) of a typhoon, and function V(x) = (V_ϕ(x), V_λ(x) ). 
w, a_0, a_1, a_2, and a_3 are constants.
The noise ϵ obeys a Gaussian distribution with mean zero and variances σ^2.

We fix w=0.93, a_0=(0.792, 0.538), a_1=(0.122, 0.371), a_2=(-0.513, 0.583), a_3=(0.770, -0.387), σ=0.4.
The target region A is {(x_ϕ, x_λ); 138.5 ≤ x_ϕ≤ 139.5, 34.5 ≤ x_λ≤ 35.5 }.
Since there is no range constraint on the distribution of the final speed v_f at the target, 
we adopt a uniform distribution with a suitably wide range U_f; here, U_f is defined as the region 
{
(v_ϕ, v_λ); 
V_ϕ(x_A)    -3 ≤ v_ϕ≤ V_ϕ(x_A)+3,  
V_λ(x_A) -3 ≤ v_λ≤ V_λ(x_A)+3 
},
where x_A is the center of target region A.

We also assume that the initial condition is uniformly distributed in
D
=
{
(x, v); 
111 ≤ x_ϕ≤ 129, 
-4 ≤ x_λ≤ 14, 
v ∈ U_0
},
where U_0 is defined as the region 
{
(v_ϕ, v_λ); 
V_ϕ(x_0)    -1.5 ≤ v_ϕ≤ V_ϕ(x_0) + 1.5,  
V_λ(x_0) -1.5 ≤ v_λ≤ V_λ(x_0) + 1.5 
}
and 
x_0 = (120, 5).
This corresponds to the case wherein typhoons that occurred in the Philippines travel to the Tokyo area exactly with a small probability. 
We set M to 10^8 and N to 16.
Examples of Monte Carlo paths for both simulations are given in Figs. <ref> and <ref>.



This simulation was carried out on a laptop computer with 2.3 Ghz Intel core i5 and 8 GBytes memory. 
The computational time of TRMC in this simulation for generating 10^8 Monte Carlo paths is around 4.0 × 10^3 seconds.

Table <ref> shows the result of computational experiments for the stochastic typhoon model.
It shows that the probabilities of FS and TRMC agree within the error bars.
If we ignore the factor defined by Eq. (<ref>), it does not reproduce the unbiased probability as in the case of the stochastic difference equation.
Furthermore, it shows that TRMC is 4.2 times in terms of ρ_2 and 7.3 times in ρ1) more efficient than FS. 
Fig. <ref> shows the convergence of TRMC when the number of Monte Carlo paths M increases. 
It reveals that our algorithm converges correctly on increasing the number of Monte Carlo paths M.


To simulate events with smaller probabilities, we make the target region A smaller as 
{(x_ϕ, x_λ); 138.75 ≤ x_ϕ≤ 139.25, 34.75 ≤ x_λ≤ 35.25 }.
It shows that the smaller the probability, the more efficient our algorithm becomes as compared with the FS.

In Fig. <ref>, a few Monte Carlo paths are shown to have moved northward.
To prevent this from happening and improve its efficiency, we restrict the velocity distribution of Monte Carlo paths to the tendency to move southward.
We change the range U_f of the final speed v_f to 
{
(v_ϕ, v_λ); 
V_ϕ(x_A)    -3 ≤ v_ϕ≤ V_ϕ(x_A)+3,  
V_λ(x_A) -2 ≤ v_λ≤ V_λ(x_A)+2 
}.
We call this simulation TRMC (restricted) in Fig. <ref>.
Table <ref> shows that the probabilities of TRMC and TRMC (restricted) agree within error bars.
Because the number of unnecessary Monte Carlo paths moving northward decreases, TRMC (restricted) is more efficient than TRMC. 
More severe constraint v_λ≥ 0, however, causes a small bias in the estimated probabilities: 
see TRMC (v_λ≥ 0) in Table <ref>.
Fig. <ref> also shows that our algorithm converges correctly on increasing the number of Monte Carlo paths M.

So far, we consider the case that the number of time steps from the initial position to Tokyo (N=16) is precisely known. 
We can relax this assumption, but we should be careful with a limitation of a discrete time model.
A typhoon can pass nearby Tokyo, for example, between N=15 and N=16, which causes an underestimation of the actual risk when we only consider hit at integral time steps. 
A way to reduce this effect is to develop models with smaller steps, while it is also possible to
introduce some initialization or interpolation method into a backward simulation.
However, we leave this as a future problem, because this study aims to check whether the concept of backward simulation is mathematically valid.



	








 §.§ Lorenz 96 Model


As a higher-dimensional example, we evaluate the efficiency of our algorithm for the Lorenz 96 model <cit.>.
The Lorenz 96 model is an atmospheric model and was introduced by Edward Lorenz in 1996.
It is defined as the set of coupled ordinary differential equations


    dx_k/dt   = f_k(x) + ϵ_k, 
    
    	f_k(x)    = -x_k-2x_k-1 + x_k-1x_k+1 - x_k + F,  
    
    	k    = 1… K,

where x={x_k; k=1… K } is the state of the system and F is a constant. 
We set K=9 and introduce Gaussian noise ϵ_k with mean zero and variance σ^2.
Here, we choose F=8, a value known to cause weak chaotic behavior and often used as a benchmark in data assimilation <cit.>. 

To simulate Eq. (<ref>), 
we have to discretize it.
While many discretization schemes are available, we focus on the simplest and most common scheme, the Euler scheme.
The time-discretized version of Eq. (<ref>) by the Euler scheme is 


    x_k, i+1 = x_k, i + f(x_i)Δ t + ϵ_kΔ t, k=1 … K,

where we set Δ t to 0.001 and σ to 0.1/√(Δ t).  

The target region A is 
{(x_1, …, x_K)| -5.0 ≤ x_i ≤ 7.0; i=1 … K }. 
We also assume that the initial state for x_i is uniformly distributed in D={(x_1, …, x_K)| 1.5 ≤ x_i ≤ 8.5; i=1 … K }. 
We set M to 10^7 and N to 100. 

We conduct this simulation on the environment described in Sect. <ref>.
The computational time of TRMC in this simulation for generating 10^7 Monte Carlo paths is around 3.0 × 10^3 seconds.

Table <ref> shows the result of computational experiments for the Lorenz 96 model. 
It shows that the probabilities of TRMC and FS agree within the error bars.
The case where we ignore the factor defined by Eq. (<ref>) does not reproduce the same unbiased probability as the other computational experiments.
The result shows that TRMC can perform better for estimating the probabilities in the  high-dimensional case.
TRMC is 5.18 times more efficient than FS in terms of ρ_2, and 8.23 times in ρ1 more efficient in Table <ref>. 
Fig. <ref> shows the convergence of TRMC when the number of Monte Carlo paths M increases. 
It reveals that our algorithm converges correctly on increasing the number of Monte Carlo paths M.












§ IMPROVED SCHEME WITH RESAMPLING


Let us consider cases with a larger number of time steps.
The proposed algorithm may not always work efficiently in this situation.
For example, we consider the case where N is equal to 500 in the Lorenz 96 model (Fig. <ref>); these weights are normalized such that their sum is 1, i.e., ∑ _j=1^MW^(j) = 1.
The inset located at the top right of the figure shows the graph with a logarithmic scale on the x-axis.
This style is also used in Fig. <ref>.
The weight distribution corresponding to  N=500  in Fig. <ref> has a heavy-tailed distribution.
This phenomenon is referred to as degeneracy, and it means that the weights become unbalanced, and a few weights dominate all the others.
This consequently causes a decrease in computational efficiency <cit.>. 

We introduce an improved scheme to solve this problem, which is realized by resampling <cit.>.
Hereafter, we denote it as TRMC (RS).
This algorithm is effective when both the number of time steps and the amount of noise are large.

Note that our algorithm is based on time-reversed dynamics 
and uses SMC differently from the previous studies 
 <cit.> 
on rare event sampling.




We assume that the resampling procedure modifies the weight at s time step

    ∏_i=0^s-1 W_i

of each Monte Carlo path to an unweighted one by eliminating Monte Carlo paths having small weights and by multiplying
Monte Carlo paths having large weights.

We denote the jth Monte Carlo path as x^(j)= {x^(j)_0, …, x^(j)_s}.
The procedure of resampling is as follows:

	
  * Define normalized weights 
	
    W̃^(j) = ∏_i=0^s-1W_i^(j)/∑_j=1^M∏_i=0^s-1 W_i^(j).

	
  * Resample M times with replacement from set { x^(j)}_j=1^M of Monte Carlo paths, where the probability of sampling set of x^(j) is proportional to W̃^(j).

After a resampling step, 
Monte Carlo paths {x^(j)}_j=1^M 
and associated weights
{W^(j)}_j=1^M 
are replaced by the set of replicated Monte Carlo paths with an equal importance weight
W^(j)=1/M∑_j=1^M∏_i=0^s-1W_i^(j). 
Degeneracy is estimated using the effective sample size <cit.>:

    M_eff = 1/∑_j=1^M (W̃^(j))^2.

A small value of M_eff corresponds to high degeneracy. 
Hence, a resampling procedure is performed when this value is lower than a certain threshold Θ=α M, where
α is a relative threshold.
That is, a resampling procedure is performed when M_eff/M < α.

We can use the Eq. (<ref>) to evaluate probabilities in this case.
Figure <ref> shows a graphical scheme of resampling.	


Using this resampling, 
we simulate the Lorenz 96 model with σ=0.3/√(Δ t), which is larger than that in Sect. <ref>.
We set the threshold α to 0.05, 0.5, and 0.9.
The simulations with these threshold values of α are denoted by α=5%,50%, and 90% respectively.

Table <ref> shows the result of computational experiments for the Lorenz 96 model.
It shows that the probabilities of FS, TRMC, and TRMC (RS, α=50%) agree within the error bars.



On the other hand, Fig. <ref> shows that TRMC (RS) is more efficient than TRMC in a wide range of threshold values. We also show the weight distributions of TRMC and TRMC (RS) 
in Fig. <ref>.
The variance of the distribution is much smaller for TRMC (RS) than for TRMC.














§ DISCUSSION



The examples provided in the preceding sections show that backward simulations using TRMC provide correct averages and can be more efficient than forward simulations. In these examples, the computational efficiencies of TRMC are 3–16 times higher than those obtained by forward simulation, when the calculated probability of hitting the target is 2 × 10^-3–10^-5. 
Note that TRMC can be used to calculate the probability for an arbitrarily small target region; this would be impossible by using forward simulation.

There are, however, cases in which TRMC is inefficient. First, TRMC is not advantageous if the time-reversed paths rarely encounter a region in which the initial density p(x_0) is high; this can occur when the initial density is not broad. Another case in which TRMC can be inefficient is when the weight in Eq. (<ref>) (or, in the continuous time version, Eq. (<ref>)) is highly time dependent. If paths with smaller weights in the initial stage of backward simulation acquire larger weights in the latter stage, resampling of the path (particle splitting) in SMC may not be effective. In this case, if TRMC with SIS is ineffective, TRMC with SMC also shows poor performance. 

To discuss the possible improvement of the algorithm, it is useful to introduce optimal backward dynamics. Although it is not easy to obtain these dynamics a priori, the formal definition is derived as follows. First, the marginal probability at step n obtained from the forward simulation is defined as

    p(x_n) =∫ dx_0:n-1{∏_i=0^n-1  p( x_i+1 | x_i) } p(x_0) _,

which satisfies the relation

    p(x_i+1) =∫  p(x_i+1|x_i)p(x_i) dx_i _.

Using Eqs. (<ref>) and (<ref>), the transition probability q^* of the optimal backward dynamics is defined as

    q^*(x_i+1→ x_i)= 
    	p(x_i+1|x_i)p(x_i)/∫  p(x_i+1|x_i)p(x_i) dx_i
    	=  p(x_i+1|x_i)p(x_i)/p(x_i+1)_.
    	.

Note that Eq. (<ref>) appears similar to the formulas used in Bayesian inference when the probability p(x_i) obtained by forward simulations is regarded as an analog of the prior distribution of x_i. In terms of the selection of Eq. (<ref>) for backward dynamics, the following relation holds:

    {∏_i=0^n-1  p( x_i+1 | x_i) } p(x_0)
    	=
    	p(x_N) ∏_i=N-1^0  q^*(x_i+1→  x_i),

Eq. (<ref>)  means that the combined probability of time-reversed paths defined by forward simulation is recovered by the backward dynamics Eq. (<ref>). Specifically, the time-reversed paths initialized by p(x_N) automatically converge to their initial density p(x_0) using the backward dynamics Eq. (<ref>).
In this sense, q^*(x_i+1→  x_i) in Eq. (<ref>) is considered as the optimal backward dynamics.
Implementation of these dynamics, however, requires the probabilities p(x_i), i=1, … N, which are usually not available prior to the simulations. Note that the backward dynamics defined by the Langevin equation in previous studies <cit.> can be derived from Eq. (<ref>) as a continuous-time limit.

Equations (<ref>)  and  (<ref>) were previously discussed <cit.> in the field of time-series data analysis, where approximations of the marginal probabilities  p(x_i), i=1, … N  are used to define the backward dynamics  q(x_i+1→  x_i). In these studies, the observed data were available at many of the time steps i=0, …, N, whereas the target is given only at i=N in our problem. Then, approximations of probabilities p(x_i), i=1, … N are derived using forward simulations constrained with the observed data (“filtering stage”).

It is, however, difficult to apply these methods to our problem. If we were to apply a similar method to our problem, we would have to run a number of forward simulations to estimate p(x_i), i=1, … N before executing backward simulations. This would be computationally expensive and seems unrealistic without a highly efficient method for the probability estimation. Methods such as those discussed previously <cit.> may be applied to optimize the backward dynamics in our problem, but this is left for future study.

On the other hand, when some observed data are available outside the equations that describe the stochastic process, we may use these data to approximate p(x_i), i=1, … N and hence use them to approximate the optimal backward dynamics. 
In this case, we avoid the use of a large amount of forward computation to construct the optimized backward dynamics. 
This seems possible for the stochastic typhoon model, where data from actual observations of real typhoons are available. 
Note that this idea is different from data assimilation (i.e., inference with simulations combined with observed data), because here we use observed data only to improve the computational efficiency; they do not cause the bias of calculated probabilities.   

In this paper, we assumed that the number of time steps is fixed.
This is not, however, always clear in advance in realistic problems.
In the case of a realistic typhoon model, we must consider the case of passing through Tokyo between two discrete time steps.
The interpolation method in this case will be developed in a future work.




§ CONCLUDING REMARKS



We discussed methods for the backward simulation of the stochastic process. These methods trace a time-reversed path from the target region to the initial configuration. A naïve approach to this problem was shown not to function as expected. To resolve the difficulties, the time reverse Monte Carlo method (TRMC) was introduced. The TRMC method is based on SIS and SMC, and is designed to provide the probabilities of events correctly. TRMC with SIS was tested for the stochastic typhoon model and the Lorenz 96 model; it converges more efficiently than forward simulations in some of these examples. For simulations with a larger number of steps, TRMC with SMC was shown to be advantageous. We also discussed the limitation and possible improvement of TRMC and its relation to the Bayes formula.






§ DEVIATION OF EQ. (<REF>)


The aim of this appendix is to prove Eq. (<ref>).
Up to the first-order Δ t, the Jacobian det(J_g(x)) is given by

    det(J_g(x))
    	   =
    	det(
    	I
    	+
    	∇ f(x) Δ t
    	) 
       =
    	1 
    	+ 
    	(
    	∇ f(x) Δ t
    	)
    	+ 
    	O( (Δ t )^2 ) 
       =
    	exp[
    	div f(x) Δ t
    	]
    	+ O( (Δ t )^2 ),

where I is the unit matrix of order D × D.
D is the dimension of stochastic process X.

Using Eq. (<ref>),
we obtain in the limit as Δ t → 0

    J(y_1, …, y_N)
    	   =
    	∏_i=0^N-1|
    	det(J_g^-1(y_i+1))
    	| 
       =
    	exp[
    	∑_i=1^N
    	- 
    	div f(x_i) Δ t
    	] 
    	+
    	O( (Δ t )^2 )
    	
       exp[
    	-
    	∫_0^Tdivf(x_t) dt
    	].

The above equation is Eq. (<ref>) in the main text.



§ DEVIATION OF EQ. (<REF>)


The aim of this appendix is to prove (<ref>).

Up to the first-order Δ t, the weight at time t_i is given by


	
    W_i   =
    		exp[
    		-( f(x_i+1) - f(x_i) )^TΣ^-1(
    		(x_i+1 - x_i)
    		-
    		Δ t/2(
    		f(x_i+1) + f(x_i)
    		)
    		)
    		] 
    		
       =
    		exp[
    		(  
    		-( f(x_i+1) - f(x_i) )^TΣ^-1(
    		(x_i+1 - x_i)
    		-
    		Δ t/2(
    		f(x_i+1) + f(x_i)
    		)
    		)
    		)
    		] 
    		
       =
    		exp[
    		-
    		(
    		( 
    		∇ f(x_i) 
    		(x_i+1 - x_i)
    		)^TΣ^-1
    		(x_i+1 - x_i)
    		)
    		+
    		o(Δ t)
    		] 
    		
       =
    		exp[
    		-
    		(
    		∇ f(x_i) ^TΣ^-1
    		(x_i+1 - x_i)
    		(x_i+1 - x_i)^T)
    		+
    		o(Δ t)
    		].




In the limit as Δ t → 0, 
Eq. (<ref>) becomes the following stochastic differential equation:

    dX_t = f(X_t) dt + σ dW_t,

where W_t is a standard Brownian motion.
Here, we used Ito's rule <cit.>, in which we substitute √(dt) for dW_t and consider up to the order of dt.  
Using Eq. (<ref>), we obtain the following relation in the limit as Δ t → 0

	
    (x_i+1 - x_i) (x_i+1 - x_i)^T    
    		dx_t dx_t^T
    		= 
    		(f(x_t) dt + σ dW_t)
    		(f(x_t) dt + σ dW_t)^T
        =
    		σ dW_t dW_t^T σ^Tdt
    		+
    		O(dt)			
    		=
    		Σ dt,


where we used the relationships dW_t dW_t^T = dt and σσ^T = Σ.

As a result, we obtain using Eq. (<ref>) and (<ref>)

	
    ∏_i=0^N-1W_i
    		   =
    		exp[
    		-
    		∑_i=0^N-1(
    		∇ f(x_i) ^TΣ^-1
    		(x_i+1 - x_i)
    		(x_i+1 - x_i)^T)
    		+
    		o(Δ t)
    		] 
       exp[
    		-
    		∫_0^T(
    		∇ f(x_t) ^T)
    		dt
    		] 
    		=
    		exp[
    		-
    		∫_0^Tdivf(x_t) dt
    		],


which is Eq. (<ref>) in the main text.



jpsj 

