text:
\begin{document}

\title[PCE for general multivariate distributions]{Polynomial Chaos Expansion for general multivariate distributions with correlated variables}
\author{Maria Navarro, Jeroen Witteveen, Joke Blom$\,^\star$}

\date{\today}

\begin{abstract}
Recently, the use of Polynomial Chaos Expansion (PCE) has been increasing to study the uncertainty
in mathematical models for a wide range of applications and several extensions of the
original PCE technique have been developed to deal with some of its limitations.
But as of to date PCE methods still have the restriction that the random variables have to be statistically
independent.
This paper presents a method to construct a basis of the probability space of orthogonal polynomials
for general multivariate distributions with correlations between the random input variables.
We show that, as for the current PCE methods, the statistics like mean, variance and Sobol' indices can be
obtained at no significant extra postprocessing costs. We study the behavior of the proposed method for a range of correlation
coefficients for an ODE with model parameters that follow a bivariate normal distribution.  In all cases
the convergence rate of the proposed
method is analogous to that for the independent case. Finally, we show, for a canonical enzymatic reaction,
how to propagate experimental errors through the process of fitting parameters to a probabilistic
distribution of the quantities of interest, and we demonstrate the significant difference in the results assuming
independence or full correlation compared to taking into account the true correlation.

\vspace{2mm}

\noindent {\it Keywords.}
Uncertainty quantification $\star$ polynomial chaos expansion $\star$ correlated inputs

\vspace{2mm}

\noindent CWI, Science Park 123, 1098 XG Amsterdam, the Netherlands

\noindent $^\star$Corresponding author. Tel.:+31 205924263

\noindent {\it Email addresses.} {\scriptsize \tt{maria.navarro@cwi.nl}, \tt{jeroen.witteveen@cwi.nl}, \tt{joke.blom@cwi.nl}}

\end{abstract}

\maketitle

\section{Introduction}
To describe real-life phenomena we often make use of deterministic mathematical models.
Typical constituents of such models are assigned a definite value and we seek a deterministic solution
to the problem.
In reality, however, those phenomena will almost always have uncertain
components: unknown parameters, imprecise experimental data, etc.. A precise mathematical description should
reflect such uncertainties. In other words, the parameters of a mathematical model of a real-life problem
possess randomness with, most likely, some degree of correlation between them as well.
In studying the propagation  of the uncertainty through the model and its effect on the final solution
the great majority of the current techniques usually ignores  the correlations between random inputs.

The simplest approach to study the model uncertainty is to apply a Monte Carlo sampling \cite{MC},
where - once a probability density function (pdf) for the random inputs is assumed based on the
a priori knowledge about them - the mean and other characteristics
can be estimated from the output distribution, by sampling repeatedly from the assumed probability density
function and simulating the model for each sample.
Although with this approach it is possible to consider the correlations between variables,
it yields reasonable results only if the number of samples is quite large, requiring a great
computational effort. Moreover, the order of convergence is merely $\sqrt{N}$, where $N$ is the number of
samples. To decrease the computational effort, several modifications leading to new methods were introduced:
Latin hypercube sampling \cite{LHS}, the Quasi-Monte Carlo (QMC) method \cite{QMC},
the Markov Chain Monte Carlo method (MCMC) \cite{MCMC}, the Response Surface method (RSM) \cite{RSM}, etc..

Alternatively, deterministic methods to study parameter sensitivity have also been developed, e.g.,
perturbation methods\cite{KleiberHien1992}, local expansion-based methods like Taylor series, and so on.

In this paper we consider Polynomial Chaos Expansion (PCE), which experienced an increasing relevance
during the last years. The method is based on Wiener's \cite{PCE} homogeneous chaos theory published in
1938. Cameron and Martin \cite{CandM} proved convergence for the
classical Wiener-Hermite PC expansions based on the Hermite polynomial functionals  in terms
of Gaussian random variables. In the sixties of the last century some successful applications appeared
\cite{classic1,classic2}. For non-Gaussian random variables, however, PCE was very slow,
leading to a decrease of interest in the method.
To solve this problem, the PCE method was extended

to polynomials of the Askey scheme in \cite{XK}.
It sparked again new interest for this method, but application tasks demanded further
adaptation of the method to general distributions.
For this reason, several PCE extensions have been developed in the last years \cite{aPCE,g2,GSPC2008,aPCE2007,G-S}.
A theoretical framework for arbitrary multivariate probability measures was laid out in \cite{g1}.

These polynomial chaos expansion methods have been shown to be effective in a
large amount of applications in different fields as
can be seen in the literature:
air-water flows in soil \cite{water}, chemical reactions \cite{react},
fluid dynamics \cite{fluid,FD}, stability and control \cite{stcon},

etc.. New and
promising fields could be (bio)chemistry and artificial intelligence, e.g., uncertain case-based reasoning
\cite{CBRGolosnoy2008,CBRMaria2009}.
However, the method still has its limitations: {\em (i)} for large numbers of random variables, PCE
becomes - unacceptably - computationally expensive, and {\em (ii)} the random variables have to be
statistically independent.
Therefore, usually the uncertain inputs are assumed to be independent or fully correlated.
To apply PCE when the inputs are linearly correlated, some methods have been
proposed: linear transformations \cite{Nataf,Rosenblatt}, K-L expansion \cite{KL} or proper
orthogonal decomposition \cite{OrD}. All of these fix the problem by applying transformations to remove the
correlations, which increases the complexity of the problem and degrades the convergence rate
because of increased nonlinearity \cite{EldredEtAl2008}.

\medskip
In this paper, we propose a - more fundamental - solution for the latter of those limitations, since
correlations can have strong dynamical effects on the final solution as we will also show.
We present a method to construct an orthogonal polynomial basis for any general multivariate
distribution, including those with correlated random variables.

The outline of this paper is as follows: In Section 2 the extension is presented of the
polynomial chaos expansion to general multivariate distributions and the
derivation of the statistics of
the Quantities of Interest (QoI) - like mean, variance, and Sobol' indices - is discussed.
In Section 3 we show for an ODE dependent on a bivariate normal distribution the effect of the strength of
the correlation - from uncorrelated to fully correlated - on the solution. We also show that in all these cases
the convergence rate for increasing expansion order is the same, i.e.,
the proposed method shows the same favourable convergence rate as the
original PC method for an uncorrelated multivariate normal distribution.
We then apply the  method on a simple but realistic example from
biochemistry, viz. the canonical enzymatic reaction \cite{BrownMM02}, that has been used in a previous
paper \cite{FEBS} to display
the effect of noisy data on the reliability of the estimated parameters. We show that {\em true} propagation
of the uncertainty in the parameters - i.e., including the correlation - effects the uncertainty in the QoI,
the concentration of the product, significantly.
Section 4 contains a discussion and concluding remarks.

\section{PCE for multivariate arbitrarily distributed input variables}
In this section we extend the polynomial chaos expansion for arbitrarily distributed
{\em independent} random variables to the
{\em multivariate} arbitrarily distributed PC expansion.

\subsection{Polynomial chaos expansion}
Consider the following stochastic equation in the probability space
$(\Omega,\mathcal{A},P)$ where
$\Omega$ is the event space, $\mathcal{A}\subseteq 2^{\Omega}$ its $\sigma$-algebra and $P$ its probability measure
\begin{equation}
\mathcal{L}(\bsx,\bsxi(\bsom); \bsu )=\bff(\bsx,\bsxi(\bsom)),\;\;\bsx\in\text{X and }\bsom\in\Omega,
\label{equ:esto}
\end{equation}
where $\bsu(\bsx,\bsxi(\bsom))$ is the stochastic solution vector,
$\mathcal{L}$ is a
differential operator, and $\bff$ a
source function;
$\bsx$ is the vector of {\em deterministic} input variables describing, e.g., time or space, and
$\bsxi(\bsom)$ is the $n$-dimensional vector of {\em random} input variables, with joint probability density
function $\rho(\bsxi(\bsom))$; $\bsxi$ typically contains the uncertainties in model parameters
or initial and boundary conditions.
In order to make the notation less cumbersome, we denote the realization
of a random vector $\bsxi(\bsom)$, for $\bsom \in \Omega$, by $\bsxi \in \Xi$, with $\Xi$ the support of the pdf.

In the polynomial chaos method the random input $\bsxi$ and the solution $\bsu$ are expanded into a
series of polynomials
\begin{eqnarray}
\bsxi&=&\sum_{i=0}^\infty \bsxi_i \Phi_i(\bsxi)\nonumber\\
\bsu(\bsx,\bsxi) &=&
\sum_{i=0}^\infty \bsu_i(\bsx)\Phi_i(\bsxi),
\label{equ:expansion}
\end{eqnarray}
separating the {\em deterministic} and the {\em random} variables. The
$\Phi_i(\bsxi)$ are $n$-dimensional polynomials that are mutually orthogonal
with respect to the probability density function $\rho(\bsxi)$
\begin{equation}
\left<\Phi_i(\bsxi),\Phi_j(\bsxi)\right> = \int_\Xi \Phi_i(\bsxi) \Phi_j(\bsxi) \rho(\bsxi)\text{d}\bsxi =
||\Phi_i||^2 \delta_{ij},
\label{equ:ortrelation}
\end{equation}
with $\delta_{ij}$ the Kronecker delta
and $||\Phi_i||^2=\left<\Phi_i,\Phi_i\right>$. For independent random variables the $\Phi_i$'s are
tensor products of one-dimensional polynomials, $\Phi_i(\bsxi)=\prod_{j=1}^n \tilde\Phi_{j_i}(\xi_j)$.

To determine the polynomial chaos expansion coefficients $\bsxi_i$ and $\bsu_i(\bsx)$
there are two main-stream methods, viz. {\em Spectral Projection} and {\em Galerkin}. Both project onto
the polynomial space, but
whereas in the spectral projection approach the {\em expansion} (\ref{equ:expansion}) is projected,
in the Galerkin approach the {\em governing equation} (\ref{equ:esto}), with the expansion substituted,
is projected onto the polynomial space. Typically, for the random input variables the first approach
is chosen, resulting in
\begin{equation}
 \bsxi_j = \frac{\left<\bsxi,\Phi_j(\bsxi)\right>}{\left<\Phi_j(\bsxi),\Phi_j(\bsxi)\right>}=
 \frac{1}{||\Phi_j||^2}
 \int_\Xi \bsxi \, \Phi_j(\bsxi) \, \rho(\bsxi)\, \text{d}\bsxi,\qquad j=0,1, ... \;.
\label{eq:SPtheta}
\end{equation}

If the function is truely nonlinear aliasing is a threat for the accuracy \cite{Xiu2007},

therefore we will use the Galerkin method to obtain the solution.
The expansions (\ref{equ:expansion})
are substituted into the stochastic equation (\ref{equ:esto}) before being projected
\begin{eqnarray}
&&\left<\mathcal{L}\left(\bsx,\sum_{i=0}^\infty \bsxi_i \Phi_i(\bsxi); \sum_{i=0}^\infty \bsu_i(\bsx)\Phi_i(\bsxi)\right),\Phi_j(\bsxi)\right>=\nonumber\\
&&\qquad\qquad\qquad\quad=\left<\bff\left(\bsx,\sum_{i=0}^\infty \bsxi_i \Phi_i(\bsxi)\right),\Phi_j(\bsxi)\right>,\qquad j=0,1,... \;.
\label{equ:Gu}
\end{eqnarray}
In practice, the number of expansion terms is truncated to
\begin{equation}
 (N+1)=\frac{(n+p)!}{n!\,p!},
 \label{equ:dim}
\end{equation}
where $p$ is the highest order of the polynomials and the integrals are either computed exactly or
approximated using quadrature rules or Monte Carlo sampling. E.g., for a system of ODEs the differential operator in the left-hand side of (\ref{equ:esto})
reduces to $\dot{\bsu}(t,\bsxi)$ and its projection onto $\Phi_j$ to $\dot{\bsu}_j(t)\;||\Phi_j||^2$.

\subsection{Multivariate arbitrarily distributed input variables}
In \cite{G-S} a method has been proposed to extend the polynomial chaos method to an arbitrarily distributed
univariate distribution. The one-dimensional polynomials $\Phi_i(\xi)$ are constructed to be mutual
orthogonal with respect to an arbitrarily pdf $\rho(\xi)$ using the well-known and robust Gram-Schmidt
orthogonalization method (see e.g. \cite{GL3}). This method can be extended to
multivariate {\em independent} random variables, where the orthogonal multidimensional polynomials are
the product of the constructed one-dimensional orthogonal polynomials.
Here we use the same orthogonalization approach
to construct directly a multidimensional orthogonal polynomial basis, $\{\Phi_{j}(\bsxi)\}_{j=0}^{N}$,
for {\em correlated} multivariate random input variables.
Since the Gram-Schmidt method constructs from an arbitrary basis a basis that is orthogonal with respect to
a given innerproduct, the first step is the choice of a suitable set of linearly independent polynomials.
It should be noticed that any linearly independent set of polynomials can be used in this method,
but for simplicity we will use the set of monic polynomials $\{e_j(\bsxi)\}_{j=0}^{N}$
given by
\begin{equation}
e_j(\bsxi)=\prod_{l=1}^n\xi_l^{j_l},
\quad j=0, \ldots,N,\;j_l \in \{0,\ldots,p\},\text{ and }\sum_{l=1}^n j_l \leq p,
\label{equ:pol}
\end{equation}

where $n$ is the dimension of the random vector $\bsxi$, $p$ the highest polynomial
degree chosen, and $N+1$ the dimension of the polynomial basis given by Equation (\ref{equ:dim}).
E.g., if $n=p=2$ then $N+1=6$ and the set of linearly independent polynomials
$\{e_{j}(\xi_{1},\xi_{2})\}_{j=0}^5$ equals $\{1, \xi_{1},\xi_{2},\xi_{1}^{2},\xi_{2}^{2}, \xi_{1}\xi_{2} \}$.
Next, the orthogonal polynomial basis $\{\Phi_{j}(\mbox{\boldmath$\xi$})\}_{j=0}^{N}$
is constructed  from $\{e_{j}(\mbox{\boldmath$\xi$})\}_{j=0}^{N}$ using the Gram Schmidt algorithm
\begin{eqnarray}
\Phi_{0}(\bsxi)&=&1, \nonumber\\
\Phi_{j}(\bsxi)&=&e_{j}(\bsxi)-\sum_{k=0}^{j-1}c_{jk}\Phi_{k}(\bsxi)\quad\text{for }1\leq j \leq N,
\label{equ:recur-formu}
\end{eqnarray}
where the coefficients $c_{jk}$ are given by
\begin{equation}
c_{jk}=\frac{\left<e_{j}(\mbox{\boldmath$\xi$}),\Phi_{k}(\mbox{\boldmath$\xi$})\right>}{\left<\Phi_{k}(\mbox{\boldmath$\xi$}),\Phi_{k}(\mbox{\boldmath$\xi$})\right>}
\label{equ:c_coeff}
\end{equation}
and the innerproduct is taken with respect to the pdf $\rho(\bsxi)$. Note, that the basis is not
unique, it is dependent on the choice and on the ordering of the set of polynomials $\{e_j\}$.

Since each polynomial of the orthogonal basis
$\{\Phi_{j}(\mbox{\boldmath$\xi$})\}_{j=0}^{N}$
can be written as a sum of the monic
polynomials $e_j(\bsxi)$, the inner products in Equation (\ref{equ:c_coeff}) can be calculated as sum of
raw moments
\begin{equation}
\mu^{r_{1}\cdots r_{n}}_{\bsxi}=\int_\Xi\prod_{l=1}^n\xi_{l}^{r_l}\rho(\bsxi)d\bsxi,\quad
r_l \in \{0,\ldots,2p\},\text{ and }\sum_{l=1}^n r_l \leq 2p.
\end{equation}
With this procedure we immediately get a set of multidimensional orthogonal polynomials. Note, that if the
random variables are independent these polynomials are the same as the ones obtained by taking the tensor-products
of the one-dimensional orthogonal polynomials.
Once these multidimensional orthogonal polynomials have been computed the PCE method can be applied analogously
to the independent case, i.e., projection onto the polynomial space of either the expansion or
the equation with the expansion substituted.

\subsubsection{Statistics}
Part of the ease of use of PCE is the simplicity with which one obtains the most used statistics of the
QoI: mean, variance, and Sobol' indices can be directly expressed using the expansion coefficients \cite{Sudret2008}.
As a consequence of the orthogonality of the basis $\{\Phi_{j}(\mbox{\boldmath$\xi$})\}_{j=0}^{N}$
this favorable feature still holds
for the mean and variance , i.e., mean and variance of the solution vector $\bsu$ in
Equation (\ref{equ:expansion}) are given by
\begin{eqnarray}
 \bsmu_\bsxi(\bsx) &=& \int_\Xi \bsu(\bsx,\bsxi) \rho(\bsxi) d\bsxi = \bsu_0(\bsx),\text{ and}\\
 \bssi^2_\bsxi(\bsx) &=& \int_\Xi \left(\bsu(\bsx,\bsxi)-\bsmu_\bsxi(\bsx)\right)^2 \rho(\bsxi) d\bsxi =
 \sum_{i=1}^\infty \bsu^2_i(\bsx)||\Phi_i(\bsxi) ||^2
 .
\end{eqnarray}
The Sobol' indices \cite{Sobol1990,Li-globalsa_2010, Sudret2008, Owen2013} measure the influence of varying
{\em only the specified combination of} variables on the total variance.
They are based on the terms of the Sobol' decomposition of the QoI
\begin{equation}
 \bsu(\bsx,\bsxi) = \bsM_0(\bsx) + \sum_{\bsl \subseteq \{1,...,n\}} \bsM_\bsl(\bsx,\bsxi_\bsl),
 \label{eq:Soboldec}
\end{equation}
where $\bsM_0(\bsx) = \bsmu_\bsxi(\bsx)$, and the terms
$\bsM_\bsl$ are recursively given by
\begin{equation}
 \bsM_\bsl(\bsx,\bsxi_\bsl) = \E_{\bsxi_{- \bsl}}\left[\bsu|\bsxi_\bsl\right] -
 \left(\bsM_0(\bsx) + \sum_{\bsk \subset \bsl} \bsM_\bsk(\bsx,\bsxi_\bsk)\right),
\end{equation}
with  the marginal expectation given by
\begin{equation}
\E_{\bsxi_{-\bsl}}\left[\bsu(\bsx,\bsxi)|\bsxi_{\bsl}\right] =
 \int_{\Xi_{-\bsl}} \bsu(\bsx,\bsxi)
 \left(\int_{\Xi_\bsl} \rho(\bsxi) d\bsxi_{\bsl}\right)
 d\bsxi_{-\bsl},
\end{equation}
where $\bsxi_{-\bsl}$ indicates all elements of $\bsxi$ except $\bsxi_\bsl$. It is easy to see that this
expansion holds, since $\bsM_{1,...,n}(\bsx,\bsxi)$ equals $\bsu(\bsx,\bsxi)$ minus all previous terms.
Note, that this is an extension to general distributions of the original definition of Sobol' \cite{Sobol1990},
which was for $\mathcal{U}[0,1]^n$.

Using decomposition (\ref{eq:Soboldec}) the variance of the QoI is decomposed
\begin{eqnarray}\label{eq:vardec}
&& \Var\left[\bsu\right] =\Cov\left[\bsu,\bsu\right]= \Cov[\bsM_0 + \sum_{\bsl \subseteq \{1,...,n\}} \bsM_\bsl,\bsu] = \nonumber\\
&&  \sum_{\bsl \subseteq \{1,...,n\}} \Cov\left[\bsM_\bsl,\bsu\right]
= \sum_{\bsl \subseteq \{1,...,n\}} \left(\Var\left[\bsM_\bsl \right] + \Cov\left[\bsM_\bsl,\bsu - \bsM_\bsl \right]\right).
\end{eqnarray}

The Sobol' indices are now defined as the - normalized - covariance between the respective terms in the expansion
and the QoI
\begin{equation}
 \bsS_\bsl(\bsx) = \frac{\Cov\left[\bsM_\bsl(\bsx,\bsxi_\bsl),\bsu\right]}
 {\Var\left[\bsu\right]},\quad\text{for }\bsl \subseteq \{1,...,n\}.
\end{equation}

Following \cite{Li-globalsa_2010} (cf. also Eq. (\ref{eq:vardec})),
we decompose the Sobol' indices in a part that is dependent on the respective
variables and a part that is due to the correlation with all other variables
\begin{eqnarray}
 \bsS_\bsl(\bsx) &=&  \bsS_\bsl^u(\bsx) + \bsS_\bsl^c(\bsx), \text{ with}\label{eq:Si}\\
 \bsS_\bsl^u(\bsx) &=& \frac{\Var\left[\bsM_\bsl(\bsx,\bsxi_\bsl)\right]}
 {\Var\left[\bsu\right]}, \text{ and}\label{eq:Si_u}\\
 \bsS_\bsl^c(\bsx) &=& \frac{\Cov\left[\bsM_\bsl(\bsx,\bsxi_\bsl),\bsu-\bsM_\bsl(\bsx,\bsxi_\bsl)\right]}
 {\Var\left[\bsu\right]};\label{eq:Si_c}
\end{eqnarray}
where the index $\bsS_\bsl(\bsx)$ represents that part of
the variance of the QoI that is due to the variance of the set of variables $\bsxi_{\bsl}$,
$\bsS_\bsl^u(\bsx)$ represents the uncorrelated
share of it,
that is the contribution to the variance that comes from the set of variables $\bsxi_{\bsl}$
by themselves, and $\bsS_\bsl^c(\bsx)$ represents the correlated share,
the contribution to the variance of the QoI that comes from
the correlation of the set of the variables $\bsxi_{\bsl}$ with
the set of variables $\bsxi_{-_\bsl}$.
The {\em total influence} of a single variable, $\xi_{j}$, including that part of the variance due to
variable $\xi_{j}$ alone
and the fraction due to any combination of $\xi_{j}$ with the remaining variables is given by
\begin{equation}
 \bsS_{\bs T_{j}}(\bsx) = \sum_{\bsl\subseteq \{1,...,n\} \wedge j\in\bsl} \bsS_\bsl(\bsx),
\end{equation}
and analogous definitions for $\bsS_{\bs T_{j}}^{u}(\bsx)$, and $\bsS_{\bs T_{j}}^{c}(\bsx)$.

By definition $\bsS_\bsl^u(\bsx)$ is always positive, but $\bsS_\bsl^c(\bsx)$ can be either
positive or negative and therefore $\bsS_\bsl(\bsx)$ can also be either positive or negative.
This makes the interpretation and more specific the influence-based ranking of the variables less easy
than in the case of independent variables. For the variables with small $\bsS^c$ values the ``normal ranking''
can still be used, but when the Sobol' indices for a variable have a large $\bsS^c$ part, the interpretation
is still an open issue.

\medskip
For the computation of the Sobol' indices it is useful to notice that
a Sobol' index can be written as a
sum of Sobol' indices of the $\Phi$-terms with the appropriate coefficient,
since $\bsu(\bsx,\bsxi)$ is a linear combination of $\Phi$'s (see also \cite{Sudret2008}).
For independent random variables the polynomials $\Phi$ (and the distribution $\rho$) are products of
univariate contributions which implies that

\begin{eqnarray}
 \E_{\bsxi_{- \bsl}}\left[\Phi_j(\bsxi) | \bsxi_\bsl\right] &=&
 \int_{\Xi_{- \bsl}}  \prod_{k=1}^n \tilde\Phi_{j_k}(\xi_k)
 \left(\int_{\Xi_\bsl} \rho(\bsxi) d\bsxi_\bsl\right)
 d\bsxi_{- \bsl} =\nonumber\\
 &=& \prod_{k \in \bsl} \tilde\Phi_{j_k}(\xi_k)
 \int_\Xi \prod_{k \in -\bsl} \tilde\Phi_{j_k}(\xi_k) \rho(\bsxi)d\bsxi =\nonumber\\
 &=& \left\{ \begin{array}{ll}
              \prod_{k \in \bsl} \tilde\Phi_{j_k}(\xi_k) &
              \text{if }\prod_{k \in -\bsl} \tilde\Phi_{j_k}(\xi_k)=1(=\tilde\Phi_0),\\
              0 & \text{otherwise.}
             \end{array}
      \right.
 \label{eq:CondE_indep}
\end{eqnarray}
For a non-zero marginal expectation, the contribution to the numerator of the Sobol' index, consisting only of the
uncorrelated part (\ref{eq:Si_u}), is then given by
\begin{equation}
 \Var\left[\E_{\bsxi_{-_\bsl}}\left[\Phi_j(\bsxi)|\bsxi_\bsl\right]\right] =
 \int_\Xi  \left(\prod_{k \in \bsl} \tilde\Phi_{j_k}(\xi_k)\right)^2 \rho(\bsxi)d\bsxi
 =||\prod_{k \in \bsl} \tilde\Phi_{j_k}||^2.
\end{equation}
So, if the polynomial
\begin{equation*}
 \Phi_j(\bsxi) = \prod_{k \in \bsl} \tilde\Phi_{j_k}(\xi_k) \prod_{k \in -\bsl} \tilde\Phi_0(\xi_k),
\end{equation*}
the contribution to the Sobol' index is given by the expansion coefficient $u_j(\bsx)$
and otherwise it is zero.

For correlated random variables the orthogonality property used in (\ref{eq:CondE_indep}) no longer holds
and a Sobol' index is no longer given by a simple combination of expansion coefficients.
Fortunately, the same trick of separating the variables can still be used to compute the marginal
expectation and the covariance thereof, albeit at a lower level:
each $\bsu(\bsx,\bsxi)$ is a linear combination of $\Phi$'s which themselves are linear combinations
of the monic polynomials $e_j(\bsxi)$, so each Sobol' index can be calculated as a linear combination of
already computed raw moments with respect to the full pdf
\begin{eqnarray}
 \E_{\bsxi_{- \bsl}}\left[e_j(\bsxi) | \bsxi_\bsl\right] &=&
 \int_{\Xi_{- \bsl}}  \prod_{k=1}^n \xi_k^{j_k}
 \left(\int_{\Xi_\bsl} \rho(\bsxi) d\bsxi_\bsl\right)
 d\bsxi_{- \bsl}\nonumber\\
 &=& \prod_{k \in \bsl} \xi_k^{j_k}
 \int_\Xi \prod_{k \in -\bsl} \xi_k^{j_k} \rho(\bsxi)d\bsxi\nonumber\\
 &=& \mu^{j_{-\bsl}}_{\bsxi} \prod_{k \in \bsl} \xi_k^{j_k}.
 \label{eq:CondE_corr}
\end{eqnarray}
Note, that in this short-hand notation
$\mu^{i_{\bsl}}_{\bsxi}\;\mu^{j_{\bsk}}_{\bsxi} \neq \mu^{i_{\bsl}+j_{\bsk}}_{\bsxi}$.\\
The $e_j$-contribution to the numerator of the Sobol index $\bsS_\bsl^u$ (\ref{eq:Si_u}) is then given by
\begin{eqnarray}
 \Var\left[\E_{\bsxi_{-_\bsl}}\left[e_j(\bsxi)|\bsxi_\bsl\right]\right] &=&
 (\mu^{j_{-\bsl}}_{\bsxi})^2 \left(
 \int_\Xi \prod_{k \in \bsl} \xi_k^{2 j_k} \rho(\bsxi)d\bsxi
 -\left(\int_\Xi \prod_{k \in \bsl} \xi_k^{j_k} \rho(\bsxi)d\bsxi\right)^2
 \right)\nonumber\\
 &=& (\mu^{j_{-\bsl}}_{\bsxi})^2 [\mu^{2 j_{\bsl}}_{\bsxi}-(\mu^{j_{\bsl}}_{\bsxi})^2].
\end{eqnarray}
Analogously, the numerator of the Sobol' index $\bsS_\bsl^c$ (\ref{eq:Si_c}) can be written as a combination of
moments, since
\begin{eqnarray}
 \lefteqn{\Cov\left[\E_{\bsxi_{-_\bsl}}\left[e_j(\bsxi)|\bsxi_\bsl\right],e_i(\bsxi)\right]}\nonumber\\
 \qquad&=&
 \mu^{j_{-\bsl}}_{\bsxi}
 \int_{\Xi} \prod_{k \in \bsl} \xi_k^{j_k} \prod_{k=1}^n \xi_k^{i_k}
 \rho(\bsxi) d\bsxi - (\mu^{j_{-\bsl}}_{\bsxi} \mu^{j_{\bsl}}_{\bsxi}) \mu^{i_{1...n}}_{\bsxi}
 \nonumber\\
 &=& \mu^{j_{-\bsl}}_{\bsxi} \mu^{j_{\bsl}+i_{1...n}}_{\bsxi}-(\mu^{j_{-\bsl}}_{\bsxi} \mu^{j_{\bsl}}_{\bsxi}) \mu^{i_{1...n}}_{\bsxi}.
\end{eqnarray}
All these moments have already been computed in the construction of the polynomial basis
$\{\Phi_{j}(\mbox{\boldmath$\xi$})\}_{j=0}^{N}$, which
implies that also the Sobol' indices can be obtained without additional cost, although some bookkeeping is
required.

\bigskip
{\bf Remark.} Note, that the method is applicable to any type of continuous or discrete input probability distribution including
an experimentally obtained one. The only requirement is that the innerproducts of monic polynomials can
be calculated, which are moreover only dependent on the distribution and not on the problem at hand.

\section{Examples}
We illustrate the proposed method with two examples to demonstrate the
significant influence of correlations between stochastic inputs on the distribution of the QoI.
The first example is a scalar ODE in two-dimensional random space; for this example we study also its
numerical behavior.
The second example is a set of four ODEs in three-dimensional random space, where
we show the propagation of the resulting uncertainty in the parameters into the uncertainty in the QoI,
the concentration of the product.

\subsection{Decay equation}
Consider the scalar ODE
\begin{equation}
y'(t; \alpha,\beta) = -\alpha(y(t; \alpha,\beta)-\beta), \qquad 0 \leq t\leq 1, \qquad y(0) = 0,
\label{equ:Example1}
\end{equation}
with two random jointly distributed input variables $(\alpha, \beta) \sim \mathcal{N}(\bsmu,\bsSi)$,
where the mean $\bsmu$ and the covariance matrix $\bsSi$ are given by
\begin{equation}
 \bsmu=\left(\begin{array}{c}1\\1\end{array}\right) \text{ and }
 \bsSi=\left(\begin{array}{c c} \sigma_\alpha^2 & \varrho\sigma_\alpha\sigma_\beta\\
                                \varrho\sigma_\alpha\sigma_\beta & \sigma_\beta^2
             \end{array}\right) \text{, with }
 \sigma_\alpha=\sigma_\beta=0.25.
\end{equation}
The correlation coefficient $\varrho$ varies from uncorrelated to fully correlated, i.e.,
$\varrho = \{0,\pm 0.5,\pm 0.9,\pm 1\}$.
Note, that the random input to the problem is univariate when the correlation coefficient $\varrho=\pm1$
\begin{equation}
y'(t; \alpha)=\left\{
\begin{array}{ll}
-\alpha(y(t; \alpha)-\alpha),&\varrho=+1\\
-\alpha(y(t; \alpha)+\alpha-2),&\varrho=-1
\end{array}
\right\}\quad
\begin{array}{l}
0 \leq t\leq 1,\; y(0)=0,\\
\alpha \sim \mathcal{N}(1,0.0625).
\end{array}

\label{equ:Example1DPN}
\end{equation}

We want to show the effects of correlation and the convergence rate of the method for increasing expansion
order uncontaminated with errors in the computation of the polynomials or the projections.
Since the analytic expression of the  pdf of the correlated Gaussian distribution is known,
all moments, and the integrals thereof - necessary to compute the polynomial basis, the projection of the
right hand-side, and the statistics - are computed exactly with the moment-generating function using the Symbolic toolbox of Matlab \cite{Matlab}.
The projection of the truncated Eq. (\ref{equ:Gu}) onto the polynomial basis results in a system of ODEs
of the size of the expansion order. This system is solved with the accurate Matlab ODE solver ODE45 with
AbsTol 1e-6 (default) and RelTol changed to 1e-6.

\begin{figure}[ht!]
   \centering
	 \subfloat[order 0, $\varrho=-0.5$]{
        \label{fig:Neg0}
        \includegraphics[width=0.31\textwidth]{Negative_0}}
   \hspace{0.001\linewidth}

   \subfloat[order 1, $\varrho=-0.5$]{
        \label{fig:Neg12}
        \includegraphics[width=0.31\textwidth]{Negative_1_1}
         \includegraphics[width=0.31\textwidth]{Negative_1_2}}\\

  	 \subfloat[order 0, $\varrho=0$]{
        \label{fig:U0}
        \includegraphics[width=0.31\textwidth]{Uncorrelated_0}}
   \hspace{0.001\linewidth}

   \subfloat[order 1, $\varrho=0$]{
        \label{fig:U12}
        \includegraphics[width=0.31\textwidth]{Uncorrelated_1_1}
         \includegraphics[width=0.31\textwidth]{Uncorrelated_1_2}}\\

  \subfloat[order 0, $\varrho=0.5$]{
        \label{fig:Pos0}
        \includegraphics[width=0.31\textwidth]{Positive_0}}
   \hspace{0.001\linewidth}

   \subfloat[order 1, $\varrho=0.5$]{
        \label{fig:Pos12}
        \includegraphics[width=0.31\textwidth]{Positive_1_1}
         \includegraphics[width=0.31\textwidth]{Positive_1_2}}\\

   \caption{Polynomial basis. Order 0 and 1, for $\varrho = \{0,\pm 0.5\}$}
  \label{fig:Pol}
\end{figure}

\begin{figure}[ht!]
   \centering
	 \subfloat{
        \label{fig:Neg23}
        \includegraphics[width=0.31\textwidth]{Negative_2_1}}
   \hspace{0.001\linewidth}

        \subfloat[order 2, $\varrho=-0.5$]{
        \label{fig:Neg24}
        \includegraphics[width=0.31\textwidth]{Negative_2_2}}
   \hspace{0.001\linewidth}
   \subfloat{
        \label{fig:Neg25}
        \includegraphics[width=0.31\textwidth]{Negative_2_3}}\\
        	 \subfloat{
        \label{fig:U23}
        \includegraphics[width=0.31\textwidth]{Uncorrelated_2_1}}
   \hspace{0.001\linewidth}

        \subfloat[order 2, $\varrho=0$]{
        \label{fig:U24}
        \includegraphics[width=0.31\textwidth]{Uncorrelated_2_2}}
   \hspace{0.001\linewidth}
   \subfloat{
        \label{fig:U25}
        \includegraphics[width=0.31\textwidth]{Uncorrelated_2_3}}\\
        	 \subfloat{
        \label{fig:Pos23}
        \includegraphics[width=0.31\textwidth]{Positive_2_1}}
   \hspace{0.001\linewidth}

        \subfloat[order 2, $\varrho=+0.5$]{
        \label{fig:Pos24}
        \includegraphics[width=0.31\textwidth]{Positive_2_2}}
   \hspace{0.001\linewidth}
   \subfloat{
        \label{fig:Pos25}
        \includegraphics[width=0.31\textwidth]{Positive_2_3}}\\
   \caption{Polynomial basis. Order 2, for $\varrho = \{0,\pm 0.5\}$}
  \label{fig:Pol2}
\end{figure}

\begin{figure}[htbp!]
   \centering
   \subfloat[Mean]{
        \includegraphics[width=6.5cm,height=6.5cm]{MeansDecay.pdf}}
   \hspace{0.01\linewidth}
   \subfloat[Standard deviation]{
        \includegraphics[width=6.5cm,height=6.5cm]{StdsDecay.pdf}}
   \hspace{0.01\linewidth}
	\caption{Effect of correlation on mean and standard deviation of decay equation (\ref{equ:Example1})}
   \label{fig:Ex1-MSd}
\end{figure}

Figures \ref{fig:Pol}-\ref{fig:Pol2} show all polynomials of order 0, 1 and 2 for correlation coefficients $\varrho = \{0,\pm 0.5\}$.
The plots show clearly the non-uniqueness of the orthogonal basis, viz., the influence of the ordering  of the monic polynomials (\ref{equ:pol}): the first polynomials of a new order, $\Phi_{\{0,1,3\}}$,  are independent of the correlation coefficient and only determined by the first monic polynomial $e_{\{0,1,3\}}$ of that order, the next ones of the same order, $\Phi_{\{2,4,5\}}$, are different and dependent on the correlation.
In Figure \ref{fig:Ex1-MSd} we show the stochastic mean and standard deviation of the solution of
Eq. (\ref{equ:Example1}) for the seven values of the correlation coefficient
 using a PCE expansion  of 45 terms, which corresponds to a polynomial order of 8
and the polynomial basis calculated with the proposed method. The results are plot accurate,  in the cases of no or full correlation with the ones obtained with the Hermite polynomial basis and same order, and
with Monte Carlo simulations for the correlated distributions.

The importance of propagating the {\em true} input distribution is especially seen
in the standard deviation, where neither uncorrelated nor fully correlated give a reasonable approximation
for a distribution with correlation $\pm 0.5$. This is also reflected in the Sobol' indices.

Figure \ref{fig:Sobol_Ex1} shows for all correlation coefficients
the evolution of the Sobol' indices (\ref{eq:Si}) over time and the decomposition into the uncorrelated
(\ref{eq:Si_u}) and the correlated (\ref{eq:Si_c}) part. These plots and Table~\ref{tab:Sobol_tf_Ex1} show how nicely the
contributions of the variables itself and of the correlated variables follow the strength and sign of the
correlation coefficient $\varrho$.\\
{\bf Remark.} The fully correlated and uncorrelated results do not necessarily bound the correlated ones.
E.g., for a delta function response which is located just off the diagonal, full correlation results in a
zero variance, and no correlation gives a lower variance than correlation.

\begin{figure}[htbp!]
\begin{tabular}{cccc}
\includegraphics[scale=0.21]{SNeg09.pdf}
\includegraphics[scale=0.21]{SNeg09U.pdf}
\includegraphics[scale=0.21]{SNeg09C.pdf}\parbox{2cm}{\vspace{-3cm} $\varrho=-0.9$}\\
\includegraphics[scale=0.21]{SNeg05.pdf}
\includegraphics[scale=0.21]{SNeg05U.pdf}
\includegraphics[scale=0.21]{SNeg05C.pdf}\parbox{2cm}{\vspace{-3cm} $\varrho=-0.5$}\\
\includegraphics[scale=0.21]{SU.pdf}
\includegraphics[scale=0.21]{SUU.pdf}
\includegraphics[scale=0.21]{SUC.pdf}\parbox{2cm}{\vspace{-3cm} $\varrho=0$}\\
\includegraphics[scale=0.21]{SPos05.pdf}
\includegraphics[scale=0.21]{SPos05U.pdf}
\includegraphics[scale=0.21]{SPos05C.pdf}\parbox{2cm}{\vspace{-3cm} $\varrho=0.5$}\\
\includegraphics[scale=0.21]{SPos09.pdf}
\includegraphics[scale=0.21]{SPos09U.pdf}
\includegraphics[scale=0.21]{SPos09C.pdf}\parbox{2cm}{\vspace{-3cm} $\varrho=0.9$}\\
\end{tabular}
\caption{Dynamic plots of the Sobol' indices of decay equation (\ref{equ:Example1}).
Left: Sobol' index, middle: contribution from the respective variables, right: contribution due to correlation with
other variables.}
\label{fig:Sobol_Ex1}
\end{figure}

\begin{table}[htbp!]
\centering
\begin{tabular}{|r|l|ccc|}
\hline\hline
 $\varrho$& & $S_i$ & $S_i^u$ & $S_i^c$ \\
\hline
    &$S_1$&              -0.337996430555719&   1.079010380744548 &   -1.417006811300311\\
-0.9&$S_2$&                1.169072254592321 &   2.794020096839133 &    -1.624947842246789\\
    &$S_{12}$&     0.168924175966794 &    0.134133434465047&    0.034790741501840\\
\hline
&Sum&                             1.000000000003396&   4.007163912048728&  -3.007163912045260\\
\hline
    &$S_1$&               0.125070219450471&    0.462450913298287&  -0.337380693847799\\
-0.5&$S_2$&                0.819760985060438&   1.197491238203163&   -0.377730253142723\\
    &$S_{12}$&    0.055168795489075&    0.037711983172792&    0.017456812316329\\
\hline
&Sum&     0.999999999999984&  1.697654134674242&   -0.697654134674193\\
\hline
    &$S_1$&             0.273826682991019 & 0.273826682991019&  0.000000000000004 \\
  0 &$S_2$&              0.709059149322051 &  0.709059149322055& -0.000000000000000\\
    &$S_{12}$&   0.017114167686936&   0.017114167686939&  0.000000000000015\\
\hline
&Sum&     1.000000000000006&   1.000000000000013&   0.000000000000019\\
\hline
    &$S_1$&             0.340423008768462&   0.196827595099748 &  0.143595413668727\\
 0.5&$S_2$&              0.660699595355148 &  0.509674242072905&   0.151025353282240 \\
    &$S_{12}$&  -0.001122604123609 &  0.016050907606319&  -0.017173511729914\\
\hline
&Sum&    1.000000000000001&  0.722552744778972&  0.311794278680881\\
\hline
    &$S_1$&              0.374323408284474 &   0.161822550653440 &   0.212500857631028\\
 0.9&$S_2$&                0.636757390356394 &   0.419027897208577&   0.217729493147815\\
    &$S_{12}$&    -0.011080798640891&  0.020103679953051&    -0.031184478593928\\
\hline
&Sum&    0.999999999999977&  0.600954127815068&   0.399045872184915\\
\hline\hline
\end{tabular}
\caption{Sobol' indices and the contributions from variables itself and from correlated variables in the
final timepoint $t_f=1$. The respective contributions follow nicely the strength and sign of the correlation
coefficient.}
\label{tab:Sobol_tf_Ex1}
\end{table}

To answer the question whether the method has the same favourable convergence behavior as the original expansion
for independent normal random variables (cf. \cite{CandM}),
we have performed an error convergence study using the following error measures
\begin{equation}
\varepsilon_{\mu}(p)=\left|\mu_{p}(t_{f})-\mu_{exact}(t_{f})\right|, \qquad
\varepsilon_{\sigma}(p)=\left|\sigma_{p}(t_{f})-\sigma_{exact}(t_{f})\right|,
\end{equation}
where $t_{f}=1$ is the final time point,
$\mu_{p}$  and  $\sigma_{p}$  correspond to the mean and the standard deviation obtained with
PC expansion of order $p$, and
$\mu_{exact}$  and $\sigma_{exact}$ are the reference solutions, in this case obtained with $p=8$.

\begin{figure}[ht!]
   \centering
   \subfloat[Mean]{
        \includegraphics[width=6.5cm,height=6.5cm]{errorM.pdf}}
   \hspace{0.01\linewidth}
   \subfloat[Standard Deviation]{
        \includegraphics[width=6.5cm,height=6.5cm]{errorS.pdf}}
   \hspace{0.01\linewidth}
	\caption{Convergence of the error in mean and standard deviation for the
	decay equation (\ref{equ:Example1})}
   \label{fig:Decay-Error-MS}
\end{figure}
As we can see in Figure \ref{fig:Decay-Error-MS} the convergence rate of the error for increasing expansion
order is for all correlation coefficients $\varrho$ the same, showing that  the convergence rate of
the original PCE method for independently normally distributed input variables is conserved.

\subsection{Enzymatic reaction}
Consider the canonical enzymatic reaction
\begin{eqnarray*}
\E + \S &\equil{1}{2}& \C\\
\C &\product{3}& \E +\P
\end{eqnarray*}
where the concentrations of the substrate \S,
the enzyme \E, the complex \C, and the product \P\ are the state variables and
$k_1$, $k_2$, and $k_3$ are the kinetic rate constants.
This example has been used in a previous paper \cite{FEBS} to show the effect of noisy data on the
reliability of the estimated parameters using linear regression. However, in most cases biologists
are not so much interested in accurate model parameters, but in an accurate solution of the QoI, in this
case the concentration of the product \P. We now can propagate the
uncertainty in the parameters, given by the Fisher information matrix, so including correlation
between the kinetic parameters, through the model to obtain the pdf for the QoI \c\P.
The mathematical model for this enzymatic reaction is given by the following system of ODEs:
\begin{eqnarray}
\frac{d\c\S(t)}{dt} &=& -k_1\c\E(t)\c\S(t)+k_2\c\C(t) \nonumber\\
\frac{d\c\C(t)}{dt} &=& k_1\c\E(t)\c\S(t)-k_2\c\C(t)-k_3\c\C(t) \nonumber\\
\frac{d\c\E(t)}{dt} &=& -k_1\c\E(t)\c\S(t)+k_2\c\C(t)+k_3\c\C(t) \\
\frac{d\c\P(t)}{dt} &=& k_3\c\C(t)\nonumber
\label{equ:FEBS}
\end{eqnarray}
Data for \c\C\ are available at regular intervals during the time interval [0,20] and
the initial concentrations are known
$\c\S(0)=1$, $\c\C(0)=0$, $\c\E(0)=1$, and $\c\P(0)=0$.
The parameters $k_1$, $k_2$, and $k_3$ have been obtained by linear regression using the noisy data for \c\C.
The optimal value of the parameters in \cite{FEBS} was  $\bsk=(0.683, 0.312, 0.212)$, with (independent)
confidence intervals $\Delta=[0.033, 0.028, 0.005]$.
To study the propagation of the uncertainty through the model, the parameters are assumed to be random and their pdf
is built based on the results in \cite{FEBS}. Thus, we assumed that the parameters follow a uniform distribution with mean
the optimal value of the parameters $\bsmu=(0.683, 0.312, 0.212)$, standard deviation $\bssi=(0.206, 0.175, 0.031)$
based on the amplitude of the confidence intervals (with multiplication factor limited to 6.25 to avoid
entering the non-physical negative parameter space),
and $C$ the correlation matrix also obtained from the
Fisher information matrix; as in the previous example, the fully correlated case, $C_f$,
and the uncorrelated case, $C_u$, are also studied
\begin{eqnarray*}
&& C = \begin{pmatrix}
1 & 0.9& -0.37 \\
0.9 & 1& -0.45 \\
-0.37 & -0.45& 1 \\
\end{pmatrix}\;
C_f = \begin{pmatrix}
1 & 1& -1\\
1 & 1& -1 \\
-1 & -1& 1 \\
\end{pmatrix}\;
C_u = \begin{pmatrix}
1 & 0& 0\\
0 & 1& 0 \\
0& 0& 1 \\
\end{pmatrix}\
\end{eqnarray*}

As the analytic expression for the pdf of the joint distribution of two or more correlated
uniform variables is unknown, we cannot use the moment-generating function to  calculate all the required moments.
Therefore, Monte Carlo integration will be used in this example.
To generate sampling points from a correlated multidimensional uniform  distribution we used the standard
approach of computing them from the correlated normal distribution (see, e.g., \cite{Ucorr}).

Again, the problem collapses to one-dimensional in the fully correlated case, correlation matrix  $C_f$.

The only parameter in this case follows the standard uniform distribution $\mathcal{U}(0,1)$.

\begin{figure}[ht!]

   \subfloat[\c{S}]{

        \includegraphics[width=0.33\linewidth]{meansS}}

   \subfloat[\c{E}]{

        \includegraphics[width=0.33\linewidth]{meansE}}

   \subfloat[\c{P}]{

        \includegraphics[width=0.33\linewidth]{meansP}}

  \caption{Mean of the stochastic state variables of Eq.~(\ref{equ:FEBS}).}
     \label{fig:Means}
\end{figure}

\begin{figure}[ht!]

   \subfloat[\c{S}]{

       \includegraphics[width=0.33\linewidth]{stdS}}

   \subfloat[\c{C} and \c{E}]{

        \includegraphics[width=0.33\linewidth]{stdC}}

   \subfloat[\c{P}]{

        \includegraphics[width=0.33\linewidth]{stdP}}
				\hspace{0.01\linewidth}

 \caption{Standard deviation of the stochastic state variables of Eq.~(\ref{equ:FEBS})
 for the uncorrelated case (grey), the fully correlated case (red), and the correlation matrix resulting from
 the experimental errors (blue).}
   \label{fig:Stds}
\end{figure}

For all three cases - correlated, uncorrelated, and fully correlated - we applied PCE
with Galerkin projection to Eq. (\ref{equ:FEBS}).
\begin{figure}[htbp!]
\begin{tabular}{ccc}
\centerline{\includegraphics[width=0.7\linewidth]{S_legendu}}\\
\hspace{-0.7cm}
\subfloat[\c{S}]{
\includegraphics[width=0.33\linewidth]{S_S.pdf}}
\subfloat[\c{C} and \c{E}]{
\includegraphics[width=0.33\linewidth]{S_C.pdf}}
\subfloat[\c{P}]{
\includegraphics[width=0.33\linewidth]{S_P.pdf}}
\end{tabular}
\caption{Dynamic plots of the Sobol' indices of the enzymatic reaction (\ref{equ:FEBS}) for the uncorrelated
case.}
\label{fig:sobolfebsU}
\end{figure}
\begin{figure}[htbp!]
\begin{tabular}{ccc}
\includegraphics[width=0.33\linewidth]{SC_S}
\includegraphics[width=0.33\linewidth]{SC_Su}
\includegraphics[width=0.33\linewidth]{SC_Sc}\\
\c{S}\\
\includegraphics[width=0.33\linewidth]{SC_C}
\includegraphics[width=0.33\linewidth]{SC_Cu}
\includegraphics[width=0.33\linewidth]{SC_Cc}\\
\c{C} and \c{E}\\
\includegraphics[width=0.33\linewidth]{SC_P}
\includegraphics[width=0.33\linewidth]{SC_Pu}
\includegraphics[width=0.33\linewidth]{SC_Pc}\\
\c{P}
\end{tabular}
\caption{Dynamic plots of the Sobol' indices of enzymatic reaction (\ref{equ:FEBS}), for the correlation
matrix resulting from the experimental errors.
Left: Sobol' index $S$, middle: contribution from the respective variables, $S^u$,
right: contribution due to correlation with other variables, $S^c$.}
\label{fig:sobolfebsC}
\end{figure}
Figures \ref{fig:Means}  and \ref{fig:Stds} show the obtained mean and standard deviation of all concentrations,
respectively. Note, that the concentration of the complex \C\ has an equal dynamic
behaviour - with opposite sign - as that of the enzyme \E. Again, it can be clearly seen from Fig.~\ref{fig:Stds}
that the amount of correlation between the random parameters has a significant impact on the variance of the
stochastic state variables. Without correlation the pre-equilibrium dynamics are seen only in the complex
\C, and as a consequence in the enzyme \E, whereas the dynamics of the variance in the product \P\ simply
follow those of the substrate \S. With correlation the pre-equilibrium dynamics show up also in the substrate \S.
This difference is displayed too in the evolution of the Sobol' indices over time, see Figs.~\ref{fig:sobolfebsU}
and \ref{fig:sobolfebsC} (left column). From the plots in Fig. \ref{fig:sobolfebsC} it is again obvious
that the interpretation of the Sobol' indices - and thus the ranking of the random variables - is less clear
in case of correlated random variables. First of all, the Sobol' indices are no longer positive, and secondly, the
contribution due to correlation can completely cancel the contribution from the variable itself, resulting
in a small Sobol' value, but fixing such a variable would have a large impact on the outcome.

\begin{table}[h!]
\centering
\begin{tabular}{|c|ccc|ccc|}
\hline
\c{S}& $S_i$& $S_i^\mathrm{u}$& $S_i^\mathrm{c}$& $S_i$&  $S_i^\mathrm{u}$& $S_i^\mathrm{c}$ \\
\hline
$S_1$&       0.35&  0.36& -0.01& -0.52&  1.57& -2.09\\
$S_2$&       0.45&  0.46& -0.01&  0.93&  1.85& -0.92\\
$S_3$&       0.08&  0.08& -0.00&  0.42&  0.28&  0.13\\
$S_{23}$&    0.02&  0.02& -0.00&  0.10&  0.10&  0.00\\
$S_{13}$&    0.00&  0.00& -0.00&  0.05&  0.04&  0.00\\
$S_{12}$&    0.09&  0.09& -0.00&  0.15&  0.39& -0.24\\
$S_{123}$&   0.00&  0.00&  0.00& -0.13&  0.03& -0.16\\
\hline
Sum&         0.99&  1.01& -0.02&  1.00&  4.26& -3.28\\
\hline
\hline
\c{C}\ \& \c{E}& $S_i$& $S_i^\mathrm{u}$& $S_i^\mathrm{c}$& $S_i$&  $S_i^\mathrm{u}$& $S_i^\mathrm{c}$ \\
\hline
$S_1$&       0.17&  0.17& -0.00& -0.14&  0.19& -0.33\\
$S_2$&       0.24&  0.24& -0.00&  0.28&  0.31& -0.03\\
$S_3$&       0.58&  0.58& -0.00&  0.87&  0.80&  0.06\\
$S_{23}$&    0.00&  0.00& -0.00& -0.08&  0.01& -0.08\\
$S_{13}$&    0.00&  0.00&  0.00& -0.06&  0.00& -0.06\\
$S_{12}$&    0.00&  0.00&  0.00& -0.01&  0.00& -0.01\\
$S_{123}$&   0.00&  0.00&  0.00&  0.12&  0.02&  0.10\\
\hline
Sum&         0.99&  0.99&  0.00&  0.98&  1.33& -0.35\\
\hline
\hline
\c{P}& $S_i$& $S_i^\mathrm{u}$& $S_i^\mathrm{c}$& $S_i$&  $S_i^\mathrm{u}$& $S_i^\mathrm{c}$ \\
\hline
$S_1$&       0.31&  0.31& -0.00& -0.33&  0.75& -1.08\\
$S_2$&       0.41&  0.41& -0.00&  0.61&  0.95& -0.34\\
$S_3$&       0.23&  0.23& -0.00&  0.66&  0.54&  0.12 \\
$S_{23}$&    0.01&  0.01& -0.00&  0.01&  0.03& -0.02\\
$S_{13}$&    0.00&  0.00&  0.00& -0.00&  0.01& -0.01\\
$S_{12}$&    0.04&  0.04& -0.00&  0.05&  0.09& -0.04 \\
$S_{123}$&  -0.00&  0.00& -0.00& -0.00&  0.00& -0.00\\
\hline
Sum&         1.00&  1.00& -0.00&  1.00&  2.37& -1.37\\
\hline
\end{tabular}
\caption{Sobol' indices in the final time point $t_f=20$. The left part of the table corresponds
to the uncorrelated case, correlation matrix $C_{u}$;
the right part to the correlation matrix $C$ obtained from the FIM.
The columns $S_{i}^{u}$ show  the contribution from the variables itself and
the columns $S_{i}^{c}$ the contribution from correlated variables; with exact computations the $S_i^c$-column
in the left part of the table should be zero.}
\label{tb:Sobolfebs}
\end{table}

\begin{table}[h!]
\centering
\begin{tabular}{|c|ccc|ccc|}
\hline
\c{S}& $S_{T_{i}}$& $S_{T_{i}}^\mathrm{u}$& $S_{T_{i}}^\mathrm{c}$& $S_{T_{i}}$& $S_{T_{i}}^\mathrm{u}$& $S_{T_{i}}^\mathrm{c}$ \\
\hline
$S_{T_1}$&  0.44&   0.45& -0.01& -0.45&  2.03& -2.49\\
$S_{T_2}$&  0.56&   0.57& -0.01&  1.05&  2.37& -1.32\\
$S_{T_3}$&  0.10&   0.10&  0.00&  0.44&  0.45& -0.03\\
\hline
\hline
\c{C}\ \& \c{E}& $S_{T_{i}}$& $S_{T_{i}}^\mathrm{u}$& $S_{T_{i}}^\mathrm{c}$& $S_{T_{i}}$& $S_{T_{i}}^\mathrm{u}$& $S_{T_{i}}^\mathrm{c}$ \\
\hline
$S_{T_1}$&  0.17&   0.17& -0.00& -0.09&  0.21& -0.30\\
$S_{T_2}$&  0.24&   0.24& -0.00&  0.31&  0.34& -0.02\\
$S_{T_3}$&  0.58&   0.58&  0.00&  0.85&  0.83&  0.02\\
\hline
\hline
\c{P}& $S_{T_{i}}$& $S_{T_{i}}^\mathrm{u}$& $S_{T_{i}}^\mathrm{c}$& $S_{T_{i}}$& $S_{T_{i}}^\mathrm{u}$& $S_{T_{i}}^\mathrm{c}$ \\
\hline
$S_{T_1}$&  0.35&   0.35& -0.00& -0.28&  0.85& -1.13\\
$S_{T_2}$&  0.46&   0.46& -0.00&  0.67&  1.07& -0.38\\
$S_{T_3}$&  0.24&   0.24&  0.00&  0.67&  0.58&  0.09 \\
\hline
\hline
\end{tabular}
\caption{Total Sobol' indices in the final time point $t_f=20$.
The left part of the table corresponds to the uncorrelated case, correlation matrix $C_{u}$;
the right part to the correlation matrix $C$ obtained from the FIM.}
\label{tb:SobolfebsT}
\end{table}

Ranking variables is often done based on the Sobol' indices. Depending on the time-frame of interest, one can
use as measure, e.g., the integral over time of the Sobol' indices, or the Sobol' indices in a specific point.
As an example of the latter we study at the final time-point, $t_f=20$, which part of the variance of the
concentrations is due to the variance of respective input variables.
Tables~\ref{tb:Sobolfebs} and \ref{tb:SobolfebsT} give the Sobol' indices
and the total Sobol' indices, respectively, for the random variables $k_1$, $k_2$, and $k_3$. Note that the
sum of the Sobol' indices equals one and for the uncorrelated case the $S^c$ values should be zero, so the
deviation of these values gives the error in the approximation, in this case mostly due to the approximation of
the high-order moments needed for the Sobol' indices. From Table~\ref{tb:SobolfebsT} one can deduce that, in
the uncorrelated case, decomposing
the variance of the concentration of the substrate \S\  would lead to the conclusion that
the random variable $k_3$ can be fixed to the nominal deterministic value, since its influence on the
variance of \S\ is negligible. In the correlated case, however, this is no longer true. One can also see
that especially for $k_1$ and $k_2$ the interpretation of the Sobol' indices is not trivial, e.g., for \C\ the
influence of $k_1$ is very small but this is due to cancellation of the $S^u$ and $S^c$ contribution.
As stated before, the interpretation of the Sobol' indices for correlated random variables is largely an open
issue, but it is clear that in general less random input variables can be fixed to their nominal deterministic
values.

\section{Discussion and concluding remarks}
In this paper we extended the PCE method to general multivariate distributions, including correlations
between the random input variables, by constructing an orthogonal polynomial basis for the probability
space. We showed that the usual statistics like mean, variance and Sobol' indices can be obtained with
no extra cost, as is the case for the current PCE implementations. This method makes it possible to compute
for a set of random input variables with any multivariate distribution, including an experimentally
determined one, the stochastic distribution of a Quantity of Interest. An application we studied is
the true propagation of experimental errors through the parameter-fitting process onto the state variables,
so that a QoI based on these state variables is described as a distribution and, e.g.,
optimal experimental design can be applied to reduce the variance in such a QoI.
In this paper we used either exact integration or Monte Carlo to compute the moments and the projection onto
the polynomial basis, but for the latter Gauss quadrature is under development.

\medskip
There are a few remarks to make with respect to the method
\begin{itemize}
 \item As for other PCE methods, the polynomial base is only dependent on the multivariate distribution
 of the random input variables. Once computed, all problems with random input variables described by
 this pdf can be solved using the same base.
 \item To obtain the polynomial basis raw moments are needed which have to be computed accurately enough
 not to pollute the results. This is not a trivial task for high dimensional problems or a high polynomial order. E.g., in our first example, the size of
 the moments range in order of magnitude from 1 for the first ten orders to $10^{15}$ for order 75 for
 an uncorrelated Normal distribution. If the correlation coefficient $\varrho$ equals $0.9$ the latter
 order of magnitude is even much higher: $10^{21}$ for moments of order 75.
 \item The orthogonal polynomial basis is not unique, it is dependent on the choice and ordering of
 the original set of linearly independent polynomials. It is an open question at the moment whether there
 is an optimal choice for this original set of polynomials, e.g., to reduce the magnitude of the moments
 required or to simplify computations.

 \item Whereas in the uncorrelated case interpretation of Sobol' indices and ranking of the influence
 of the variation of the random input variables on the variance in the QoI is more or less straighforward,
 this is no longer true for correlated random variables.
\end{itemize}

\section*{Acknowledgements}
MN and JB acknowledge support by European Commission's 7th Framework Program, project BioPreDyn grant number 289434.

\newpage

config: arxiv


TITLE: Does the UKF work? - or - how would a correct implementation of UKF look like?
QUESTION [0 upvotes]: I am a 5th year engineering student, and I am currently struggling to implement the Unscented Kalman Filter. I tried the regular version and a square root version and they both failed. Research showed many people failed to implement it, and only one succesful though simple implementation un the matlab central website.
This brings me to the question - does it even work? Is the UKF a reasonable solution to any random process, or maybe just to very simple cases? From its formulation, it seems to be able to handle very complicated systems, but I cant seem to be able to make it work, nor could i find too many evidence of succesful work by anyone. Currently, I am assuming my implementation is just wrong, even though I am working on it over a year (!). If this really is the case, I'd be greatfull if anyone could enlight me, as my projects deadline is this weekend and I am helpless.
The system i work on tries to estimate the orientation of an object, and contains 3 angles. The equation:$$$$
$[ \dot x_1] =\matrix [1, \space sin(x_1)tan(x_2), \space cos(x_1)tan(x_2)]*\omega $
$ [\dot x_2] =\space[0, \space cos(x_1), \space -sin(x_1)]*\omega$
$ [\dot x_3] =\space[0, \space sin(x_1)sec(x_2), \space cos(x_1)sec(x_2)]*\omega$
Represents the relation between the state and (given) angular velocity $\omega$ to the state derivation $\dot x$, so my state equation would be adding this times the sampling time dt to the state provided by last iteration.
The measurement function is given by: original picture of matrix for comparison
$$R_i^b=\left[\matrix{\cos\theta\cos\psi & \cos\psi\sin\theta\sin\phi - \cos\phi\sin\psi & \cos\psi\sin\theta\cos\phi + \sin\phi\sin\psi \\
\cos\theta\sin\psi & \sin\psi\sin\theta\sin\phi + \cos\phi\cos\psi & \sin\psi\sin\theta\cos\phi - \sin\phi\cos\psi \\
-\sin\theta        & \cos\theta\sin\phi                            & \cos\theta\cos\phi}\right]$$
$$
\hat y = \left[\matrix{R_i^b\\ R_i^b}\right]
\cdot\left[\begin{smallmatrix}0\\0\\ -g \\ M_h\\0\\M_z\end{smallmatrix}\right], \qquad\text{where $g$: gravity const; $M_h$, $M_z$: magnetic const} $$
here, note: $ x=[\phi, \space \theta, \space \psi ]^T$ - the 3 angles to estimate.
As mentioned, 2 versions of the UKF were tested, but none gave the right results. The original form of the filter kept failing with the error of covariance matrix not being positive definite, and after fixing it with various approaches, just returned wrong results, showing more like noise than anything. the square root form, given in http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.80.1421&rep=rep1&type=pdf , just kept failing at the last step of the iteration, where a downdate needs to be preformed, saying the matrix cannot be downdate due to - you guessed it - the matrix not being positive definite.
as long as the  question may be, I would appreciate any help you can give me, wether consists on your experiment with the filter or just with the math or matlab in general.
with respect, Joseph.

REPLY [1 votes]: I'm answering 4 months later but if you have been working for more than a year maybe it is still useful.
Yes, the UKF filter does work, if used correctly.
Getting better results than measurements is not just about having a right state matrix - which is still important. It is mainly about
 - correct error estimation - through measurement error (usually R matrix) and confidence in your model (usually named Q)
 - correct means, especially when talking about angles, since the mean of two angles is not just (a+b)/2 : (359° + 1°)/2 = 0°, not 180°
Example about the model correctness : is the acceleration is not just due to gravitation but also due to angular speed, the model is wrong. That must be indicated in Q.
Now, about implementation : When talking about orientation measurement, quaternions are more used that euler angles, because of euler singularities https://en.wikipedia.org/wiki/Gimbal_lock.
To get orientation from IMU and You will find very widely used library here http://x-io.co.uk/open-source-imu-and-ahrs-algorithms/.
Hope this helps.


TITLE: Schanuel conjecture and its set theoretic status
QUESTION [5 upvotes]: I have once been told that Schanuel's conjecture (which is a transcendence statement about tuples of complex numbers and their exponentials) is set theoretically absolute, which means that its truth (or falsity) is independent of the set theoretical model (or assumptions), because it is a $\Pi^1_2$ statement. Is this correct?

REPLY [7 votes]: This is an extremely sloppy reference to the Shoenfield absoluteness theorem.
This theorem says that given an arithmetical statement about the natural numbers which is not "too complicated" (first-order statements are not too complicated, but not just them), then the statement is true in every two transitive models of $\sf ZFC$ which have the same ordinals.
This doesn't mean that the truth value is independent of $\sf ZFC$. It just means that we cannot prove independence using the two "usual" tools of set theory:

Forcing, which does not change the ordinals of a model, nor destroy its transitivity; and
Inner models, which are by definition a subclass of the model which is transitive and contains all the ordinals.

There might be other ways to prove independence, e.g. compactness arguments of some sort. But things like compactness usually destroy nice properties of models of $\sf ZFC$, and the absoluteness theorem no longer applies. On the other hand, we are less interested in not-nice models of $\sf ZFC$ in general, which is why set theorists focus on the two methods above, rather than other methods.
So the conclusion is that if we can force the truth value of the Schanuel conjecture, then we actually proved it.


text:
\begin{document}

\maketitle

\begin{abstract}
We study $\mathrm{wsat}\big({G}(n, p), K_s\big)$, the weak $K_s$-saturation number of the Erd\H{o}s--R\'enyi random graph $G(n, p)$, where $K_s$ is the complete graph on $s$ vertices. Kor\'{a}ndi and Sudakov proved that the weak $K_s$-saturation number of $K_n$ is stable, i.e. it remains the same after removing edges with constant probability. We prove that there exists a threshold function for this stability property and give upper and lower bounds on the threshold. Having this, we simplify the proof of Kor\'andi and Sudakov and make their result more general. A general upper bound for  $\mathrm{wsat}\big({G}(n, p), K_s\big)$ is also provided.
\\[.3cm]
\noindent {\bf Keywords:}    Weak saturation,  Random graph, Stability. \\[1mm]
\noindent {\bf AMS Mathematics Subject Classification\,(2020):}    05C80, 05C35, 60K35. \\[6mm]
\end{abstract}

\section{Introduction}

Given a graph $F$, an {\it $F$-bootstrap percolation process} is a sequence of graphs $H_0\subset H_1\subset\cdots\subset H_m$ such that, for every $i\in\{1,\ldots,m\}$, $H_{i}$ is obtained from $H_{i-1}$ by adding an edge that belongs to a copy of $F$ in $H_i$. The $F$-bootstrap processes were introduced by Bollobas 50 years ago~\cite{bol} and can be seen as special cases of the `cellular automata' introduced by von Neumann~\cite{Neumann} after a suggestion of Ulam~\cite{Ulam}. The $F$-bootstrap percolation is also similar to the $r$-neighborhood bootstrap percolation model having applications in physics (see, e.g.,~\cite{Adler,Fontes,Morris}).

A fundamental question about bootstrap percolation is, how large is the smallest possible percolating set in the $F$-bootstrap process? Formally, given two graphs $G$ and $F$,  a spanning subgraph  $H$  of $G$  is  said to be a  {\it weakly $F$-saturated subgraph of $G$} if $H$ contains no subgraph  isomorphic to $F$ and there exists an $F$-bootstrap percolation process $H=H_0\subset H_1\subset\cdots\subset H_m=G$. The  minimum number of edges  in  a weakly $F$-saturated  subgraph of $G$ is called {\it weak $F$-saturation number} of $G$, denoted by $\mathrm{wsat}(G, F)$.

As usual, we denote by $G(n,p)$, the binomial random graph on vertex set $[n]:=\{1,\ldots,n\}$ with the edge probability $p$ (recall that this graph contains every edge with probability $p$ independently of all the others).
Kor\'andi and   Sudakov   initiated  the study of  weak saturation number of random graphs in \cite{kor}. They   proved that, for every fixed $p\in(0, 1)$ and integer  $s\geqslant 3$, $\mathrm{wsat}\big({G}(n, p), K_s\big)=\mathrm{wsat}\big(K_n, K_s\big)$ whp (this is the abbreviation for the common notion `with high probability' meaning with a probability approaching 1 as $n\to\infty$), which is in turn ${s-2 \choose 2}+(n-s+2)(s-2)$ by a classic result proved by Lov\'{a}sz \cite{L77}.
They also noticed that the same is true for $p\geqslant n^{-\varepsilon}$ for small
enough $\varepsilon> 0$ and asked about smaller $p$ and about possible {\it threshold probability}~\cite[Chapter 1.5]{RG} for the
property of having the weak $K_s$-saturation number exactly ${s-2 \choose 2}+(n-s+2)(s-2)$, which we denote by $\mathcal{A}_s$.
In this paper, we prove that this threshold exists and present lower and upper bounds on that.\\

Before proceeding with the results, let us list the notation we use in the paper.
For   a graph $G$, we denote  the vertex set and the edge set of   $G$  by $V(G)$ and $E(G)$, respectively.  For  $S\subset V(G)$, we denote the induced subgraph of $G$ on $S$  by $G[S]$. For a vertex $v\in V(G)$, set $N_G(v)=\{x\in V(G) \, | \, v \text{   is adjacent to } x\}$. Also for a set $U \subset V(G)$ let $N_G(U)=\cap_{u\in U} N_G(u)$, and $N_G[U]=U \cup N_G(U)$.

\section{A lower bound on $p$ and the existence of threshold}

In this section, we prove the existence of threshold and present a lower bound.

\begin{theorem}
	Let $s\geq 3$ be fixed,
	$$
	q_s(n)=n^{-\frac2{s+1}}(\ln n)^\frac{2}{(s-2)(s+1)},\quad
	c_s=\left[2\left(1-\frac{1}{s+1}\right)(s-2)!\right]^{\frac{2}{(s+1)(s-2)}}.
	$$
	Then the following hold:
	\begin{enumerate}
		\item[(i)] There exists a threshold probability for the property $\mathcal{A}_s$.
		\item[(ii)] If $p(n)\leq c_sq_s(n)$, then whp the property $\mathcal{A}_s$ does not hold in ${G}(n, p)$.
	\end{enumerate}
	\label{lower}
\end{theorem}

{\it Proof.} Below, we assume that $p>\frac{1}{n\ln n}$ since, otherwise, whp there is no $K_s$ in $G(n,p)$ and the number of edges is asymptotically smaller than $n$ (it follows from Markov's inequality, see, e.g., Chapter 3.1 of \cite{RG}).

Let $\varepsilon\in(0,1)$. Then, whp the number of edges of $G(n,p)$, hereinafter denoted by $e(G(n,p))$, belongs to $\left(\frac{n^2p}{2(1+\varepsilon)},\frac{n^2p}{2(1-\varepsilon)}\right)$ since it has binomial distribution with parameters ${n\choose 2}$ and $p$ (it follows, for example, from Chebyshev's inequality).

Let $X_s$ be the number of $s$-cliques in $G(n,p)$ and let $\varepsilon>0$. It is well known (see, e.g., Chapter 3 of \cite{RG}) that,
for $p\ll n^{-\frac{2}{s+1}}$, whp $X_s=o(n^2p)$, implying $\mathrm{wsat}\big({G}(n, p), K_s\big)=e({G}(n, p))(1+o(1))$ whp since $e({G}(n, p))-X_s\leq\mathrm{wsat}\big({G}(n, p), K_s\big)\leq e({G}(n, p))$. From this, we immediately get that whp $\mathrm{wsat}\big({G}(n, p), K_s\big)\neq\mathrm{wsat}\big(K_n, K_s\big)$ for $p\ll n^{-\frac{2}{s+1}}$ such that $p\notin I:=\left((1-\varepsilon)\frac{2s-4}{n},(1+\varepsilon)\frac{2s-4}{n}\right)$. For $p\in I$ and $s\geq 4$, whp $X_s=0$ and, therefore, the weak saturation number is exactly $e({G}(n, p))$ that is not concentrated on a single value. Finally, if $p\in\left((1-\varepsilon)\frac{2}{n},(1+\varepsilon)\frac{2}{n}\right)$, then whp the weak $K_3$-saturation number equals the difference between the number of edges and the number of triangles (whp all the triangles in $G(n,p)$ are disjoint since whp there are no subgraphs with at most 5 vertices and at least two cycles; see Theorem 3.4 of \cite{RG}). The latter random variable is also not concentrated in a unit set since the number of triangles is bounded from above by asymptotically Poisson random variable (for a constant $L$, the property of having at most $L$ triangles is {\it decreasing}~\cite[Chapter 1.3]{RG} while the asymptotical distribution of the number of  triangles in $G(n,2(1+\varepsilon)/n)$ is Poisson; see Lemma 1.10 and Theorem 3.19 of \cite{RG}).
Summing up, for $p\ll n^{-\frac{2}{s+1}}$, whp $G(n,p)$ does not have $\mathcal{A}_s$.

Now, let $\gamma>0$ be fixed and small enough, and $p\geq n^{-\frac{2}{s+1}-\gamma}$. To finish the proof, we need the following technical lemma about the property $\mathcal{B}_s$ saying that every edge belongs to $K_s$.

\begin{lemma}
	Let
	$$
	p(n)=c_sq_s(n)\left[1+\frac{2\ln\ln n}{(s-2)^2(s^2+s)\ln n}+\frac{w_n}{\ln n}\right].
	$$
	Then the following hold:
	\begin{enumerate}
		\item[(i)] If $w_n\to\infty$, then whp $G(n,p)$ has $\mathcal{B}_s$.
		\item[(ii)] If $w_n\to-\infty$, then whp $G(n,p)$ does not have $\mathcal{B}_s$.
	\end{enumerate}
	\label{edge_in_clique}
\end{lemma}

{\it Proof.} In~\cite{SP90}, it is proven that there exists $R$ such that, for $p>Rq_s(n)$, whp every pair of vertices have $s-2$ adjacent common neighbors. Therefore, we may assume that $w_n\leq R\ln n$.\\

Below, we, as usual, apply, so called, second moment method and the Janson inequality~\cite{Janson}. Let $N_s$ be the number of edges that do not belong to any $K_s$.

Let $u,v\in[n]$ be distinct, and $W:=\{w_1,\ldots,w_{s-2}\}\subset[n]\setminus\{u,v\}$. Consider an event $K[W]$ saying that each  $w_i$ is adjacent to both $u,v$ and every other $w_j$ in $G(n,p)$. Let $\mu(u,v)$ counts the number of sets $W$ as above such that $K[W]$ happens. Then,
\begin{equation}
{\sf E}N_s={n\choose 2}p{\sf P}(\mu(1,2)=0).
\label{expectation-mu}
\end{equation}
If $s=3$, then ${\sf P}(\mu(1,2)=0)=(1-p^2)^{n-2}$. For $s\geq 4$, compute
$$
\sum_W {\sf P}(K[W])={n-2\choose s-2}p^{\frac{(s+1)(s-2)}{2}}=:\lambda,
$$
$$
\sum_{W_1,W_2}{\sf P}(K[W_1]\cap K[W_2])=\sum_{\ell=1}^{s-3}{n-2\choose s-2}{s-2\choose\ell}{n-s\choose s-2-\ell}p^{(s+1)(s-2)-\frac{(\ell+3)\ell}{2}}=\Delta,
$$
where the first summation is over $W\in{[n]\setminus\{u,v\}\choose s-2}$, and the second one is over distinct $W_1,W_2\in{[n]\setminus\{u,v\}\choose s-2}$ having a non-empty intersection. Then, the Janson inequality implies that
$$
\mathrm{exp}\left[-\lambda\left(1+O\left(p^{\frac{(s+1)(s-2)}{2}}\right)\right)\right]\leq{\sf P}(\mu(1,2)=0)\leq
\mathrm{exp}\left[-\lambda+\Delta/2\right].
$$
For the considered values of $p$, $\Delta=o(\lambda\cdot n^{-1/(s+1)}\ln n)$. Therefore,
\begin{equation}
{\sf P}(\mu(1,2)=0)=\exp\left[-\lambda\left(1+o\left(n^{-1/(s+1)}\ln n\right)\right)\right]=e^{-\lambda}\left(1+o\left(n^{-1/(s+1)}\ln^2 n\right)\right).
\label{mu}
\end{equation}

If $w_n\to\infty$, then
$$
{\sf E}N_s\leq{n\choose 2}p e^{-\lambda+\Delta/2}\leq\mathrm{exp}[2\ln n+\ln p-\lambda+O(1)]\to 0.
$$
Part (i) follows.\\

If $w_n\to-\infty$, then
$$
{\sf E}N_s\geq{n\choose 2}p \exp\left[-\lambda\left(1+O(n^{-(s-2)}\ln n)\right)\right]\geq\exp[2\ln n+\ln p-\lambda+o(1)]\to\infty.
$$
To finish the proof of Part (ii), we need the second moment. Similarly, using the Janson inequality, we estimate ${\sf P}(\mu(1,2)=0,\mu(3,4)=0)$ and ${\sf P}(\mu(1,2)=0,\mu(2,3)=0)$.

For $W:=\{w_1,\ldots,w_{s-2}\}\subset[n]\setminus[4]$, consider an event $K'[W]$ saying that each  $w_i$ is either adjacent to both $1,2$, or adjacent to both $3,4$ and every other $w_j$ in $G(n,p)$. Let $\mu'$ counts the number of sets $W$ as above such that $K'[W]$ happens. Then,
\begin{equation}
{\sf P}(\mu(1,2)=0,\mu(3,4)=0)\leq{\sf P}(\mu'=0).
\label{mu-one}
\end{equation}
If $s=3$, then ${\sf P}(\mu'=0)=(1-p^2)^{2(n-4)}$. For $s\geq 4$, compute
$$
\sum_W {\sf P}(K'[W])={n-4\choose s-2}p^{\frac{(s+1)(s-2)}{2}}(2-p^{2(s-2)})=:\lambda'=2\lambda(1+O(1/n)),
$$
$$
\sum_{W_1,W_2}{\sf P}(K'[W_1]\cap K'[W_2])=
\sum_{\ell=1}^{s-3}{n-4\choose s-2}{s-2\choose\ell}{n-s-2\choose s-2-\ell}\times
$$
$$
p^{(s-2)(s-3)-\frac{(\ell-1)\ell}{2}}2p^{2\ell}\left(2p^{2(s-2)-2\ell}\right)^2(1+O(p))=:\Delta'=8\Delta(1+O(q_s(n))),
$$
where the first summation is over $W'\in{[n]\setminus[4]\choose s-2}$, and the second one is over distinct $W_1,W_2\in{[n]\setminus[4]\choose s-2}$ having a nonempty intersection. Then, the Janson inequality implies that
$$
{\sf P}(\mu'=0)\leq
\mathrm{exp}\left[-\lambda'+\Delta'/2\right]=\mathrm{exp}\left[-2\lambda(1+O(1/n))+4\Delta(1+O(q_s(n)))\right]=
$$
\begin{equation}
\exp\left[-2\lambda\left(1+o\left(n^{-1/(s+1)}\ln n\right)\right)\right]=e^{-2\lambda}\left(1+o\left(n^{-1/(s+1)}\ln^2 n\right)\right).
\label{mu-one-from-above}
\end{equation}

In the same way, let us estimate ${\sf P}(\mu(1,2)=0,\mu(2,3)=0)$. For $W:=\{w_1,\ldots,w_{s-2}\}\subset[n]\setminus[3]$, consider an event $K''[W]$ saying that each  $w_i$ is either adjacent to both $1,2$, or adjacent to both $2,3$ and every other $w_j$ in $G(n,p)$. Let $\mu''$ counts the number of sets $W$ as above such that $K''[W]$ happens. Then,
\begin{equation}
{\sf P}(\mu(1,2)=0,\mu(2,3)=0)\leq{\sf P}(\mu''=0).
\label{mu-two}
\end{equation}
If $s=3$, then ${\sf P}(\mu''=0)=(1-p+p(1-p)'')^{n-3}$. For $s\geq 4$, compute
$$
\sum_W {\sf P}(K''[W])={n-3\choose s-2}p^{\frac{(s+1)(s-2)}{2}}(2-p^{s-2})=:\lambda''=2\lambda(1+o(1/\sqrt{n})),
$$
$$
\sum_{W_1,W_2}{\sf P}(K''[W_1]\cap K''[W_2])=
\sum_{\ell=1}^{s-3}{n-3\choose s-2}{s-2\choose\ell}{n-s-1\choose s-2-\ell}\times
$$
$$
p^{(s-2)(s-3)-\frac{(\ell-1)\ell}{2}}p^{2(s-2)-\ell}8p^{\ell}p^{2((s-2)-\ell)}(1+O(p))=:\Delta''=8\Delta(1+O(q_s(n))),
$$
where the first summation is over $W'\in{[n]\setminus[3]\choose s-2}$, and the second one is over distinct $W_1,W_2\in{[n]\setminus[3]\choose s-2}$ having a non-empty intersection. Then, the Janson inequality implies that
$$
{\sf P}(\mu''=0)\leq
\mathrm{exp}\left[-\lambda''+\Delta''/2\right]=\mathrm{exp}\left[-2\lambda(1+O(1/\sqrt{n}))+4\Delta(1+O(q_s(n)))\right]=
$$
\begin{equation}
\exp\left[-2\lambda\left(1+o\left(n^{-1/(s+1)}\ln n\right)\right)\right].
\label{mu-two-from-above}
\end{equation}
Then, relations (\ref{expectation-mu})--(\ref{mu-two-from-above}) imply
$$
\mathrm{Var}N_s={\sf E}N_s^2-({\sf E}N_s)^2={\sf E}N_s+{n\choose 2}{n-2\choose 2}p^2{\sf P}(\mu(1,2)=\mu(3,4)=0)+
$$
$$
n(n-1)(n-2)p^2{\sf P}(\mu(1,2)=\mu(2,3)=0)-\left({n\choose 2}p{\sf P}(\mu(1,2)=0)\right)^2\leq
$$
$$
{\sf E}N_s+{n\choose 2}^2 p^2 e^{-2\lambda}\left(1+o\left(n^{-1/(s+1)}\ln^2 n\right)\right)+
$$
$$
n^3p^2 e^{-2\lambda}\left(1+o\left(n^{-1/(s+1)}\ln^2 n\right)\right) - {n\choose 2}^2 p^2 e^{-2\lambda}\left(1+o\left(n^{-1/(s+1)}\ln^2 n\right)\right)=o({\sf E}N_s)^2.
$$
Part (ii) follows from Chebyshev's inequality. $\Box$\\

Let $p(n)$ be as in Lemma~\ref{edge_in_clique} and $w_n\to-\infty$. As far as $G(n, p)$ contains an edge which is not contained in any copy of $K_s$, we cannot hope for $\mathcal{A}_s$ to be held. Indeed, assuming that there exists an edge $e$ in $G(n, p)$ that is not contained in any $K_s$, we immediately get that $e$ should belong to a weakly $K_s$-saturated graph. Assume that $F$ is a weakly $K_s$-saturated subgraph of $G(n,p)$ with ${s-2 \choose 2}+(n-s+2)(s-2)$ edges. Then $F\setminus e$ has ${s-2 \choose 2}+(n-s+2)(s-2)-1$ edges and is weakly $K_s$-saturated in $G(n,p)\setminus e$. But this is impossible since whp ${G}(n, p)\setminus \{e\}$ is weakly $K_s$-saturated in $K_n$ (the threshold for ${G}(n, p)$ to be weakly $K_s$-saturated in $K_n$ is at most $n^{-2/(s+1)-\gamma}$ \cite{bal}, and the same argument works for $G(n,p)\setminus \{e\}$). Part (ii) of Theorem~\ref{lower} follows.\\

To prove Part (i), consider $\mathcal{B}^*_s$, the graph property that every pair of vertices have $s-2$ adjacent common neighbors. This property is {\it increasing}~\cite[Chapter 1.3]{RG} with the {\it sharp threshold probability}~\cite[Chapter 1.6]{RG} $q^*_s(n)=(2(s-2)!)^{\frac{2}{(s+1)(s-2)}}q_s(n)$ (in the same way as in the proof of Lemma~\ref{edge_in_clique}, this fact follows from the second moment method and the Janson inequality, see details in~\cite{SP90}). Clearly, $\mathcal{A}_s\cap\mathcal{B}^*_s$ increases as well. Therefore, it has a threshold probability $r_s(n)\geq q^*_s(n)$. The following lemma concludes the proof of Theorem~\ref{lower}. $\Box$

\begin{lemma}
	$r_s(n)$ is threshold for $\mathcal{A}_s$.
\end{lemma}

{\it Proof}. Clearly, whp $G(n,p)$ has the property $\mathcal{A}_s$ if $p\gg r_s(n)$. It remains to prove the opposite for $p\ll r_s(n)$.

Let us first assume that there exists a sequence $n_k$ and a constant $C$ such that $r_s(n_k)<Cq^*_s(n_k)$ for all large enough $k$. Then, for this sequence $n_k$, the result immediately follows from Part (ii) of Theorem~\ref{lower}.

Therefore, we may assume that $r_s(n)>w_n q^*_s(n)$ for all $n$ and some $w_n\to\infty$ as $n\to\infty$. If $p(n) < c_sq_s(n)$ for infinitely many $n$, then whp (for these $n$) $G(n,p)$ does not have the property $\mathcal{A}_s$ by Part (ii) of Theorem~\ref{lower}. If $(1+\varepsilon)q^*_s(n_k)\leq p(n_k)\ll r_s(n_k)$ for some $\varepsilon>0$ and an infinite sequence $n_k$, then
$$
{\sf P}(G(n_k,p)\in\mathcal{A}_s)\leq{\sf P}(G(n_k,p)\notin\mathcal{B}^*_s)+{\sf P}(G(n_k,p)\in\mathcal{A}_s\cap\mathcal{B}^*_s)\to 0,\quad k\to\infty.
$$

Finally, assume that there exists an $\varepsilon>0$ and an infinite sequence $n_k$ (below, we omit $k$ for simplicity of notations) such that
$$
p_*(n_k):=c_s q_s(n_k)\leq p(n_k)\leq (1+\varepsilon)q^*_s(n_k)
$$
and ${\sf P}(G(n_k,p(n_k))\in\mathcal{A}_s)\geq\delta>0$.

Let $p_0(n)=p_*(n)\left[1+\frac{1}{\sqrt{\ln n}\ln\ln n}\right]$. If $p<p_0$, then $G(n,p_0)$ can be obtained from $G(n,p)$ by drawing every missing edge with probability $\frac{p_0-p}{1-p}$. Then, by~(\ref{expectation-mu}) and ~(\ref{mu}), the probability that there exists an edge $xy$ in $G(n,p_0)\setminus G(n,p)$ such that $x,y$ do not have $s-2$ adjacent common neighbors in $G(n,p)$ is at most
$$
{n\choose 2}\frac{p_0-p}{1-p}\exp\left[-{n-2\choose s-2}p_*^{\frac{(s+1)(s-2)}{2}}+o(1)\right]=O\left(\frac{1}{\ln\ln n}\right).
$$
Therefore,
$$
{\sf P}(G(n,p_0)\in\mathcal{A}_s)\geq{\sf P}(G(n,p)\in\mathcal{A}_s)-O\left(\frac{1}{\ln\ln n}\right)\geq\delta-o(1).
$$

So we may assume that $p\geq p_0$. In the same manner, $G(n,(1+\varepsilon)q^*_s(n))$ can be obtained from $G(n,p)$ by drawing every missing edge. Then, the probability that there exists an edge $xy$ in $G(n,(1+\varepsilon)q^*_s(n))\setminus G(n,p)$ such that $x,y$ do not have $s-2$ adjacent common neighbors in $G(n,p)$ is at most
$$
{n\choose 2}\frac{(1+\varepsilon)q^*_s(n)-p}{1-p}\exp\left[-{n-2\choose s-2}p_0^{\frac{(s+1)(s-2)}{2}}+o(1)\right]=O\left([\ln n]^{\frac{2}{(s-2)(s+1)}}\exp\left[-\frac{\sqrt{\ln n}}{\ln\ln n}\right]\right).
$$
Therefore,
$$
{\sf P}(G(n,(1+\varepsilon)q^*_s(n))\in\mathcal{A}_s)\geq{\sf P}(G(n,p)\in\mathcal{A}_s)-o(1)\geq\delta-o(1),
$$
a contradiction. $\Box$

\section{Upper bound on $p$}

\begin{theorem}
	Let $s\geqslant 3$ be an integer. If $p\geq n^{-\frac{1}{2s-3}}(\ln n)^2$, then $G(n,p)$ has the property $\mathcal{A}_s$ whp.
	\label{th_upper}
\end{theorem}

{\it Proof.} Since for $p$ in the range, whp $G(n,p)$ is weakly $K_s$-saturated in $K_n$ (see~\cite{bal}), we get that whp $\mathrm{wsat}(G(n,p),K_s)\geq\mathrm{wsat}(K_n,K_s)$.

It remains to prove that whp there exists a weakly $K_s$-saturated subgraph of $G(n,p)$ with ${s-2\choose 2}+(n-s+2)(s-2)$ edges. Below, we define two graph properties that imply such weak saturation `stability'. To define these properties, let us recall the following definition. {\it The $k$-th power of a graph} $G$ is the graph $G^k$ on the same vertex set such that $x,y$ are adjacent in $G^k$ if and only if the distance between them in $G$ is at most $k$.

Let $\mathrm{EXT}$ be the graph property `every set of $s$ vertices has $s-2$ common neighbors inducing a complete graph' and $\mathrm{HAM}$ be the graph property `for every set of $s-1$ vertices, its common neighborhood induces a graph containing $(s-2)$-th power of a Hamiltonian path'.

The following two lemmas conclude the proof of Theorem~\ref{th_upper}. $\Box$

\begin{lemma}\label{l1}
	If $\mathrm{EXT}$ and $\mathrm{HAM}$ hold for a graph $G$ with $n$ vertices, then $\mathrm{wsat}(G,K_s)\leq\mathrm{wsat}(K_n,K_s)$.
\end{lemma}

\begin{lemma}\label{l2}
	If $p\geq n^{-\frac1{2s-3}}(\ln n)^{2}$, then ${\sf P}(G(n, p)\in\mathrm{EXT}\cap\mathrm{HAM})\to 1$ as $n\to\infty$.
\end{lemma}

{\it Proof of Lemma~\ref{l1}.} Assume that $G$ is a graph satisfying $\mathrm{EXT}$ and $\mathrm{HAM}$. Let $H_0$ be a copy of $K_{s-2}$ in $G$. Define a weakly $K_s$-saturated spanning subgraph $H \subset G$ as follows. $H$ contains all edges of $H_0$, and also all edges of $G$ between $H_0$ and $N_G(H_0)$. We still have to add some other edges going outside $N_G[H_0]$. For every $v \in V(G)\setminus N_G[H_0]$, we add $s-2$ edges adjacent to $v$ described below. By the property $\mathrm{HAM}$, the graph $F_v$ induced on $N_G\big(V(H_0) \cup \{v\}\big)$ contains a copy of $(s-2)$-th power of a Hamiltonian path. Starting from an arbitrary vertex, denote the vertices of $F_v$ going in the natural order induced by the Hamiltonian path, by $x^v_1,x^v_2,\ldots,x^v_{f_v}$, where $f_v=|V(F_v)|$. We add edges $vx^v_1,vx^v_2,\ldots,vx^v_{s-2}$ to $H$. It is easy to see that $H$ is of size ${s-2 \choose 2}+(n-s+2)(s-2)$, so it suffices to prove that $H$ is weakly $K_s$-saturated in $G$.

First, all edges between the vertices of $N_G(H_0)$ can be `infected' since they belong to $K_s$ containing $H_0$. Second, for each $v \in V(G)\setminus N_G[H_0] $, we may `infect' edges $vx^v_{s-1},vx^v_s,\ldots,xv^v_{f_v}$ one by one since every such edge belongs to the $s$-clique containing previous $s-2$ vertices of the $(s-2)$-th power of the Hamiltonian path. Finally, each edge $xy$ between vertices of $V(G) \setminus N_G[H_0]$ can be `infected', since by the property $\mathrm{EXT}$, $N_G(\{x,y\}\cup V(H_0))$ contains a copy of $K_{s-2}$, say $K^*$. But note that $V(K^*)\subset N_G(V(H_0))$, so $xy$ is the last edge of the $s$-clique $G[V(K^*)\cup \{x,y\}]$ and hence can be `infected' as well. $\Box$\\

{\it Proof of Lemma~\ref{l2}.} It is known~\cite{SP90} that, for $p\gg n^{-\frac{2}{3(s-1)}}[\ln n]^{\frac{2}{3(s-2)(s-1)}}$, whp $G(n,p)$ has the property $\mathrm{EXT}$.

To prove the second property, notice that, for $k\in\mathbb{N}$, $r\in\mathbb{N}$, if $p(m)\gg m^{-\frac 1k}(\ln m)^4$, then $G(m, p)$ contains the $k$-th power of a Hamiltonian path with probability at least $1-m^{-r}$. This follows from the fact (see  \cite{posa}, Theorem 1.6 of \cite{fischer}, and Theorem 1.1 of \cite{kuhn}) that
there exists $C$ such that if $p^*(m) > C m^{-\frac 1k}(\ln m)^3$, then ${G}(m, p^*)$ contains the $k$-th power of a Hamiltonian path with probability at least $1-1/e$. To boost this probability to $1-m^{-r}$ it suffices to take the union of $r\ln m $ independent copies of ${G}(m, p^*)$ on $[m]$ and set $p=1-(1-p^*)^{r\ln m}$.

Fix $W_0\subset[n]$ of size $s-1$ and set $m=\frac{np^{s-1}}{2}$. Then, by the Chernoff bound,
the probability that the common neighborhood of $W_0$ does not contain the $(s-2)$-th power of a Hamiltonian path is at most
$$
{\sf P}(|N_{G(n,p)}(W_0)|<m)+m^{-3s}{\sf P}(|N_{G(n,p)}(W_0)|\geq m)\leq e^{-\frac{m}{4}}+n^{-s}.
$$
From the union bound, the result follows. $\Box$

\begin{remark}
	The power in the log factor in Theorem~\ref{th_upper} is not the best possible. Using best known estimations of threshold probability for containing the $k$-th power of a Hamiltonian path, we may replace $[\ln n]^2$ with $[\ln n]^{1/3+\varepsilon}$ for $s=3$, $[\ln n]^{8/5}$ for $s=4$ and $[\ln n]^{1/2}$ for $s\geqslant 5$.
	\label{rmrk}
\end{remark}

\section{Upper bound on $\mathrm{wsat}\big({G}(n, p), K_s\big)$}

From the previous arguments, we know that whp $\mathrm{wsat}\big({G}(n, p), K_s\big)=e(G(n,p))(1+o(1))$ when $p\ll n^{-\frac{2}{s+1}}$ and $\mathrm{wsat}\big({G}(n, p), K_s\big)={s-2\choose 2}+(n-s+2)(s-2)$  when $p\geq n^{-\frac{1}{2s-3}}(\ln n)^{\beta}$, where $\beta=1/3+\varepsilon$ for $s=3$, $\beta=8/5$ for $s=4$ and $\beta=1/2$ for $s\geqslant 5$. In this section, we prove an upper bound for $\mathrm{wsat}\big({G}(n, p), K_s\big)$ for all the remaining values of $p$ (the trivial lower bound equals ${s-2\choose 2}+(n-s+2)(s-2)$ whp).

\begin{thm}\label{upperbd}
	Let $\gamma=1$ for $s=3$, $\gamma=6$ for $s=4$ and $\gamma=0$ for $s\geq 5$. Let $w_n\to\infty$ as $n\to\infty$. Then  whp $$\mathrm{wsat}\big(G(n, p), K_s\big)\leqslant n(s-2)+ \frac{w_n(\ln n)^{2(\gamma+s-2)}}{p^{2s-3}}.$$
\end{thm}
\begin{proof}

For $p\leq\left((\ln n)^{\gamma+s-2}/n\right)^{1/(s-1)}$, the result is immediate since whp $|E(G(n,p))|\ll\frac{w_n(\ln n)^{2(\gamma+s-2)}}{p^{2s-3}}$. For $p\geq n^{-\frac{1}{2s-3}}(\ln n)^{\beta}$ the result follows from Remark~\ref{rmrk}. Let
$$
\left((\ln n)^{\gamma+s-2}/n\right)^{1/(s-1)}<p<n^{-\frac{1}{2s-3}}(\ln n)^{\beta}.
$$

Let $G$ be a graph and let  $X$ be a subset of  $V(G)$ having  two properties:
	\begin{itemize}
		\item [(i)] For  every vertex $v\in \overline{X}=V(G)\setminus X$ the subgraph of $G$ induced on $N_G(v)\cap X$  contains an $(s-2)$-th power of a Hamiltonian path.
		\item [(ii)] Any two distinct vertices $v,w\in \overline{X}$ have $s-2$ pairwise adjacent common neighbors in $X$.
	\end{itemize}
 Let	$H$ be a spanning subgraph of $G$ containing all edges with both endpoints in $X$ and $(s-2)|\overline{X}|$ more edges such that every vertex $v\in\overline{X}$ belongs to exactly $s-2$ edges $\{v,w_1\},\ldots,\{v,w_{s-2}\}$ where $w_1,\ldots,w_{s-2}$ are initial vertices of an $(s-2)$-th power of a Hamiltonian path of $G[N_G(v)\cap X]$.

Let us show that $H$ is weakly $K_s$-saturated in $G$. In the same way, as in the proof of Lemma~\ref{l1}, all edges having one vertex in $X$ can be `infected'. Let $\tilde H$ have all edges of $H$ and all edges between $X$ and $\overline{X}$ that are presented in $G$. Let $v,w\in\overline{X}$ be adjacent in $G$. Then, the property (ii) ensures that adding the edge $\{v,w\}$ to $\tilde H$ creates a copy of $K_s$. Thus, $\mathrm{wsat}(G, K_s)\leqslant (s-2)|\overline{X}|+e(G[X])$.

Let $G=G(n, p)$. Let $m=\left\lfloor\frac{(\ln n)^{\gamma+s-2}}{p^{s-1}}\sqrt{w_n}\right\rfloor$, $X=[m]$. Since whp $e(G[X])\leq m^2p $, it remains to show that whp (i) and (ii) hold for $G$ with our choice of $X$.

As we noticed in the proof of Lemma~\ref{l2}, for any  $\ell\geq\left\lfloor\frac{(\ln n)^{\gamma+s-2}}{p^{s-1}}\right\rfloor$, the relation $p\gg \ell^{-\frac {1}{s-2}}(\ln n)^{\gamma/(s-2)+1}$ implies that $G(n,p)[\ell]$ contains the $(s-2)$-th power of a Hamiltonian path with probability at least $1-o(1/n)$ (actually, the power of log-factor in the proof of Lemma~\ref{l2} is 4, but it can be improved in this way due to best known estimations of threshold probability for containing the $k$-th power of a Hamiltonian path, see~\cite{Frieze_Survey}). Given a vertex $v\in[n]\setminus[m]$, with probability $o(1/n)$, it has less than $\frac{(\ln n)^{\gamma+s-2}}{p^{s-1}}$ neighbors in $[m]$ by the Chernoff bound. Property (i) follows.

Finally, by the Janson inequality~\cite{Janson}, the probability that, there exist distinct  $v,w\in \overline{X}$ without $s-2$ pairwise adjacent common neighbors in $[\tilde m]$, $\tilde m=\left\lfloor\frac{\ln n}{p^{(s+1)/2}}\right\rfloor<m$, is at most
\begin{align*}
 {n\choose 2}&\exp\left(-{\tilde m\choose s-2}p^{{s\choose 2}-1}+\sum_{t=1}^{s-3}{\tilde m\choose s-2}{s-2\choose t}{\tilde m-s+2\choose s-2-t}p^{(s+1)(s-2)-2t-{t\choose 2}}\right)\leq\\
 &\exp\left(2\ln n-\frac{1}{(s-2)!}\tilde m^{s-2}p^{(s+1)(s-2)/2}(1-o(1))\right)\to 0,\quad n\to\infty. \qedhere
\end{align*}
\end{proof}

\section*{Acknowledgments}

Maksim Zhukovskii gratefully acknowledge the financial support from the Ministry of Educational and Science of the Russian Federation in the framework of MegaGrant no 075-15-2019-1926


TITLE: Is there an analogue of Mordell-Weil theorem for other fields?
QUESTION [0 upvotes]: The Mordell-Weil theorem states that for an abelian variety $A$ over a number field $K$ the group of $K$-rational points of $A$ is finitely generated and abelian.
What if $K$ is not a number field, e. g. Q$_p$ or a transcendental extension of $\mathbb Q$?

REPLY [2 votes]: For $K=\Bbb Q_p$ there will be uncountably many $K$-rational points.
The exact structure of $K$-rational points will depend on the
reduction type of the curve, but points close to the base-point $O$ are
in essence parameterised by the formal group of the curve. See Silverman's book for much more detail.