# A Semiparametric Generalized Exponential Regression Model with a Principled
Distance-based Prior for Analyzing Trends in Rainfall

Arijit Dey and Arnab Hazra CONTACT Arnab Hazra. Email<EMAIL_ADDRESS>Department of Mathematics and Statistics, Indian Institute of Technology
Kanpur, Kanpur, 208016, India.

###### Abstract

The Western Ghats mountain range holds critical importance in regulating
monsoon rainfall across Southern India, with a profound impact on regional
agriculture. Here, we analyze daily wet-day rainfall data for the monsoon
months between 1901–2022 for the Northern, Middle, and Southern Western Ghats
regions. Motivated by an exploratory data analysis, we introduce a
semiparametric Bayesian generalized exponential (GE) regression model; despite
the underlying GE distribution assumption being well-known in the literature,
including in the context of rainfall analysis, no research explored it in a
regression setting, as of our knowledge. Our proposed approach involves
modeling the GE rate parameter within a generalized additive model framework.
An important feature is the integration of a principled distance-based prior
for the GE shape parameter; this allows the model to shrink to an exponential
regression model that retains the advantages of the exponential family. We
draw inferences using the Markov chain Monte Carlo algorithm. Extensive
simulations demonstrate that the proposed model outperforms simpler
alternatives. Applying the model to analyze the rainfall data over 122 years
provides insights into model parameters, temporal patterns, and the impact of
climate change. We observe a significant decreasing trend in wet-day rainfall
for the Southern Western Ghats region.

###### keywords:

Climate change; Generalized exponential distribution; Markov chain Monte
Carlo; Penalized complexity prior; Semiparametric Bayesian regression; Western
Ghats; Wet-day precipitation modeling

††articletype: Research Article

##  1 Introduction

The Western Ghats, a prominent mountain range along the western coast of
India, plays a crucial role in shaping the climatic patterns and hydrological
dynamics of Southern India. Known for its exceptional biodiversity, lush
forests, and vital water resources, the Western Ghats has long captured the
attention of researchers and environmentalists [32, 50]. Among the various
climatic parameters that influence this ecologically significant region,
rainfall is a crucial driver of its diverse ecosystems, water availability,
and overall environmental health. The Western Ghats, characterized by its
rugged terrain and proximity to the Arabian Sea, experiences a unique and
intricate rainfall pattern heavily influenced by monsoon dynamics [48]. Over
the last century, this region has experienced notable climatic shifts due to
global climate change and local human activities [1, 49]. Analyzing wet-day
rainfall in this region during monsoon months over an extended period of 122
years using a flexible statistical model offers a valuable opportunity to gain
insights into long-term trends, variability, and potential shifts in the
monsoonal regime.

Researchers have widely employed the exponential distribution to model
rainfall data [46, 16] in literature; its simplicity integrates seamlessly
into hydrological and climatological frameworks. However, contemporary
research increasingly recognizes the need for innovative probability
distribution models for better encompassing complex real-world data patterns.
This realization has prompted the introduction of novel probability classes
with far-reaching implications across diverse research domains. [2] provides
an excellent overview of newly developed distributions. A notable collection
of models are generalized distributions, gaining attention from both practical
and theoretical statisticians for their adaptability to various datasets. Some
examples include the Marshall-Olkin generalized exponential distribution [31],
generalized inverse Gaussian distribution [21], generalized Rayleigh
distribution [26], etc. A more comprehensive examination of these
distributions is available in [45].

[10] introduced another crucial generalized distribution called the
generalized exponential (GE) distribution, which emerges as a specific case
within the three-parameter exponentiated-Weibull model. The GE distribution
has two parameters- a shape and a rate (or scale, defined as the inverse of
rate) parameter. This distribution boils down to an exponential distribution
when the shape parameter is one. Thus, with an additional shape parameter, it
expands the capabilities of the exponential distribution, making it more
adaptable to various datasets. Since its introduction, many researchers have
integrated substantial advancements in exploring different properties,
estimation strategies, extensions, and applications of this distribution. For
instance, [11] found the efficacy of GE distribution compared to gamma or
Weibull distributions, whereas [12] discussed different methods of estimating
the parameters of GE distribution. [20], [38], and [25] explored Bayesian
estimation and prediction methods in this context. [13] reviewed the existing
results and discussed some new estimation methods and their characteristics.
Numerous researchers have modeled the experimental data using GE distribution
across several disciplines like meteorological studies [29, 15]; flood
frequency analysis [30]; reliability analysis [3]; lifetime analysis [4]; risk
analysis [42]. But, as of our knowledge, any GE regression-type model has
never been proposed in the literature. In this study, taking a step beyond
exponential regression, we employ the GE regression framework to explore
rainfall patterns.

Rainfall data collected over a century are inherently nonstationary. Here,
modeling the temporal trend using traditional parametric regression would
struggle to capture the intricate and evolving short-term temporal patterns.
In this context, a semiparametric regression setup emerges as a promising
approach. In the existing literature, many researchers have delved into
applying semiparametric regression techniques for analyzing rainfall patterns
[53, 34]. While a generalized linear model (GLM) assumes the link function to
be a linear combination of the covariates, the more flexible generalized
additive models [GAM, 14] allow the link function to be a sum of nonlinear
smooth functional forms of the underlying covariates. We generally model each
smooth function in GAMs as a linear combination of basis functions like cubic
B-splines. Instead of estimating the entirely unknown function, we draw
inferences based on basis function coefficients [8]. Henceforth, instead of
GAM, we use the term ‘semiparametric regression’, which is common in Bayesian
nonparametrics. The rate parameter of the GE distribution is always positive,
and hence, it would be reasonable to model the log-rate in a semiparametric
regression framework.

Within the Bayesian methodology, priors hold a pivotal role in inference, and
the literature provides a diverse spectrum of prior distributions utilized for
regression coefficients in semiparametric regression frameworks. For instance,
Gaussian prior was employed by [7], while [27] opted for Laplace prior. [28]
utilized Zellner’s $g$-prior, while [6] considered flat priors, and [23] used
the Normal-Gamma prior. On the other hand, the gamma distribution has
consistently been considered the most natural prior choice for the shape
parameter of the GE distribution; authors who introduced the GE distribution
chose a gamma prior for the shape parameter in [25] as well. Besides [38] and
[22] also employed gamma prior for the shape parameter. However, the
literature demonstrates that a handful of alternative prior choices have also
been utilized. For example, [33] employed a Jeffrey’s prior, indicating their
preference for an objective prior, and [5] opted for a non-informative prior
in their study.

The Penalized Complexity (PC) prior, introduced by [43], has emerged in recent
literature, which mitigates the model complexity through penalization. In
cases where a model extends from a simpler foundational model by incorporating
an additional parameter, this type of prior becomes applicable; it penalizes
the escalation in model complexity that arises when favoring the extended
model over its more straightforward counterpart. Existing literature
encompasses instances of this approach across various models [47]. [51]
developed PC priors for estimating the effective degrees of freedom in
Bayesian penalized splines (P-splines), while [35] discussed a PC prior for
the skewness parameter of the power links family, and [44] proposed
interpretable and comprehensive PC priors for the coefficients of a stationary
autoregressive process.

In this paper, along with modeling the wet-day rainfall of the Western Ghats
region for the last century by a semiparametric Bayesian GE regression model,
we employ the PC prior for the GE shape parameter, which allows the GE
regression to shrink towards an exponential regression. Thus, the exponential
distribution is considered the base model for the GE distribution. In several
practical examples [16, 17], the exponential distribution is already a
reasonable model and enjoys several benefits of being a member of the
exponential family, thus shrinking the GE distribution to its base model
through shrinking the shape parameter to one is justified. On the other hand,
we opt for the independent Gaussian priors for the regression coefficients. We
draw inferences using the Markov chain Monte Carlo (MCMC) algorithm; here,
conjugate priors are not available for the model parameters, and thus we
update them using Metropolis-Hastings steps. We conduct a thorough simulation
study by simulating 1000 datasets from each combination of the model
generating and model fitting scenarios, and we compare the performances of
parametric and semiparametric Bayesian GE regression models under the
conventional gamma prior choices for the GE shape parameter along with our
proposed one. We study the coverage probabilities for the shape parameter and
the rate functions and compare these two models using widely applicable
information criterion [WAIC, 52]. We implement the proposed methodology to the
daily wet-day precipitation spanning from 1901 to 2022 in different regions of
the Western Ghats mountain range, using the year as a covariate and wet-day
precipitation as a response. We study the convergence and mixing of the MCMC
chains and compare different model fits in terms of WAIC.

The paper is structured as follows: Section [2] delves into the GE
distribution, thoroughly examining its properties. Section [3] discusses an
exploratory data analysis that justifies our semiparametric GE model
assumption for the wet-day precipitation data. In Section [4], we introduce
the GE regression model. Proceeding to Section [5], we concentrate on
delineating the prior specifications for the regression model, including
introducing a principled distance-based prior for the shape parameter of the
GE distribution. Bayesian parameter inference is addressed in Section [6].
Section [7] presents the outcomes of the simulation study, while Section [8]
discusses the results obtained from our proposed model and some simpler
alternatives. Finally, Section [9] summarizes our findings and contributions.

##  2 Background: Generalized Exponential (GE) Distribution

We say a random variable $X$ follows GE distribution if its cumulative
distribution function (CDF) is given by

$F(x;\alpha,\lambda)=\left(1-e^{-\lambda x}\right)^{\alpha};\hskip
10.0ptx,\alpha,\lambda>0,$

where $\alpha$ is the shape parameter and $\lambda$ is the rate parameter. The
corresponding probability density function (PDF) is given by

$f(x;\alpha,\lambda)=\alpha\lambda\left(1-e^{-\lambda
x}\right)^{\alpha-1}e^{-\lambda x};\hskip 10.0ptx,\alpha,\lambda>0.$ (1)

The GE distribution is a more complex model than the exponential distribution,
as it incorporates an extra shape parameter. Both models coincide when
$\alpha$ = 1.

###  2.1 Properties of GE

The hazard function of the GE distribution is given by

$h(x;\alpha,\lambda)=\frac{f(x;\alpha,\lambda)}{1-F(x;\alpha,\lambda)}=\frac{\alpha\lambda\left(1-e^{-\lambda
x}\right)^{\alpha-1}e^{-\lambda x}}{1-\left(1-e^{-\lambda
x}\right)^{\alpha}};\hskip 10.0ptx>0.$

The GE distribution has an increasing or decreasing hazard rate depending on
the value of the shape parameter. The hazard function is decreasing for
$\alpha<1$, constant for $\alpha=1$, and increasing for $\alpha>1$. The moment
generating function (MGF) of the GE distribution is given by

$M_{X}(t)=\frac{\Gamma(\alpha+1)\Gamma(1-\frac{t}{\lambda})}{\Gamma(1+\alpha-\frac{t}{\lambda})};\hskip
5.0pt0\leq t<\lambda,$

and differentiating the log of the MGF with respect to $t$ repeatedly and then
setting $t=0$, we get the expectation, variance, and skewness of GE
distribution as

$\displaystyle\textrm{E}(X)$ $\displaystyle=$
$\displaystyle\lambda^{-1}\left[\psi(\alpha+1)-\psi(1)\right],$
$\displaystyle\textrm{V}(X)$ $\displaystyle=$
$\displaystyle\lambda^{-2}\left[\psi^{(1)}(1)-\psi^{(1)}(\alpha+1)\right],$
$\displaystyle\text{Skewness}(X)$ $\displaystyle=$
$\displaystyle\left[\psi^{(2)}(\alpha+1)-\psi^{(2)}(1)\right]\bigg{/}\left[\psi^{(1)}(1)-\psi^{(1)}(\alpha+1)\right]^{\frac{3}{2}},$

where $\psi^{(m)}(z)=\dfrac{\partial^{m}}{\partial
z^{m}}\psi(z)=\dfrac{\partial^{m+1}}{\partial z^{m+1}}\ln\Gamma(z)$ is the
polygamma function of order $m$; for $m=0$, it denotes the digamma function.

Figure 1: Generalized Exponential probability density function (top-left),
hazard function (top-right), mean (bottom-left), variance (bottom-middle), and
skewness (bottom-right) functions. Both top panels share the same legend.

Figure 1 sheds light on different aspects of the GE distribution, e.g., PDF,
hazard function, mean, variance, and skewness. The top-left panel of Figure 1
shows that for $\alpha<1$, the curve depicting the PDF of the GE distribution
has an asymptote at the Y-axis and then decreases exponentially and
monotonically as we move across the positive real line. With $\alpha=1$, GE
coincides with the exponential distribution, thus having mode at zero (with
value $=\lambda$) and gradually decreasing similarly as the previous case.
When $\alpha>1$, the curve initiates at zero, then increases over a range of
values, and eventually decreases monotonically, having a unique mode at
$\log(\alpha)/\lambda$. As mentioned earlier, the top-right panel of Figure 1
shows that the hazard function is monotonically increasing when $\alpha<0$,
monotonically decreasing when $\alpha>0$, and constant (the value being
$\lambda=1$) when $\alpha=1$. The mean and variance of GE behave somewhat in a
similar manner. From the bottom-left and the bottom-middle panel of Figure 1,
we see for a fixed value of $\alpha$, both mean and variance decrease with
increasing $\lambda$ and for a fixed value of $\lambda$, both increase as
$\alpha$ increases. On the other hand, the skewness of the GE distribution
depends only on the shape parameter and decreases exponentially with
increasing $\alpha$ (bottom-right panel of Figure 1).

##  3 Data Preprocessing and Exploratory Data Analysis

In this section, we describe the preprocessing steps involved in obtaining the
dataset comprising average daily precipitation data for rainy days in the
Northern, Middle, and Southern Western Ghats regions during the monsoon months
between 1901–2022. Besides, we discuss pieces of evidence based on exploratory
data analysis that confirm the suitability of a GE distribution for fitting
the data and to determine whether a semiparametric mean structure is necessary
or a linear trend structure would suffice.

We obtain daily gridded rainfall (in mm) data over the Western Ghats region
with a spatial resolution of $1.0\degree\times 1.0\degree$, covering the
period from 1901–2022. The data was sourced from the official website of the
Indian Meteorological Department, Pune
(https://www.imdpune.gov.in/cmpg/Griddata/Rainfall_1_NetCDF.html). The gridded
data product was obtained through spatial interpolation of ground-station data
following the procedure described in [36]. We extract the daily rainfall
information for June, July, August, and September (JJAS) throughout all the
years. Additionally, we exclude days within the JJAS months where recorded
rainfall amounts were zero. Out of the pixels representing the Western Ghats
area, we group them into three distinct significant regions: the Northern,
Middle, and Southern regions (the regions are shown in the supplementary
material). We compute the daily rainfall values for each region by calculating
the average of the corresponding pixel values within that region. Afterward,
we conduct a further analysis based on these regions.

Figure 2: Bar diagrams of the annual average wet-day rainfall during June
through September along with fitted mean curves based on twelve cubic
B-splines with equidistant knots, for Northern, Middle, and Southern Western
Ghats regions (top panels). Histograms of the detrended residuals from the
daily rainfall overlapped with the fitted GE densities (bottom panels).

Given our dataset (after preprocessing) spans over a century, our initial
focus involves performing necessary analyses to address any potential trends
within the data. In the top panels of Figure 2, we present a bar diagram
depicting the average yearly rainfall for each year. No clear long-term linear
trend is observable for any of the three regions. However, several short-term
upward and downward patterns are noticeable. We use a basis spline regression
approach to explore such short-term trends, which treats daily rainfall values
as response variables and corresponding years as covariates. Considering the
residuals from this regression, we can effectively eliminate any potential
trends embedded in the data. We overlap the estimated means with the bar
diagrams in the top panels of Figure 2. Firstly, the estimated mean curve
aligns well with the visualized bar diagram. Moreover, both components
highlight the presence of a nonstationary rainfall pattern. This pattern, in
turn, underscores the suitability of employing a semiparametric regression
model, which can effectively accommodate and incorporate these nonstationary
patterns within the data.

Subsequently, after conducting a fundamental analysis to identify and remove
outliers via the popular adjusted-boxplot method developed by [19], we present
two important visualizations in the bottom panels of Figure 2, where the
panels correspond to three regions of interest. Firstly, we showcase a
histogram illustrating the distribution of the detrended residuals obtained by
exponentiating the residuals obtained by fitting a semiparametric regression
curve to the log-transformed wet-day rainfall observations, which aligns with
the standard link function formulations for generalized additive models.
Additionally, a red line denotes the fitted density of the GE distribution,
with parameters estimated from the detrended residuals. We observe a strong
alignment between the estimated density and the associated histograms,
indicating a favorable fit. This visual representation significantly supports
the rationale behind the GE regression model proposed in this paper.
Additionally, the plots highlight a marked convergence of the GE distributions
towards their foundational model, the exponential distribution. This
convergence also reinforces our second consideration: using a novel distance-
based prior for the shape parameter of the GE distribution.

##  4 Generalized Exponential (GE) Regression

The GE regression model is a statistical model that we can use for modeling
continuous positive-valued response variables. It can be considered an
extension of the standard linear regression model that allows for non-Gaussian
and asymmetric distributions, accommodating heteroscedasticity and skewness in
the data. The regression framework assumes that the response variable $Y$
follows a GE distribution, and it models the relationship between the $Y$ and
the covariates $\bm{X}=(X_{1},\ldots,X_{P})^{\prime}$ through a linear
predictor $\eta$. The linear predictor is a linear combination of the
covariates with associated regression coefficients, given as

$\eta(\bm{X})=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\dots+\beta_{P}X_{P},$
(2)

where $\bm{\beta}=(\beta_{0},\beta_{1},\ldots,\beta_{p})^{\prime}$ is the
vector of regression coefficients.

In the GE regression model, the shape parameter $\alpha$ is considered an
inherent property of the distribution that characterizes the shape and
asymmetry of the distribution, allowing for a more flexible modeling approach
compared to a standard linear regression with a Gaussian error component. On
the other hand, the rate parameter $\lambda$ is the parameter of interest in
the regression model, which captures the association between the covariates
and the response variable. By incorporating the covariates through the rate
parameter, the model captures how changes in the covariate values influence
the rate or dispersion of the response variable. Moreover, given that the rate
parameter of the GE distribution is positive, we relate it to the linear
predictor from (2) using a link function, which is designed to ensure the rate
parameter always stays positive. Thus, we can conceptualize the GE Regression
model as

$Y_{i}|\bm{X_{i}}=\bm{x_{i}}\sim GE(\alpha,\lambda_{i}),$

where $g(\lambda_{i})=\eta(\bm{x_{i}})$, with $g(\cdot)$ representing the
appropriate link function and
$\bm{x_{i}}=(x_{i1},x_{i2},\dots,x_{iP})^{\prime}$. The parameters in the
Generalized Exponential regression model, including the regression
coefficients and the shape parameter, are typically estimated using maximum
likelihood estimation or other suitable estimation methods.

The above-explained theory introduces a parametric framework for the GE
regression model, aiming to capture the relationship between the response
variables and the covariates by utilizing rate parameters. However, this paper
adopts a semiparametric approach to achieve the same objective, which allows
us to incorporate both parametric and nonparametric components into the model,
providing flexibility and accommodating potential complexities in the
relationship between the covariates and the response variables. We introduce
the specific form of the predictor function and other essential formulations
following a general introduction to semiparametric regression models.

###  4.1 Semiparametric Regression

Parametric regression assumes a known distribution for the response variable
and estimates a finite-dimensional set of parameters. It is straightforward,
provides easy estimation, and offers interpretability. However, it may
struggle to capture nuanced aspects of the data, such as non-linearity or
variable interactions. In contrast, nonparametric regression assumes no
specific form for the relationship between the response and explanatory
variables, allowing flexibility based on data-derived information. While this
approach allows for more flexible modeling, it is computationally intensive,
less interpretable, and can be affected by the curse of dimensionality.
Semiparametric regression integrates the above two approaches, allowing us to
have the best of both regimes. It incorporates the interpretability of the
parametric setup and flexibility of the nonparametric setup.

In linear or generalized linear models (GLM) that fall under the parametric
setup, we assume the conditional mean of the distribution of the response
variable is linked with the linear predictor through a linear combination of
the covariates or their functions. Semiparametric setup extends this domain of
regression models by introducing the nonparametric component in the linear
predictors. In this setup, the most general formulation of the linear
predictor can be given as

$\eta(\bm{x})=\sum_{p=1}^{P}f_{p}(x_{p}),$ (3)

where $f_{j}$’s are smoothing functions of continuous covariates and
$\bm{x}=(x_{1},x_{2},\dots,x_{P})^{\prime}$. Regarding smoothing functions,
most semiparametric methods assume they can be expressed as a linear
combination of finite basis functions, often denoted as

$f_{p}(z)=\sum_{k=1}^{K_{p}}\beta_{p,k}B_{p,k}(z);\hskip 8.0ptp=1,2,\dots,P,$
(4)

where $B_{p,k}(\cdot)$’s are known basis functions and $\beta_{p,k}$ are
unknown basis function coefficients that determine the shape of the smoothing
function $f_{p}(z)$. A basis expansion of $M$ terms can match the true curve
$f_{p}(\cdot)$ at any $M$ points $X_{1},\dots,X_{M}$ in the range of
covariates. Hence, increasing $M$ gives us an arbitrarily flexible model.

In this study, we employ a semiparametric model akin to (3) for the rate
parameter of the GE distribution. With a covariate vector comprising $P$
components and the appropriate logarithmic link function, the regression model
takes the form:

$Y_{i}|\bm{X_{i}}=\bm{x_{i}}\sim
GE\big{(}\alpha,\lambda(\bm{x_{i}})\big{)}\text{ with
}\log\left\\{\lambda(\bm{x_{i}})\right\\}=\sum_{p=1}^{P}\sum_{k=1}^{K_{p}}\beta_{p,k}B_{p,k}(x_{ip}),$
(5)

where $B_{p,k}(\cdot)$’s are cubic B-splines and $\beta_{p,k}$’s are the
spline coefficients representing the weights assigned to the corresponding
spline functions. To provide a brief overview, a cubic B-spline is a
piecewise-defined polynomial cubic function that is defined on a set of knots
or control points, taking the form $B_{s}(x)=(x-v_{s})_{+}^{3}$ where
$v_{s}$’s are the fixed knots that span the range of $x$ and ‘+’ denotes the
positive part.

##  5 Prior Specification

Selecting an appropriate prior is one of the most crucial parts of a Bayesian
analysis. While there is no universal rule for constructing optimal prior
distribution, choosing an appropriate prior can significantly enhance the
quality of the study. A well-chosen proper prior can stabilize the posterior
distribution and yield better results compared to an improper prior [Chapter
4, 39]. This section defines the prior distribution parameters of the
regression model presented in (5).

In semiparametric Bayesian regression, instead of formulating an explicit
prior distribution for $\lambda$, independent prior distributions are
explicitly specified for the spline coefficients
$\bm{\beta}=(\beta_{1},\beta_{2},\dots,\beta_{K})^{\prime}$. This paper
considers independent weakly-informative Gaussian priors for $\beta_{k}$’s. As
for the shape parameter of the GE distribution, we employ a newly developed
class of priors. In cases where a model is constructed based on a simpler base
model, the chosen prior should accurately reflect the characteristics of the
model considered and capture its departure from the base model. This type of
prior construction is founded upon the work of [43], who introduced the
Penalized Complexity (PC) prior. The PC prior is a model-based prior that
imposes a penalty on the deviation of the model of consideration from its
simpler base version at a logarithmic constant rate. The following subsections
discuss the PC prior for $\alpha$.

###  5.1 Penalized Complexity (PC) prior

The PC prior is an informative proper prior that exhibits robustness
properties of high quality and invariance under reparameterization. It aims to
penalize the complexity that arises when we move from a simpler base model to
a more complex one, thereby preventing overfitting and adhering to Occam’s
razor principle [43]. By using the PC prior, we uphold the principle of
parsimony, which suggests a preference for simpler models until sufficient
evidence supports more complex alternatives.

The PC prior is established based on the statistical difference between the
proposed complex model and its base model. We quantify this distance using the
Kullback-Leibler divergence (KLD) [24]. KLD is an information-based measure
that essentially measures how much information we lose when we substitute a
complex model having PDF $f$ with its simpler version having PDF $g$. For the
GE distribution, the exponential distribution is commonly selected as the
appropriate base model. Hence, for our purposes, we consider $f$ as the GE PDF
and $g$ as the exponential PDF.

For two continuous distributions with PDF $f$ and $g$ defined over the same
support, KLD is defined as

$\text{KLD}(f\parallel
g)=\int_{-\infty}^{\infty}f(y)\log\left(\frac{f(y)}{g(y)}\right)\,dy\,.$ (6)

We define the distance between the two models by the ‘unidirectional’ distance
function $d(f\parallel g)=\sqrt{2\text{KLD}(f\parallel g)}$ [43]. The absence
of symmetry in KLD is not a concern within this context. Our main focus is on
quantifying the additional complexity that arises from employing the intricate
model rather than the other way around.

The main idea of the PC prior involves assigning priors to the distance
between two models, $d(\alpha)$, rather than directly on the model parameters,
and then by employing a change-of-variables approach, one can obtain a prior
distribution for the parameter of interest. While constructing the PC prior
for the shape parameter $\alpha$, we take this distance as a function of
$\alpha$, i.e.,
$d(\alpha)=\sqrt{2\text{KLD}(\alpha)}=\sqrt{2\text{KLD}(f\parallel g)}$ with
$f$ and $g$ being GE and exponential PDFs respectively.

To incorporate the fact that the prior should have a decaying nature as a
function of the distance between the two models, we take the constant rate
penalization assumption and construct the PC prior by assigning an exponential
prior to the distance, i.e., $d(\alpha)\sim\text{Exp}(\theta)$ with
$\theta>0$. This gives us the PC prior for $\alpha$ as

$\pi(\alpha)=\theta e^{-\theta d(\alpha)}\left|\frac{\partial
d(\alpha)}{\partial\alpha}\right|,$ (7)

where $\theta>0$ is a user-defined quantity that controls the prior mass at
the tail. We impose this quantity in such a way that this characterizes how
informative the PC prior we want. This is achieved by imposing a condition
$\Pr[d(\alpha)>U]=\xi$ where $U$ is the upper bound of the tail-event and
$\xi$ is the weight of the event [43].

###  5.2 PC prior for the shape parameter of GE distribution

This section contains the KLD between our complex model, the GE model $(f)$
and its natural base model, Exponential $(g)$, and the PC prior of the shape
parameter, $\alpha$.

######  Theorem 5.1.

The KLD, with $f$ being the PDF of the GE distribution given in (1) and $g$
being the PDF of the exponential distribution with rate $\lambda$, is given by

$\text{KLD}(\alpha)=\log(\alpha)+1/\alpha-1.$

######  Theorem 5.2.

The PC prior for the shape parameter ($\alpha$) of the Generalized Exponential
(GE) distribution is supported over the positive real line and given as

$\pi(\alpha)=\frac{\theta}{2}\cdot\exp\left(-\theta\,\sqrt{2\log(\alpha)+\frac{2(1-\alpha)}{\alpha}}\right)\cdot\left(2\log(\alpha)+\frac{2(1-\alpha)}{\alpha}\right)^{-\frac{1}{2}}\cdot\left|\frac{1}{\alpha}-\frac{1}{\alpha^{2}}\right|.$

###### Proof.

Proof for Theorem 5.1 is given in the Appendix A and Theorem 5.2 follows
directly from (7) and the expression $d(\alpha)=\sqrt{2\text{KLD}(\alpha)}$.
Moreover, Theorem 5.2 includes a scaling factor which ensures that
$\pi(\cdot)$ integrates to one. ∎

Figure 3: The PC prior for the shape parameter of the GE distribution for
different choices of the hyperparameter $\theta$.

In Figure 3 we illustrate $\pi(\alpha)$ for different hyperparameter
specifications $\theta$. We notice a proportional relationship between the
value of $\theta$ and the extent of contraction to the base model. As the
value of $\theta$ decreases, the tails become heavier, resulting in reduced
contraction towards the base model. We also observe that for $\theta\leq 4/3$,
the mode of the PDF occurs at a value of $\alpha$ less than one, but for
$\theta\geq 4/3$, the mode is at $\alpha=1$ (Figure 3). One might expect the
prior would have mode at $\alpha=1$ irrespective of the value of $\theta$.
Here, we do not necessarily need the mode at $\alpha$ = 1, and rather we
should rely on a prior that is consistent with the principles of PC prior.

##  6 Inference

This paper employs Bayesian estimation methods to infer and quantify the
uncertainty surrounding the parameters of interest. Fundamentally, Bayesian
inference encompasses amalgamating prior knowledge or beliefs (prior
probability) regarding a particular event or hypothesis with evidence or data
(likelihood) via the Bayes theorem and derives a revised or updated belief
(posterior probability), which provides a more comprehensive understanding of
the values of the parameters in question and their inherent uncertainty.

In our context, the likelihood function based on $n$ observations from the GE
distribution under the regression setting from (5) is given by

$L(\alpha,\bm{\beta}|\bm{y})=\prod_{i=1}^{n}f\big{(}y_{i};\alpha,\lambda(\bm{x_{i}})\big{)},$
(8)

where $\bm{y}=(y_{1},y_{2},\dots,y_{n})^{\prime}$ is the observed data,
$f(y;\alpha,\lambda)$ is the PDF of the GE distribution from (1) and
$\lambda(\bm{x_{i}})$ taking form as given in (5). Also, let $\pi(\alpha)$ and
$\pi(\bm{\beta})$ denote the specified mutually independent priors for the
parameters $\alpha$ and $\bm{\beta}$.

Combining the priors for shape parameter $\alpha$ and the regression
coefficients $\bm{\beta}$, and the likelihood function as given in (8), we
obtain the joint posterior distribution as
$\pi(\alpha,\bm{\beta}|\bm{y})\propto
L(\alpha,\bm{\beta}|\bm{y})\cdot\pi(\bm{\beta})\cdot\pi(\alpha)$ from which
Bayesian inference is facilitated. However, the explicit form of the marginal
posterior density of the parameters is not analytically tractable, leading to
employing simulation-based techniques such as Markov chain Monte Carlo (MCMC)
methods or numerical approximation methods like Integrated Nested Laplace
Approximations (INLA), introduced by [41]. In this paper, we employ MCMC
techniques for parameter inference, specifically utilizing the adaptive
Metropolis-Hastings algorithm within Gibbs sampling. We iteratively adjust the
variance of the proposal distribution within the chain so that the acceptance
rate remains between 0.3 and 0.5. We initiate the MCMC chains with an initial
value of 1 for $\alpha$ and the maximum likelihood estimate for $\bm{\beta}$
as calculated under $\alpha$ = 1. Our approach employs an algorithm that
updates one parameter at a time. Regarding computational resources, conducting
the simulations outlined in Section 7, involving 144 different configurations,
each with 1000 datasets, and employing an MCMC output of 10000 iterations,
required nearly 8 hours to finalize. This computational task was carried out
on a desktop system equipped with a Ryzen 9 5900x 12-core 24-threaded
processor and 64GB of RAM. On average, the computation time for a single MCMC
chain was around 1 second for the parametric setup (assuming $\log(\lambda)$
to be a simple linear function of year) using a simple linear regression model
and approximately 3 seconds for our proposed semiparametric model.

Furthermore, it is feasible to describe the ultimate distribution of the
posterior estimates for the parameters $\alpha$ and $\bm{\beta}$. In this
setting, the Bernstein-von Mises theorem outlines the shape of this asymptotic
distribution, delineating how the parameters $\alpha$ and $\bm{\beta}$ behave
as the sample size tends towards infinity. To gauge the level of uncertainty
linked with these parameter estimations, we investigate the asymptotic
variance of the parameters, which is encapsulated by the inverse of the
information matrix. Additional elaboration can be found in the supplementary
materials.

##  7 Simulation Study

We conducted an extensive simulation study to demonstrate the effectiveness of
the PC prior and the proposed semiparametric model over a parametric model,
where the GE rate parameter (in log scale) is modeled as a simple linear
function of the covariate(s). We designed two separate simulation setups for
this purpose. In the first setup, we compared the PC prior to the conventional
gamma prior for the shape parameter of the GE distribution under four
different scenarios. These scenarios included– Setting 1: generating data from
a linear setup and fitting a parametric model, Setting 2: generating data from
a nonlinear setup and fitting a parametric model, Setting 3: generating data
from a linear setup and fitting a semiparametric model, and Setting 4:
generating data from a nonlinear setup and fitting a semiparametric model. For
each setting, we used four different prior specifications: the PC prior with
parameters 2.5 and 5 and the gamma prior with parameters $(0.01,0.01)$ and
$(1,1)$. We compared their effectiveness in estimation. Additionally, we
considered two cases for the number of samples, namely $n$ = 24 and $n$ = 99,
to gain insights into parameter estimation in scenarios with small and large
sample sizes, respectively. In each case of this simulation setup, we
calculated the coverage probability of $\alpha$ based on $95\%$ credible
intervals and the absolute bias in estimating $\alpha$ to facilitate the
comparison.

Figure 4: Coverage probabilities (top illustration) and absolute bias values
(bottom illustration) computed based on imposing a PC prior. For each
illustration, from the left-most to the right-most columns: (i) Setting 1 with
$n=24,99$, (ii) Setting 2 with $n=24,99$, (iii) Setting 3 with $n=24,99$, (iv)
Setting 4 with $n=24,99$.

In the second simulation setup, we compared the efficiency of the
semiparametric model with that of the parametric model. We considered four
settings for this comparison– Setting 5: generating data from a linear setup
and using a gamma prior for $\alpha$, Setting 6: generating data from a linear
setup and employing a PC prior for $\alpha$, Setting 7: generating data from a
nonlinear setup and using a gamma prior for $\alpha$ and finally, Setting 8:
generating data from a nonlinear setup and employing a PC prior for $\alpha$.
We fitted the parametric and semiparametric models for each setting and
compared their goodness of fit. We also examined variations in hyperparameters
for each case. Specifically, in Settings 5 and 7, we considered the gamma
prior with parameters $(0.01,0.01)$ and $(1,1)$, respectively. In Settings 6
and 8, we employed the PC prior with parameters 2.5 and 5, respectively.
Additionally, we explored two cases for the number of samples, namely $n$ = 24
and $n$ = 99, with the same objective as in the previous setup. In every case
within this simulation setup, we have computed the absolute fitting error and
the WAIC of the model fit with the estimated parameters as part of the
comparison process.

All the simulations were conducted using three different values of $\alpha$,
the shape parameter: 0.5, 1, and 2. The choice of these values allowed us to
explore various scenarios; the value of $\alpha$ being 1 represented the true
scenario, resembling an exponential distribution. In contrast, the values of
$\alpha$ being 0.5 and 2 indicated deviations from the exponential-like
behavior. We utilized two different covariate sequences depending on the
sample size. In (5), we assume $P=1$ under all the simulation settings. When
the sample size is $n=24$, the covariate sequence we considered was
$\bm{X}=(0.04,0.08,\dots,0.96)^{\prime}$. For a larger sample size of $n=99$,
the covariate sequence was $\bm{X}=(0.01,0.02,\dots,0.99)^{\prime}$. To
generate data from the linear model, we constructed a design matrix with two
columns: one for the intercept and the other for $\bm{X}$. However, for
nonlinear data generation, we modified the second column of the design matrix
to include $\sin(2\pi\bm{X})$ instead of $\bm{X}$. During fitting the
parametric model, we used the same design matrix as the one employed during
the linear data generation process when fitting the data. However, we
introduced basis splines with ten basis functions for the semiparametric
setup. In total, we generated 1,000 datasets from the GE distribution. We
employed MCMC techniques to deduce the parameters, utilizing 4,000 MCMC
samples in total. The initial 2,000 samples were considered burn-in samples
and excluded from analysis, while a thinning interval of 5 was applied.

Figure 4 corresponds to the first simulation setup. The coverage probabilities
are the proportion of times the true value of $\alpha$ falls under the 95%
credible intervals of $\alpha$ obtained from MCMC. The top illustration of
Figure 4 presents the coverage probabilities of different simulation setups.
Four columns represent Settings 1, 2, 3, and 4, with two rows representing
sample sizes of 24 and 99, respectively. Each panel showcases four prior
specifications represented by different lines, while the X-axis represents the
different $\alpha$ values considered. In the bottom illustration of Figure 4,
the identical simulation setups are depicted, but instead of coverage
probabilities, it focuses on illustrating the absolute bias in the estimation
of $\alpha$.

Figure 5: Absolute fitting error (top illustration) and WAIC values (bottom
illustration) based on fitting a semiparametric GE regression model. For each
illustration, from the left-most to the right-most columns: (i) Setting 5 with
$n=24,99$, (ii) Setting 6 with $n=24,99$, (iii) Setting 7 with $n=24,99$, (iv)
Setting 8 with $n=24,99$. (For the bottom illustration, the lines representing
the gamma prior with parameters $(0.01,0.01)$ and the PC prior with $\theta=5$
are excluded since they closely resemble the demonstrated prior
specification).

Figure 5 focuses on the second simulation setup. In the top illustration of
Figure 5, we compare the goodness of fit between the parametric and
semiparametric models using the absolute fitting error. Similarly to previous
figures, there are four columns representing Settings 5, 6, 7, and 8, with two
rows representing sample sizes of 24 and 99. Each panel displays
specifications of different values of $\alpha$ at the X-axis, and lines
correspond to either the parametric setup or the semiparametric setup, with
varying hyperparameter specifications. Moving to the bottom illustration of
Figure 5, it provides a similar comparison but focuses on the WAIC of the
model fit instead of the absolute fitting error.

The two aforementioned simulation steps offer substantial validation for the
efficacy of the semiparametric regression setup and the utilization of the PC
prior. The top illustration in Figure 4 clearly illustrates that when the true
value of $\alpha$ approximates one, the PC prior exhibits superior coverage
probability compared to the conventional gamma prior. This pattern holds
across all four configurations in the initial simulation setup, except Setting
2, where both the gamma and PC priors yield undesirable outcomes due to
attempts to fit a parametric linear model to highly non-linear data.
Furthermore, the bottom illustration in Figure 4 highlights a reduction in
estimation bias when $\alpha$ equals one, aligning with the inherent
characteristic of the PC prior to shrink the estimate towards the base model.
Additionally, the lower row within the same figure empirically confirms the
well-established hypothesis that as the sample size increases, the influence
of the prior gradually diminishes. This pattern is evident as the lines
representing absolute bias nearly overlap, regardless of the chosen prior.

On the other hand, the top illustration in Figure 5 highlights a substantial
increase in the absolute fitting error as we transition from a dataset
demonstrating a linear trend to one illustrating a nonlinear trend. However,
the plots associated with settings involving the generation of linear data and
fitting a parametric model (Settings 5 and 6) exhibit minimal differences and
occasionally yield superior results. When examining the settings where data
originates from a nonlinear context (Settings 7 and 8), a distinct and
considerable gap emerges between the lines representing the parametric and
semiparametric models, where the latter consistently performs less favorably
than the former, occupying a lower position on the graph. The lower
illustration in Figure 5 conveys a similar observation. Despite Settings 5 and
6 pertaining to data generated from a linear setup, we observe that the
parametric model consistently exhibits lower WAIC scores. In Settings 7 and 8,
this ordering of WAIC scores is entirely reversed, with the lines representing
the semiparametric model falling consistently below those of the parametric
model. This reversal underscores the clear advantage of employing
semiparametric modeling, particularly in scenarios characterized by high
nonlinearity.

##  8 Data Application

We utilized the dataset introduced in Section 3 to drive our objective and
conducted data analysis. The dataset comprises daily average rainfall for the
wet days of the monsoon months during 1901–2022, focusing on Northern, Middle,
and Southern Western Ghats regions. As outlined in Section 3, to account for
the nonstationary nature of the dataset, we employed semiparametric regression
techniques, using rainfall as the response variable and the corresponding year
as the covariate. This analysis was performed separately for all three
regions. With the validation from Section 3, we assumed the daily rainfall
follows the GE distribution. Within a regression framework, we examined how
the rate parameter of the GE distribution correlates with the covariate.

###  8.1 Model Description

We formulate our model such that the rate parameter of the GE distribution in
the regression is influenced by the covariate ‘Year’ $(T)$. Then, if $Y_{t}$
represents our response variable of daily rainfall, our model is given as

$Y_{T}|T=t\sim GE(\alpha,\lambda(t))\text{ where }\lambda(t)\text{ is
dependent on the covariate }T=t.$

We conducted the analysis using two distinct models. The first model employs a
parametric approach to model the rate parameter, while the second model is our
proposed semiparametric formulation. We employed a simple linear regression
model for the rate parameter in the parametric setting, given as
$\lambda_{(L)}(t)$ in (9). On the other hand, for the semiparametric
regression, we adopted the suggested basis spline regression form presented in
(4). While Section 4.1 outlines the theory in a broader context, it is
important to note that in practical application, we lack access to multiple
covariate data. We have only one covariate (Year), and consequently, in this
section, we take $P=1$. With this, our model takes the form given as
$\lambda_{(NL)}(t)$ in (9) given by

$\lambda_{L}(t)=\exp(\beta_{0}+\beta_{1}~{}t),~{}\lambda_{NL}(t)=\exp\left[\sum_{k=1}^{K}\beta_{k}B_{k}(t)\right].$
(9)

We employed Bayesian methods to estimate the model parameters. As the
exploratory analysis shows a strong contraction of the considered GE
distribution towards its base model, the exponential one, we used the proposed
PC prior for the shape parameter $\alpha$ and chose independent flat Gaussian
priors for the regression parameters. As discussed in Section [6], we resorted
to the MCMC techniques to draw inferences about the model parameters.

We chose to utilize $K=12$ basis functions for our analysis. With data
spanning over 122 years, adopting 12 splines enabled us to capture decadal
patterns using each spline effectively. To optimize the hyperparameter
$\theta$ for the PC prior and achieve the most precise fit in our
semiparametric regression, we computed WAIC values across a range of $\theta$
values, spanning within the interval {0.5, 1, …, 5} for $\theta$. After
examining the northern, middle, and southern regions individually, we
identified the optimal values for $\theta$ that yielded the lowest WAIC
values. Specifically, for the Northern Western Ghats region, $\theta$ = 4.5
demonstrated the most favorable outcome. Similarly, for the middle Western
Ghats region, $\theta$ = 3.5 was identified as optimal, while for the southern
Western Ghats region, $\theta$ = 1.5 exhibited superior performance. As a
result, these chosen $\theta$ values were employed for their respective
regions during the final model fitting stage.

We fitted six models corresponding to the two models and all three regions.
For each of the model fits, we generated 10,000 MCMC samples for each model
parameter. The initial 3,000 samples were removed as burn-in and subsequently
excluded from the analysis. Additionally, we employed a thinning interval of
5. To evaluate convergence, mixing, and the behavior of the chains derived
from the MCMC process, we visualize the trace plots of the parameters
associated with both model fits in the supplementary materials. Specifically,
we present the trace plots of the shape parameter $\alpha$ for each region.
The regression parameters also exhibited similar satisfactory mixing behavior.

###  8.2 Model Comparison

In this section, we carry out a detailed comparison between the two model fits
mentioned in (9). These model fits were obtained using a Bayesian approach
through MCMC simulations, as elaborated upon in Section 8.1. The central focus
of this comparison is visually represented in Figure 6. In this figure, we
visually present the estimated mean daily rainfall on wet days for each year
spanning the period from 1901 to 2022. These estimates are provided for the
three regions separately, and both model fits – one characterized as
semiparametric (depicted in red) and the other as parametric (depicted in
blue). Furthermore, distinct panels are employed to illustrate various regions
(the top panel represents the Northern region, the middle panel corresponds to
the Middle region, and the bottom panel depicts the Southern region). In each
panel, we display the estimated trajectory alongside a bar diagram, which
offers a clear view of the annual average of daily rainfall on wet days for
each year.

Figure 6: Estimated mean of daily wet-day rainfall (in mm) with Semiparametric
(red line) and Parametric (blue line) models given by (9), along with
corresponding point-wise 95% credible intervals (ribbons). The top, middle,
and bottom panels show the results for the Northern, Middle, and Southern
Western Ghats regions, respectively.

Across all three regions, a noticeable trend emerges: the semiparametric
models exhibit a notably superior fit. This distinction becomes evident as we
observe multiple abrupt fluctuations in the bars representing the annual
averages of wet-day precipitation. Remarkably, the semiparametric model
effectively captures these fluctuations. Particularly noteworthy is the
ability of the semiparametric model to accurately capture the nonstationarity
present in the precipitation patterns. This heightened ability to encapsulate
the dynamic variations in precipitation is a notable strength of the
semiparametric model fitting. In Figure 6, we also provide a visual
representation of a 95% credible interval for the trajectory, estimated from
the MCMC samples. This interval allows us to understand the uncertainty
associated with the estimation process.

###  8.3 Inferences about Western Ghats Rainfall

The fitted model offers valuable insights into the intricate short-term trends
and features of rainfall in the Western Ghats region. Similar to the numerous
instances in the literature where the authors modeled rainfall data using an
exponential distribution [46, 16], our study also echoes this trend, aptly
captured by the PC prior. The estimated GE shape parameters (posterior means)
for the Northern, Middle, and Southern Western Ghats regions are 0.859, 0.949,
and 0.873, respectively. Corresponding posterior standard deviations for these
regions are 0.096, 0.100, and 0.097, respectively. These shape parameter
values indicate a pronounced alignment with the exponential distribution of
wet-day rainfall (the shape parameter is one). Consistent with the fluctuating
pattern in the annual average of daily wet-day rainfall, the fitted mean lines
for each region also demonstrate periodic crests and troughs. Besides, a
consistent and stable mean rainfall trend is noticeable across the Northern
and Middle Western Ghats regions. In the Southern Western Ghats region, the
fitted parametric and semiparametric models distinctly reveal a decaying
pattern in the annual averages of daily wet-day rainfall.

We present two significant insights into the rainfall patterns within these
regions: the overarching decade-long shifts and individual region-specific
probability rainfall plots. The calculation of the decadal change involves
determining the overall rainfall shift and dividing it by the number of
decades, resulting in $\\{\mu(2022)-\mu(1901)\\}/12.1$, where
$\mu(t)=\lambda(t)^{-1}\left[\psi(\alpha+1)-\psi(1)\right]$, and $t$
representing the corresponding year. In this equation, $\alpha$ represents the
estimated value of the shape parameter specific to the region, while
$\lambda(t)$ denotes the fitted rate parameter values for the given year, $t$.
Subsequently, the calculated decadal shifts in rainfall amount to 0.458 mm,
0.078 mm, and -0.367 mm for the northern, middle, and southern regions,
respectively. In Figure 7, the probability rainfall graphs are displayed for
three distinct probabilities: 0.3 (red line), 0.5 (blue line), and 0.7 (green
line). In agrometeorology, $100p$% probability rainfall means the $(1-p)^{th}$
quantile of the probability distribution of rainfall. These plots hold
significant implications in agriculture as they empower farmers to formulate
their harvesting strategies based on the anticipated likelihood of rainfall,
allowing them to align their plans with the rainfall patterns to fulfill their
specific requirements. Figure 7 showcases the estimated probability rainfall
graphs, along with pointwise 95% credible intervals for the estimated
rainfall. We derived these intervals from the MCMC samples; they illustrate
the uncertainty associated with the estimation process.

Figure 7: 30% (red line), 50% (blue line), and 70% (green line) probability
rainfall (in mm) with corresponding point-wise 95% credible intervals
(ribbons).

As a crucial component of our comprehensive analysis, we further discuss the
dynamic nature of annual average rainfall for the Western Ghats region over
the past century by exploring the plot of its rate of change. Interpreting
this quantity unveils insights into trends, variations, and shifts in mean
values over time, offering glimpses into rainfall behavior. A higher magnitude
implies swift changes, while a lower one indicates gradual shifts. A positive
rate denotes an increasing fitted mean rainfall over time, potentially
signaling rising annual average rainfall. Conversely, a negative rate
signifies a decreasing fitted mean, indicating a declining pattern and drier
conditions. A rate of change near zero indicates a stable fitted mean
rainfall. Fluctuations around zero imply short-term variations within a steady
range. We compute this quantity by taking the derivative of the fitted mean
from our semiparametric model with respect to the time component ($t$). This
process involves differentiating the cubic B-splines from $\lambda_{NL}(t)$ in
(9), given by

$\frac{\partial\mu(t)}{\partial
t}=-\frac{\psi(\alpha+1)-\psi(1)}{[\lambda(t)]^{2}}\sum_{k=1}^{K}\beta_{k}\frac{\partial
B_{k}(t)}{\partial t}$ (10)

where we computed the derivatives of the cubic B-splines using fda package
[37] in R. Figure 8 illustrates the plots depicting the rate of change over
the years for each of the three regions (the left panel for the Northern
region, the middle panel for the Middle region, and the right panel for the
Southern region).

Figure 8: Rate of change in annual average of daily wet-day Rainfall in the
monsoon months across the year 1920–2022, given by
$\frac{\partial\mu(t)}{\partial t}$ in (10). The black line represents zero
value.

For different regions of the Western Ghats mountain range, the data uncovers
varying trends in rainfall patterns. Initially, the Northern and Middle
regions exhibit more pronounced fluctuations in the rate-of-change graphs than
the Southern region. This pattern suggests more rapid variations in rainfall
trends in the Northern and Middle regions, while a more stable rainfall
pattern is visible for the Southern area. Moreover, the small-scale positive
and negative rate-of-change instances are well-balanced for the Northern and
Middle regions. This pattern implies that over the past century, changes in
rainfall have been relatively symmetric in terms of increase and decrease,
with no significant alterations in long-term patterns. In contrast, the
Southern region displays a substantial portion of years with graphs below the
zero line, signifying a prevalent decreasing trend in rainfall. The rate-of-
change in mean for the last 30 years shows consistent negative values in the
Southern sector, indicating the declining rainfall trend, while the graphs for
the other two regions consistently exhibit positive values, indicating an
increasing trend in rainfall over the past three decades in those areas. The
pointwise 95% credible intervals include the zero line for the Northern and
Middle regions; hence, the positive values for the last years are not
significant. On the other hand, while the posterior mean rate-of-change
remains negative for the Southern region in general, the credible intervals
indicate that the negative values of rate-of-change are significant for
several timestamps; however, the positive values are generally insignificant.

##  9 Discussions and Conclusions

The Western Ghats, a formidable mountain range running parallel to the western
coast of the Indian subcontinent, have a significant role in shaping
precipitation patterns in Southern India. This impact is especially notable
during the monsoon season, responsible for a substantial portion of the yearly
rainfall of the region that is essential for ecosystem vitality and
agricultural sustenance. The Western Ghats can be divided into Northern,
Middle, and Southern regions. The proposed semiparametric generalized
exponential (GE) regression model provides a reasonable fit for the wet-day
rainfall data for all three regions. The model allows the marginal
distributions to be the popular GE distribution.

With its shape and rate parameters, the GE distribution felicitates more
rigorous skewness attributes than several other distributions. Thus, it is a
better choice as a potential flexible model to incorporate high positive
skewness in the data. Additionally, depending on the shape parameter, the
varying nature of the hazard function provides GE distribution more
compatibility in fitting into complex data structures. In the regression
arena, semiparametric regression is a powerful statistical method that
combines the flexibility of the nonparametric models with the interpretability
and efficiency of the parametric models. The superiority of our proposed model
in capturing nonlinearity compared to the corresponding parametric model is
depicted in Sections 7 and 8. On the other hand, PC prior is a principled
distance-based prior that penalizes departure from a base model and is used
for specifying priors on parameters where it is difficult to elicit from
expert knowledge directly. This paper introduces a PC prior for the shape
parameter in the GE distribution, with the motivation of driving the GE
distribution closer to the characteristics of the exponential distribution, a
well-known probability distribution model for classical rainfall modeling.

There are several directions for extending this research. In addition to
modeling the rate parameter, we can consider treating the shape parameter as a
time-dependent variable. Instead of utilizing splines for the rate parameter,
an alternative approach could involve employing a Gaussian process prior.
Moreover, to ensure the applicability of our comparisons to large datasets, we
may explore various approximation techniques like Gaussian Markov random
fields [40]. While this paper has primarily focused on the temporal analysis
of rainfall data, further enhancements can be made by incorporating spatial
components [9]. This extension involves investigating the variability in
rainfall patterns across diverse geographical regions or watersheds [54].
Additionally, there is potential for developing a real-time rainfall
prediction system, offering timely information for tasks such as flood
forecasting, reservoir management, and emergency response, based on the
foundation provided by this model. For the high-dimensional spatial problems,
our model can be implemented as a two-stage model where the GE parameters can
be estimated at each spatial location ignoring the spatial structure and those
estimates can be smoothed using a Gaussian process [18].

From the application perspective, we observed a consistent overall trend with
periodic fluctuations in the Northern and Middle Western Ghats regions.
However, a clear declining trend was evident in the Southern Western Ghats
region. This observation is further supported by the decadal analysis of
rainfall changes in these three regions, where only the Southern region
exhibited a clear and significant negative value indicating the effects of
climate change. This research not only enhances our comprehension of the
intricate climatic dynamics within the Western Ghats but also emphasizes the
critical role of precise predictive models in anticipating seasonal rainfall
variations.

## Data availability statement

The dataset used in this paper can be downloaded (in a gridded data format)
from https://www.imdpune.gov.in/cmpg/Griddata/Rainfall_1_NetCDF.html.

## Disclosure statement

No potential conflict of interest was reported by the authors.

## References

  * Abe and James [2013] Abe, G. and James, E. (2013) Impacts of anthropogenic regulation on streamflow in the humid tropics of Western Ghat regions of Kerala state.  _International Journal of Advances in Engineering & Technology_ 6(4), 1895.
  * Ahmad _et al._ [2019] Ahmad, Z., Hamedani, G. G. and Butt, N. S. (2019) Recent developments in distribution theory: a brief survey and some new generalized classes of distributions.  _Pakistan Journal of Statistics and Operation Research_ pp. 87–110.
  * Aslam and Shahbaz [2007] Aslam, M. and Shahbaz, M. Q. (2007) Economic reliability test plans using the generalized exponential distribution.  _Journal of Statistics_ 14(1), 53–60.
  * Cota-Felix _et al._ [2009] Cota-Felix, J. E., Rivas-Davalos, F. and Maximov, S. (2009) An alternative method for estimating mean life of power system equipment with limited end-of-life failure data.  In _2009 IEEE Bucharest PowerTech_ , pp. 1–4.
  * Dey [2010] Dey, S. (2010) Bayesian estimation of the shape parameter of the generalised exponential distribution under different loss functions.  _Pakistan Journal of Statistics and Operation Research_ pp. 163–174.
  * Fahrmeir and Lang [2001] Fahrmeir, L. and Lang, S. (2001) Bayesian semiparametric regression analysis of multicategorical time-space data.  _Annals of the Institute of Statistical Mathematics_ 53, 11–30.
  * Gelfand and Kottas [2003] Gelfand, A. E. and Kottas, A. (2003) Bayesian semiparametric regression for median residual life.  _Scandinavian Journal of Statistics_ 30(4), 651–665.
  * Ghosal and Van der Vaart [2017] Ghosal, S. and Van der Vaart, A. (2017) _Fundamentals of nonparametric Bayesian inference_. Volume 44.  Cambridge: Cambridge University Press.
  * Gotway and Stroup [1997] Gotway, C. A. and Stroup, W. W. (1997) A generalized linear model approach to spatial data analysis and prediction.  _Journal of Agricultural, Biological, and Environmental Statistics_ pp. 157–178.
  * Gupta and Kundu [1999] Gupta, R. D. and Kundu, D. (1999) Theory & methods: Generalized exponential distributions.  _Australian & New Zealand Journal of Statistics_ 41(2), 173–188.
  * Gupta and Kundu [2001a] Gupta, R. D. and Kundu, D. (2001a) Exponentiated exponential family: an alternative to gamma and Weibull distributions.  _Biometrical Journal: Journal of Mathematical Methods in Biosciences_ 43(1), 117–130.
  * Gupta and Kundu [2001b] Gupta, R. D. and Kundu, D. (2001b) Generalized exponential distribution: different method of estimations.  _Journal of Statistical Computation and Simulation_ 69(4), 315–337.
  * Gupta and Kundu [2007] Gupta, R. D. and Kundu, D. (2007) Generalized exponential distribution: Existing results and some recent developments.  _Journal of Statistical Planning and Inference_ 137(11), 3537–3547.
  * Hastie and Tibshirani [1990] Hastie, T. J. and Tibshirani, R. J. (1990) _Generalized additive models_. Volume 43.  Boca Raton: CRC press.
  * Hazra [2022] Hazra, A. (2022) Minimum density power divergence estimation for the generalized exponential distribution.  _arXiv preprint arXiv:2206.08216_ .
  * Hazra _et al._ [2018] Hazra, A., Bhattacharya, S. and Banik, P. (2018) A Bayesian zero-inflated exponential distribution model for the analysis of weekly rainfall of the Eastern Plateau Region of India.  _Mausam_ 69(1), 19–28.
  * Hazra and Ghosh [2019] Hazra, A. and Ghosh, A. (2019) Robust statistical modeling of monthly rainfall: The minimum density power divergence approach.  _arXiv preprint arXiv:1909.08035_ .
  * Hazra _et al._ [2021] Hazra, A., Huser, R. and Jóhannesson, Á. V. (2021) Latent gaussian models for high-dimensional spatial extremes.  _arXiv preprint arXiv:2110.02680_ .
  * Hubert and Vandervieren [2008] Hubert, M. and Vandervieren, E. (2008) An adjusted boxplot for skewed distributions.  _Computational Statistics & Data Analysis_ 52(12), 5186–5201.
  * Jaheen [2004] Jaheen, Z. F. (2004) Empirical Bayes inference for generalized exponential distribution based on records.  _Communications in Statistics-Theory and Methods_ 33(8), 1851–1861.
  * Jorgensen [1982] Jorgensen, B. (1982) Statistical properties of the generalized inverse Gaussian distribution.  _Lecture Notes in Statistics_ 9.
  * Kim and Song [2010] Kim, C. and Song, S. (2010) Bayesian estimation of the parameters of the generalized exponential distribution from doubly censored samples.  _Statistical Papers_ 51, 583–597.
  * Koop and Poirier [2004] Koop, G. and Poirier, D. J. (2004) Bayesian variants of some classical semiparametric regression techniques.  _Journal of Econometrics_ 123(2), 259–282.
  * Kullback and Leibler [1951] Kullback, S. and Leibler, R. A. (1951) On information and sufficiency.  _The Annals of Mathematical Statistics_ 22(1), 79–86.
  * Kundu and Gupta [2008] Kundu, D. and Gupta, R. D. (2008) Generalized exponential distribution: Bayesian estimations.  _Computational Statistics & Data Analysis_ 52(4), 1873–1883.
  * Kundu and Raqab [2005] Kundu, D. and Raqab, M. Z. (2005) Generalized Rayleigh distribution: different methods of estimations.  _Computational Statistics & Data Analysis_ 49(1), 187–200.
  * Lee and Sison-Mangus [2018] Lee, J. and Sison-Mangus, M. (2018) A Bayesian semiparametric regression model for joint analysis of microbiome data.  _Frontiers in Microbiology_ 9, 522.
  * Li and Hanson [2014] Li, L. and Hanson, T. E. (2014) A Bayesian semiparametric regression model for reliability data using effective age.  _Computational Statistics & Data Analysis_ 73, 177–188.
  * Madi and Raqab [2007] Madi, M. T. and Raqab, M. Z. (2007) Bayesian prediction of rainfall records using the generalized exponential distribution.  _Environmetrics_ 18(5), 541–549.
  * Markiewicz _et al._ [2015] Markiewicz, I., Strupczewski, W. G., Bogdanowicz, E. and Kochanek, K. (2015) Generalized exponential distribution in flood frequency analysis for Polish rivers.  _PloS One_ 10(12), 1–26.
  * Marshall and Olkin [1997] Marshall, A. W. and Olkin, I. (1997) A new method for adding a parameter to a family of distributions with application to the exponential and Weibull families.  _Biometrika_ 84(3), 641–652.
  * Mathew _et al._ [2021] Mathew, M. M., Sreelash, K., Mathew, M., Arulbalaji, P. and Padmalal, D. (2021) Spatiotemporal variability of rainfall and its effect on hydrological regime in a tropical monsoon-dominated domain of Western Ghats, India.  _Journal of Hydrology: Regional Studies_ 36, 100861.
  * Naqash _et al._ [2016] Naqash, S., Ahmad, S. and Ahmed, A. (2016) Bayesian analysis of generalized exponential distribution.  _Journal of Modern Applied Statistical Methods_ 15(2), 38.
  * Nguyen-Huy _et al._ [2020] Nguyen-Huy, T., Deo, R. C., Mushtaq, S. and Khan, S. (2020) Probabilistic seasonal rainfall forecasts using semiparametric d-vine copula-based quantile regression.  In _Handbook of Probabilistic Models_ , pp. 203–227. Elsevier.
  * Ordoñez _et al._ [2023] Ordoñez, J. A., Prates, M. O., Bazán, J. L. and Lachos, V. H. (2023) Penalized complexity priors for the skewness parameter of power links.  _Canadian Journal of Statistics (Published online)_ .
  * Rajeevan _et al._ [2008] Rajeevan, M., Bhate, J. and Jaswal, A. K. (2008) Analysis of variability and trends of extreme rainfall events over India using 104 years of gridded daily rainfall data.  _Geophysical Research Letters_ 35(18).
  * Ramsay [2023] Ramsay, J. (2023) _fda: Functional Data Analysis_.  R package version 6.1.4.
  * Raqab and Madi [2005] Raqab, M. Z. and Madi, M. T. (2005) Bayesian inference for the generalized exponential distribution.  _Journal of Statistical Computation and Simulation_ 75(10), 841–852.
  * Reich and Ghosh [2019] Reich, B. J. and Ghosh, S. K. (2019) _Bayesian statistical methods_.  Boca Raton, USA: CRC Press.
  * Rue and Held [2005] Rue, H. and Held, L. (2005) _Gaussian Markov random fields: theory and applications_.  CRC press.
  * Rue _et al._ [2009] Rue, H., Martino, S. and Chopin, N. (2009) Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 71(2), 319–392.
  * Sarhan [2007] Sarhan, A. M. (2007) Analysis of incomplete, censored data in competing risks models with generalized exponential distributions.  _IEEE Transactions on Reliability_ 56(1), 132–138.
  * Simpson _et al._ [2017] Simpson, D., Rue, H., Riebler, A., Martins, T. G. and Sørbye, S. H. (2017) Penalising model component complexity: A principled, practical approach to constructing priors.  _Statistical Science_ 32(1), 1–28.
  * Sørbye and Rue [2017] Sørbye, S. H. and Rue, H. (2017) Penalised complexity priors for stationary autoregressive processes.  _Journal of Time Series Analysis_ 38(6), 923–935.
  * Tahir and Cordeiro [2016] Tahir, M. H. and Cordeiro, G. M. (2016) Compounding of distributions: a survey and new generalized classes.  _Journal of Statistical Distributions and Applications_ 3, 1–35.
  * Todorovic and Woolhiser [1974] Todorovic, P. and Woolhiser, D. (1974) Stochastic model of daily rainfall.  _Miscellaneous Publication, US Department of Agriculture_ 1275, 232 – 246.
  * Van Niekerk _et al._ [2021] Van Niekerk, J., Bakka, H. and Rue, H. (2021) A principled distance-based prior for the shape of the Weibull model.  _Statistics & Probability Letters_ 174, 109098.
  * Varikoden _et al._ [2019] Varikoden, H., Revadekar, J., Kuttippurath, J. and Babu, C. (2019) Contrasting trends in southwest monsoon rainfall over the Western Ghats region of India.  _Climate Dynamics_ 52, 4557–4566.
  * Veerabhadrannavar SA [2022] Veerabhadrannavar SA, V. B. (2022) Assessment of Impact of Climate Change in the Western Ghats Region, India.  _Indian Journal of Science and Technology_ pp. 1466–1472.
  * Venkatesh _et al._ [2021] Venkatesh, B., Nayak, P., Thomas, T., Jain, S. K. and Tyagi, J. (2021) Spatio-temporal analysis of rainfall pattern in the Western Ghats region of India.  _Meteorology and Atmospheric Physics_ 133, 1089–1109.
  * Ventrucci and Rue [2016] Ventrucci, M. and Rue, H. (2016) Penalized complexity priors for degrees of freedom in Bayesian p-splines.  _Statistical Modelling_ 16(6), 429–453.
  * Watanabe and Opper [2010] Watanabe, S. and Opper, M. (2010) Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory.  _Journal of Machine Learning Research_ 11, 3571–3594.
  * Wigena _et al._ [2015] Wigena, A. H., Djuraidah, A. and Rizki, A. (2015) Semiparametric modeling in statistical downscaling to predict rainfall.  _Applied Mathematical Sciences_ 9(88), 4371–4382.
  * Yang _et al._ [2005] Yang, C., Chandler, R., Isham, V. and Wheater, H. (2005) Spatial-temporal rainfall simulation using generalized linear models.  _Water Resources Research_ 41(11).

##  Appendix A Kullback-Leibler divergence between generalized exponential
(GE) and exponential distribution

The density function of a random variable following a GE distribution with
shape parameter $\alpha$ and rate parameter $\lambda$ is given by
$f(y)=\alpha\lambda\left(1-e^{-\lambda y}\right)^{\alpha-1}e^{-\lambda
y},y>0$, and that of a random variable following an exponential distribution
with rate parameter $\lambda$ is given as $g(y)=\lambda e^{-\lambda y},y>0.$
Thus, the KLD between the GE distribution and the exponential distribution can
be obtained as

$\displaystyle\text{KLD}(f\parallel g)$ $\displaystyle:=\textrm{KLD}(\alpha)$
$\displaystyle=\int_{0}^{\infty}\log\left(\frac{f(y)}{g(y)}\right)\cdot
f(y)\,dy$
$\displaystyle=\log(\alpha)\int_{0}^{\infty}f(y)\,dy+(\alpha-1)\int_{0}^{\infty}\log\left(1-e^{-\lambda
y}\right)\cdot f(y)\,dy$
$\displaystyle=\log(\alpha)+(\alpha-1)\int_{0}^{\infty}\log\left(1-e^{-\lambda
y}\right)\cdot\alpha\lambda\left(1-e^{-\lambda y}\right)^{\alpha-1}e^{-\lambda
y}\,dy$ $\displaystyle=\log(\alpha)-(\alpha-1)\int_{0}^{\infty}\alpha x\cdot
e^{-x(\alpha-1)}\cdot e^{-x}\,dx\hskip
5.0pt\left[\textrm{replacing}~{}e^{-x}=1-e^{-\lambda y}\right]$
$\displaystyle=\log(\alpha)+\frac{(1-\alpha)}{\alpha}.$


$\displaystyle\delta_{N}(A,\epsilon,\phi):=\sum_{j=1}^{\log_{16}N-1}(CA\epsilon)^{j}\|S_{N/16^{j}}^{\perp}\phi\|_{H_{x}^{1}}+(CA\epsilon)^{\log_{16}N-1}\|\phi\|_{H_{x}^{1}}.$
(A.13)

Since $\|S_{N}^{\perp}\phi\|_{H_{x}^{1}}\rightarrow 0$, uniformly on any
compact set of $H^{1}$, we deduce that $\delta_{N}(A,\epsilon,\phi)$ converges
to $0$, uniformly on a compact set of $H^{1}(\mathbb{T}^{3})$. Plugging into
(A.11), we obtain that

$\displaystyle\|u_{N}-u\|_{X^{1}(I)}\leq
A_{0}\|\phi-\widetilde{\phi}\|_{H_{x}^{1}}+A_{0}\delta_{N}(A,\epsilon,\phi).$
(A.14)

By the embedding $X^{1}(I)\hookrightarrow L^{\infty}(I;H_{x}^{1})$ and the
fact that $u_{N}(t),u(t)\in C(\overline{I};H_{x}^{1})$, we obtain that

$\sup_{t\in\overline{I}}\|u_{N}(t)-u(t)\|_{H_{x}^{1}}\leq
C_{0}\|\phi-\widetilde{\phi}\|_{H_{x}^{1}}+C_{0}\delta_{N}(A,\epsilon,\phi),$

for some absolute constant $C_{0}>0$. This completes the proof of Lemma A.4
when $\sigma=1$.

The general case $\sigma>1$ follows from similar analysis. Here we only
indicate the necessary modification in the proof. Indeed, in Step 1, we
replace $X^{1},H^{1}$ norms by $X^{\sigma},H^{\sigma}$ norms in (A.2), thanks
to (A.5). In Step 2, all inequalities remain unchanged when replacing all
$X^{1},H^{1},N^{1}$ norms by $X^{\sigma},H^{\sigma},N^{\sigma}$ norms (up to
change the numerical constants $C$ in front of each inequality), thanks to
(A.5) and the trivial embedding $X^{\sigma}\hookrightarrow X^{1}$. This
completes the proof Lemma A.4. ∎

######  Definition A.8 ($(A,\delta)$-partition).

Let $A>\delta>0$. Given an interval $[-T,T]$ and $\phi\in
H^{1}(\mathbb{T}^{3})$, we define an $(A,\delta)$-partition with respect to
$\phi$ a collection of finite intervals $([\tau_{j-1},\tau_{j}])_{j=1}^{m}$
such that

$-T=\tau_{0}<\tau_{2}<\cdots<\tau_{m}=T,\quad\|\Phi(\tau_{j-1})\phi\|_{H_{x}^{1}}\leq
A,\;\|\Phi(t)\phi\|_{Z([\tau_{j-1},\tau_{j}])}<\delta,\;\forall j=1,\cdots,m.$

We collect some basic properties. The following property is immediate:

######  Proposition A.9 (Refinement of $(A,\delta)$-partition).

Any refinement of an $(A,\delta)$-partition with respect to $\phi$ is an
$(A,\delta)$-partition.

######  Proposition A.10.

Assume that $([\tau_{j-1},\tau_{j}])_{j=1}^{m}$ is an $(A,\delta)$-partition
with respect to $\phi$. Then for sufficiently small $\delta=\delta(A)>0$,

$\|{\mathrm{e}}^{i(t-\tau_{j-1})\Delta}\Phi(\tau_{j-1})\phi\|_{Z([\tau_{j-1},\tau_{j}])}<2\delta,\;\forall
j=1,\cdots,m.$

###### Proof.

Denote $u(t)=\Phi(t)\phi$, then

${\mathrm{e}}^{i(t-\tau_{j-1})\Delta}u(\tau_{j-1})=u(t)-\frac{1}{i}\int_{\tau_{j-1}}^{t}{\mathrm{e}}^{i(t-t^{\prime})\Delta}(|u(t^{\prime})|^{4}u(t^{\prime}))dt^{\prime}.$

The desired consequence follows from the Strichartz inequality as in the proof
of Lemma A.4. Hence we omit the detail. ∎

######  Proposition A.11 ($H^{1}$-Stability of an $(A,\delta)$-partition).

Assume that $([\tau_{j-1},\tau_{j}])_{j=1}^{m}$ is an $(A,\delta)$-partition
with respect to $\phi$. There exists $\delta_{1}>0$ such that for any
$\delta<\delta_{1}$, there exists $\epsilon_{0}=\epsilon_{0}(m,A,\delta)>0$,
such that for any $\phi_{1}$ in the $\epsilon_{0}$-neighborhood of $\phi$
(with respect to the $H^{1}$-topology), $([\tau_{j-1},\tau_{j}])_{j=1}^{m}$ is
an $(2A,2\delta)$-partition for $\phi_{1}$.

###### Proof.

Denote $u=\Phi(t)\phi$ and $u_{1}=\Phi(t)\phi_{1}$. Fix an interval
$I_{j}=[\tau_{j-1},\tau_{j}]$. By the Strichartz inequality,

$\|{\mathrm{e}}^{i(t-\tau_{0})\Delta}(\phi-\phi_{1})\|_{Z(I_{1})}\leq
C_{1}\|{\mathrm{e}}^{i(t-\tau_{0})\Delta}(\phi-\phi_{1})\|_{X^{1}(I_{1})}\leq
C_{1}\|\phi-\phi_{1}\|_{H_{x}^{1}}\leq C_{1}\epsilon_{0}$

for some absolute constant $C_{1}>1$. Pick $\delta<\delta_{1}$ as in Lemma
A.4, and for $\epsilon_{0}\ll\delta$ (so that $C_{1}\epsilon_{0}<\delta$) and
for any $\phi_{1}$ in an $\epsilon_{0}$-neighborhood of $\phi$ (with respect
to the $\dot{H}^{1}$-topology), the solution $u_{1}$ (which exists thanks to
Lemma A.4) with initial data $\phi_{1}$ satisfies

$\|u_{1}\|_{X^{1}(I_{1})}\leq\frac{3}{2}A,\;\|u_{1}\|_{Z(I_{1})}<\frac{3}{2}\delta.$

Applying Lemma A.2, we have

$\displaystyle\|u-u_{1}\|_{X^{1}(I_{1})}\leq$
$\displaystyle\|{\mathrm{e}}^{it\Delta}(\phi-\phi_{1})\|_{X^{1}(I_{1})}+\big{\|}|u|^{4}u-|u_{1}|^{4}u_{1}\big{\|}_{N^{1}(I_{1})}$
$\displaystyle\leq$
$\displaystyle\|\phi-\phi_{1}\|_{H_{x}^{1}}+C\|u-u_{1}\|_{X^{1}(I_{1})}(\|u\|_{Z(I_{1})}^{2}+\|u_{1}\|_{Z(I_{1})}^{2})(\|u\|_{X^{1}(I_{1})}^{2}+\|u_{1}\|_{X^{1}(I_{1})}^{2})$
$\displaystyle\leq$
$\displaystyle\|\phi-\phi_{1}\|_{H_{x}^{1}}+CA^{2}\delta^{2}\|u-u_{1}\|_{X^{1}(I_{1})}.$

Taking $\delta>0$ small enough such that $1-CA^{2}\delta^{2}>\frac{1}{2}$, we
deduce that

$\|u-u_{1}\|_{Z(I_{1})}\leq C_{1}\|u-u_{1}\|_{X^{1}(I_{1})}\leq
2C_{1}\|\phi-\phi_{1}\|_{H_{x}^{1}}<2C_{2}\epsilon_{0}.$

In particular, by the embedding property $X^{1}(I)\hookleftarrow
L^{\infty}(I;H_{x}^{1})$ , for almost every
$\tau_{1}^{*}\in(\tau_{0},\tau_{1})$,

$\|u(\tau_{1}^{*})-u_{1}(\tau_{1}^{*})\|_{H_{x}^{1}}\leq
C_{2}\|u-u_{1}\|_{X^{1}(I_{1})}\leq
2C_{2}\|\phi-\widetilde{\phi}\|_{H_{x}^{1}}\leq 2C_{2}\epsilon_{0},$

where $C_{2}>0$ is another absolute constant. By the continuity of the flows
$t\mapsto u(t),u(t_{1})$, we have

$\|u(\tau_{1})-u_{1}(\tau_{1})\|_{H_{x}^{1}}\leq
2C_{2}\epsilon_{0},\quad\|u_{1}\|_{Z(I_{1})}\leq\|u\|_{Z(I_{1})}+\|u-u_{1}\|_{Z(I_{1})}<\delta+2C_{1}\epsilon_{0},$

hence

$\|{\mathrm{e}}^{i(t-\tau_{1})\Delta}(u(\tau_{1})-u_{1}(\tau_{1}))\|_{Z(I_{2})}\leq
C_{1}\|u(\tau_{1})-u_{1}(\tau_{1})\|_{H_{x}^{1}}\leq 2C_{1}C_{2}\epsilon_{0}.$

By choosing $\epsilon_{0}$ small enough such that
$\epsilon_{0}\sum_{j=1}^{m}(2C_{1}C_{2})^{j}<\delta$, we can repeat the
argument above until $I_{m}$. In particular,
$\|u_{1}\|_{Z(I_{j})}<2\delta,\|u_{1}\|_{X^{1}(I_{j})}<2A$ for all
$j=1,2,\cdots,m$ This implies that $([\tau_{j-1},\tau_{j}])_{j=1}^{m}$ is an
$(2A,2\delta)$-partition with respect to $\phi_{1}$. The proof of Proposition
A.11 is complete. ∎

Now we are ready to prove:

######  Proposition A.12 (Long-time approximation).

Given $T\geq 1$ and $\phi\in H^{1}(\mathbb{T}^{3})$.

$\lim_{N\rightarrow\infty}\|\Phi_{N}(t)\phi-\Phi(t)\phi\|_{H^{1}(\mathbb{T}^{3})}=0,\forall
t\in[-T,T],$

uniformly for $\phi$ on a compact set $K\subset H^{1}(\mathbb{T}^{3})$.
Moreover, for any $|t|\leq T$ and $N\in\mathbb{N}$, the sets
$\Phi(t)(K),\Phi_{N}(t)(K)$ are compact in $H^{1}(\mathbb{T}^{3})$.

###### Proof.

Fix $A>0,T>0$ and $K\subset B^{H^{1}}_{A/2}(0):=\\{\phi:\;\|\phi\|_{H^{1}}\leq
A/2\\}$ a compact set of $H^{1}(\mathbb{T}^{3})$. Note that for any $\phi\in
B_{A/2}^{H^{1}}$, $H[\phi]\leq C_{0}A^{2}$. We divide the proof into several
steps. Set $A_{1}=4\sqrt{C_{0}}A$.

Step 1: Existence of a uniform $(A_{1},\delta)$-partition:

Thanks to Theorem A.1, for any $\phi\in K$ in $H^{1}(\mathbb{T}^{3})$,
$\|\Phi(t)\phi\|_{Z([-T,T])}\leq\Lambda(C_{0}A^{2},T)<\infty$. In particular,
there exists an $(\frac{A_{1}}{2},\frac{\delta}{2})$-partition
$([\tau_{j-1},\tau_{j}])_{j=1}^{m}$ with respect to $\phi$. By stability
(Proposition A.11), there exists
$\epsilon_{0}=\epsilon_{0}(m,A_{1},\delta)>0$, such that
$([\tau_{j-1},\tau_{j}])_{j=1}^{m}$ is an $(A_{1},\delta)$-partition with
respect to all $\phi_{1}\in B_{\epsilon_{0}}^{H^{1}}(\phi)$. Since $K$ is
compact, there exist finitely many $\phi_{1},\cdots,\phi_{n}\in K,$
$\epsilon_{i}>0,i=1,\cdots,n$ and $(A,\delta)$-partitions
$([\tau^{(i)}_{j-1},\tau^{(i)}_{j}])_{j=1}^{m_{i}},i=1,\cdots,n$, such that

  * (1)

$K\subset\bigcup_{i=1}^{n}B_{\epsilon_{u}}^{H^{1}}(\phi_{i})$;

  * (2)

$([\tau^{(i)}_{j-1},\tau^{(i)}_{j}])_{j=1}^{m_{i}}$ is an
$(A_{1},\delta)$-partition for all $\phi\in
B_{\epsilon_{i}}^{H^{1}}(\phi_{i})$, $i=1,\cdots,n$.

Consider a refinement $([\tau_{j-1},\tau_{j}])_{j=1}^{m}$ of partitions
$([\tau^{(i)}_{j-1},\tau^{(i)}_{j}])_{j=1}^{m_{i}}$. By Proposition A.9,
$([\tau_{j-1},\tau_{j}])_{j=1}^{m}$ is a uniform $(A_{1},\delta)$-partition
with respect to all $\phi\in K$.

Step 2: Long-time convergence:

Now we are able to iterate Lemma A.4 (with the small parameter $2\delta$
instead of $\delta$ in the statement) from $I_{1}=[\tau_{0},\tau_{1}]$ to
$I_{m}=[\tau_{m-1},\tau_{m}]$. Thanks to Proposition A.10 and the energy
conservation law, we have

$\|{\mathrm{e}}^{i(t-\tau_{j-1})\Delta}\Phi(\tau_{j-1})\phi\|_{Z(I_{j})}<2\delta,\;\|\Phi(\tau_{j-1})\phi\|_{H_{x}^{1}}\leq
A_{1}.$

In order to apply Lemma A.4 on each interval $I_{j}$ with initial data
$\Phi_{N}(\tau_{j-1})\phi$ and $\Phi(\tau_{j-1})\phi$, we have to ensure that

$\displaystyle\|\Phi_{N}(\tau_{j-1})\phi-\Phi(\tau_{j-1})\phi\|_{H_{x}^{1}}<2\delta.$
(A.15)

First, since $\delta_{N}(A_{1},2\delta,\phi)\rightarrow 0$ uniformly on the
compact set $K$, we can choose $N_{0}$ large enough, such that for all $N\geq
N_{0},\phi\in K$,

$(C_{0}+C_{0}^{2}+\cdots+C_{0}^{m})\delta_{N}(A_{1},2\delta,\phi)<\delta.$

Now we argue by induction that

$\displaystyle\|\Phi_{N}(t)\phi-\Phi(t)\phi\|_{C(\overline{I_{j}};H_{x}^{1})}\leq(C_{0}+\cdots+C_{0}^{j})\delta_{N}(A_{1},2\delta,\phi).$
(A.16)

. When $j=1$, $\Phi_{N}(\tau_{0})\phi=\Phi(\tau_{0})\phi$, and from the last
assertion of Lemma A.4, we have

$\displaystyle\|\Phi_{N}(t)\phi-\Phi(t)\phi\|_{C(\overline{I_{1}};H_{x}^{1})}\leq
C_{0}\delta_{N}(A_{1},2\delta,\phi).$

Assume that (A.16) holds for some $j\geq 1$, in particular, (A.15) holds
thanks to our choice. Then we are able to apply Lemma A.4 on the time interval
$I_{j+1}$ to obtain that

$\displaystyle\|\Phi_{N}(t-\tau_{j-1})$
$\displaystyle\Phi_{N}(\tau_{j-1})\phi-\Phi(t-\tau_{j-1})\Phi(\tau_{j-1})\phi\|_{C(\overline{I_{j}};H_{x}^{1})}$
$\displaystyle\leq$ $\displaystyle
C_{0}\|\Phi_{N}(\tau_{j-1})\phi-\Phi(\tau_{j-1})\phi\|_{H_{x}^{1}}+C_{0}\delta_{N}(A_{1},2\delta,\phi)$
$\displaystyle\leq$
$\displaystyle(C_{0}+C_{0}^{2}+\cdots+C_{0}^{j+1})\delta_{N}(A_{1},2\delta,\phi).$

Hence (A.16) holds for all $j=1,2,\cdots,m$. In particular, we have

$\|\Phi_{N}(t)\phi-\Phi(t)\phi\|_{C([-T,T];H_{x}^{1})}\leq(C_{0}+\cdots+C_{0}^{m})\delta_{N}(A_{1},2\delta,\phi)$

which converges to $0$, as $N\rightarrow\infty$, uniformly in $\phi\in K$.

This completes the proof of Proposition A.12. ∎

Now we are ready to accomplish the proof of Proposition 3.1 and Proposition
3.2:

###### Proof of Proposition 3.1.

If $\sigma=1$, the $H^{1}$-uniform bound for $\Phi_{N}(t)\phi$ and
$\Phi(t)\phi$ follows from the defocusing feature of (1.1) and the
conservation of energy. Now we assume that $\sigma>1$. By the compact
embedding $H^{\sigma}(\mathbb{T}^{3})\hookrightarrow H^{1}(\mathbb{T}^{3})$,
the ball $B_{R}^{H^{\sigma}}$ is compact with respect to the $H^{1}$-topology.
By the same argument as Step 1 in the proof of Proposition A.12, there exists
a uniform $(A,\delta)$-partition (with $A=R$ here and
$\delta<\frac{1}{\sqrt{2C(R,T)}}$ as in the proof of Corollary A.2)
$([\tau_{j-1},\tau_{j}])_{j=1}^{m}$ of $[-T,T]$, where $m$ depends only on the
$R,T$ and $\sigma$. Repeating the analysis in the proof of Corollary A.2, we
obtain that for all $|t|\leq T$ and $N\in\mathbb{N}$,

$\|\Phi(t)\phi\|_{H_{x}^{\sigma}}+\|\Phi_{N}(t)\phi\|_{H_{x}^{\sigma}}\leq
C^{m}\|\phi\|_{H_{x}^{\sigma}}.$

This completes the proof of Proposition 3.1. ∎

###### Proof of Proposition 3.2.

We assume that $\sigma>1$, otherwise, the proof is completed as Proposition
A.12. Let $K$ be a compact set of $H^{\sigma}(\mathbb{T}^{3})$. In particular,
$K$ is bounded of $H^{\sigma}(\mathbb{T}^{3})$ and compact with respect to the
$H^{1}(\mathbb{T}^{3})$-topology. To prove the uniform convergence on $K$, we
follow the same scheme of analysis as in the previous section. By Proposition
3.1, there exists a constant $D(K,T)$ depending only on $T>0$ and the compact
set $K$ in $H^{\sigma}$, such that

$\displaystyle\sup_{t\in[-T,T]}\|\Phi(t)\phi\|_{H^{\sigma}}+\sup_{t\in[-T,T]}\|\Phi_{N}(t)\phi\|_{H^{\sigma}}\leq
D(K,T)$ (A.17)

for all $N\in\mathbb{N}$. At this stage, we are able to repeat the argument as
in the proof of Proposition A.12 line by line if we replace the norms
$H^{1},X^{1},N^{1}$ everywhere by $H^{\sigma},X^{\sigma},N^{\sigma}$ and the
constants $A$ by $D(K,T)$. We omit the details and conclude. ∎

## References

  * [1] M. Aizenman, H. Duminil-Copin, Marginal triviality of the scaling limits of critical $4D$ Ising and $\phi^{4}$ models, Ann. of Math. (2), 194 (2021) 163–235,.
  * [2] H. Bahouri, P. Gérard, High frequency approximation of solutions to critical nonlinear wave equations, Amer. J. Math.121 (1999)131–175.
  * [3] J. Bourgain, Invariant measures for the 2D-defocusing nonlinear Schrödinger equation, Comm. Math. Phys. 176 (1996) 421-445.
  * [4] J. Bourgain, Global wellposedness of defocusing critical nonlinear Schrödinger equation in the radial case, J. Amer. Math. Soc. 12 (1999), 145–171.
  * [5] B. Bringmann, Almost sure local well-posedness for a derivative nonlinear wave equation, Int. Math. Res. Not. (2021), 8657–8697.
  * [6] B. Bringmann, Y. Deng, A. Nahmod, H. Yue, Invariant Gibbs measures for the three dimensional cubic nonlinear wave equation, arXiv:2206.07589.
  * [7] N. Burq, L. Thomann, Almost sure scattering for the one dimensional nonlinear Schrödinger equation, to appear in Mem. Amer. Math. Soc.
  * [8] N. Burq, L. Thomann, N. Tzvetkov, Long time dynamics for the one dimensional non linear Schrödinger equation, Ann. Inst. Fourier. 63 (2013), no. 6, 2137–2198.
  * [9] J. Colliander, M. Keel, G. Staffilani, H. Takaoka, T. Tao, Global well-posedness and scattering for the energy-critical nonlinear Schrödinger equation in $\mathbb{R}^{3}$, Ann. of Math. 167 (2008), 767–865.
  * [10] C. Fan, Y. Ou, G. Staffilani, H. Wang, 2D-defocusing nonlinear Schrödinger equation with random data on irrational tori, Stoch. Partial Differ. Equ. Anal. Comput. 9 (2021), 142–206.
  * [11] A. Debussche, Y. Tsustumi, Quasi-Invariance of Gaussian Measures Transported by the Cubic NLS with Third-Order Dispersion on $\mathbb{T}$, arXiv:2002.04899 (2020).
  * [12] Y. Deng, A. Nahmod, H. Yue, Invariant Gibbs measures and global strong solutions for 2D NLS, arXiv:1910.08492.
  * [13] Y. Deng, A. Nahmod, H. Yue, Random tensors, propagation of randomness, and nonlinear dispersive equations. Invent. Math. 228 (2022), no. 2, 539–-686.
  * [14] J. Forlano, W. Trenberth, On the transport of Gaussian measures under the one-dimensional fractional nonlinear Schrödinger equation, Ann. Inst. Henri Poincare, Anal. Non Lineaire, 36 (2019), 1987-2025.
  * [15] J. Forlano, K. Seong, Transport of Gaussian measures under the flow of one-dimensional fractional nonlinear Schrödinger equations, Comm. PDE (2022)1296–1337.
  * [16] J. Forlano, L. Tolomeo, Quasi-invariance of Gaussian measures of negative regularity for fractional nonlinear Schrödinger equations , arXiv:2205.11453.
  * [17] L. Friedlander, An invariant measure for the equation $u_{tt}-u_{xx}+u^{3}=0$, Comm. Math. Phys., 98 (1985) 1–16.
  * [18] G. Genovese, R. Lucà, N. Tzvetkov, Transport of Gaussian measures with exponential cut-off for Hamiltonian PDEs, arXiv:2103.04408 [math.AP] (2021).
  * [19] G. Genovese, R. Lucà, N. Tzvetkov, Quasi-invariance of Gaussian measures for the periodic Benjamin-Ono-BBM equation, Stoch. Partial Differ. Equ. Anal. Comput. 11 (2023), no. 2, 651–684.
  * [20] J. Glimm, A. Jaffe, Quantum physics. A functional integral point of view, Second edition. Springer-Verlag, New York, 1987. xxii+535 pp.
  * [21] T. Gunaratnam, T. Oh, N. Tzvetkov, H. Weber, Quasi-invariant Gaussian measures for the nonlinear wave equation in three dimensions, Probab. Math. Phys. 3 (2022), no. 2, 343–379.
  * [22] S. Herr, D. Tataru, N. Tzvetkov, Global well-posedness of the energy-critical nonlinear Schrödinger equation with small initial data in $H^{1}(\mathbb{T}^{3})$, Duke Math. J. 159 (2011), no. 2, 329–349.
  * [23] A. Ionescu, B. Pausader, The energy-critical defocusing NLS on $\mathbb{T}^{3}$. Duke Math. J. 161 (2012), no. 8, 1581–1612.
  * [24] C. E. Kenig and F. Merle, Global well-posedness, scattering and blow-up for the energy-critical, focusing, non-linear Schrödinger equation in the radial case, Invent. Math. 166 (2006), 645–675.
  * [25] F. Planchon, N. Tzvetkov, N. Visciglia, Transport of Gaussian measures by the flow of the nonlinear Schrödinger equation, Mathematische Annalen (2020) 378:389–423.
  * [26] F. Planchon, N. Tzvetkov, N. Visciglia, Modified energies for the periodic generalized KdV equation and applications, arXiv:2107.01926, to appear in Ann. Inst. Henri Poincaré, Anal. Non Linéaire.
  * [27] P. Malliavin, Stochastic calculus of variations and hypoelliptic operators, In: Proc. Inter. Symp. on Stocha. Diff. Equations, Kyoto 1976, Wiley 1978, 195–263.
  * [28] A. Nahmod, G. Staffilani, Almost sure well-posedness for the periodic 3D quintic nonlinear Schrödinger equation below the energy space, J. Eur. Math. Soc., 17 (2015), no. 7, 1687–1759.
  * [29] D. Nualart, The Malliavin Calculus and Related Topics , 2nd edition, Probability and its Applications (New York). Springer-Verlag, Berlin, 2006. xiv+382 pp. ISBN: 978-3-540-28328-7; 3-540-28328-5.
  * [30] T. Oh, T. Robert, N. Tzvetkov, Stochastic nonlinear wave dynamics on compact surfaces, Annales Henri Lebesgue 6 (2023) 161–223.
  * [31] T. Oh, K. Seong, Quasi-invariant Gaussian measures for the cubic fourth order nonlinear Schrödinger equation in negative Sobolev spaces, J. Funct. Anal. 281 (2021), no. 9, 109150, 49 pp.
  * [32] T. Oh, P. Sosoe, N. Tzvetkov, An optimal regularity result on the quasi-invariant Gaussian measures for the cubic fourth order nonlinear Schrödinger equation, J. Ec. Polytechnique, Math., 5 (2018), 793-841.
  * [33] T. Oh, L. Thomann, Invariant Gibbs measures for the 2-d defocusing nonlinear wave equations, Ann. Fac. Sci. Toulouse Math. (6) 29 (2020), 1–26.
  * [34] T. Oh, Y. Tsutsumi and N. Tzvetkov, Quasi-invariant Gaussian measures for the cubic nonlinear Schroedinger equation with third-order dispersion, C. R. Acad. Sci. Paris, Ser. I, 357 (2019), 366-381.
  * [35] T. Oh, N. Tzvetkov, Quasi-invariant Gaussian measures for the cubic fourth order nonlinear Schrödinger equation, Probab. Theory Relat. Fields (2017) 169:1121–1168.
  * [36] T. Oh, N. Tzvetkov, Quasi-invariant Gaussian measures for the two-dimensional defocusing cubic nonlinear wave equation, J. Eur. Math. Soc. 22 (2020), no. 6, 1785–1826.
  * [37] B. Simon, The $P(\phi)_{2}$ Euclidean (quantum) field theory, Princeton Series in Physics. Princeton University Press, Princeton, N.J., 1974. xx+392 pp.
  * [38] N. Tzvetkov, Quasi-invariant Gaussian measures for one dimensional Hamiltonian PDEs, Forum Math. Sigma 3 (2015), e28, 35 pp.
  * [39] P. Sosoe, W. Trenberth, T. Xian, Quasi-invariance of fractional Gaussian fields by the nonlinear wave equation with polynomial nonlinearity, Differential Integral Equations 33 (2020), no. 7-8, 393–430.
  * [40] C. Sun, N. Tzvetkov, Gibbs measure dynamics for the fractional nonlinear Schrödinger equation, SIAM J. Math. Anal., 52(5), 4638–4704.
  * [41] C. Sun, N. Tzvetkov, Refined probabilistic well-posedness for the weakly dispersive NLS, Nonlinear Analysis, 213(11) 112530.
  * [42] N. Visciglia, The modified energy technique and applications, Boll. Unione Mat. Ital.14 (2021) 3–16.