{'fraction_non_alphanumeric': 0.040045647888785146, 'fraction_numerical': 0.026143790849673203, 'mean_word_length': 4.314222712238148, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hierarchical Bayesian Spatio Temporal Model Comparison on the Earth Trapped Particle Forecast

We compared two hierarchical Bayesian spatio temporal (HBST) results, Gaussian process (GP) and autoregressive (AR) models, on the Earth trapped particle forecast. Two models were employed on the South Atlantic Anomaly (SAA) region. Electron of >30 keV (mep0e1) from National Oceanic and Atmospheric Administration (NOAA) 15-18 satellites data was chosen as the particle modeled. We used two weeks data to perform the model fitting on a 5°x5° grid of longitude and latitude, and 31 August 2007 was set as the date of forecast. Three statistical validations were performed on the data, i.e. the root mean square error (RMSE), mean absolute percentage error (MAPE) and bias (BIAS). The statistical analysis showed that GP model performed better than AR with the average of RMSE = 0.38 and 0.63, MAPE = 11.98 and 17.30, and BIAS = 0.32 and 0.24, for GP and AR, respectively. Visual validation on both models with the NOAA map's also confirmed the superior of the GP than the AR. The variance of log flux minimum = 0.09 and 1.09, log flux maximum = 1.15 and 1.35, and in successively represents GP and AR.


Introduction
Trapped particle has a major role for low Earth orbit (LEO) satellite anomaly. The risk is increased for the LEO Equatorial inclination satellites due to their frequent pass over the South Atlantic Anomaly (SAA) region. The SAA is an area with dense distribution of trapped particle, especially proton and low energy electron [1]. The high energy proton causes the degradation of satellite solar array whereas the low energy electron causes the satellite surface charging [2]. Overcoming with this problem, we develop a LEO trapped particle distribution forecasting system based on a hierarchical Bayesian spatio temporal (HBST) model. The HBST is a statistical model that works for spatio temporal data in a Bayesian framework [3].
In our previous work [4], we have employed a simple HBST model, i.e. Gaussian Process (GP) [5]. In this work, we attempt to employ HBST-Autoregressive (AR) model [6], a modified version of the GP. The difference between the two models is the AR considering the previous day parameter in their true value calculation (Equation 11). Detail of both models is explained in Section 2.
As in our previous work [4], we still use the data taken by National Oceanic and Atmospheric Administration (NOAA) satellite. However, in this work we attempt to employ multiple satellites data instead of only using NOAA 15, i.e. NOAA 15-18. We also improved our gridding method by shrinking the size of our grid from 10x10 of longitude and latitude in our previous work to become 5x5 of longitude and latitude as in Figure 1a. The application of 5x5 grid size makes the SAA is divided in 12 layers at the latitude direction and randomly select four points of validation in each layer (as in Figure 1b). The detail about the gridding process is explained in our previous paper [4].
Since we start to employ multiple satellites data, then we reduce the data input, from 30 days data in our previous work to 14 days data in this work. This step is taken to optimize the time consumed while running computational process. We choose to employ the implementation of both models on the 30 keV of electron energy level (mep0e1). We choose mep0e1 due to the sensitivity of the medium energy proton and electron detector (MEPED) sensor installed in the NOAA satellites. The area of South Atlantic Anomaly (SAA) is still used in this work. We set 31 August 2007 as the target date of the forecasting. Kriging interpolation method still used to draw the distribution of forecasted particle flux value over the SAA region. Detail about the Kriging interpolation of our work can be found in Suparta et al. [7]. In the following section, we describe two HBST models mentioned, i.e. GP and AR models. Section 3 explains the forecasting results of the two models, performed in both statistical and visual analyses and discuss the results. The summary and future work is compiled in Section 4.

Hierarchical Bayesian Spatio Temporal Model
First, we define the notations used in this paper [8]. Let t defines time where t=1,…, T, and T is the total number of time units. In this work, we put T = 14. This is done by selecting data from 17-30 August 2007 as the model input. For the forecasting day (T+1), we put 31 August 2007. The Z (si, t) is a logarithmic value of SAA's 30 keV electron flux. We also denote E (si, t) as a true value corresponding to Z (si, t) at site si, i = 1,…, n at time t. Let Zt = (Z (si, t),…, Z (sn, t))' and Et = (E (si, t),…, E (si, t))'. We should denote all observed data by z, and z * will yield all missing data. As for E, E defines all Et, for t=1,…, T. The total number of observations modeled is defined by N, which N=nT. Generally, in Bayesian forecast modeling, there are p covariates, denote by Xt which are in the form of 1 x p vector. These covariates could be some variable that influence Z value. The regression coefficients are denoted by β, also in the form 1 x p vector. The NOAA data does not provide the covariates values for this work. Therefore, we only employ the intercept value rather than β. Finally, the θ symbol is used to denote all parameters used in our model ( = ( , , 2 , 2 , ∅)).

Forecasting using Gaussian Process (GP) Models.
The spatio-temporal Bayesian GP model [5] is expressed as is a variance-covariance matrix with have dimension n x n, i, j = 1,…, n. 2 is the site invariant common variance, while (. ; , ) is the spatial correlation matrix that have spatial decay and smoothness parameter . Both parameters, and , are independent each other.

Results and Discussion
We use R spTimer package [9] to compute MCMC Gibbs sampling procedure as defined in previous section. We set 13,000 for j value and discharge the first 3,000 iterarations to avoid the starting effect. The forecasting acuracy of both models was then be evaluated by three statistical validations, i.e. the root mean square error (RMSE), mean absolute percentage error (MAPE) and bias (BIAS) on the validation points (Figure 1b).
The results of statistical validation of GP and AR models are displayed in Figure 2. Figures 2a and  2c, mentioned for RMSE and BIAS with units in log flux (particles/cm 2 sec str), whereas for MAPE in Figure 2b, the unit is in percentage. The accuracy of forecasting result is pointed by MAPE, which defines how close the forecasting result to the real data. Our models, GP and AR, give 11.98 and 17.30 for MAPE values which indicates their high degree of forecasting accuracy. On the other hand, RMSE result indicates variance between forecasting value and real data. We find that the average values of RMSE for GP and AR is 0.38 and 0.63, respectively. We could also evaluate the nature of flux distribution variability from the RMSE. The higher of RMSE value in a layer means the higher of variability among the flux in that layer. In this work, we could find that the L3 and L4 layers have the highest values of RMSE. It means that the high variability of flux particle distribution occurred at the SAA core. The trend of flux particle distribution variability shows that the distribution of flux is quite uniform outside the SAA area (L9 to L10). We also find that at SAA boundaries area such as L1, L2, We perform a visual validation for the results as depicted in Figure 3. Figure 3a defines the plot of raw data acquired by NOAA 15-18 on 31 August 2007, whereas Figure 3b defines the NOAA map of mep0e1 on 31 August 2007. Figures 3c and 3d define the estimation of our forecast values after being interpolated by a Kriging technique, for GP and AR, respectively. The Kriging interpolation was computed by using R fields package [10]. The quality of Kriging interpolation is defined by its variance values and represented by Figures 3e and 3f, which successively represent GP and AR. The different between our results in Figures 3c-3d to the NOAA map in Figure 3b is the period of time. NOAA's map in Figure 3b is generated for 4 days distribution (28-31 August 2013) on a nowcast basis. On the other hand, our results in Figure 3c-3d showed forecast maps and designed for one day distribution, 31 August 2007.
As illustrated in Figures 3c-3d, the AR produced a significant overestimate number of a maximum and underestimate number of a minimum. The log maximum number of real data is 7.65 and 1.09 for a minimum. Thus, the GP gives 1.09 for a maximum variance and 0.09 for a minimum, whereas the AR gives 1.35 and 1.15 for a maximum and a minimum, respectively. From the comparison between Figures 3e and 3f, we could see that the GP result Kriging interpolation is better than the AR since it produced less variance than AR. From this result, we can conclude that the autoregressive parameter (in Equation 11) is not suitable for our work.

Summary and Future work
This work has been successfully compared two HBST models to forecast trapped particle flux distribution on 31 August 2007 from 14-day input of NOAA's data. The GP obtained 88.02 % for forecasting accuracy weather the AR obtained 82.7 %. We assume that the less accurate of AR due to the AR temporal correlation parameter used in this work, that needs to be modified regarding to our data. Compared to our previous work, the employment of multiple satellites data in this work omits the missing values occurrence at validation points. In the future work, we propose to shrink the grid into 2x2size of longitude and latitude to reduce the average bias effect in the gridding process.