Approximate Bayesian Computation in the estimation of the parameters of the Forbush decrease model

Realistic modeling of the complicated phenomena as Forbush decrease of the galactic cosmic ray intensity is a quite challenging task. One aspect is a numerical solution of the Fokker-Planck equation in five-dimensional space (three spatial variables, the time and particles energy). The second difficulty arises from a lack of detailed knowledge about the spatial and time profiles of the parameters responsible for the creation of the Forbush decrease. Among these parameters, the central role plays a diffusion coefficient. Assessment of the correctness of the proposed model can be done only by comparison of the model output with the experimental observations of the galactic cosmic ray intensity. We apply the Approximate Bayesian Computation (ABC) methodology to match the Forbush decrease model to experimental data. The ABC method is becoming increasing exploited for dynamic complex problems in which the likelihood function is costly to compute. The main idea of all ABC methods is to accept samples as an approximate posterior draw if its associated modeled data are close enough to the observed one. In this paper, we present application of the Sequential Monte Carlo Approximate Bayesian Computation algorithm scanning the space of the diffusion coefficient parameters. The proposed algorithm is adopted to create the model of the Forbush decrease observed by the neutron monitors at the Earth in March 2002. The model of the Forbush decrease is based on the stochastic approach to the solution of the Fokker-Planck equation.


Introduction and background information
The Forbush decrease (Fd) of galactic cosmic ray (GCR) intensity are identified as a short period decrease (5-12 days) in the stream of the GCR particles reaching the Earth [1]. The Fds occur as an outcome the substantial disturbances in the interplanetary space occurring due to the powerful coronal mass ejecta and solar flares on the Sun. It was shown [2,3,4,5] that Fd of GCR intensity observed by super neutron monitors can provide beneficial information about the structure of the interplanetary magnetic field (IMF) turbulence. In connection with this, the theoretical modeling plays an important role to explain the information obtained from Fds observed by neutron monitors. In general modeling of the GCR transport in the heliosphere is a quite challenging task due to complicated interaction of four main processes: convection by the solar wind, diffusion on irregularities of interplanetary magnetic field (IMF), particle drifts in the non-uniform magnetic field and adiabatic cooling (e.g. [6]). The interplay of these processes represents the Parker transport equation (PTE) [7] being the second order parabolic type partial To compose the complete model of the Fd is complex problem relating with the complex structure of the disturbed vicinity of the interplanetary space responsible for the Fd of the GCR intensity. To take into account the simultaneous changes of the solar wind velocity, structure of the heliospheric neutral sheet (HNS), regular and turbulence components of the IMF in the restricted disturbed vicinity of the interplanetary space is, basically impossible.
In papers [2,3,4,5] was shown that during the Fd is observed the increase of the exponent ν of the power spectral density (PSD) of the components of the IMF (P SD ∝ f −ν , where f is the frequency). From the other side according to the quasi-linear theory [8,9] the diffusion coefficient K is proportional to the rigidity R of the GCR particles as,K ∝ R γ , where γ = 2 − ν. Thus, during Fd we observe a decrease of the diffusion coefficient due to the increase if the IMF turbulence. However, direct insertion of the exponent ν calculated based on the experimental IMF data is not possible. One reason is that in the model we should consider the spatial distribution of the exponent ν while the IMF data cover only one point of space where the spacecraft recording the IMF is placed. The second is that we observe the increase of ν mainly for the IMF components perpendicular to the field line i.e B y and B z , while in the model we should include the overall ν. The third problem is that we cannot calculate the value of ν with arbitrary resolution. The reason is the requirement of employing relatively long data series to attain the values of the exponent ν in low-frequency range f ∈< 10 −6 , 10 −5 > responsible for the scattering of the GCR particles to which neutron monitor and ground meson telescopes respond. This usually was done (e.g. [2,3]) based on the comparison of the single data series with length predetermined by the duration of the Fd; i.e. there were calculated the PSDs for three periods, before, during and after the Fd.
In this paper we employ the Approximate Bayesian Computation (ABC) method to compose the model of the Fd best fitted to the Fd recorded by the neutron monitors in the period of 18 March-2 April 2002. (Fig. 1). We assume that Fd occurs due to decrease of the diffusion coefficient K in the vicinity of space where the Fd is observed. The parameter γ responsible for the K decrease will be fitted by the ABC method.

Model of the Forbush decrease
We model the transport of GCR in the heliosphere by the stochastic approach. In this approach the individual particle motion is described as a Markov stochastic process, and the system evolves probabilistically. Consequently the Parker equation is brought into a form of the backward Fokker-Planck equation of the form: where f = f ( r, R, t) is an omnidirectional distribution function of three spatial coordinates r = r(r, θ, ϕ), particles rigidity R and time t; The coefficients A 1 , ..., A 10 are dependent on the solar wind velocity U , the drift velocity components v d , and the 3D anisotropic diffusion tensor K ij of GCR particles [11]. Based on the Ito formulation the corresponding set of the stochastic differential equations (SDE) is solved (e.g., [12]). In the case of the GCR transport the set of SDEs in three dimensional heliocentric spherical space have a form where r is the trajectory of individual pseudoparticle in the phase space and dW i is the Wiener process. All coefficients and the details of the numerical solution of the set Eq. 2 are given in [13,14,15]. In the proposed model we assume that Fd takes place due to established corotating heliolongitudinal disturbances in the interplanetary space. In this region the IMF turbulence is increased resulting in a gradual diminishing of the diffusion at the Earth orbit, causing larger scattering of the GCR particles, and in effect, fewer GCR particles reach the Earth (e.g., [3,5]). We simulate this process by the gradual decrease and then the increase of the diffusion coefficient at the Earth orbit with respect the heliolongitudes. The diffusion coefficient K II of cosmic ray particles has a form: The parameters α, β and ζ values are selected by the ABC method so that the Fd model matches the observational amplitude of Fd in March 2002 (see Fig. 1).

Approximate Bayesian Computation algorithm setup
Let λ be a parameter vector, given the prior distribution π(λ). The goal of Bayesian inference is to approximate the posterior distribution, π(λ|d obs ) ∝ π(d obs |λ)π(λ), where π(d obs |λ) is the likelihood of λ given the observed data d obs . The main idea of ABC methods is to accept λ as an approximate posterior draw if its associate (modeled) data d is close enough to the observed data d obs . Accepted parameters are a sample from π(λ|ρ(d, d obs ) < ) where the ρ(d, d obs ) is the chosen measure of discrepancy, and is a threshold defining 'closeness margin'. If is sufficiently small then the distribution π(λ|ρ(d, d obs ) < ) will be a good approximation for the posterior distribution π(λ|d obs ). It is often difficult to define an adequate distance function ρ(d, d obs ) between the simulated and observed data. So, in many cases it is replaced with a distance defined by summary statistics, ρ(S(d), S(d obs )). However, as here we consider the Fd amplitudes in specific time points, we can compare those data directly without a use of summary statistics.
In ABC methods, Sequential Monte Carlo (SMC) is used to automatically, sequentially 'clean' approximation of posterior distribution to be used to generate proposals for further steps. In ABC SMC method, the set of samples with weights, called particles sampled from the prior distribution π(λ 1 ), are propagated through a sequence of intermediate posterior distributions π(λ t |ρ(d, d obs ) < t ), t = 1, ..., T , until it represents a sample from the target distribution, π(λ T |ρ(d, d obs ) < T ). These methods aim to generate draws from p(λ t |ρ(d, d obs ) < t ), at each of a series of sequential steps t, where t define a series of thresholds. In [16] the authors proposed strategies called ABC SMC with Adaptive Weights (ABC SMC AW). In the first step algorithm 1 initializes the threshold schedule 1 > 2 >, ..., > T . Then N samples are simulated based on the predefined a priori distribution π(λ 1 ) and the corresponding acceptance condition ρ(d, d obs ) < 1 . Next, the initial uniform weights are calculated. Samples, denoted by a tilde are drawn from the previous generation with probabilities w t−1 j . Using perturbation kernel K λ,t (λ t i |λ i ) new 'fresh' samples λ t i are obtained, with the veracity of the condition ρ(d, d obs ) < t . The weights are calculated according to the formula in stage (11); in stage (12) the weights are normalized and the time step is increased t = t + 1. This procedure is repeated until t ≤ T .
In this paper the scanned parameters vector λ is with the following priori distribution on particular parameters: 12. Normalize weights w t i for i = 1, ..., N end for end for a bias towards underprediction or overprediction of data by the model. We have slightly modified the classical FB measure ρ(d, d obs ) as: under additional definition, that |di−di| |di+di| = 0 when d i = 0 andd i = 0. The most commonly used adaptive scheme for threshold choice is based on the quantile of the empirical distribution of the distances between the simulated data and observations (e.g. [18]). This method determines t in the t time-iteration by sorting the measure {ρ(d t−1 θ i , d obs )} 1<i≤N and setting t such that predetermined percent α t of the simulated data fulfill the condition ρ(d, d obs ) < t . The basic problem with this approach is that we are not able to estimate the value of the measure ρ(d, d obs ) based on the values obtained in iteration t − 1. In this paper, we propose an efficient adaptive mechanism for setting the α t . The value of the quantile is judged based on of the previous stage t − 1, but with three major improvements: ..,N measure of samples set is not discrete. To keep the procedure flexible and support the interval [0, 1] the best choice for continuous distributions seems to be the family of Beta distributions.
• The set {ρ t−1 θ i } i=1,...,N is extended by adding the values of ρ(d, d obs ) calculated in step t regardless whether the samples have been accepted or not.
• We propose to define the value of α t to be dependent on the ratio between a number of currently accepted samples (AS) and the number of all hitherto generated samples (AGS) i.e.: α t = 1 − AS/AGS.
We chose transition kernel K λ,t (λ t i |λ i ) to be a three dimensional Gaussian kernel (e.g. [17]). The set of samples from t − 1 is used to construct a kernel K λ,t (λ t i |λ i ), which will describe the probability mass. Additionally only samples keeping the parameter γ in a physically acceptable interval i.e., 0 < γ < 2 are accepted.

The simulations results and discussion
We have employed the ABC SMC with Adaptive Weights method described in section 3 to find the diffusion coefficient parameters λ ≡ (α, β, ζ) matching the Fd observed in March 2002 (Fig. 3). The code for the numerical realization of the described algorithm is realized in Julia environment [19]. The model of GCR stochastic transport in the heliosphere described in section 2 is easy to parallelize versus the number of simulated pseudoparticles using Distributed Arrays Package. This approach fallouts from the assumption that any random process is independent of the other realization, accordingly each pseudoparticles trajectory is independent on another. We performed the simulations applying the ABC SMC AW method setup described in section 3. Each calculation of the FB (Eq. 5) value requires one run of the model of stochastic transport of GCR in the heliosphere represented by the SDEs (Eq. 2). In a single run, the trajectories of N = 3000 pseudoparticles were traced backward in time in the spherical heliocentric coordinate system. The pseudoparticles with the rigidity R=10 GV were initialized in the point representing the Earth's orbit and traced backward in time until crossing the heliosphere boundary assumed at 100 AU. Details of the numerical solution of Eqs.2 are given in [13]. Fig. 1 illustrates the trace plots for all searched parameters in subsequent stages of the algorithm. The red line marks the mean value within 300 iterations. One can see, that with iterations the samples are more focused on the average value. This is the best visible in the case of β and ζ parameters. These parameters influence the most the amplitude of Fd. From the other side, we can see that the estimation of the parameter λ is not so unequivocal. The variance of this parameter is the largest. The bottom panel in Fig.1 presents the values of FB for accepted samples throughout the scanning procedure. One can see, that the F B value changes significantly in subsequent stages. In the initial iterations of given time, we observe the meaningful change of F B value. This is connected with adapting the threshold to Probability density colors the plot, the more red regions, the higher is its probability. The red vertical lines in diagonal plots mark the value of the parameter with the higher probability. The green horizontal line cuts the value 50% of maximal probability. the knowledge about the Fd amplitude, i.e., the value t is modified according to the Acceptance Ratio (AR) . Starting from 1500 iteration, no significant improvement in FB value is observed. In result, ABC SMC AW algorithm accepts only samples with F B value lower than one estimated in previous stages. This reveals the strength of the proposed procedure, which keeps the Acceptance Ratio at the appropriate level in all stages. The summarized results of reconstruction procedure are presented in Fig.2  The discrepancy results from the assumption that the Fd is modeled only by the change in diffusion coefficient, while the solar wind velocity was assumed as constant. However, one can see that proposed ABC SMC AW method in conjunction with the stochastic model of GCR transport was able to confirm that as the main cause of the Fd observed in March 2002 we can consider the decrease of the exponent γ causing the change of the diffusion coefficient. These results coincide with the experimental data for Fd in March 2002 presented in [2]. The obtained results confirm that the ABC methodology can be successfully applied to judge about the spatial and temporal changes of the parameters responsible for GCR modulation in the heliosphere. The future work will include testing other FB measures and summary statistics. The energy dependence of the diffusion coefficient and the spatial dependence of the drift velocity will also be fitted. We thank the principal investigators of Apatity and Moscow neutron monitor for the ability to use their data.