% this file is called up by thesis.tex
% content in this file will be fed into the main document

%: ----------------------- name of chapter  -------------------------
\chapter{A statistical approach to DA}\label{C:da_bayesian} % top level followed by section, subsection
The deterministic approach to \gls{DA} presented in Chapter \ref{C:da_nse} provides an accurate and noise filtering estimator of the velocity of the fluid retrieved from noisy measures. A weak point of such approach is that, relying on a specific realization of noise it produces a single estimate of the state which might not be representative of the whole state distribution. Moreover, in the deterministic approach measures are considered as if they were independent. In the context of hemodynamics, which is characterized by complex phenomena, this is a limiting feature since additional a priori information regarding the noise is not taken into account.

In medical applications an evaluation of reliability of simulations and numerical results is mandatory; for our method to be a valid tool supporting medical practice it is important to quantify the credibility of predictions in statistical terms. These include computing the most likely value for the velocity or a set of values which the state is most likely to belong to. These goals are achieved when we determine the statistical distribution of the velocity and related variables; this task requires prior information about the variables involved and of the noise affecting the observations. 

In this chapter, based on \cite{md_av_uq} (in preparation), we present a Bayesian approach to \gls{DA}; the purpose is to find a probability distribution for the state variable, treated as random. Parameters of the probability density function (\gls{PDF}), such as mean value and variance, and induced confidence regions allow to quantify the uncertainty of numerical simulations.

In Section \ref{S:multivariate_normal} we introduce basic statistical notation and properties of normal random vectors which are useful in dealing with our application. In Section \ref{S:math_form} we introduce the mathematical formulation of the Bayesian method applied to the linearized \gls{NSE} and we extend the formulation to the nonlinear case using the algorithm for the nonlinear equations introduced in Chapter \ref{C:da_nse}. Then, we present point estimators, namely the maximum a posteriori (\gls{MAP}) and the maximum likelihood (\gls{ML}) estimators, and spread estimators (confidence regions) based on the \gls{PDF} of the normal stress on the inflow boundary. In Section \ref{S:interpolation} we recall the regularization technique based on interpolation developed in Chapter \ref{C:da_nse} and we show how it affects the statistical formulation. In Section \ref{S:bayes_num_res} we present numerical results in correspondence of the same test cases used to test the deterministic formulation in Chapter \ref{C:da_nse}. We compare statistical and deterministic point estimators and we compute confidence regions for the velocity ans the wall shear stress \gls{WSS}.
%: ----------------------- paths to graphics ------------------------
% change according to folder and file names
\ifpdf
    \graphicspath{{4/figures/PNG/}{4/figures/PDF/}{4/figures/}}
\else
    \graphicspath{{4/figures/EPS/}{4/figures/}}
\fi
%: ----------------------- contents from here ------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The multivariate normal distribution}\label{S:multivariate_normal}
In this section we report basic concepts regarding the normal distribution for random vectors and we introduce properties and results specifically related to our application. In particular, we focus on statistical properties of affine transformations of normal random vectors which are used in Section \ref{S:spread_est} to derive the velocity \gls{PDF}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The multivariate normal PDF}
The multivariate normal distribution is a generalization to dimensions $d\geq 2$ of the univariate normal distribution. For the univariate case, a Gaussian random variable $X$ with mean $\mu$ and variance $\sigma^2$ has \gls{PDF}
\begin{equation}
f(x) = \dfrac{1}{\sqrt{2\pi\sigma^2}} exp\left\{-\dfrac{(x-\mu)^2}{2\sigma^{2}}\right\}\;\;\forall\;x\in(-\infty,\,\infty),
\end{equation}
and it is denoted $X\sim\mathcal{N}(\mu,\,\sigma^2)$. 

In the multivariate case the formula is generalized to a $d$-dimensional vector; in particular, the exponent has the form
\begin{equation}
(\xx - \boldsymbol{\mu})^{\TT}\Sigma^{-1}(\xx-\boldsymbol{\mu})
\end{equation}
where $\boldsymbol{\mu}\in\mathbb{R}^d$ represents the expected value of the multivariate random variable $\bf X$ and $\Sigma$ is the symmetric positive definite covariance matrix of $\bf X$, i.e. $\Sigma_{i,i} = var(\XX_i)$ and $\Sigma_{i,j} = \Sigma_{j,i} = cov(\XX_i,\,\XX_j)$ for $i,\,j= 1,\,...,\,d$. The normalizing constant $1/\sqrt{2\pi \sigma^2}$ has to be generalized to a constant such that the volume under the density surface is unitary. It can be shown, see \cite{johnson_wichern} that such constant is
\begin{equation}
\dfrac{1}{\sqrt{(2\pi)^{d} det(\Sigma)}}.
\end{equation}
Thus, the multivariate normal \gls{PDF} has the form
\begin{equation}
f(\xx) = \dfrac{1}{\sqrt{(2\pi)^{d} det(\Sigma)}} exp\left\{-(\xx - \boldsymbol{\mu})^{\TT}\Sigma^{-1}(\xx-\boldsymbol{\mu})\right\}\;\;\forall\;x_i\in(-\infty,\,\infty),\; i = 1,\,...,\,d
\end{equation}
and we say ${\bf X}\sim\mathcal{N}(\boldsymbol{\mu},\,\Sigma)$. It is easy to show that $\boldsymbol{\mu}$ is the point of maximum density, mode and expected value \cite{johnson_wichern}.
%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Contour lines} Contours of constant density $c$ play a central role in determining confidence regions for multivariate normal variables. For the $d$-dimensional distribution these are ellipsoids generated by the equation
\begin{equation}
(\xx - \boldsymbol{\mu})^{\TT}\Sigma^{-1}(\xx-\boldsymbol{\mu}) = c^2;
\end{equation}
these ellipsoids are centered in $\boldsymbol{\mu}$ and have axes $\pm c\sqrt{\lambda_j}{\bf e}_j$, where $\lambda_j$ are the eigenvalues of $\Sigma$ and ${\bf e}_j$ the associated eigenvectors for $j = 1,\,...,\,d$.
%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Properties of the multivariate normal distribution}\label{S:multivariate_properties}
In this section we report a few properties which are relevant for results discussed in this chapter and we introduce the chi-squared, $\chi^2$, distribution.

For a $d$-dimensional normal random vector $\XX\sim\mathcal{N}(\boldsymbol{\mu},\,\Sigma)$ we have the following properties:
\begin{description}
\item[P1] Affine transformations of $\XX$ are normally distributed. In particular, given ${\rm A}\in\mathbb{R}^{q\times d}$ and ${\bf b}\in\mathbb{R}^q$, then $\YY = {\rm A}\XX+{\bf b}\sim \mathcal{N}_q({\rm A}\boldsymbol{\mu} + {\bf b},\, {\rm A}\Sigma {\rm A}^{\TT})$.

\item[P2] All subsets of the components of $\bf X$ have normal distribution. If we rearrange the vector as
\begin{equation}
\XX = \left[\begin{array}{c}
\XX_1\\
\XX_2 
\end{array}\right], \quad 
\boldsymbol{\mu} = \left[\begin{array}{c}
\boldsymbol{\mu}_1\\
\boldsymbol{\mu}_2 
\end{array}\right], \quad
\Sigma = \left[\begin{array}{cc}
\Sigma_{1,1} & \Sigma_{1,2}\\
\Sigma_{2,1} & \Sigma_{2,2}   
\end{array}\right],
\end{equation}  
where $\XX_1\in\mathbb{R}^{q}$ and $\XX_2\in\mathbb{R}^{d-q}$; then, $\XX_1\sim\mathcal{N}_q(\boldsymbol{\mu}_1,\,\Sigma_{1,1})$ and $\XX_2\sim\mathcal{N}_{d-q}(\boldsymbol{\mu}_2,\,\Sigma_{2,2})$.

\item[P3] We first introduce the chi-squared distribution.
%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{The chi-squared distribution}  The chi-squared distribution with $d$ degrees of freedom is the distribution of a sum of the squares of $d$ independent standard normal random variables. If $X_1,\, ...,\, X_d$ are independent, standard, normal random variables, then the sum of their squares,
\begin{equation}\label{E:chi_sq}
Z = \sum\limits_{i=1}^d X_i^2,
\end{equation}
is distributed according to the chi-squared distribution with $d$ degrees of freedom. This is usually denoted as
\begin{equation}\label{E:chi_sq_notation}
Z\sim\ \chi^2_d.
\end{equation} 
The \gls{PDF} of $Z$ is given by
\begin{equation}\label{E:chi_sq_pdf}
f(z;\,d) = \left\{\begin{array}{ll}
\dfrac{1}{2^{d/2}\Gamma(d/2)}\,z^{d/2 - 1} e^{-z/2},  	& z \geq 0;\\
[3mm]
0, 														& \text{otherwise}
\end{array}\right.
\end{equation}
where $\Gamma$ denotes the Gamma function defined as $\Gamma(n) = (n-1)!$, for an integer $n$, and as $\Gamma(\dfrac{n}{2})= \pi\dfrac{(n-2)!!}{2^{(n-1)/2}}$ for a half-interger $\dfrac{n}{2}$.\\

The third property relates normal vectors and the chi-squared random variable and allows us to draw confidence regions for Gaussian vectors.
\begin{itemize}
\item $(\XX - \boldsymbol{\mu})^{\TT}\Sigma^{-1}(\XX-\boldsymbol{\mu})$ is distributed as $\chi_d^2$, where $\chi_d^2$ denotes the chi-squared distribution \cite{johnson_wichern} with $d$ \gls{DOF}s.
\item The $\mathcal{N}(\boldsymbol{\mu},\,\Sigma)$ distribution assigns probability $(1-\alpha)$ to the ellipsoid $\{\xx:\; (\xx - \boldsymbol{\mu})^{\TT}\Sigma^{-1}(\xx-\boldsymbol{\mu})\leq \chi^2_d(\alpha)\}$ where $\chi^2_d(\alpha)$ denotes the upper $(100\,\alpha)$-th percentile of the $\chi^2_d$ distribution.
\end{itemize}
\end{description}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Mathematical formulation}\label{S:math_form}
In this section we formulate the assimilation problem as a statistical inverse problem using the Bayesian approach; an extensive introduction to statistical inverse problems can be found in the book by Kaipio and Somersalo \cite{somersalo}.
As done in Chapter \ref{C:da_nse} we present the formulation for the linearized \gls{NSE} and then we extend it to the nonlinear formulation. Then, we derive statistical estimators for the velocity $\UU$ in $\dom$. 

The purpose of statistical inversion is the prediction of stochastic features of the variables of interest in order to quantify the credibility of numerical methods and simulations. The prediction of the uncertainty affecting the variables is based on the knowledge of (i) the measurement process; (ii) deterministic models of the unknowns \cite{somersalo}.

In formulating the inverse problem we assume to deal with discretized variables, all treated as random. The randomness is in the degree of information of their realizations and such degree resides in the probability distributions. The entities involved are \gls{PDF}s; in this respect, statistical approaches are different from deterministic ones which produce a single estimate of the unknown. Here the method delivers a distribution. In the following sections the bold variables denote random vectors while the capital plain variables a specific realization. With an abuse of notation we introduce the random variable $\HH$ which describes the normal stress of the fluid at the inflow section; $\MB$ is the random variable that describes the measures and $\boldsymbol{\varepsilon}$ the noise perturbing the measurements. We let $\pi_{pr}(H)$ be the \gls{PDF} of $\HH$, usually called {\it prior}, and $\pi_{noise}(\varepsilon)$ the one of $\boldsymbol{\varepsilon}$; these distributions are assumed to be known. The {\it likelihood} function $\pi(M|H)$ denotes the \gls{PDF} of the measurement conditioned on a realization of $\HH$. The distribution of $\HH$ conditioned on a realization of $\MB$, $\pi(H|M) = \pi_{post}(H)$, is usually called {\it posterior}. The purpose of the Bayesian procedure is to estimate the posterior exploiting the Bayes formula (see e.g. \cite{somersalo}) 
\begin{equation}\label{E:bayes_thm}
\pi(H|M) = \pi_{post} (H) = \frac{\pi(M|H) \pi_{pr}(H)}{\pi(M)};
\end{equation}
where $\pi(M)$ is the \gls{PDF} of the measures. Since we are interested in finding the variable $\HH$, that maximizes the posterior distribution, the denominator does not affect the optimization process and we can simply consider the relation 
\begin{equation}\label{E:bayes_thm_prop}
\pi(H|M) = \pi_{post} (H) \propto \pi(M|H) \pi_{pr}(H).
\end{equation}
Also, let us consider the linearized model that relates $\HH$ and $\MB$,
\begin{equation}\label{E:add_noise}
{\rm Z}\HH + \boldsymbol{\varepsilon} = \MB.
\end{equation}
Here, ${\rm Z} = {\rm D}{\rm S}^{-1}{\rm R}_{in}{\rm M}_{in}$ is the matrix introduced in Section \ref{S:da_nse_DO}, which relates the random variable representing the velocity measures to the variable $\HH$ according to the deterministic model introduced in the previous section. This relation holds in the linearized case and takes into account the presence of the noise. The random variable $\boldsymbol{\varepsilon}$ accounts for the discrepancy between ${\rm Z}\HH$ and $\MB$. Models in the form \eqref{E:add_noise} are referred to as {\it additive noise} and they are used with the assumption of mutual independence of $\HH$ and $\boldsymbol{\varepsilon}$. This assumption implies that the \gls{PDF} of $\boldsymbol{\varepsilon}$ is independent of any realization of $\HH$, say $\overline{H}$.
%When we fix $\HH=\overline{h}$ (a specific realization of $\HH$) the assumption ensures that the density function of $\boldsymbol{\varepsilon}$ remains unaltered when conditioned on $\overline{h}$. Therefore, we deduce that $\pi(y|\overline{h})$ is distributed like $\boldsymbol{\varepsilon}$ with density function translated by ${\rm Z}\overline{h}$.
Hence, the likelihood function can be expressed as
\begin{equation}\label{E:likelihood}
\pi(M|\overline{H}) = \pi_{noise}(M-{\rm Z}\overline{H}).
\end{equation}
If we consider $M = \dd$, a specific realization of $\MB$, we rewrite \eqref{E:bayes_thm_prop} as 
\begin{equation}\label{E:bayes_thm_bis}
\pi_{post}(H) \propto \pi_{noise}(\dd-{\rm Z}H) \pi_{pr}(H).
\end{equation}
When we consider Gaussian random variables we can derive an explicit form for $\pi_{post}(H)$. We introduce the \gls{PDF}s of $\HH$ and $\boldsymbol{\varepsilon}$, according to the notation reported in Section \ref{S:multivariate_normal}:
\begin{equation}\label{E:pr_like_H}
\begin{array}{l} 
\pi_{pr}(H) \propto
exp\left\{ -\dfrac{1}{2}(H-H_0)^{\TT}\Sigma^{-1}_{pr}(H-H_0) \right\},\\
[4mm]
\pi_{noise}(\varepsilon) \propto exp\left\{- \dfrac{1}{2}(\varepsilon-\varepsilon_0)^{\TT}\Sigma^{-1}_{noise}(\varepsilon-\varepsilon_0)\right\};
\end{array}
\end{equation}
where $H_0$ and $\varepsilon_0$ are the expectation values and $\Sigma_{pr}$ and $\Sigma_{noise}$ are the correlation matrices for $\HH$ and $\boldsymbol{\varepsilon}$ respectively. In the Gaussian assumption, using Theorem 3.7, Chapter 3 of \cite{somersalo}, we can write \eqref{E:bayes_thm_bis} as
\begin{equation}\label{E:postPDF}
\pi_{post}(H) \propto exp\left\{-\frac{1}{2}(H-H_{post})^{\TT}\Sigma^{-1}_{post}(H-H_{post})\right\};
\end{equation}
where the correlation matrix and the expectation value read
\begin{equation}\label{E:postH}
\begin{array}{l} 
\displaystyle\Sigma_{post} = (\Sigma_{pr}^{-1} + {\rm Z}^{\TT}\Sigma_{noise}^{-1}{\rm Z})^{-1},\\
[5mm]
\displaystyle H_{post} = \Sigma_{post}({\rm Z}^{\TT}\Sigma_{noise}^{-1}(\dd-\varepsilon_0) +  \Sigma_{pr}^{-1}H_0) \Rightarrow \\
[2mm]
\displaystyle (\Sigma_{pr}^{-1} + {\rm Z}^{\TT}\Sigma_{noise}^{-1}{\rm Z}) H_{post} = {\rm Z}^{\TT}\Sigma_{noise}^{-1}(\dd-\varepsilon_0) +  \Sigma_{pr}^{-1}H_0.
\end{array}
\end{equation}
Using this result one can calculate point and spread estimates. As an example, we can answer questions like ``how likely is the unknown inside a set of values?''.
%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The formulation for the nonlinear NSE}\label{S:bayes_nnlin}
In treating the nonlinearity we do not replace \eqref{E:add_noise} with a nonlinear model; instead, we consider an iterative approach inspired by the deterministic one. Also in this case, we rely on the Newton method for the \gls{NSE}; the algorithm reads as follows. Given a guess for the velocity at iteration $k+1$, we find $H_{post,k+1}$ solving 
\begin{equation}\label{E:nonlinearhMAP}
H_{post,k+1} = (\Sigma_{pr}^{-1} + {\rm Z}_k^{\TT}\Sigma_{noise}^{-1}{\rm Z}_k)^{-1}({\rm Z}^{\TT}\Sigma_{noise}^{-1}(\dd-\varepsilon_0) +  \Sigma_{pr}^{-1}h_0).
\end{equation}
Here, ${\rm Z}_k$ is defined as in \ref{S:da_nse_nonlinear}, i.e. ${\rm Z}_k = {\rm D}{\rm S}_k^{-1}{\rm R}_{in}^{\TT}{\rm M}_{in}$ where
\begin{equation}\label{E:discr_lin_ns_BIS}
\textrm{S}_k=\left[\begin{array}{ll}
\textrm{C} + \textrm{A}_k & \textrm{B}^\TT\\
\textrm{B} & \textrm{O} 
\end{array}\right].
\end{equation}
${\rm A}_k$ is the discretization of the advection operator with advection field $U_{k}$, the velocity vector associated with the normal stress $H_{post,k}$. We perform iterations \eqref{E:nonlinearhMAP} until a convergence criterion is satisfied. Since in our application the variable of interest is the velocity, we check the convergence using $\|U_{k}-U_{k+1}\|\leq \delta$, a user defined tolerance. Note that with this formulation $\HH$ and $\UU$, at each iteration, are related by a linear model and, for this reason, $\UU$ can still be considered normally distributed.
%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Statistical point estimators}\label{S:point_estimators}
One of the most common point estimators is the maximum a posteriori (\gls{MAP}) estimator. We let $H_{MAP}$ be the \gls{MAP} estimator for $\HH$. It is defined as the most likely value of $\HH$ given $\dd$, formally
\begin{equation}\label{E:h_map_est}
H_{MAP} = \arg\max\limits_{H}\;\pi_{post}(H).
\end{equation}
In the Gaussian assumption $H_{MAP}$ corresponds to $H_{post}$, the expected value of the posterior distribution and can be found by solving 
\begin{equation}
\displaystyle (\Sigma_{pr}^{-1} + {\rm Z}^{\TT}\Sigma_{noise}^{-1}{\rm Z}) H_{post} = {\rm Z}^{\TT}\Sigma_{noise}^{-1}(\dd-\varepsilon_0) +  \Sigma_{pr}^{-1}H_0.
\end{equation}
This leads to a computational problem similar to classical deterministic regularization methods; nevertheless, the estimate comes from different motivation and relies on different information.

Another common estimator is the maximum likelihood (\gls{ML}) estimator; this is the value of $\HH$ which is most likely to produce the data $\dd$. It is defined as
\begin{equation}\label{E:hML_def}
H_{ML} = \arg\max\limits_{H} \;\pi(M|H).
\end{equation}
For the Gaussian distribution, combining \eqref{E:likelihood} and \eqref{E:pr_like_H}, this corresponds to solving
\begin{equation}\label{E:hML_opt1}
\arg\max\limits_H \;exp\left\{- \dfrac{1}{2}(\dd-{\rm Z}\,H-\varepsilon_0)^{\TT}\Sigma^{-1}_{noise}(\dd-{\rm Z}\,H-\varepsilon_0)\right\}
\end{equation}
or, equivalently
\begin{equation}\label{E:hML_opt2}
\arg\min\limits_H \; - \dfrac{1}{2}(\dd-{\rm Z}\,H-\varepsilon_0)^{\TT}\Sigma^{-1}_{noise}(\dd-{\rm Z}\,H-\varepsilon_0).
\end{equation}
This leads to the solution of the following linear system
\begin{equation}\label{E:hML}
({\rm Z}^{\TT}\Sigma_{noise}^{-1}{\rm Z}) H_{ML} = {\rm Z}^{\TT}\Sigma_{noise}^{-1}(\dd-\varepsilon_0).
\end{equation}
It is important to note that this is not a Bayesian estimator since it does not rely on the Bayes formula; it only requires the knowledge of stochastic features of the measurement process, i.e. the likelihood function. Note that, when the matrix of system \eqref{E:hML} is singular the estimator is not well-defined. This happens when either  $\Sigma_{noise}$ is singular or ${\rm Z}$ is rank-deficient. In the latter case, as proved in Chapter \ref{C:da_nse}, Proposition \ref{P:suff_cond_Qmatrix}, singularity is due to data location. When sites on $\Gamma_{in}$ do not satisfy optimality sufficient conditions we can prevent singularity of ${\rm Z}^{\TT}{\rm Z}$ using interpolation, see Section \ref{S:wellp_2}. This process alters the noise; in Section \ref{S:interpolation} we provide details on how statistical properties of the noise are affected by data interpolation.
%and it is not well-defined in case of ill-posed inverse problems. In our case this corresponds to having data not satisfying sufficient conditions for optimality (see Proposition \ref{P:suff_cond_Qmatrix}); as a consequence, the matrix of system \eqref{E:hML} is singular.

\gls{MAP} and \gls{ML} estimators are strictly related and the choice between the two depends on the level of prior knowledge on the variable of interest; $H_{ML}$ corresponds to rejecting (or not trusting) our prior belief on $\HH$, while the \gls{MAP} estimation pulls the estimate towards the prior. 

In a certain sense the $H_{ML}$ can be considered as the limit case of $H_{MAP}$ \cite{bayes_choice}. Consider the following iterative algorithm. Let $\pi_{j-1}$ be an estimate of the \gls{PDF} of $\HH$ at iteration $j-1$; then,
\begin{enumerate}
\item set $\pi_{pr}= \pi_{j-1}$;
\item compute $\pi_{post}$ and $H_{post}$ using \eqref{E:postH};
\item set $\pi_{j}= \pi_{post}$, $H_{MAP,j} = H_{post}$, $\,j=j+1$;
\end{enumerate}
It is possible to show that, regardless of how we choose $\pi_0$ and $H_0$, iterating the procedure we have $\lim\limits_{j\rightarrow\infty}H_{MAP,j}= H_{ML}$ \cite{bayes_choice}.

In using the \gls{MAP} estimator some knowledge on the distribution of $\HH$ must be assumed. For this reason, we introduce the ``Gaussian smoothness'' priors \cite{somersalo}; they correspond to prior models with structural information on $\HH$. When \eqref{E:add_noise} comes form the discretization of \gls{PDE}s, as in our case, it is a common assumption to consider $\HH$ having smoothness properties (e.g. second order differentiability). In this work we choose $\Sigma_{pr}^{-1} = \frac{1}{\sigma^2} {\rm L}^{\TT}{\rm L}$, where L corresponds to the discretization of a differential operator. As an example, when L is the identity operator and $\sigma^2=1$, we are assuming $\HH$ to be corrupted by white noise. 
%%%%%%%%%%%%%%%%%%%%%
\subsection{Statistical spread estimators}\label{S:spread_est}
Beyond point estimates, much broader information is carried by the posterior distribution. In this section we define credibility regions for the velocity vector, i.e. set of values the variable belongs to with a certain probability.

Using properties of affine transformations of random variables and the theory of confidence regions for the multivariate normal distribution, see Section \ref{S:multivariate_properties}, we introduce the \gls{PDF} of $\UU$ (the random variable associated with the discretized velocity) and credibility regions for components of the velocity vector in each DOF.

For the linearized \gls{NSE}, an affine transformation relates the normal stress $\HH$ and the velocity field $\UU$, i.e. the discretization of the state equation: ${\rm S}\VV = {\rm R}_{in}^{\TT}{\rm M}_{in}\HH + \FF$. Solving for $\VV$, where $\VV = [\UU\; \PP]$, we have
\begin{equation}
\VV = {\rm S}^{-1}{\rm R}_{in}^{\TT}{\rm M}_{in}\HH +  {\rm S}^{-1}\FF.
\end{equation}
This is equivalent to
\begin{equation}\label{E:u_map}
\UU = {\rm E}({\rm S}^{-1}{\rm R}_{in}^{\TT}{\rm M}_{in}\HH +  {\rm S}^{-1}\FF);
\end{equation}
where $\rm E$ is the matrix that extracts the velocity, i.e. such that ${\rm E}\VV = \UU$. From now on we set ${\rm T} = {\rm E}{\rm S}^{-1}{\rm R}_{in}^{\TT}{\rm M}_{in}$. This affine transformation maps the Gaussian variable $\HH\sim \mathcal{N}(H_{post}, \,\Sigma_{post})$ into the Gaussian variable $\UU\sim\mathcal{N}(U,\,\Sigma_U)$ where $U = \TT \,H_{post}+{\rm E}{\rm S}^{-1}\FF$ and $\Sigma_U = \TT\Sigma_{post}\TT^{\TT}$, since we do not consider $\FF$ affected by uncertainties. Our goal is to have a measure of how likely the velocity, or one of its components, falls within a set of values, or, conversely, how large is the set of values corresponding to probability $\alpha$ for the velocity to belong to such set. To this end, we introduce confidence regions for the velocity field; more specifically, we focus on credibility regions for the horizontal and vertical velocity in each \gls{DOF} of the computational grid, i.e. $[\UU_i \;\; \UU_{i+N_u/2}]^{\TT}\in \mathbb{R}^2$, for $i = 1,\,...,\,N_u/2$. Note that we assume horizontal velocities to be stored in the first half of the vector. This 2-dimensional vector has a bivariate normal distribution; 
\begin{equation}\label{E:u_distribution}
\left[\begin{array}{l}
\UU_i\\
\UU_{i+N_u/2}
\end{array}
\right]\sim \mathcal{N}\left(
\left[\begin{array}{l}
U_i\\
U_{i+N_u/2}
\end{array}
\right], \; 
\left[\begin{array}{ll}
\Sigma_{U;\,i,i} & \Sigma_{U;\,i,i+N_u/2}\\
\Sigma_{U;\,i+N_u/2,i} & \Sigma_{U;\,i+N_u/2,i+N_u/2}
\end{array}\right]\right)
\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.7\textwidth]{ellips1.pdf}
\end{center}
\caption{Confidence region associated with $\chi_2^2(\alpha)$.}
\label{F:ellips1}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.7\textwidth]{ellips2.pdf}
\end{center}
\caption{Confidence region for the standardized variable corresponding to percentile 1.}
\label{F:ellips2}
\end{figure} 
%%%%%%%%%%%%%%%%%%%%%%%%%%
For the sake of notation, we use the same notation of Section \ref{S:multivariate_normal} and we rewrite \eqref{E:u_distribution} as ${\bf X} \sim \mathcal{N}(\boldsymbol{\mu},\,\Sigma)$. Confidence intervals for $\UU$ can be found using property {\bf P3} (see section \ref{S:multivariate_properties}) for the 2-dimensional case. Since $(\XX - \boldsymbol{\mu})^{\TT}\Sigma^{-1}({\bf X}-\boldsymbol{\mu})$ has a chi-squared distribution with 2 degrees of freedom, probability of $1-\alpha$ can be assigned to the ellipsis $\{\xx:\; (\xx - \boldsymbol{\mu})^{\TT}\Sigma^{-1}(\xx-\boldsymbol{\mu})\leq \chi^2_2(\alpha)\}$ in the $(x,\,y)$ plane, where $\chi^2_2(\alpha)$ denotes the upper $(100\alpha)$-th percentile of the distribution. Using this result, for a fixed value of $\alpha$, we can draw regions in the 2-dimensional plane as shown in Figure \ref{F:ellips1}. Spectral properties of the correlation matrix $\Sigma$ play a fundamental role; the eigenvectors of $\Sigma$ are the directions of the main axes and the eigenvalues are a relative measure of the variability of the velocity components. On the other hand the ellipsis in Figure \ref{F:ellips2} is a rotation and translation of the previous one and corresponds to the equation
\begin{equation}\label{E:standard_ellipsis}
\dfrac{x^2}{\lambda_{max}} + \dfrac{y^2}{\lambda_{min}} = 1, \qquad {\rm or} \;\; \xx^{\TT}{\Lambda}^{-1}\xx=1,
\end{equation} 
where, $\Lambda$ is the diagonal matrix of the eigenvalues of the correlation matrix. Note that $\xx^{\TT}{\Lambda}^{-1}\xx$ still has a chi-squared distribution with two \gls{DOF}s; moreover $1\cong\chi_2^2(40\%)$. Thus, $\lambda_{max}$ represents the maximum deviation from the mean for the components of ${\bf X}$ in a 60\% confidence region. In our investigation we consider a map, over the computational grid, of the max eigenvalue of the correlation matrix for the bivariate distribution and, for selected \gls{DOF}s, an 80\% confidence region.
%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{remark}
When using the posterior distribution to infer statistical properties of the variable it is important to understand the meaning of confidence regions. In fact, the result of Bayesian inversion is not the real distribution of the state (or related variables); it is rather an estimate of such \gls{PDF} based on two assumptions: the knowledge of the noise distribution and of the prior. For this reason, finding a confidence region based on the posterior does {\it not} correspond to finding the set of values the ``true'' unknown belongs to, but just an estimate of that \cite{somersalo}. In fact, it can happen that a significant part of ``true'' values lies outside a 90\% confidence region (see \cite{somersalo} Section 3.5 for a formal explanation). 
\end{remark}
%%%%%%%%%%%%%%%%%%%%%%%%
\section{Interpolation of data}\label{S:interpolation}
In Chapter \ref{C:da_nse} on the basis of Proposition \ref{P:suff_cond_Qmatrix}, we introduced a form of regularization alternative to Tikhonov techniques. Such result states that, in order to have the warranty of a unique minimizer without using Tikhonov regularization, the selection matrix, restricted to inflow boundary site, should have rank $N_{in}$. Using available data on $\Gamma_{in}$ not satisfying sufficient conditions, we build an interpolating function to be evaluated on \gls{DOF}s on $\Gamma_{in}$, now considered part of the sites set. More specifically, if we assume data $\widetilde{\dd}$ to be given, then, after interpolation, we can write $\dd = \Pi\widetilde{\dd}$, where $\Pi\in \mathbb{R}^{N_s+N_{in},N_s}$ is the interpolation matrix associated with piece-wise linear Lagrangian polynomials (see Section \ref{S:wellp_2} for a 2-dimensional example). 

It is important to note that the interpolation process alters the noise distribution (specifically $\Sigma_{noise}$) introducing additional correlation among the data. Let us assume to know the correlation matrix of $\boldsymbol{\varepsilon}$, say $\Sigma_{\varepsilon}$, this is associated with data $\widetilde{\dd}$ on available sites. After interpolation we write \eqref{E:add_noise} as ${\rm Z}\HH + \boldsymbol{\eta} = {\bf M}$. Here, the random vector $\boldsymbol{\eta} = \Pi \boldsymbol{\varepsilon}$ is associated with a new set of sites: the given ones and the \gls{DOF}s on $\Gamma_{in}$. In the Gaussian assumption, $\boldsymbol{\eta}$ is still normally distributed with mean $\Pi\varepsilon_0$ and correlation matrix $\Sigma_{\eta} = \Pi\Sigma_{\varepsilon}\Pi^{\TT}$, see Section \ref{S:multivariate_properties}. Note that $\Sigma_{\eta}$ is singular by construction; this fact makes $H_{ML}$, in equation \ref{E:hML}, not well-defined and, as anticipates above, Gaussian smoothness priors or modified interpolation matrices should be considered; more details are given in the following section.
%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical results}\label{S:bayes_num_res}
In this section we present numerical results obtained testing the Bayesian formulation on 2-dimensional and axisymmetric 3-dimensional geometries for both the linearized and nonlinear \gls{NSE}. In order to test the effectiveness of this approach we compare the \gls{MAP} and \gls{ML} estimators with deterministic estimates obtained using the formulation of Chapter \ref{C:da_nse}. In particular, we compare accuracy and number of iterations in the solution of linear systems. The index of accuracy is related to the velocity fields retrieved from $H_{MAP}$, $H_{ML}$ and $H_{det}$ (the deterministic estimate); it is defined as $E_{U}=\frac{\|\UU-\UU_{anl}\|_2}{\|\UU_{anl}\|_2}$, where $\UU_{anl}$ is the discretized analytic solution, see equation \eqref{E:exact_ns}, assumed known. We also define an average error over a set of noise realizations, $\overline{E}_{U} = \frac{1}{n}\sum_{i=1}^n\overline{E}_{U,i}$ where $E_{U,i}$ is associated with the $i$-th realization. In addition, we consider a measure of the gain, $\gamma$, in using statistical estimators as opposed to deterministic ones: $\gamma = 1-\frac{\overline{E}_{U, stat}}{\overline{E}_{U,det}}$ where $stat$ stands for either \gls{MAP} or \gls{ML}.

Statistical spread estimators for the velocity and for the \gls{WSS} are derived and possible interpretations are discussed.

We test two possible sites configurations on $\Gamma_{in}$ and we discuss numerical results in terms of accuracy and statistical interpretation. We consider the solution of the linearized problem first; then, we treat the nonlinear case.
%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Implementation details}\label{p:implem_details}We implement the \gls{FE} method with choice of compatible \gls{FE} spaces P1bubble-P1 for velocity and pressure respectively. Numerical results are obtained with the C++ \gls{FE} library {\tt lifeV} and post-processed with {\tt ParaView}. For the solution of the linear system \eqref{E:postH} we refer to the \gls{GMRES}R method \cite{fgmres_gmresr,gmresr}; for the matrix $\rm S$ we consider the \gls{PCD} preconditioner proposed by Kay in \cite{kay} and by Silvester et al. in \cite{silvester}. 

As smoothness prior we use $\Sigma^{-1}_{pr} = \alpha{\rm L}^{\TT}{\rm L}$, $\alpha\geq 0$, where $\rm L$ is the discrete gradient. This is a common choice for \gls{PDE} constrained problems and it is not based on any physical consideration. Nevertheless, we note that in hemodynamics applications some prior knowledge of the velocity and pressure distribution can be derived from the anatomy of the human body; this might be used to infer statistical properties of the control variable $\HH$ using the relation between flow and boundary stress.
%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Data generation}\label{P:data_generation}
We assume to have sites on $\Gamma_{in}$ and in $\dom$ either on selected layers (featuring the configuration of Figure \ref{F:domain}) or uniformly distributed in the domain. At this preliminary stage data are synthetic, i.e. generated adding random noise to a known analytic solution (see equation \eqref{E:exact_ns}). Moreover, we set the available data (the specific realization of $\MB$) $\dd= {\rm Q}\UU_{anl} + \boldsymbol{\varepsilon}$, $\QQ$ being the selection matrix introduced in Chapter \ref{C:da_nse}. Here, the Gaussian random vector $\boldsymbol{\varepsilon}$ is generated coloring white noise: $\boldsymbol{\varepsilon}= {\rm G}\ww$, where $\ww \sim \mathcal{N}(0,\,{\rm I})$ and ${\rm G}$ is such that ${\rm GG}^{\TT} = \Sigma_{\varepsilon}$, the correlation matrix for the noise on sites, assumed known. We compute ${\rm G}$ as the Cholesky factor of $\Sigma_{\varepsilon}$.
%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{The likelihood function}\label{P:posterior}
In the choice of the noise expected value it is reasonable to consider $\varepsilon_0=\mathbf{0}$ (personal communication of Dr. Brummer, Emory CHOA). We assume that the noise correlation depends on the mutual distance between sites. In particular, we assume that close data highly affect each other, and that data significantly far from each other are almost independent. For this reasons, we consider an exponential decay with respect to the mutual distance. Formally,
\begin{equation}\label{E:sigma_noise}
[\Sigma_{noise}]_{ij} = [\Sigma_{\varepsilon}]_{ij} = exp\left\{-\dfrac{1}{l^2}\|\xx_i^m - \xx_j^m\|_2^2 \right\};
\end{equation}
$l$ being a reference distance for the domain of interest. In this way $\Sigma_{noise}$ is a symmetric positive definite matrix. As mentioned in the previous section, when we perform interpolation the correlation matrix is of the form $\Sigma_{noise} = \Sigma_{\eta} = \Pi\Sigma_{\varepsilon}\Pi^\TT$ and, in general, it is singular. 
%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Point estimators}
\paragraph{Test case I}\label{S:test1} In the square domain $\dom = [-0.5,\; 1.5]\times [0,\; 2]$ we consider the analytic solution \eqref{E:exact_ns}. Data on the inflow boundary satisfy the conditions for optimality (specifically, one measurement per grid node, which ensures well-posedness using the \gls{FE} pair P1bubble-P1), \gls{IS} internal data are sparse in $\dom$ and the reference distance is $l$ = 1. 

Table \ref{T:err_cov_comp} reports $E_{U}$ and $\alpha$, fixed a priori, for $H_{det}$, $H_{ML}$ and $H_{MAP}$; with \gls{SNR} set to 20. 
%Due to the ill-conditioned nature of the problem caused by the presence of the noise, the tolerance on the residual for the numerical solver is set to $\delta = 0.005$; this critical value is determined empirically. In fact, in our numerical tests, we observe the typical L-curve behavior \cite{Hansen}, for the approximation error; future work includes an estimate of the optimal $\delta$.
%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{ c | l l | l l l | l l l | }
\cline{2-9}
	&	\multicolumn{2}{|c|}{$H_{det}$} & \multicolumn{3}{|c|}{$H_{ML}$} & \multicolumn{3}{|c|}{$H_{MAP}$}\\ \hline
\multicolumn{1}{|c|}{\gls{IS}}	&	$\overline{E}_{U}$	&	$\alpha$				
							&	$\overline{E}_{U}$	&	$\alpha$		& 	$\gamma$
							&	$\overline{E}_{U}$	&	$\alpha$		&	$\gamma$\\
\hline
\multicolumn{1}{|c|}{0}		&	0.08077	&	0.1		&	0.06290		&	0 & 22\%	&	0.06296		&	0.001	&	22\%\\	
\hline
\multicolumn{1}{|c|}{100}	&	0.03080	&	0		&	0.03033		&	0 & 2\%		&	0.04170		&	0.15	&	-35\%\\	
\hline
\multicolumn{1}{|c|}{200}	&	0.06868	&	0.15	&	0.03751		&	0 & 41\%	&	0.04853		&	0.15	&	29\%\\	
\hline
\multicolumn{1}{|c|}{300}	&	0.08837	&	0.15	&	0.02049		&	0 & 77\%	&	0.02662		&	0.2		&	70\%\\	
\hline
\end{tabular}
\caption{Accuracy results for statistical and deterministic solutions.}
\label{T:err_cov_comp}
\end{center}
\end{table}
Note that the most accurate solution corresponds to the \gls{ML} estimator. In fact, system \eqref{E:hML} is non-singular and we do not make any assumptions on the prior probability. In this case, we can observe that the gain increases as we increase the number of internal data; the more noisy information we have, the better we can get from the statistical solution. Furthermore, the number of iterations required by the \gls{GMRES}R for convergence is approximately 8 for both the approaches.
%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Test case II}\label{S:test2}
In this case, using the same analytical solution and reference distance $l$ = 1, we consider data on the inflow boundary not satisfying sufficient conditions for optimality and additional internal data located on internal slices parallel to $\Gamma_{in}$ (with the purpose of imitating Figure \ref{F:domain}). We want to compare statistical and deterministic estimators and to investigate the role of interpolation.

We first focus on the regularized case {\it without} interpolation and we compare $H_{MAP}$ and $H_{det}$. The regularization parameter $\alpha$ is strongly data dependent; since we want to test the estimators over a set of noise realizations, finding the optimal $\alpha$ would be highly time consuming. For this reason, we compare deterministic and statistical approaches in correspondence of three fixed values of $\alpha$.

In Table \ref{T:hDEThSTAT} (top) we report relative errors $\overline{E}_{U}$, for $H_{det}$ and $H_{MAP}$ in correspondence of $\alpha$ = 0.5, 0.05, 0.005. In this case we consider 300 internal data distributed on 10 slices (30 sites per slice, not in correspondence of grid nodes); also, \gls{SNR} = 20 and 10, the number of noise realizations is $n$ = 20, $\Delta$ = 0.076, and the number of sites on $\Gamma_{in}$ is $N_{s,in}$ = 14. 
%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline
\gls{SNR}		&	$\alpha$		&	$\overline{E}_{U,det}$	&	$\overline{E}_{U,MAP}$	&	$\gamma$\\	\hline
20		&	0.5		&	0.0665						&	0.0530						&	24\%	\\	\hline
20		&	0.05	&	0.0666						&	0.0550						&	17\%	\\	\hline
20		&	0.005	&	0.0706						&	0.0579						&	18\%	\\	\hline
10		&	0.5		&	0.1272						&	0.0946						&	26\%	\\	\hline
10		&	0.05	&	0.1514						&	0.1032						&	32\%	\\	\hline
10		&	0.005	&	0.1256						&	0.1059						&	28\%	\\	\hline
\end{tabular}

\vspace{.4cm}

\begin{tabular}{|c|c|c|c|}
\hline
\gls{SNR}		&	$\overline{E}_{U,det}$	&	$\overline{E}_{U,ML},\; (mod)$	&	$\gamma$\\	\hline
20		&	0.0709						&	0.0552								&	22\%	\\	\hline
10		&	0.1518						&	0.1256								&	17\%	\\	\hline
\end{tabular}
\caption{Accuracy results for statistical and deterministic solutions for the linearized NSE.}
\label{T:hDEThSTAT}
\end{center}
\end{table}

From these results we infer the following facts. Compared to the deterministic estimator, $H_{MAP}$ is always more accurate since it takes into account additional information brought by statistical properties of the data. We stress that, since synthetic data are generated using the matrix $\rm G$, the noise correlation matrix is exact. 

In average the number of \gls{GMRES} iterations is higher for $H_{MAP}$ (1.3 times bigger than the one required by $H_{det}$); this is due to the presence of $\Sigma_{noise}^{-1}$ in the formulation. 

As mentioned above, since using regularization is computationally expensive, we also want to compare $H_{ML}$ and $H_{det}$ without the additional regularizing term (i.e. $\alpha$ = 0 for both of them) and using interpolation. We recall that when interpolating the induced correlation matrix $\Sigma_{\eta}$ is singular by construction making the problem potentially ill-posed. In our numerical results we consider two options: the first consists in solving the correlation matrix $\Sigma_{\eta}$ in a least-squares sense, i.e. $\Sigma_{\eta}^\TT\Sigma_{\eta} \xx = \Sigma_{\eta}^\TT\yy$ where $\xx$ is the unknown and $\yy$ is a right hand side. The second option consists in using a modified correlation matrix constructed with the assumption of independent and identically distributed interpolated data. This is equivalent to setting the diagonal block of $\Sigma_{\eta}$ corresponding to \gls{DOF}s on $\Gamma_{in}$ equal to $\sigma^2{\rm I}$, where $\rm I$ is the identity matrix of dimension $N_{in}$, which avoids singularity. Since the statistical information associated with interpolated data is redundant, the modification to the correlation matrix does not compromise the accuracy of the statistical estimator; in fact, the statistical surplus value is brought only by statistical properties of measurements. In Table \ref{T:hDEThSTAT} (bottom) we report results obtained in correspondence of the same 20 noise realizations used in the regularized case.
%As expected, in correspondence of $\Sigma_{\eta}$, results are less accurate; also, a higher number of iterations is required. 
Using the modified matrix $\Sigma_{\eta}$, system \eqref{E:hML} is well-posed and the additional knowledge brought by the noise correlation matrix makes the statistical estimator more accurate than $H_{det}$. On the other hand, the number of iterations in solving for $H_{ML}$ is 1.5 times bigger than the one required by $H_{det}$. As mentioned above, the higher number of iterations is expected and it is due to the presence of $\Sigma_{noise}^{-1}$; nevertheless, we believe that the gain in accuracy makes the computational effort worthwhile.
%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Test case III}\label{S:test3}
We present numerical results obtained solving the \gls{NSE} using the same analytic solution and reference distance of test case I and II. We consider data on $\Gamma_{in}$ not satisfying sufficient conditions for optimality and internal data located on 10 internal slices (with 30 data/slice). The numerical problem is solved using algorithm \eqref{E:algorithm}. In Table \ref{T:hDEThSTAT_NS} we report results obtained in correspondence of $\Delta$ = 0.076 and \gls{SNR} = 20 and 10. In the computation of $H_{MAP}$ the regularization parameter is chosen relying on results obtained in the linearized case; we choose $\alpha$ = 0.5 since it corresponds to the most accurate numerical solution. The poor gain in correspondence of \gls{SNR} = 20 means that statistical information associated with a low amount of noise is not significant enough to make a considerable difference with respect to deterministic estimates in terms of accuracy.
%nevertheless, we still get a more accurate result than the deterministic estimator. 
In the other cases statistical estimators are significantly more accurate and yield a gain up to the 25\% in correspondence of low values of \gls{SNR}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
\gls{SNR}		&	$\overline{E}_{U, det}$	&	$\overline{E}_{U, MAP}$	&	$\gamma$	\\	\hline
20		&	0.0822							&	0.07371							&	10\%	\\	\hline
10		&	0.1394							&	0.1041							&	25\%	\\	\hline
\end{tabular}
\hspace{.4cm}
\begin{tabular}{|c|c|c|c|}
\hline
\gls{SNR}		&	$\overline{E}_{U, det}$	&	$\overline{E}_{U, ML}$	&	$\gamma$	\\	\hline
20		&	0.0855							&	0.0579							&	6\%		\\	\hline
10		&	0.1675							&	0.1363							&	18\%	\\	\hline
\end{tabular}
\caption{Accuracy results for statistical and deterministic solutions for the NSE.}
\label{T:hDEThSTAT_NS}
\end{center}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Test case IV}\label{S:test4}
We consider an axisymmetric formulation for the numerical solution of the 3-dimensional Poiseuille flow in a cylinder of length $L$ =6, diameter $D$ =3 and viscosity $\nu$ = 1 in a 2-dimensional structured grid of dimension 80$\times$20. Data on the inflow boundary do not satisfy sufficient conditions and internal data are located on 5 internal slices with 10 data per slice. Here, the reference distance is $l$ = 1.5. We perform interpolation and, in order to prevent potential ill-conditioning, we add a regularization term (or smoothness prior) with $\alpha$ = 1e-7. Results for \gls{SNR} = 20 and 10 are reported in Table \ref{T:hDEThSTAT_axi}; in both cases we have a significant gain in accuracy. Also, the number of iterations, in average, is the same.
%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
\gls{SNR}		&	$\overline{E}_{U, det}$	&	$\overline{E}_{U, MAP}$	&	$\gamma$\\	\hline
20		&	0.0396						&	0.0308						&	22\%	\\	\hline
10		&	0.1423						&	0.0978						&	31\%	\\	\hline
\end{tabular}
\caption{Accuracy results for statistical and deterministic solutions for the axisymmetric case.}
\label{T:hDEThSTAT_axi}
\end{center}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
\multicolumn{4}{|c|}{test case I}								\\  \hline
\gls{SNR}		&	$\alpha$		&	$E_{U, MAP}$		&	$\lambda_{MAX}$	\\	\hline
20		&	0.5				&	0.06003				&	0.09759			\\	
10		&	0.5				&	0.07754				&	0.1744			\\	\hline
\end{tabular}
\hspace{.2cm}
\begin{tabular}{|c|c|c|c|}
\hline
\multicolumn{4}{|c|}{test case II}								\\  \hline
\gls{SNR}		&	$\alpha$		&	$E_{U, MAP}$		&	$\lambda_{MAX}$	\\	\hline
20		&	1.e-5			&	0.01720				&	0.3759			\\	
10		&	1.e-5			&	0.02565				&	0.7024			\\	\hline
\end{tabular}
\caption{Accuracy results and maximum deviation for different values of SNR.}
\label{T:lambdaMAX}
\end{center}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.41\textwidth]{sd_SQUARE_20.png} &
\includegraphics[width=0.5\textwidth]{node1_20_SQUARE.pdf} \\
\includegraphics[width=0.5\textwidth]{node2_20_SQUARE.pdf} &
\includegraphics[width=0.5\textwidth]{node3_20_SQUARE.pdf}
\end{tabular}
\end{center}
\caption{Standard deviation map (top left) and ellipses corresponding to node 1 (top right), 2 (bottom left) and 3 (bottom right).}
\label{F:ellip_square}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Spread estimators}
\paragraph{Test case I} We consider the square geometry $[0.5,\;1.5]\times[0,\;2]$ and the analytic solution \eqref{E:exact_ns}; data on $\Gamma_{in}$ satisfy conditions for optimality in Proposition \ref{P:suff_cond_Qmatrix}. For this test case the computational grid is sufficiently small to allow explicit computation and storage of correlation matrices; matrix operations are performed in {\tt Matlab}. In particular, for the computation of $\Sigma_{post}^{-1} = \Sigma_{pr} + {\rm Z}^{\TT}\Sigma_{noise}^{-1}{\rm Z}$ and $\Sigma_U = {\TT}\Sigma_{post}^{-1}{\TT}^{\TT}$ we rely on the conjugate gradient method for matrix inversion. For each sub matrix of $\Sigma_U$ defined as in \eqref{E:u_distribution}, corresponding to horizontal and vertical velocities on one \gls{DOF}, say the $j$-th \gls{DOF}, we compute the square root of the maximum eigenvalue, $\lambda_{max,\,j}$; as described in Section \ref{S:spread_est} this is the maximum deviation from the mean in a 60\% confidence region. In Figure \ref{F:ellip_square} (top left) we report the map of $\sqrt{\lambda_{max,\,j}}$ over the computational grid for \gls{SNR} = 20. We observe a peak on the inflow boundary, this is due to the fact that the noise of data on $\Gamma_{in}$ has high influence on the solution, the result is a higher uncertainty for the velocity components. In correspondence of three areas featuring different values of deviation we compute a 80\% confidence region corresponding to the set 
\begin{equation}
(\UU_j - U_j)^{\TT}\Sigma_U^{-1}(\UU_j - U_j)\leq \chi_2^2(20\%) \cong 3.219,
\end{equation}
where $\UU_j$ corresponds to the vector of horizontal and vertical velocities on the $j$-th \gls{DOF} and $U_j$ is the corresponding mean value. 
Such ellipses are reported in Figure \ref{F:ellip_square}. Results for \gls{SNR} = 20 and \gls{SNR} = 10 are reported in Table \ref{T:lambdaMAX} (left). Here $\lambda_{MAX} = max_{j} (\lambda_{max,\,j})$ is the maximum deviation over the computational grid. 

Far from being realistic, this test case highlights the fact that Bayesian inversion is capable of reducing the standard deviation of the noise affecting observed measures, 0.1467 for \gls{SNR} = 20 and 0.2953 for \gls{SNR} = 10.
%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.95\textwidth]{sd_AXI_20.png}
\end{center}
\caption{Standard deviation map.}
\label{F:nodes_axi}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.48\textwidth]{node1_20_AXI.pdf} &
\includegraphics[width=0.46\textwidth]{node2_20_AXI.pdf}
\end{tabular}
\end{center}
\caption{Ellipses corresponding to node 1 (left) and 2 (right).}
\label{F:ellip_axi}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Test case II} We consider the 3-dimensional axisymmetric formulation and data on $\Gamma_{in}$ not satisfying conditions for optimality; thus, we perform interpolation. Recall that, in this case, $\Sigma_{noise} = \Pi\Sigma_{\varepsilon} \Pi^{\TT}$, where $\Sigma_{\varepsilon}$ is the correlation matrix for the measures. Matrices operations are treated as in the previous test case. In Figure \ref{F:nodes_axi} we report the map of $\sqrt{\lambda_{max,\,j}}$ in each \gls{DOF} for \gls{SNR} = 20. Here, the peak is located in the lower left corner. This result is reasonable and it is related to boundary conditions. On the wall and on the outflow section we prescribe exact boundary data; for this reason, we do not have uncertainty. As we move far from such boundaries, i.e. towards the axis of symmetry and the inflow section we have a higher deviation from the mean. This behavior is enhanced by the fact that the boundary condition on $\Gamma_{in}$ is the control variable which is highly affected by the measurement noise. Nevertheless, we observe that except for the bottom left corner the deviation is lower that the one of the data, 0.325 for \gls{SNR} = 20 and 0.65 for \gls{SNR} = 10. As in Test case I, we compute the 80\% confidence region in areas featuring different behavior. Results for \gls{SNR} = 20 are reported in Figure \ref{F:ellip_axi}; note that in this case we have a horizontal flow, this explains the low deviation from the mean for vertical components of the velocity.
%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.5\textwidth]{carotid_data.png} &
\includegraphics[width=0.5\textwidth]{carotid_bayesian.png}
\end{tabular}
\end{center}
\caption{Data configuration and magnitude of the velocity field.}
\label{F:carotid_data_velocity}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline
\gls{SNR}		&	$n$		&	$\overline{E}_{U, det}$	&	$\overline{E}_{U, ML}$	& $\gamma$	\\	\hline
20				&	20		&	0.05273					&	0.03617					& 31\%		\\	\hline
\end{tabular}
\caption{Comparison results between deterministic and ML estimator: accuracy with respect to the reference solution.}
\label{T:carotid_bayes}
\end{center}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.5\textwidth]{carotid_FE_ML.png} &
\includegraphics[width=0.5 \textwidth]{carotid_FE_DET.png}
\end{tabular}
\end{center}
\caption{Comparison results between deterministic and ML estimator. Map of the difference between statistical and reference solution, and between deterministic and reference solution}
\label{F:carotid_U_comparison}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.5\textwidth]{carotid_bayesian_estimator.png}
\end{center}
\caption{Standard deviation map.}
\label{F:carotid_lambda}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Towards real geometries}\label{S:bayesian_real_geom}
In this section we consider a more complex geometry as done in Chapter \ref{C:da_nse}, Section \ref{S:da_nse_real_geometries}. We reproduce a 2-dimensional section of the carotid. We compare statistical and deterministic estimators and we compute confidence intervals for the velocity and for the \gls{WSS}.

Data are located as reported in Figure \ref{F:carotid_data_velocity} (left), internal data lie on four layers in correspondence of grid nodes. These are generated adding Gaussian noise to a reference solution $\UU_{\rm FE}$, the numerical solution obtained solving, with \gls{FE}, the forward problem on a very fine grid (see Section \ref{S:da_nse_real_geometries}, test case I, for implementation details). The correlation matrix for the Gaussian noise is defined as in the previous test cases, i.e. as in equation \eqref{E:sigma_noise}. The indexes of accuracy are the average relative error $\overline{E}_{U} = \frac{1}{n}\sum_{i=1}^n E_{U,i}$ where $E_{U,i}$ is associated with the $i$-th realization and it is defined as $E_{U}=\frac{\|\UU-\UU_{\rm FE}\|_2}{\|\UU_{\rm FE}\|_2}$, and the gain $\gamma = 1-\frac{\overline{E}_{U, stat}}{\overline{E}_{U,det}}$.

In Table \ref{T:carotid_bayes} we report comparison results between statistical and deterministic estimators; since data satisfy conditions for optimality we set the regularization parameter $\alpha = 0$, which corresponds to computing the \gls{ML} estimator in the statistical formulation. These results show that also for more complex geometries statistical estimators are robust with respect to the noise and they are significantly more accurate than the deterministic ones. As en example, in Figure \ref{F:carotid_data_velocity} we report the magnitude of the assimilated velocity field; in Figure \ref{F:carotid_U_comparison} the difference between statistical and reference solution, and between deterministic and reference solution. It is interesting to note that far from the inflow boundary the two estimators are very similar, while close to $\Gamma_{in}$ the statistical solution is more accurate. This shows that the \gls{ML} estimator is more robust with respect to the noise that significantly affects the numerical solution on $\Gamma_{in}$.

In view of the application to real cardiovascular problems we want to estimate the confidence intervals for the velocity. In Figure \ref{F:carotid_lambda} we report the map of $\sqrt{\lambda_{max}}$, see Section \ref{S:spread_est}, over the computational grid; consistently with previous results we observe that the standard deviation increases far from the wall since we work in a rigid assumption and we prescribe homogeneous Dirichlet conditions on the walls. Furthermore, since close to $\Gamma_{in}$ the solution is more sensitive to the noise in this area we have the highest deviation from the mean.
%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{The wall shear stress}
In medical practice other than estimating the flow and assessing its statistical properties it is fundamental to estimate other flow related variables as the \gls{WSS} since, as a matter of fact, they are crucial for predicting the occurrence of diseases or even the rupture of a blood vessel. At the moment \cite{brummer} the \gls{WSS} is computed using the observed velocity data (either projected onto the whole grid or collected on each node); these estimators are not accurate enough and are more sensitive to the presence of the noise. Using point and statistical estimators for the velocity we compute statistical parameters of the \gls{WSS} distributions; also in this case we assume this random vector to be normally distributed, ${\bf WSS} \sim \mathcal{N}(WSS,\,\Sigma_{WSS})$, mean and correlation matrix are determined by the PDF of the velocity. In fact, velocity and WSS are related by the linear transformation
\begin{equation}
{\bf WSS} = \TT_w{\UU};
\end{equation}
where $\TT_w\in \mathbb{R}^{N_u,N_w}$ maps the discretized velocity into the discretized \gls{WSS}; $N_w$ being the number of \gls{DOF}s of the \gls{WSS}.
%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.5\textwidth]{carotid_wss_bayesian.png} & 
\includegraphics[width=0.5\textwidth]{wss_magn20.pdf} 
\end{tabular}
\end{center}
\caption{WSS on a selected wall; on the left the vector field, on the right the magnitude as a function of the $x$-coordinate.}
\label{F:bayes_WSS}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%
According to property {\bf P1} in section \ref{S:multivariate_properties}, $\bf WSS$ has a Gaussian distribution with mean $WSS = \TT_w U$ and covariance $\Sigma_w = \TT_w\Sigma_U\TT_w^{\TT}$. Furthermore, according to {\bf P2} the vector of horizontal and vertical components of the \gls{WSS} on each \gls{DOF} is distributed as follows 
\begin{equation}\label{E:w_distribution}
\left[\begin{array}{l}
{\bf WSS}_i\\
{\bf WSS}_{i+N_w/2}
\end{array}
\right]\sim \mathcal{N}\left(
\left[\begin{array}{l}
WSS_i\\
WSS_{i+N_w/2}
\end{array}
\right], \; 
\left[\begin{array}{ll}
\Sigma_{WSS;\,i,i} & \Sigma_{WSS;\,i,i+N_u/2}\\
\Sigma_{WSS;\,i+N_u/2,i} & \Sigma_{WSS;\,i+N_u/2,i+N_u/2}
\end{array}\right]\right).
\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.5\textwidth]{wss_sd20PRES.pdf} &
\includegraphics[width=0.5 \textwidth]{wss_sd_all_data20PRES.pdf}
\end{tabular}
\end{center}
\caption{Standard deviation of the assimilated and non-assimilated WSS as a function of the $x$-coordinate.}
\label{F:wss_sd}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%
In Figure \ref{F:bayes_WSS} (left) we report the \gls{WSS} computed on the left inner wall of the bifurcation; the relative error, in norm $l_2$ with respect to the reference solution is of the order of the 15\% with \gls{SNR} = 20. On the right, the magnitude of the vector in each \gls{DOF} of the selected boundary is displayed as a function of the horizontal coordinate.
As we did for the velocity, we can compute in each \gls{DOF} $\sqrt{\lambda_{max}}$, where $\lambda_{max}$ is the maximum eigenvalue of the 2-dimensional correlation matrix in \eqref{E:w_distribution}, this is a measure of the maximum variation from the mean in a confidence region of the 60\%. In Figure \ref{F:wss_sd} (left) we report on such selected wall the value of $\sqrt{\lambda_{max}}$ for each \gls{DOF}; on the right we report the same quantity for a ``non-assimilated'' \gls{WSS}. More specifically we assume ${\bf WSS} = \TT_w \UU_{obs}$, where $\UU_{obs}$ is the vector of observed velocities. Thus, it features the correlation matrix 
\begin{equation}
[\Sigma_{obs}]_{ij} = exp\left\{-\dfrac{1}{l^2}\|\xx_i^m - \xx_j^m\|_2^2 \right\};
\end{equation} 
where $\xx_i^m$, $i$ = 1, ... $N_u$, are observation sites (here corresponding to \gls{DOF}s of the velocity). According to these relations, ${\bf WSS} \sim \mathcal{N}(\TT_w U_{obs},\,\Sigma_{obs})$. From the results we note that the maximum deviation from the mean is much higher than the assimilated one; note that the flat configuration of the non-assimilated \gls{WSS} in the interval $(-6,\,-1)$ is due to the grid-structure, very regular in correspondence of that area.  Even if the result is preliminary and depends on the data configuration, the huge reduction in the standard deviation obtained using \gls{DA}, lead us to conclude that we can effectively filter the noise.

% ---------------------------------------------------------------------------
%: ----------------------- end of thesis sub-document ------------------------
% ---------------------------------------------------------------------------

