text:
\section{Introduction}
In most system identification problems, the input signal---that is, the
independent variable---is perfectly known~\cite{ljung1999system}. Often, the
input signal is the result of an identification experiment, where a signal with
certain characteristics is designed and applied to the system to measure its
response. However, in some applications, the hypothesis that the input signal is
known may be too restrictive. In this work, we propose a new model structure
that accounts for partial knowledge about the input signal and we show how many
classical system identification problems can be seen as problems of identifying
instances of this model structure.

The proposed model structure, which we call \emph{uncertain-input model}, is composed of a linear
time-invariant dynamical system (the \emph{linear system}) and of a signal
 of which partial information is available (the \emph{unknown input}). In the
next section, we characterize formally the unknown input; before that, we give
some examples of classical models that can be seen as uncertain-input models.
The \emph{Hammerstein model} is a cascade composition of a static nonlinear
function followed by a linear time-invariant dynamical
system~\cite{bai1998optimal,giri2010block,risuleo2015kernel}. In the Hammerstein
model, the (perfectly known) input signal passes through the unknown static
nonlinear function. After the nonlinear transformation, the signal which
is fed to the linear system, is completely unknown.
However, some characteristics of the signal may be known; for instance, we may
known that the nonlinear function is smooth or we may have a set of candidate
basis functions to choose among. Another instance of a
model where the input is not perfectly known is the \emph{errors-in-variables}
model~\cite{soederstroem2010system}. In the errors-in-variables formulation, the
input in known up to noisy measurements. The noise in the input introduces many difficulties and
special techniques have been developed to deal with
it~\cite{soederstroem2003why,soederstroem2007errors}. Closely related to
errors-in-variables models, \emph{blind system identification} methods are
used when the input signal is completely unknown~\cite{abedmeraim1997blind}.
These are particularly useful in telecommunications, image reconstruction, and
biomedical
applications~\cite{nakajima1993blind,moulines1995subspace,mccombie2005laguerre}.
Blind system identification problems are generally ill posed, and certain
assumptions on the input signal are needed to recover a
solution~\cite{ahmed2014blind}. Similar to blind problems are the problems of
\emph{system identification with missing data}. In these cases, the missing
data are estimated, together with a description of the system, by making
hypotheses on the mechanism that generated the missing
data~\cite{wallin2014maximum,markovsky2013structured,pillonetto2009bayesian,risuleo2016kernel,linder2017identification}.

In all the applications we have outlined, we can identify the common thread of
a linear system fed by a signal about which we have limited \emph{prior
information}. This leads us naturally to consider a Bayesian
framework where we can use prior distributions to encode beliefs about the
unknown quantities~\cite[Section~2.4]{bernardo2000bayesian}. Within the vast
framework of Bayesian methods, we concentrate on Gaussian
processes~\cite{rasmussen2006gaussian}. These enable us to
compute many quantities in closed form and to reason about identification in
terms of a limited number of sufficient statistics. For these reasons,
Gaussian-process modeling has become a popular approach in system
identification~\cite{pillonetto2011new,frigola2014identification,svensson2017flexible}. Although Gaussian processes
are typically analytically convenient, the structure of the uncertain-input problem leads
to an intractable inference problem: even though we model the system and the
input as Gaussian processes, the output of the system depends on their
convolution and therefore does not admit a Gaussian description. To perform the
inference---that is find the posterior distribution of the unknowns given the
observations---we need approximation methods. We propose two different
approximation methods for the posterior distribution of the unknowns: one Markov
Chain Monte Carlo~(MCMC, see~\cite{gilks1996markov}) method and one variational
 approximation~\cite{beal2003variational} method. In the MCMC
 method, we use the Gibbs sampler~\cite{geman1984stochastic} to draw
 particles from the posterior distribution and we approximate expectations as
 averages computed with the particles. In the variational method, we find the
 factorized distribution that best approximates the posterior distribution in
 Kullback-Leibler distance.

To give flexibility to the model, we allow the Gaussian priors to depend on
certain parameters (called \emph{hyperparameters}) that need to be estimated
from data together with the measurement noise variances. To estimate these
parameters, we use the \emph{empirical Bayes method} which requires maximizing
the marginal distribution of the data (sometimes called \emph{evidence},
see~\cite{maritz1989empirical}). To
this end, we propose an iterative algorithm based on the
Expectation-Maximization
(EM) method~\cite{dempster1977maximum}. The EM method alternates between the
computation of the expected value of the joint likelihood of the data, of the
unknown system, and of the input (E-step), and the maximization of this expected
value with respect to the unknown parameters (M-step). We show that the E-step
can be computed using the same approximations of the posterior distributions
that are used in the inference and that the M-step consists in a series of
simple and independent optimization problems that can be solved easily.

As mentioned above, the uncertain-input model encompasses several classical
model structures that have been object of research in the system-identification
community for decades. Two important contributions of this work are as follows.
\begin{enumerate}
  \item We unify the problems of identifying systems that are usually regarded
    as belonging to different model classes into a single identification
    framework.
  \item We formalize a method to apply the new tools of Gaussian processes and
    Bayesian inference to classical system identification problems.
\end{enumerate}

To support the validity of the proposed methods, we present identification
experiments on synthetic datasets of cascaded linear systems and of
Hammerstein systems.

\subsection{Notation}
The notation $[A]{}_{i,j}$ indicates the element
of matrix $A$ in position $i,j$ (single subscripts are used for vectors).
``$\mathrm{T}_{N\times n}(v)$'' denotes the $N$ by
$n$ lower-triangular Toeplitz matrix of the $m$ dimensional vector $v$:
\begin{equation}
\sbr{\mathrm{T}_{N\times n}(v)}_{i,j} = \begin{cases} v_{i-j+1} & 0\leq i-j+1\leq m\\0
                                                        & \text{otherwise}\end{cases}
\end{equation}
If $v$ is a vector, then $V$ is the $N$ by $N$ Toeplitz matrix whose
elements are given by $v$. The notation ``$\|a\|^2_{M}$'' is shorthand for $a^T Ma$. The
notation ``$\mathcal{N}(\alpha,\Sigma)$'' indicates the Gaussian distribution with mean
vector $\alpha$ and covariance matrix $\Sigma$. The notation
``$\mathcal{GP}(\mu,\Sigma)$'' indicate a Gaussian process with mean function
$\mu$ and covariance function $\Sigma$. Random variables and their realizations
have the same symbol. The notation ``$x;\theta$'' indicates that the random variable
$x$ depends on the parameter $\theta$. If $x$ is a random variable, $\mathrm{p}(x)$
denotes its density. The symbol ``$\cong$'' indicates equality up to an
additive constant and ``$\delta$'' is the Dirac density.

\section{Uncertain-input systems}\label{sec:ui_systems}
In this work, we propose a new model structure called the \emph{uncertain-input
model}. Consider the block scheme in Figure~\ref{fig:ui_block_scheme}. Many system
identification tasks can be formulated as the identification of a linear system $S$,
subject to an input sequence $\{w_t\}$. In this work, we consider problems in
which we have partial information about the input sequence, and this partial
information depends on the specific problem at hand.

\begin{figure}[htb]
  \centering
  \includegraphics{blockscheme.pdf}
  \caption{A block scheme of the general uncertain input system.}\label{fig:ui_block_scheme}
\end{figure}

We assume that the linear system $S$ is time invariant, stable, and causal. Therefore,
it is uniquely described by the sequence $\{g_t\}$ of its impulse response
samples, and the output of the system generated by an input $\{w_t\}$ can be
represented as the discrete convolution of the system impulse response with the
input signal---that is, at time $t$, the measurements of the output can be
written as the noise-corrupted discrete convolution
\begin{equation}\label{eq:model_y}
  y_t = \del{w\ast g}_t + \varepsilon_t,
\end{equation}
where $\{\varepsilon_t\}$ is a stochastic process that describes additive
measurement noise, and where ``$\ast$'' denotes the discrete time convolution
\begin{equation}\label{eq:convolution}
  \del{w\ast g}_t = \sum_{k=1}^\infty g_{k}w_{t-k}.
\end{equation}

In the uncertain-input model, we consider that the input signal is measured with
additive white noise described by a stochastic processes $\{\eta_t\}$. This
assumption allows us to write, for the input measurements, the model
\begin{equation}\label{eq:model_v}
  v_t = w_t + \eta_t.
\end{equation}

We assume that the noise processes $\{\eta_t\}$ and $\{\varepsilon_t\}$ are independent
Gaussian white-noise processes. This means that every noise sample has a
Gaussian distribution,
\begin{equation}\label{eq:model_noise}
    \eta_t \sim \mathcal{N}(0,\sigma_v^2),\quad
    \varepsilon_t \sim \mathcal{N}(0,\sigma_y^2),
\end{equation}
and that $\varepsilon_t$ is independent of $\varepsilon_s$, for
$s\neq t$, and of $\eta_s$ for any $s$. To allow for models where some
observations are missing, we assume infinite variance for those noise
components that correspond to the missing samples.

To encode the prior information we have about the input signal and about the linear
system, we use Gaussian process models. We model the unknown input signal
and the impulse response of the linear system as a realization of a joint Gaussian
processes with suitable mean and covariance functions,
\begin{equation}\label{eq:gp_joint}
\sbr[4]{\begin{matrix}w\\g\end{matrix}}\!\!\sim \!\mathcal{GP} \!\del{\!\!
    \sbr[4]{\begin{matrix}\mu_w(\cdot;\theta)\\
    \mu_g(\cdot;\rho)\end{matrix}}\!\!,\!
    \sbr[4]{\begin{matrix}K_w(\cdot;\theta) &
        {K_{gw}(\cdot,\cdot;\rho,\theta)}^T\\
      K_{gw}(\cdot,\cdot;\rho,\theta) &
      K_{g}(\cdot,\cdot;\rho)
  \end{matrix}}
  }\!.
\end{equation}
The mean functions of the Gaussian processes, $\mu_g(\cdot\,;\theta)$ and
$\mu_w(\cdot\,;\rho)$, may depend on the parameter vectors $\theta$ and
$\rho$, called \emph{hyperparameter vectors}, which can be used to shape the prior
information to the specific application. The same goes for the covariance
functions $K_w(\,\cdot\,,\,\cdot\,;\theta)$,
$K_g(\,\cdot\,,\,\cdot\,;\rho)$, and $K_{gw}(\,\cdot\,,\,\cdot\,;\rho,\theta)$
which may depend on (possibly different) hyperparameters.

For notational convenience, we present the explicit computations in the case of
independent Gaussian process models for $g$ and $w$---that is, we consider the
case where
\begin{equation}\label{eq:gp_models}
  \begin{split}
    w &\sim
    \mathcal{GP}\big(\mu_w(\,\cdot\,\,;\theta),K_w(\,\cdot\,,\,\cdot\,;\theta)\big),\\
    g &\sim \mathcal{GP}\big(\mu_g(\,\cdot\,;\rho),K_g(\,\cdot\,,\,\cdot\,;\rho)\big),
  \end{split}
\end{equation}
and the cross-covariance of processes
is zero. However, all results we show hold also in the more general case.

We assume that we have collected $N$ measurements of the processes $\{v_t\}$
and $\{y_t\}$ and, for sake of simplicity, we also assume that $w_t=0$ for
$t<0$ (see~\cite{risuleo2015estimation} for a way to extend the proposed
framework to unknown initial conditions).
From~\eqref{eq:model_y}, we see that the output measurements only depend on the
values of the impulse response at the discrete time instants $t=1,2,\ldots,N$; therefore,
we can consider the joint distribution of the samples $g_t$ for $t=1,2,\ldots,N$.
From the Gaussian process model~\eqref{eq:gp_models}, we have that, if we
collect the samples of $\{g_t\}$ into an $N$-dimensional column vector $g$,
this vector has a joint Gaussian distribution given by
\begin{equation}\label{eq:model_g}
    g \sim \mathcal{N}\big(\mu_g(\rho),K_g(\rho)\big),
\end{equation}
where we have defined the mean vector and the
covariance matrix induced by~\eqref{eq:gp_models} as
\begin{equation}
  {\big[\mu_g(\rho)\big]}_j := \mu_g(j\,;\rho), \quad {\big[K_g(\rho)\big]}_{i,j} :=
  K_g(i,j\,;\rho).
\end{equation}

From~\eqref{eq:model_v} and~\eqref{eq:model_y}, we have that the $N$
measurements of the input and output only depend on the samples $w_t$ for
$t=1,\ldots,N$; therefore, we can consider the joint distribution of these
samples, collected in an $N$-dimensional vector $w$.
This distribution is Gaussian, and it is given by
\begin{equation}\label{eq:model_w}
    w \sim \mathcal{N}\big(\mu_w(\theta),K_w(\theta)\big),
\end{equation}
where we have defined the mean vector and the
covariance matrix induced by~\eqref{eq:gp_models} as
\begin{equation}
  {\big[\mu_w(\theta)\big]}_j := \mu_w(j\,;\theta), \quad {\big[K_w(\theta)\big]}_{i,j} :=
  K_w(i,j\,;\theta).
\end{equation}

Assembling the models for the different components, given
by~\eqref{eq:model_y},~\eqref{eq:model_v},~\eqref{eq:model_noise},~\eqref{eq:model_g},
and~\eqref{eq:model_w}, we arrive at the following definition
of the uncertain-input model:
\begin{equation}\label{eq:ui_model}
  \left\lbrace
  \begin{aligned}
    y &= Wg + \varepsilon\\
    v &= w + \eta\\
    g & \sim \mathcal{N}\big(\mu_g(\rho),K_g(\rho)\big)\\
    w & \sim \mathcal{N}\big(\mu_w(\theta),K_w(\theta)\big)\\
    \varepsilon & \sim \mathcal{N}\big(0,\sigma_y^2 I_N\big)\\
    \eta & \sim \mathcal{N}\big(0,\sigma_v^2 I_N\big)\\
    g,\,&w,\,\varepsilon,\,\eta\text{ mutually independent}
  \end{aligned}
  \right.
\end{equation}
where we have collected the output measurements $\{y_t\}$ in a vector $y$ and
where $\varepsilon$ and $\eta$ are the vectors of the first $N$ input and output
noise samples.  The matrix
$W$ is the $N\times N$ Toeplitz matrix of the input,
$W := \mathrm{T}_{N\times N}(w)$,
 which represents the discrete-time convolution~\eqref{eq:convolution} as the
product $Wg$. If we define the $N\times N$ Toeplitz matrix of the impulse
response samples, $G := \mathrm{T}_{N\times N}(g)$, then we have the property
\begin{equation}\label{eq:toeplitz_property}
  Wg = Gw.
\end{equation}
In the next section, we give examples of some classical system identification
problems that can be cast as uncertain-input identification problems.

\section{Examples of uncertain-input models}
The uncertain-input framework is a generalization of many classical
system-identification problems. All these classical problems can be analyzed
using the tools of uncertain-input models; furthermore, under the right
conditions, the identification approach that we propose for uncertain-input
models reduces to classical system-identification approaches.
\subsection{Linear predictor model}
Consider the \emph{output-error} transfer-function model~\cite{ljung1999system},
\begin{equation}
  y_t = \frac{B(q;\rho)}{F(q;\rho)} u_t + \varepsilon_t,
\end{equation}
where $B(q;\rho)$ and $F(q;\rho)$ are polynomials in the one-step shift operator
$q$ and $\varepsilon_t$ is Gaussian white noise. If we consider the parametric
predictor of the output-error model, we can write it as
\begin{equation}
  \hat y_{t|t-1} = \frac{B(q;\rho)}{F(q;\rho)} u_t = \del{g(\rho)\ast u}_t,
\end{equation}
where $g_t(\rho)$ is the impulse response of the predictor transfer function.
We can see this model as a degenerate uncertain-input model with
\begin{equation}
  \begin{aligned}
    {[\mu_g(\rho)]}_i &= g_i(\rho), & {[K_g(\rho)]}_{i,j} &= 0,\\
    {[\mu_w(\theta)]}_i &= u_i, & {[K_w(\theta)]}_{i,j} &= 0.
  \end{aligned}
\end{equation}
We can also incorporate the framework of Bayesian identification of finite
impulse-response models with first order stable-spline kernels (for a survey,
see~\cite{pillonetto2014kernel}) with the choice
\begin{equation}\label{eq:stable-spline}
  \begin{aligned}
    {[\mu_g(\rho)]}_i &= 0, & {[K_g(\rho)]}_{i,j} &= \rho_1\,\rho_2^{\max(i,j)},\\
    {[\mu_w(\theta)]}_i &= u_i, & {[K_w(\theta)]}_{i,j} &= 0,
  \end{aligned}
\end{equation}
where $\rho_1\geq 0$ is a scaling parameter and $\rho_2\in \intcc{0,1}$
regulates the decay rate of $g$ (see,~\cite{pillonetto2010new}).
Note that, in this formulation, any kernel can be used to model $g$ (see, for
instance,~\cite{chen2014constructive,dinuzzo2015kernels}).

\subsection{Errors-in-variables system identification}
Errors-in-variables models are often described by the set of
equations~\cite{soederstroem2007errors},
\begin{align}
  y_t &= \del{g \ast w}_t + \eta_t,\\
    v_t &= w_t + \varepsilon_t.
\end{align}
It is clear that this type of models naturally fit into the
uncertain-input framework of~\eqref{eq:ui_model}.
In particular, we can consider the classical errors-in-variables
problem of identifying a parametric model of $S$ when $w_t$ is the realization of
a stationary stochastic signal with a rational
spectrum~\cite{castaldi1996identification}. In this case, we can
write $\{w_t\}$ as the filtered white noise process
\begin{equation}
  w_t = \frac{C(q;\theta)}{D(q;\theta)}e_t,
\end{equation}
where $e_t$ is unitary variance Gaussian white noise,
and $C(q;\theta)$ and $D(q;\theta)$ are complex polynomials in the one-step
shift operator $q$.

From this expression, we see that $w$ is a Gaussian process with
zero mean and covariance matrix $\Sigma_w(\rho)$ that depends on the
parameterization of the input filter. Using a parametric model for the
system, we obtain an uncertain-input system with
\begin{equation}
  \begin{aligned}
    {[\mu_g(\rho)]}_i &= g_i(\rho), & {[K_g(\rho)]}_{i,j} &= 0,\\
    {[\mu_w(\rho)]}_i &= 0, & {[K_w(\rho)]}_{i,j} &= \Sigma_w(\rho).
  \end{aligned}
\end{equation}
Alternatively, we could estimate all samples of the input signal
with the choice ${[\mu_w(\rho)]}_i = \theta_i$ and
$[K_w(\rho)]=0$, even though this may lead to nonidentifiability of the
model~\cite{soederstroem2003why,zhang2015errors,risuleo2016kernel}.

\subsection{Blind system identification}
Blind system identification can also be cast as the problem of identifying an
uncertain-input model by setting the input noise variance to $\sigma_v^2 =
\infty$ (this indicates that no input measurements are available). In this
case, different parameterizations of the input
lead to different models for the input process. For instance, we can consider
the parameterization of the input as a switching signal with known switching
instants $T_0<T_1<\cdots<T_p$; in this case we can choose
${[\mu_w(\rho)]}_i = h_i^t\theta$,
where $h_i$ is a selection vector that is nonzero in the $i$th interval:
\begin{equation}
  \begin{cases}{}
    {[h_i]}_j = 1 & \text{if}\;\;T_{i-1}< t \leq T_i,\\
    {[h_i]}_j = 0 & \text{otherwise}.
  \end{cases}
\end{equation}
Models similar to this one were used, for instance, in~\cite{ohlsson2014blind}
and~\cite{bottegal2015blind}.

\subsection{Cascaded system identification}\label{sec:cascaded_systems}
In cascaded linear systems, the output of one linear system is used as the input
to a second linear system (see Figure~\ref{fig:cascade}).
\begin{figure}[htb]
  \centering
  \includegraphics{cascade_blockscheme.pdf}
  \caption{Cascaded linear systems.}\label{fig:cascade}
\end{figure}

For sake of argument, we consider nonparametric models for both linear systems
(the reasoning also holds for parametric models):
\begin{equation}
  g_1 \sim \mathcal{N}\big(0,K_1(\theta)\big),\qquad g_2 \sim \mathcal{N}\big(0,K_2(\rho)\big).
\end{equation}
Because $g_1$ is a Gaussian vector, the intermediate variable $w$ is also a
Gaussian vector, with zero mean and covariance matrix given by
\begin{equation}\label{eq:cascaded_input}
  K_w(\theta) = UK_1(\theta) U^T,
\end{equation}
where $U:=\mathrm{T}_{N\times N}(u)$ is the Toeplitz matrix of the input signal $u_t$. Therefore, we
can model the linear cascade as an uncertain-input model with input modeled as a
zero-mean process with covariance matrix given by~\eqref{eq:cascaded_input}
where, for instance, we use the first-order stable spline kernel introduced
in~\eqref{eq:stable-spline}. The same choice of kernel can be made for $K_2(\rho)$.

\subsection{Hammerstein model identification}\label{sec:hammersteins}
The Hammerstein model is a cascade of a static nonlinear function followed by a
linear dynamical system (see Figure~\ref{fig:hammerstein}).
\begin{figure}[htb]
  \centering
\includegraphics{hammerstein_blockscheme.pdf}
  \caption{The Hammerstein model.}\label{fig:hammerstein}
\end{figure}

In the Hammerstein model, the intermediate variable $w_t$ is not observed (which,
symbolically, corresponds to an infinite $\sigma_v^2$). If we consider models
for the input block  that are  combinations of known basis
functions~\cite{bai1998optimal}, according to
\begin{equation}
  w_t = \sum_{j=1}^p \theta_j\, \varphi_j(u_t)
\end{equation}
we can collect the values of the unknown input $w$ and the parameters in a vector such that
\begin{equation}
  w = \Phi \theta,\qquad {[\Phi]}_{i,j} = \varphi_j(u_i).
\end{equation}
This can be modeled as an uncertain-input model with
$\mu_w(\theta) = \Phi\theta$ and $K_w(\theta) = 0$.

The uncertain-input framework also encompasses nonparametric models for the
input nonlinearity. For instance, we can model the Hammerstein cascade as the
uncertain-input model with the Gaussian radial-basis-function kernel as input model:
\begin{equation}\label{eq:rbf}
  {[K_w(\theta)]}_{i,j} = \theta_1
  \exp\left\{-\frac{1}{\theta_2}{(u_i-u_j)}^2\right\}.
\end{equation}
As for the linear system, we can use either parametric or nonparametirc
modeling approaches (see~\cite{bai1998optimal,risuleo2015new}).

\section{Estimation of uncertain-input models}\label{sec:estimation}
As discussed in Section~\ref{sec:ui_systems}, we suppose that we have collected
$N$ samples of the output $y_t$ and, possibly, $N$ samples of the noisy input
signal $v_t$ (in some applications, such as Hammerstein models and blind system
identification, these samples are not available). Whenever present, we assume
that the external input $u_t$ is completely known. We consider the following
identification problem.

\begin{prob}\label{prob:sysid}
  Given the $N$-dimensional vectors of measurements $y$ and $v$,  generated
  according to~\eqref{eq:ui_model}, estimate the impulse response $g$, the
  unknown input $w$, and the hyperparameters
  $\tau = \{\rho,\theta,\sigma_y^2,\sigma_v^2\}$.
\end{prob}
Because we are using the Gaussian process model~\eqref{eq:gp_models}, we have
natural candidates for the estimates of $g$ and $w$. Interpreting~\eqref{eq:model_g} and~\eqref{eq:model_w} as prior distributions of the
unknowns, we know that the best estimates given the data (in the minimum
mean-square error sense) are the conditional expectations
\begin{equation}
  g^\star = \mathbf{E}\big[ g\big| y,v\big],\quad w^\star = \mathbf{E}\big[w\big|y,v\big].
\end{equation}
However, these conditional expectations depend on the value of the
hyperparameter vector $\tau$. Because this value is not available, we follow an
empirical Bayes approach~\cite{maritz1989empirical} and we approximate
the true conditional expectations---that correspond to the true values of the
hyperparameters $\tau$---with the conditional expectations
\begin{equation}\label{eq:solution_posteriormean}
    \hat g := \int g\, \mathrm{p}(g|y,v\,;\hat\tau) \dif g,\qquad
    \hat w := \int w\, \mathrm{p}(w|y,v\,;\hat\tau) \dif w,
\end{equation}
where we are using estimated values $\hat \tau$ of the hyperparameters. In the
empirical Bayes approach, the estimates of the hyperparameters are chosen
by maximizing the marginal likelihood of the data,
\begin{equation}\label{eq:solution_marginal}
  \hat \tau := \arg \max_{\tau} \log \mathrm{p}(y,v\,;\tau),
\end{equation}
where $\mathrm{p}(y,v\,;\tau)$ is the marginal distribution of the measurements according
to the model in~\eqref{eq:ui_model}.

Solving~\eqref{eq:solution_marginal} yields the marginal likelihood estimate of
the hyperparameters that can be used to find the empirical Bayesian estimates of
$g$ and $w$ in~\eqref{eq:solution_posteriormean}. However, this approach
requires distributions that, in general, are not available in closed form.
Furthermore,~\eqref{eq:solution_marginal} is possibly a
high-dimensional optimization problem that does not admit an analytical
expression. To address this last problem, we use the
EM method to derive an iterative algorithm that
solves~\eqref{eq:solution_marginal}. We start by rewriting the marginal
likelihood as
\begin{equation}
  \mathrm{p}(y,v\,;\tau) = \int \mathrm{p}(y,v,g,w\,;\tau) \dif g \dif w.
\end{equation}
With this observation, we can see~\eqref{eq:solution_marginal} as a
maximum likelihood problem with latent variables, where the latent variables are
$g$ and $w$. Appealing to the theory of the EM method, we have that iterating
the two steps
\begin{description}
  \item[E-step:] Given an estimate $\hat \tau^{(k)}$ of $\tau$, construct the
    following lower bound of the marginal likelihood
    \begin{equation}\label{eq:Qtrue}
      \hspace{-1em} Q(\tau,\hat \tau^{(k)})\! = \!\int \!\log \mathrm{p}(y,v,g,w\,;\tau)
      \mathrm{p}(g,w|y,v\,;\hat \tau^{(k)})\dif g \dif w;
    \end{equation}
  \item[M-step:] Update the hyperparameter estimates as
    \begin{equation}
      \hat \tau^{(k+1)} = \arg\max_{\tau} Q(\tau,\hat \tau^{(k)});
    \end{equation}
\end{description}
from an arbitrary initial condition $\hat \tau^{(0)}$,
we obtain a sequence of estimates $\{\hat \tau^{(k)}\}$ of increasing
likelihood, which converges to a stationary point of the marginal likelihood of
the data. In practice, this stationary point will always be a local maximum:
saddle points are numerically unstable and minimal perturbations will drive the
sequence of updates away from them~\cite{mclachlan2007algorithm}.

Using the EM method, we have transformed the problem of maximizing the marginal
likelihood into a sequence of optimization problems. The whole point of the EM method is that
these problems should be simpler to solve than the original optimization
problem.

In addition to using the EM method to solve the marginal likelihood problem, we
can rewrite~\eqref{eq:solution_posteriormean} as
\begin{equation}\label{eq:posteriormeans}
   \hat g := \int g\, \mathrm{p}(g,w|y,v\,;\hat\tau) \dif w\dif g,\qquad
   \hat w := \int w\, \mathrm{p}(g,w|y,v\,;\hat\tau) \dif g\dif w.
\end{equation}
Comparing~\eqref{eq:posteriormeans} and the $Q$ function in the E-step, we see
that the solution of Problem~\ref{prob:sysid} using the procedure we have
described depends on expectations with respect to
the distribution $\mathrm{p}(g,w\,| y,v\,;\tau)$. This
distribution is, in general, not available in closed form. In the next section
we present three special cases when this distribution can be computed in closed form
and we present the resulting estimation algorithms. In
Section~\ref{sec:approximations}, we show two different ways to approximate
this joint posterior distribution in the general case.

\section{Cases with degenerate prior distributions}
There are cases where the integrals~\eqref{eq:Qtrue}
and~\eqref{eq:posteriormeans}, required to estimate uncertain-input
systems, admit closed-form solutions. This happens when either the prior for $g$
or for $w$ (or both) are degenerate distributions. This means that,
symbolically, we let the covariances $K_g(\rho)$ and $K_w(\theta)$ go to zero
and, respectively,
\begin{equation}
  \mathrm{p}(g;\rho) \to \delta\del{g-\mu_g(\rho)}, \quad
  \mathrm{p}(w;\theta) \to \delta\del{w-\mu_w(\theta)}.
\end{equation}
From these expressions, we see that the models of the unknown quantities $g$ and
$w$ are uniquely determined by the parameter vector $\tau$ (there is no
uncertainty or variability): therefore, we refer to these kind of models as
\emph{parametric models}. We now present three cases of parametric models that
admit closed form expressions for the EM algorithm.

\subsection{Semiparametric model}
The first model is called \emph{semiparametric}. It is obtained when
$K_w(\theta) \to 0$. This effectively means that the prior
density~\eqref{eq:model_w} collapses into the Dirac density centered around the
mean function, and the posterior distributions of the unknowns admit closed
form expressions:

\begin{lem}\label{lem:pigs_posterior_w}
  Consider the uncertain-input system~\eqref{eq:ui_model}. In the
  limit when $K_w(\theta)\to 0$, we have that
$\mathrm{p}(w|y,v;\tau) = \delta\del{w\!-\!\mu_w(\theta)}$.
\end{lem}
\begin{proof}\label{pf:pigs_posterior_w}
  When $K_w(\theta)\to 0$, the prior density becomes the degenerate normal
  distribution $\mathrm{p}(w;\theta) = \delta\big(w-\mu_w(\theta)\big)$.
  From the law of conditional expectation, we have
  \begin{equation}\label{eq:pf_pigs_posterior_w}
    \mathrm{p}(w|y,v;\tau) = \frac{\mathrm{p}(y,v|w;\tau)\delta(w - \mu_w(\theta))}{\mathrm{p}(y,v;\tau)};
\end{equation}
in addition, the evidence becomes
\begin{equation}
  \mathrm{p}(y,v;\tau) = \int \mathrm{p}(y,v|w;\tau)\mathrm{p}(w;\theta)\dif w =
  p(y,v|\mu_w(\theta);\tau).
\end{equation}
Plugging this expression into~\eqref{eq:pf_pigs_posterior_w} we have the
result.
\end{proof}

\begin{lem}\label{lem:pigs_posterior_g}
  Consider the uncertain-input system~\eqref{eq:ui_model}. In the limit when
  $K_w(\theta)\to 0$, the posterior distribution $\mathrm{p}(g|y,v;\tau)$ is Gaussian with
  covariance matrix and mean vector given by
  \begin{equation}\label{eq:pigs_posterior_g}
    P_g   = \del{\frac{1}{\sigma_y^2}{M_w(\theta)}^T M_w(\theta) +
    {K_g(\rho)}^{-1}}^{-1},\qquad
    m_g = P_g\del{\frac{1}{\sigma_y^2}{M_w(\theta)}^T
  y+{K_g(\rho)}^{-1}\mu_g(\rho)},
\end{equation}
where $M_w(\theta):=\mathrm{T}_{N\times N}\big(\mu_w(\theta)\big)$.
\begin{equation}
\end{equation}
\end{lem}
\begin{proof}
  Note that $y|g,w,v;\tau$ is an affine transformation of the
  Gaussian random variable $\varepsilon$; hence, it is Gaussian. By
  the law of conditional expectation and ignoring terms independent of $g$, we
  have that
  \begin{equation}
    \begin{aligned}
      \log\, \mathrm{p}(g|w,y,v;\tau) &\cong \log p(y|g,w;\tau) + \log p(g;\rho)
                            \cong -\frac{1}{2\sigma_y^2}\enVert{y - M_w(\theta)g}^2 - \frac{1}{2}\enVert{g}^2_{{K_g(\rho)}^{-1}}\\
                            &\cong -\frac{1}{2} \enVert{g}_{P_g^{-1}} + g^T m_g\cong - \frac{1}{2}\enVert{g - m_g}^2_{{P_g}^{-1}}
    \end{aligned}
  \end{equation}
  where $P_g$ and $m_g$ are defined in~\eqref{eq:pigs_posterior_g}.
  Because it is quadratic, the posterior distribution of $g$ is Gaussian,
  with the indicated covariance matrix and mean vector.
\end{proof}

Thanks to Lemma~\ref{lem:pigs_posterior_w} and
Lemma~\ref{lem:pigs_posterior_g}, the E-step can be computed analytically when
$K_w(\theta) = 0$, and
the function $Q(\tau, \hat \tau^{(k)})$ admits a closed-form expression. To
this end,
let $\Delta^{k}$ be the $N$ by $N$ matrix given by
\begin{equation}
  {[\Delta^k]}_{i,j} =
  \begin{cases}{}
    1 & \text{if } i+j-1 = k\\
    0 & \text{otherwise}
  \end{cases}
\end{equation}
and let
\begin{equation}\label{eq:R}
  R = \begin{bmatrix}
    \Delta^1 & \Delta^2 & \Delta^3 & \cdots & \Delta^N
  \end{bmatrix}.
\end{equation}
Then, we have the following result.

\begin{thm}\label{thm:Qpigs}
  Consider a semiparametric uncertain-input model with $K_w(\theta) = 0$.
  Let $\hat \tau^{(k)}$ be estimates of the hyperparameters at the $k$th
  iteration of the EM method and let $\hat g^{(k)}$ and $\hat P_g^{(k)}$ be the moments
  in~\eqref{eq:pigs_posterior_g} when $\tau=\hat \tau^{(k)}$. Define
  \begin{equation}
    \begin{aligned}
    \hat R_y(\theta) &= \enVert[1]{ y - \hat G^{(k)}\mu_w(\theta) }^2,\\
    \hat R_v(\theta) &= \enVert[1]{v - \mu_w(\theta)}^2,\\
    \hat S_g^{(k)} &=  R^T \del[1]{I_N \otimes \hat P_g^{(k)}} R;
    \end{aligned}
  \end{equation}
  where $\hat G^{(k)} = T_{N\times N}(\hat g^{(k)})$.
  Then, the function $Q(\tau,\hat \tau^{(k)})$ is
  given by
   \begin{equation}\label{eq:Qpigs}
     \begin{aligned}
       Q(\tau,\hat \tau^{(k)}) &= -\frac{1}{2\sigma_v^2}\hat R_v(\theta) - \frac{N}{2}\log \sigma_v^2
     -\frac{1}{2\sigma_y^2}\del{\hat R_y(\theta) + \enVert{\mu_w(\theta)}_{
     \hat S_g^{(k)}}^2}
     - \frac{N}{2}\log
     \sigma_y^2\\
     &- \frac{1}{2} \enVert{\hat g^{(k)} - \mu_g(\rho)}^2_{{K_g(\rho)}^{-1}} - \frac{1}{2} \trace\cbr{{K_g(\rho)}^{-1}\hat P_g^{(k)}} -
     \frac{1}{2}\log\det K_g(\rho).
     \end{aligned}
   \end{equation}
\end{thm}
\begin{proof}
  See Appendix~\ref{pf:Qpigs}.
\end{proof}

From~\eqref{eq:Qpigs}, we see that the optimization with respect to $\theta$ is
not independent of $\sigma_y^2$ and $\sigma_v^2$. Therefore, to update the
hyperparameter $\theta$, we use a conditional-maximization
step~\cite{meng1993maximum}, where we keep the noise variances fixed
to their values at the previous iterations. The use of the conditional
maximization step allows us to write the updates of the EM method in closed
form:

\begin{cor}\label{cor:Qpigs}
  At the $k$th iteration of the EM method, the parameters can be updated as
  \begin{equation}
  \begin{aligned}
    \hat \rho^{(k+1)} &= \arg\min_{\rho}
  \trace\cbr[1]{{K_g(\rho)}^{-1}\hat P_g^{(k)}} + \log \det K_g(\rho)\\
  &+ \enVert[1]{\hat g^{(k)}-\mu_g(\rho)}^2_{K_g^{-1}(\rho)},\\
  \hat\theta^{(k+1)}& = \arg \min_\theta\!
  \frac{\hat R_v(\theta)}{2\hat \sigma_v^{(k)\,2}} +
  \frac{1}{2\hat\sigma_y^{(k)\,2}}\!\del[2]{\!
    \hat R_y(\theta)\!+\! \enVert{\mu_w(\theta)}^2_{
    \hat S_g^{(k)}}\!}\!,\\
    \hat \sigma_y^{(k+1)} &= \frac{1}{N}\del{\hat R_y(\hat \theta^{(k+1)}) +
    \enVert[1]{\mu_w(\hat \theta^{(k+1)})}^2_{S_g^{(k)}}},\\
    \hat \sigma_v^{(k+1)} &= \frac{1}{N}\hat R_v(\hat \theta^{(k+1)})\,.
\end{aligned}
\end{equation}
\end{cor}
\begin{proof}
  Follows from the two-step maximization of~\eqref{eq:Qpigs}: first, maximize
  with respect to $\theta$ and $\rho$ keeping $\sigma_y^2$ and $\sigma_v^2$
  fixed to their values at the previous iteration; then, maximize with respect
  to $\sigma_y^2$ and $\sigma_v^2$ using the updated values of the
  hyperparameters.
\end{proof}

Thanks to Corollary~\ref{cor:Qpigs}, we have a simple way to compute the EM
estimates of the kernel hyperparameters and of the noise variances for
semiparametric models with $K_w(\theta) = 0$: starting
from an initial value of the unknown parameters, we first update the
hyperparameters $\rho$ and
$\theta$, then we use the new values to update the noise variances $\sigma_y^2$
and $\sigma_v^2$. Under mild regularity conditions, this procedure yields a
sequence of estimates that converges to a local maximum of the marginal
likelihood~(it is a Generalized EM sequence, see~\cite{wu1983convergence}).

\begin{rem}
  In this section, we have presented the case when $K_w(\theta)\to 0$. However,
  thanks to the symmetry of the model assured by~\eqref{eq:toeplitz_property},
  the same kind of algorithm works when $K_g(\rho)\to 0$ (by exchanging the
  roles of $g$ and $w$).
\end{rem}

\subsection{Parametric model}
In case we let both the input and the system covariance matrices go to zero,
all the variability in the model is removed, and we are left with classical
\emph{parametric} models. In this case, the marginal likelihood of the data collapses
into the likelihood where the impulse response and the input are
replaced with the parametric models $\mu_g(\rho)$ and $\mu_w(\theta)$
\begin{equation}
  \mathrm{p}(y,v;\rho,\theta,\sigma^2) = \int\!\!
  \mathrm{p}(y,v|g,w;\sigma^2)\mathrm{p}(g;\rho)\mathrm{p}(w;\theta)\dif g\dif w
  =\mathrm{p}(y,v|\,g\!=\!\mu_g(\rho),\,w\!=\!\mu_w(\theta);\,\sigma^2).
\end{equation}
In other words, the marginal likelihood of the data is the distribution of the
data conditioned on the events $g=\mu_g(\rho)$ and $w=\mu_w(\theta)$. This
distribution is given in closed form by
\begin{equation}\label{eq:Qparametric}
  \log \mathrm{p}(y,v|\mu_g(\rho),\mu_w(\theta);\sigma^2) =
  -\frac{1}{2\sigma_y^2}\enVert{y-M_w(\theta)\mu_g(\rho)}^2
  -\frac{N}{2}\log\sigma_y^2
  -\frac{1}{2\sigma_v^2}\enVert{v-\mu_w(\theta)}^2-\frac{N}{2}\log \sigma_v^2.
\end{equation}
where $M_w(\theta)$ is the Toeplitz matrix of $\mu_w(\theta)$.

In this parametric-model case, we have that the posterior means reduce to the
prior means and the maximum marginal-likelihood criterion collapses into the
classical maximum-likelihood or prediction-error estimation method. To
estimate the system, we first maximize~\eqref{eq:Qparametric} to find the
parameter values $\hat \tau$; then, we estimate the system with
\begin{equation}
  \hat g = \mu_g(\hat\rho),\quad \hat w=\mu_w(\hat \theta).
\end{equation}

The strategy to maximize~\eqref{eq:Qparametric} depends on the specific
structure of the problem. In some applications, concentrated-likelihood or
integrated-likelihood approaches have been proposed (for a review,
see~\cite{berger1999integrated}). An interesting consistent approach, for the
parametric EIV case, has been proposed in~\cite{zhang2015errors}.
In~\cite{bai2004convergence}, the authors show that if $g$ and $w$ are linearly
parameterized, alternating between estimation of $g$ and of $w$ leads to the
minimum of~\eqref{eq:Qparametric}.

\begin{rem}
  The EM based algorithm presented in Section~\ref{sec:estimation} cannot be
  used in the parametric model case because of the impulsive posterior
  distributions: during the M-step, the method is overconfident in the current
  value of the parameters and no update occurs. However, the EM method can be
  used in the parametric case by considering a covariance matrix that shrinks
  toward zero at every iteration.
\end{rem}

\section{Approximations of the joint posterior
distribution}\label{sec:approximations}
In the previous section, we have shown three cases in which the collapse of the prior
distribution allows us to express the marginal likelihood of the data and the
posterior distributions in closed form. In general, however, these distributions do not
have a closed form expression. Therefore, in this section, we present two ways to
approximate the joint posterior distribution
$\mathrm{p}(g,w\,| y,v\,;\tau)$.
In the first, we make a particle approximation. The particles are drawn
from the joint posterior using an MCMC method. In the second, we make a
variational approximation of the joint posterior.

\subsection{Markov Chain Monte Carlo integration}\label{sec:mcmc}
Monte Carlo methods are built around the concept of \emph{particle
approximation}. In a particle approximation method, a density with a complicated
functional form is approximated with a set of point probabilities---that is, we
approximate a density $\mathrm{p}(x)$ according to
\begin{equation}
  \mathrm{p}(x) \approx \frac{1}{M}\sum_{j=1}^M \delta(x-x_j).
\end{equation}
If the particle locations $x_j$ are drawn from $\mathrm{p}(x)$, and the number of
particles $M$ is large enough, the expectation of any measurable function $f(x)$
over any set can be approximated as
\begin{equation}\label{eq:monte_carlo_integral}
  \mathbf{E}\{f(x)\} = \int\!f(x)\,\mathrm{p}(x) \dif x \approx \frac{1}{M}\sum_{j=1}^M f(x_j),
\end{equation}
where $\{x_j\}$ are drawn from  $\mathrm{p}(x)$. This result comes directly from the sampling
property of the Dirac density $\delta(\,\cdot\,)$. From a different
perspective, we can see~\eqref{eq:monte_carlo_integral} as an estimation of the
true expectation. With this interpretation, we have that this estimator is
unbiased,
\begin{equation}
  \mathbf{E}\Bigg\{ \frac{1}{M}\sum_{j=1}^M
   f(x_j)\Bigg\} = \mathbf{E}\{f(x)\},
\end{equation}
and its covariance is inversely proportional to the number of samples used,
\begin{equation}
  \mathbf{cov}\Bigg\{ \frac{1}{M}\sum_{j=1}^M
  f(x_j)\Bigg\} = \frac{1}{M}\mathbf{cov}\{f(x)\}.
\end{equation}
In practice, the number of samples needed depends on the specific application: in
certain applications, few particles (say 10 or 20) may suffice; in other
applications, we might need a much larger number of particles (in the order of
thousands; for a complete treatment, see~\cite[Chapter~11]{bishop2006pattern}).

When implementing Monte Carlo integrations, a common approach is MCMC\@. In
these methods, we set up a Markov chain whose stationary distribution is the
distribution we want to approximate and we run it to
collect samples~\cite{gilks1996markov}.

One convenient way to create a Markov chain is \emph{Gibbs sampling}. Using
this method, we obtain a particle approximation of a joint distribution (called
the \emph{target distribution}) by
sampling from all the \emph{full conditional} distributions---the distribution of
one random variable conditioned on all other variables---in sequence. This
procedure results in a Markov chain that has the target distribution as its
stationary distribution. Contrary to many other sampling methods, Gibbs
sampling does not include a rejection step; this means that the samples
proposed at every step are accepted as samples from the chain. This may lead to
faster mixing and decorrelation of the chain compared to other MCMC
methods~\cite[Chapter~11]{bishop2006pattern}.

The main drawback with Gibbs sampling is that we must sample the full
conditional distributions of all variables. Therefore, it is only applicable if
these distributions have a functionally convenient form. In the case at hand,
we have the following results.
\begin{lem}\label{lem:cond_g} Consider the uncertain-input model~\eqref{eq:ui_model}. The
  density $\mathrm{p}(g|y,w;\tau)$ is Gaussian with covariance
  matrix and mean vector given by
\begin{equation}\label{eq:conditional_g_pars}
    P_g = {\left(\frac{W^T W}{\sigma_y^2} +
    {K_g(\rho)}^{-1}\right)}^{-1},\qquad
    m_g =  P_g\left(\frac{W^T y}{\sigma_y^2} +
  {K_g(\rho)}^{-1}\mu_g(\rho)\right).
\end{equation}
\end{lem}
\begin{proof}
  The proof follows along the same line of reasoning as the proof of Lemma~\ref{lem:pigs_posterior_g}.
\end{proof}

\begin{lem}\label{lem:cond_w} Consider the uncertain-input model~\eqref{eq:ui_model}. The
  density $\mathrm{p}(w|y,v,g;\tau)$ is Gaussian with covariance matrix
  and mean given by
\begin{equation}\label{eq:conditional_w_pars}
    P_w = {\left(\frac{G^T G}{\sigma_y^2} + \frac{I_N}{\sigma^2_v} +
  {K_w(\theta)}^{-1}\right)}^{-1},\qquad
 m_w =P_w\left(\frac{G^T y}{\sigma_y^2} + \frac{v}{\sigma_v^2} +
 {K_w(\theta)}^{-1}\mu_w(\theta)\right).
\end{equation}
\end{lem}
\begin{proof}
  Because $y$ and $v$ are conditionally independent given $w$ and $g$, we have
  that
  \begin{equation}
    \begin{aligned}
      &\log\mathrm{p}(w|y,v,g;\tau) \cong
      \log\del{\mathrm{p}(y|g,w;\sigma_y^2)\mathrm{p}(v|w;\sigma_v^2)\mathrm{p}(w;\theta)}\\
          &\cong -\frac{1}{2\sigma_y^2}\enVert{y \!-\! Gw}^2
          \!\!-\!\frac{1}{2\sigma_v^2}\enVert{v \!-\! w}^2\! \!-\! \frac
      {1}{2} \enVert{w \!-\! \mu_w(\theta)}^2_{{K_w(\theta)}^{-1}}\\
      &\cong -\frac{1}{2}\enVert{w}^2_{P_w^{-1}}
      + w^T m_w \cong -\frac{1}{2} \enVert{w - m_w}_{P_w^{-1}}^2
    \end{aligned}
  \end{equation}
  where $P_w$ and $m_w$ are given in~\eqref{eq:conditional_w_pars}. The
  log-density of $w|y,v,g,w;\tau$ is quadratic and, hence, it is Gaussian with
  the indicated mean vector and covariance matrix.
\end{proof}

\begin{rem}\label{rem:general_mc}
  In case we consider the more general Gaussian process
  model~\eqref{eq:gp_joint}, where $g$ and $w$ are a priori dependent,
  Lemma~\ref{lem:cond_g} and Lemma~\ref{lem:cond_w} still hold with slightly
  modified expressions for the mean vectors and covariance matrices (to account
  for the prior correlation). For instance, the conditional density of $g$ is
  Gaussian with covariance matrix and mean vector given by
  \begin{equation}
    P_g \!=\! {\left(\frac{W^T W}{\sigma_y^2} +
    \Lambda_g(\rho,\theta)\right)}^{-1},\qquad
    m_g \!=\!  P_g\!\left( \!\frac{W^T y }{\sigma_y^2} +
    \Lambda_g(\rho,\theta)\mu_g(\rho) + \Lambda_{gw}(\rho,\theta)(w \!-\!
  \mu_w(\theta))\!\right)\!,
\end{equation}
where $\Lambda_{gw}(\rho,\theta)$ and $\Lambda_{g}(\rho,\theta)$ are,
respectively, the lower left and right blocks of the inverse of the prior
covariance matrix.
\end{rem}

In view of Lemma~\ref{lem:cond_g} and Lemma~\ref{lem:cond_w}, we can easily
set up the Gibbs sampler to draw from the joint posterior distribution: from any
initialization of the impulse response $g^{(0)}$ and of the input signal $w^{(0)}$, we sample
\begin{equation}\label{eq:gibbs}
  \begin{aligned}
    g^{(j+1)}&|w^{(j)},y,v;\tau \sim \mathcal{N}(m_g^{(j)},P_g^{(j)}),\\
    w^{(j+1)}&|g^{(j+i)},y,v;\tau  \sim \mathcal{N}(m_w^{(j)},P_w^{(j)}).
  \end{aligned}
\end{equation}
where $m_g^{(j)}$ and $P_g^{(j)}$ are the mean and covariance
in~\eqref{eq:conditional_g_pars} when $w=w^{(j)}$, and where $m_w^{(j)}$ and
$P_w^{(j)}$ are the mean and covariance in~\eqref{eq:conditional_w_pars} when
$g=g^{(j+1)}$.

Because it is a Markov chain, the samples drawn using~\eqref{eq:gibbs} are
correlated, and subsequent samples have memory about the initial conditions and
are far away from the stationary distribution (which is equal to the target
distribution). Therefore, we discard the first samples of the Markov chain, and
we only retain the $M$ samples after a \emph{burn-in} of $B$ samples:
\begin{equation}\label{eq:burnin}
  \bar g^{(j)} = g^{(j+B)}, \quad
  \bar w^{(j)} = w^{(j+B)}, \quad j = 1,\ldots,M.
\end{equation}
If the burn-in is large enough, the Markov chain has lost its memory
about the initial conditions and is producing samples that come form the
stationary distribution. The choice of the length of the burn-in is a
difficult problem, and some heuristic algorithms have been
proposed (see~\cite[Section~1.4.6]{gilks1996markov}).

When we have drawn enough samples from
the Markov chain, we compute the Monte Carlo estimate of the function $Q$; in
other words, we replace the E-step in the EM method with a Monte Carlo E-step
(this is sometimes known as the MCEM method; see~\cite{wei1990monte}). We create the approximate lower
bound (at the $k$th iteration of the EM method) by setting
\begin{equation}
  Q^{\text{\textsc{mc}}}(\tau,\hat \tau^{(k)}) = \frac{1}{M_k} \sum_{j=1}^{M_k} \log \mathrm{p}(y,v,\bar g^{(j,k)},
    \bar w^{(j,k)};\tau)
\end{equation}
where $\bar g^{(j,k)}$ and $\bar w^{(j,k)}$ are samples from the stationary
distribution of~\eqref{eq:gibbs} at the $k$th iteration of the EM method. In
the uncertain-input case, the function $Q^{\text{\textsc{mc}}}$ is available in closed form as
a function of the sample moments of $g$ and $w$.
\begin{thm}\label{thm:Qmc}
  Let $\cbr[1]{\bar g^{(j,k)}}_{j=1}^{M_k}$ and
  $\cbr[1]{\bar w^{(j,k)}}_{j=1}^{M_k}$ be
  samples from the stationary distribution of the Gibbs
  sampler~\eqref{eq:gibbs} at the $k$th
  iteration of the EM method and define
  \begin{equation}\label{eq:mcem_moments}
    \begin{aligned}
      \hat g^{(k)} &= \frac{1}{M_k}\sum_{j=1}^{M_k} \bar g^{(j,k)},\quad \hat w^{(k)}
      = \frac{1}{M_k}\sum_{j=1}^{M_k} \bar w^{(j,k)},\\
      \hat P_g^{(k)} &= \frac{1}{M_k}\sum_{j=1}^{M_k} \del{\bar g^{(j,k)} -\hat g^{(k)}}\del{\bar g^{(j,k)} - \hat g^{(k)}}^T\\
      \hat P_w^{(k)}  &= \frac{1}{M_k}\sum_{j=1}^{M_k} \del{\bar w^{(j,k)} - \hat w^{(k)}}\del{
      \bar w^{(j,k)} - \hat w^{(k)}}^T\\
      \hat R_v^{(k)} &= \frac{1}{M_k}\sum_{j=1}^{M_k} \enVert{v - \bar w^{(j,k)}}^2,\\
      \hat R_y^{(k)} &= \frac{1}{M_k}\sum_{j=1}^{M_k} \enVert{y - \bar G^{(j,k)}\bar w^{(j,k)}}^2.
\end{aligned}
  \end{equation}
  Then, the function $Q^{\text{\textsc{mc}}}(\tau, \hat \tau^{(k)})$ is given by
  \begin{equation}\label{eq:Qmc}
    \mathmakebox[0.88\columnwidth][l]{\begin{aligned}
        \!Q^{\text{\textsc{mc}}}(\tau,\hat \tau^{(k)}) &= \!- \frac{\hat R_v^{(k)}}{2\sigma_v^2}
      \!-\! \frac{N}{2} \log
      \sigma_v^2 \!-\!\frac{\hat R_y^{(k)}}{2\sigma_y^2} \!-\! \frac{N}{2}\log \sigma_y^2
      \!-\! \frac{1}{2}\trace\cbr{{K_g(\rho)}^{-1} \hat P_g^{(k)}}\\
      \!&-\!\frac{1}{2} \enVert{\hat g^{(k)} \!-\!
    \mu_g(\rho)}^{2}_{{K_g(\rho)}^{-1}}
  \!-\!\frac{1}{2} \trace\cbr{{K_w(\theta)}^{-1} \hat P_w^{(k)}}
    \!-\!\frac{1}{2} \enVert{\hat w^{(k)} \!-\!
  \mu_w(\theta)}^{2}_{{K_w(\theta)}^{-1}}\\
  &\!-\!\frac{1}{2} \log \det K_g(\rho) \!-\!\frac{1}{2} \log \det K_w(\theta).
    \end{aligned}}
  \end{equation}
\end{thm}
\begin{proof}
  See Appendix~\ref{pf:Qmc}.
\end{proof}

In the M-step, we update the hyperparameters $\hat \tau^{(k)}$ by maximizing
the approximate lower bound of the marginal likelihood, $Q^{\text{\textsc{mc}}}$. Because of the closed
form expression in Theorem~\ref{thm:Qmc}, the M-step splits into the decoupled
optimization problems for the kernel hyperparameters and the noise variances
according to the following:
\begin{cor}\label{cor:Qmc}
  At the $k$th iteration of the EM method, the kernel hyperparameters can be
  updated as
  \begin{equation}
    \begin{aligned}
      \hat \rho^{(k+1)} &= \arg\min_\rho
      \enVert{\hat g^{(k)}-\mu_g(\rho)}^{2}_{{K_g(\rho)}^{-1}}
      + \trace\cbr{{K_g(\rho)}^{-1} \hat P^{(k)}_g} + \log \det K_g(\rho),\\
      \hat \theta^{(k+1)} &= \arg\min_\theta
      \enVert{\hat w^{(k)}-\mu_w(\theta)}^{2}_{{K_w(\theta)}^{-1}}
      + \trace\cbr{{K_w(\theta)}^{-1} \hat P^{(k)}_w} + \log \det K_w(\theta),
    \end{aligned}
  \end{equation}
  and  the noise variances can be updated as
  \begin{equation}
    \hat \sigma_v^{2\,(k+1)} = \frac{\hat R_v^{(k)}}{N},\qquad
    \hat \sigma_y^{2\,(k+1)} = \frac{\hat R_y^{(k)}}{N}.
  \end{equation}
\end{cor}
\begin{proof}
  Follows from direct maximization of~\eqref{eq:Qmc}.
\end{proof}

Thanks to Theorem~\ref{thm:Qmc} and Corollary~\ref{cor:Qmc}, we have a simple
way to compute the MCEM estimates of the kernel hyperparameters and of the
noise variances; starting from an initial value of the hyperparameters, we iterate
the following three steps:
\begin{enumerate}
  \item Run a Gibbs sampler according to~\eqref{eq:gibbs}.
  \item Collect the samples according to~\eqref{eq:burnin} and compute the moments according
    to~\eqref{eq:mcem_moments}.
  \item Update the parameters according to Corollary~\ref{cor:Qmc}.
\end{enumerate}
Under mild regularity conditions, these iterations yield a sequence of
parameter estimates that converges to a stationary point of the marginal likelihood of the
data (under the condition that the number of particles $M_k$ at iteration $k$ is
such that that $\sum_{k=1}^\infty M_k^{-1} = \infty$;
see~\cite{neath2013convergence}). Then, using the estimated
hyperparameters, we can run a new Gibbs sampler and approximate the integrals
in~\eqref{eq:posteriormeans} with averages over the samples:
\begin{equation}\label{eq:particle_means}
  \hat g \approx \frac{1}{M} \sum_{j=1}^M \bar g^{(j)},\quad
  \hat w \approx \frac{1}{M} \sum_{j=1}^M \bar w^{(j)}.
\end{equation}

\subsection{Variational Bayes approximation}\label{sec:vb}
The second  method we present is a variational approximation method.
Instead of approximating the unknown joint posterior density using sampling, we
propose an analytically tractable family of distributions and we look for the best
approximation of the unknown posterior density within that family.

The variational Bayes method hinges on the fact that
\begin{equation}
  \log\mathrm{p}(y,v,g,w;\tau) = \log(g,w|y,v;\tau) + \log\mathrm{p}(y,v;\tau).
\end{equation}
Hence, for any proposal distribution $\mathrm{q}$ in some family of
distributions $\mathcal{Q}$, we can write
\begin{equation}
  \log \mathrm{p}(y,v;\tau) = \log \frac{\mathrm{p}(y,v,g,w;\tau)}{\mathrm{q}(g,w)} - \log
  \frac{\mathrm{p}(g,w|y,v;\tau)}{\mathrm{q}(g,w)}.
\end{equation}
Taking the expectation with respect to $\mathrm{q}$ and observing that the left hand
side is independent of $g$ and $w$, we get that
\begin{equation}\label{eq:marginal_split}
  \log \mathrm{p}(y,v;\tau) = L(\mathrm{q}) + KL(\mathrm{q}),
\end{equation}
where we have defined the functional
\begin{equation}
  L(\mathrm{q}) = \int \log\left(\frac{\mathrm{p}(y,v,g,w;\tau)}{\mathrm{q}(g,w)}\right)\mathrm{q}(g,w)\dif
  g \dif w,
\end{equation}
and the \emph{Kullback-Leibler (KL) distance}~\cite{kullback1951information}
\begin{equation}
  KL(\mathrm{q}) = \int \log\left(\frac{\mathrm{q}(g,w)}{\mathrm{p}(g,w|y,v;\tau)}\right)\mathrm{q}(g,w)\dif
  g \dif w.
\end{equation}
Although the KL distance is not a metric---it is not symmetric
and it does not satisfy the triangle inequality---it is a useful measure of
similarity between probability
distributions (see~\cite[Section~1.6.1]{bishop2006pattern}).

Because the left hand side of~\eqref{eq:marginal_split} is independent of $\mathrm{q}$,
we can find the distribution $\mathrm{q}^\star$ with minimum distance
(in the KL sense) to the target distribution by maximizing
the functional $L(\mathrm{q})$ with respect to $\mathrm{q}\in \mathcal{Q}$,
\begin{equation}\label{eq:kl_optimization}
  \mathrm{q}^\star(g,w) =\arg\min_{\mathrm{q} \in \mathcal{Q}} KL(\mathrm{q}) = \arg\max_{\mathrm{q}\in \mathcal{Q}} L(\mathrm{q}).
\end{equation}
This technique allows us to use the known functional $L(\mathrm{q})$ to find the $\mathrm{q}$
with minimum KL distance to the unknown joint posterior distribution.

To use the variational approximation, we need to fix a family of distributions
$\mathcal{Q}$ among which to look for $\mathrm{q}^\star$. In this work, we use a
\emph{mean-field approximation}, meaning that we look for an approximation of
the posterior distribution where $g$ and $w$ are independent given the data; in
other words we consider proposal distributions that factorize into two
independent factors according to
\begin{equation}
  \mathrm{q}(g,w) = \mathrm{q}_g(g)\mathrm{q}_w(w).
\end{equation}

After choosing the family of proposal distributions, we need to find the best
approximation $\mathrm{q}^\star$ in terms of KL distance to the unknown
posterior distribution; in view of~\eqref{eq:kl_optimization}, the solution is
given by
\begin{equation}
  \mathrm{q}^\star(g,w) = \arg\max_{\mathrm{q}_g,\,\mathrm{q}_w} L(\mathrm{q}_g\mathrm{q}_w).
\end{equation}
Consider first the factor $\mathrm{q}_g$. We have that
\begin{equation}
  \begin{aligned}
    L&(\mathrm{q}_g\mathrm{q}_w) = \int
  \log\left(\frac{\mathrm{p}(y,v,g,w;\tau)}{\mathrm{q}_g(g)\mathrm{q}_w(w)}\right)
  \mathrm{q}_g(g) \mathrm{q}_w(w)\dif g \dif w,\\
  &\!\cong\!\int \!\sbr{\log\mathrm{p}(y,v,g,w;\tau)\mathrm{q}_w(w)\dif w - \log \mathrm{q}_g(g)
}\mathrm{q}_g(g)\dif g,
  \end{aligned}
\end{equation}
ignoring terms independent of $\mathrm{q}_g$.
If we define the distribution $\mathrm{p}_w(y,v,g;\tau)$ such that
\begin{equation}
  \log \mathrm{p}_w(y,v,g;\tau) = \int \log \mathrm{p}(y,v,g,w;\tau) \mathrm{q}_w(w)\dif w,
\end{equation}
we have that, again ignoring terms independent of $\mathrm{q}_g(g)$,
\begin{equation}
  L(\mathrm{q}_g\mathrm{q}_w) \cong - \int \log\del{\frac{\mathrm{p}_w(y,v,g;\tau)}{\mathrm{q}_g(g)}}\mathrm{q}_g(g)
  \dif g,
\end{equation}
which is the negative KL distance between the factor $\mathrm{q}_g$ and the
density $\mathrm{p}_w(y,v,g;\tau)$. Because the KL distance is nonnegative,
by choosing $\mathrm{q}_g^\star(g)=\mathrm{p}_w(y,v,g;\tau)$ (where the KL distance is zero)
we are maximizing the functional $L$ with respect to
$\mathrm{q}_g$. Considering now $\mathrm{q}_w(w)$, we can trace the same argument and find that
the optimal choice is
\begin{equation}\label{eq:q_w_star}
  \log \mathrm{q}^\star_w(w) = \int \log p(y,v,g,w;\tau)\mathrm{q}^\star_g(g)\dif g,
\end{equation}
where $\mathrm{q}^\star_g(g)$ is the solution of
\begin{equation}\label{eq:q_g_star}
  \log \mathrm{q}^\star_g(g) = \int \log p(y,v,g,w;\tau)\mathrm{q}^\star_w(w)\dif w.
\end{equation}
The maximum of $L(\mathrm{q}_g\mathrm{q}_w)$ is, therefore, the simultaneous solution
of~\eqref{eq:q_w_star} and~\eqref{eq:q_g_star}. The solution can be found with
the following iterative procedure: from an initialization
$\mathrm{q}^{(0)}_g$ and $\mathrm{q}_{w}^{(0)}$ of the densities, compute
\begin{equation}\label{eq:variational_bayes}
  \begin{aligned}
    \!\log \mathrm{q}^{(j+1)}_w(w) &= \int \log p(y,v,g,w;\tau)\mathrm{q}^{(j)}_g(g)\dif g,\\
    \!\log \mathrm{q}^{(j+1)}_g(g) &= \int \log p(y,v,g,w;\tau)\mathrm{q}^{(j+1)}_w(w)\dif w.
  \end{aligned}
\end{equation}
This iterative procedure will converge to the
simultaneous solution of~\eqref{eq:q_g_star} and~\eqref{eq:q_w_star}
(see~\cite[Chapter~10]{bishop2006pattern}; see also~\cite{boyd2004convex}).

As was the case for the Gibbs sampler, which can be used only if it easy to
sample from the full conditional distributions, the variational approximation of
the joint posterior is only useful if it is possible to compute the expectations
in~\eqref{eq:q_w_star} and~\eqref{eq:q_g_star}. In the uncertain-input case, we
have the following result.
\begin{thm}\label{thm:vb_gaussian}
  Let $\mathrm{q}_g^\star\mathrm{q}^\star_w$ be the factorized density with minimum KL distance
  to posterior density $\mathrm{p}(g,w|y,v;\tau)$, for a fixed value of the
  hyperparameters. Then, $\mathrm{q}^\star_g$ and $\mathrm{q}^\star_w$ are Gaussian distributions.
\end{thm}
\begin{proof}
  See Appendix~\ref{pf:vb_gaussian}.
\end{proof}

Theorem~\ref{thm:vb_gaussian} allows us to compute expectations with respect to
$\mathrm{q}^\star_g$ and $\mathrm{q}^\star_w$ easily. In addition, at every iteration
of~\eqref{eq:variational_bayes} the approximating densities remain Gaussian.
This allows us to write the update~\eqref{eq:variational_bayes} in terms
of the first and second moments of the approximating densities:
\begin{cor}\label{cor:iterative_vb}
  Let $ w^{(j)}$ and $ g^{(j)}$ be the mean vectors of $\mathrm{q}_w^{(j)}$
  and $\mathrm{q}_g^{(j)}$ at the $j$th iteration of~\eqref{eq:variational_bayes}
  and let $ P_w^{(j)}$ and $ P_g^{(j)}$ be the covariance matrices. Let
  $ g^{(j+1)}$, $ w^{(j+1)}$, $ P^{(j+1)}_g$, and
  $ P_w^{(j+1)}$ be the mean vectors and covariance matrices at the
  $(j+1)$th iteration. Let
  \begin{equation}
    \begin{aligned}
      T_g^{(j)}&= R\del{I_N \otimes \sbr{ P_g^{(j)}+  g^{(j)}
       g^{(j)\,^T}}}R^T,\\
      T_w^{(j+1)}&= R\del{I_N \otimes \sbr{ P_w^{(j+1)}+  w^{(j+1)}
           w^{(j+1)\,T}}}R^T.
    \end{aligned}
  \end{equation}
  where the matrix $R$ is defined in~\eqref{eq:R}. Then,
  \begin{equation}\label{eq:vb_moments}
    \mathmakebox[0.88\columnwidth][l]{
  \begin{aligned}
   \!P_w^{(j+1)}&\!=\!\del[3]{\frac{1}{\sigma_y^2}T_g^{(j)} + \frac{1}{\sigma^2_v}I_n +
    {K_w(\theta)}^{-1}}^{-1},\\
   \!w^{(j+1)} &\!=\! P_w^{(j+1)}\!\del[3]{\frac{ G^{(j)\,T}}{\sigma_y^2}y + \frac{1}{\sigma_v^2}v +
  {K_w(\theta)}^{-1}\mu_w(\theta)}\!,\\
     \!P_g^{(j+1)} &= \del[3]{\frac{1}{\sigma_y^2}T_w^{(j+1)} +
      {K_g(\rho)}^{-1}}^{-1},\\
     \!g^{(j+1)} &=  P_g^{(j+1)}\del[3]{\frac{ W^{(j+1)\,T}}{\sigma_y^2}y +
  {K_g(\rho)}^{-1}\mu_g(\rho)}.
  \end{aligned}}
\end{equation}
\end{cor}
\begin{proof}
  See Appendix~\ref{pf:iterative_vb}.
\end{proof}

Thanks to Corollary~\ref{cor:iterative_vb}, we can iteratively update the
moments of the Gaussian factors, and the iterations will converge to the moments
of optimal variational approximation of the joint posterior distribution.

\begin{rem}
  In case we consider the more general Gaussian process
  model~\eqref{eq:gp_joint}, the results of
  Theorem~\ref{thm:vb_gaussian}  and of Corollary~\ref{cor:iterative_vb} still
  hold with minor modifications (similarly to what is presented in
  Remark~\ref{rem:general_mc}). However, the approximation of posterior
  independence may not make sense when using a-priori dependent Gaussian
  process models.
\end{rem}

Using the factorized approximation of the joint distribution,
we can approximate the E-step in the EM method with a variational
E-step (this is sometimes known as the VBEM method,
see~\cite{beal2003variational}). We
create the variational approximation of the lower bound (at the $k$th iteration
of the EM method) by setting
\begin{equation}
  Q^{\text{\textsc{vb}}}(\tau,\hat\tau^{(k)}) \!:=\! \int\!
  \log\mathrm{p}(y,v,g,w;\tau)\hat \mathrm{q}_g^{(k)}(g)\hat \mathrm{q}_w^{(k)}(w)\dif w\dif g,
\end{equation}
where $\hat \mathrm{q}_g^{(k)}$ and $\hat\mathrm{q}_w^{(k)}$ are the limits of the variational
Bayes iterations with the hyperparameters set to $\hat \tau^{(k)}$.

Because the complete-data likelihood is quadratic in $g$ and $w$, the
approximation $Q^\text{\textsc{vb}}$ admits the closed form expression
in function of the moments of $g$ and $w$.
\begin{thm}\label{thm:Qvb}
  Let $\hat g^{(k)}$ and $\hat w^{(k)}$ be the mean vectors of $\hat\mathrm{q}_g^{(k)}$
  and of $\hat \mathrm{q}_w^{(k)}$, respectively, and let $\hat P^{(k)}$ and
  $\hat P^{(k)}$  be their covariance matrices. Define
  \begin{equation}
    \begin{aligned}
      \hat S_w^{(k)} &= R \del[1]{I_n \otimes \hat P_g^{(k)}} R^T,
       &\hat T_w^{(k)} &= \hat S_w^{(k)} \!+\! \hat W^{(k)T}\hat W^{(k)},\\
      \hat R_v^{(k)}  &= \enVert[1]{v - \hat w^{(k)}}^2,
      &\hat R_y^{(k)} &= \enVert[1]{y - \hat W^{(k)} \hat g^{(k)}}^2,
    \end{aligned}
  \end{equation}
  where $R$ is defined in~\eqref{eq:R}.
  Then,
  \begin{equation}\label{eq:Qvb}
    \begin{aligned}
      &Q^\text{\textsc{vb}}(\tau,\hat \tau^{(k)}) = -\frac{\hat R_v^{(k)}}{2\sigma_v^2}
      -\frac{N}{2}\log \sigma_v^2 - \frac{N}{2}\log \sigma_y^2- \frac{1}{2\sigma_y^2}\del{\hat R_y^{(k)} -
    \enVert[1]{\hat g^{(k)}}^2_{\hat S_w^{(k)}} - \trace\cbr{\hat T_w^{(k)} P_g^{(k)}}}\\
      &- \frac{1}{2}\trace\cbr{{K_g(\rho)}^{-1} \hat P_g^{(k)}} -\frac{1}{2} \enVert{\hat g^{(k)} -
    \mu_g(\rho)}^{2}_{{K_g(\rho)}^{-1}}-\frac{1}{2}
    \trace\cbr{{K_w(\theta)}^{-1} \hat P_w^{(k)}}\\
    &-\frac{1}{2} \enVert{\hat w^{(k)} -
  \mu_w(\theta)}^{2}_{{K_w(\theta)}^{-1}}-\frac{1}{2} \log \det K_g(\rho)
  -\frac{1}{2} \log \det K_w(\theta).
    \end{aligned}
  \end{equation}
\end{thm}
\begin{proof}
  See Appendix~\ref{pf:Qvb}.
\end{proof}

Thanks to the structure of the function $Q^\text{\textsc{vb}}(\tau,\hat \tau^{(k)})$,
the M-step splits into
decoupled optimization problems for the kernel hyperparameters and for the
noise variances.
\begin{cor}\label{cor:Qvb}
  At the $k$th iteration of the EM method, the kernel hyperparameters can be
  updated as
\begin{equation}\label{eq:vbem_update_rho}
    \begin{aligned}
      \hat \rho^{(k+1)} &= \arg\min_\rho
      \enVert{\hat g^{(k)}-\mu_g(\rho)}^{2}_{{K_g(\rho)}^{-1}}+
      \trace\cbr{{K_g(\rho)}^{-1} \hat P^{(k)}_g} + \log \det K_g(\rho),\\
      \hat \theta^{(k+1)} &= \arg\min_\theta
      \enVert{\hat w^{(k)}-\mu_w(\theta)}^{2}_{{K_w(\theta)}^{-1}}+ \trace\cbr{{K_w(\theta)}^{-1} \hat P^{(k)}_w} + \log \det K_w(\theta),
    \end{aligned}
\end{equation}
and the noise variances can be updated as
\begin{equation}
  \begin{aligned}
    \hat \sigma_v^{(k+1)} &= \frac{\hat R^{(k)}_v}{N},\\
    \hat \sigma_y^{(k+1)} &=
  \frac{\hat R_y^{(k)}\!+\!
    \enVert[1]{\hat g^{(k)}}^2_{\hat S_w^{(k)}} \!+\! \trace\cbr{\hat T_w^{(k)}
  P_g^{(k)}}}{N}.
  \end{aligned}
\end{equation}
\end{cor}
\begin{proof}
  Follows from direct maximization of~\eqref{eq:Qvb}.
\end{proof}

Thanks to Theorem~\ref{thm:Qvb} and Corollary~\ref{cor:Qvb}, we have a simple
iterative proceduce to compute the VBEM estimates of the kernel hyperparameters
and of the noise variances; starting from an inital value of the
hyperparameters, we iterate the following two steps:
\begin{enumerate}
  \item Compute the moments of the variational approximation according to
    Corollary~\ref{cor:iterative_vb}.
  \item Update the hyperparameters according to Corollary~\ref{cor:Qvb}.
\end{enumerate}

Under mild regularity conditions, these iterations yield a sequence of
parameter estimates that converges to a stationary point of the marginal
likelihood of the data (see~\cite[Section~2.2]{beal2003variational}).
Then, we can run the iterations in Corollary~\ref{cor:iterative_vb} again to find the
posterior mean estimates of $g$ and $w$.

\section{Simulations}
In this section, we evaluate the methods proposed on some problems that can be
cast as problems of identifying uncertain-input systems.

\subsection{Cascaded linear systems}
In this numerical experiment, we estimate cascaded systems with the structure
presented in Section~\ref{sec:cascaded_systems}. We perform a Monte Carlo
experiment consisting of 500 runs. In each run, we generate two systems by
randomly sampling 40 poles and 40 zeros, in complex conjugate pairs, using the
following technique. We sample the poles randomly, with magnitudes uniformly between 0.4
and 0.8 and phases uniformly between 0 and $\pi$. We sample the zeros randomly,
with magnitudes uniformly between 0 and 0.92 and phases uniformly between 0 and
$\pi$. All systems are generated with unitary static gain. The noise variances on the
input and output measurements are $1$, respectively $1/100$, times the variance
of the corresponding noiseless signals; this means that the
sensor at the output of $S_2$ is considerably more accurate than the sensor at
the output of $S_1$.

We simulate the responses of the systems with a Gaussian
white-noise input with variance 1. We collect $N=200$ samples of the output, from zero
initial conditions, and we estimate the samples of the impulse responses of the
two systems.

As described in Section~\ref{sec:cascaded_systems}, the systems are modeled as
zero-mean Gaussian processes with first order stable-spline kernels.
All the methods are initialized with the choices $\rho_1=\theta_1=1$ and
$\rho_2=\theta_2 = 0.6$. The noise variances are initialized from the
sample variances of the errors of the linear least squares estimates of $g_1$
and $g_2$ from the noisy data.

In the experiment, we compare the following estimators.
\begin{description}
  \item[C-MCEM] The method described in Section~\ref{sec:mcmc}. It uses an
    MCMC approximation of the joint posterior with $B=400$ and $M=2000$. The
    EM iterations are stopped once the relative change in the parameter values
  is below $10^{-2}$.
  \item[C-VBEM] The method described in Section~\ref{sec:vb}. It uses a
    variational approximation of the joint posterior. The
    EM iterations are stopped once the relative change in the parameter values
  is below $10^{-2}$.
\item[C-2Stage] A kernel-based two-stage method. First, it estimates the
  first system in the cascade from $u$ and $v$. Then, it simulates the
  intermediate signal $\hat w$ as the response of the estimated system to  $u$
  and uses $\hat w$ and $y$ to estimate the second system in the cascade.
\item[C-Naive] A naive kernel-based estimation method. It estimates the
  first system in the cascade from $u$ and $v_t$ and the second system from
  $v_t$ and $y_t$. It corresponds to using the noisy signal $v_t$ as if it were
  the noiseless input to the second system in the cascade.
\end{description}

To evaluate the performance of the estimators, we use the following goodness-of-fit metric
\begin{equation}\label{eq:goodnessoffit}
  \mathrm{Fit}^g_j = 1- \frac{\enVert{g_j - \hat g_j}}{\enVert{g_j -
  \mathrm{mean}(g_j)}}
\end{equation}
where $g_j$ is the impulse response of the system at the $j$th Monte Carlo run,
and  $\hat g_j$ is an estimate of the same impulse response.

The results of the experiment are presented in
Figure~\ref{fig:boxplot_cascade}. The figure shows the boxplots of the fit of
the estimated impulse responses of the two blocks in the cascade over the
systems in the dataset.

\begin{figure}[htb]
  \centering
\includegraphics{boxplot_cascade.pdf}
  \caption{Results of the estimation of cascaded linear systems.}\label{fig:boxplot_cascade}
\end{figure}

From the figure, it appears that the proposed approximation methods are able to
reconstruct the cascaded model with higher accuracy than the alternative
approaches we have considered. Furthermore, there seems to be no clear
disadvantage in using the variational Bayes approximation as compared to the,
more correct, sampling-based approximation. Regarding the performance of the
methods in estimating $g_1$, we see that the methods C-MCEM and C-VBEM perform
better than the other methods (which give the same result). Both C-2Stage and
C-Naive only use the information in $v$ to estimate $g_1$, whereas C-MCEM and
C-VBEM use the full joint distribution of $v$ and $y$ to estimate $g_1$. Given
that in our setting the noise on $y$ is much lower than the noise on $v$,
there is information in $y$ that the joint methods are able to
leverage to improve the estimate of $g_1$ (similar phenomena were already
observed in~\cite{hjalmarsson2009system}, and in~\cite{everitt2013geometric}).
This allows C-MCEM and C-VBEM to better estimate $g_1$.

\subsection{Hammerstein systems}
In this numerical experiment, we estimate Hammerstein systems with the
structure presented in Section~\ref{sec:hammersteins}. We perform four Monte
Carlo experiments consisting of 500 runs. In each run, we generate a stable
transfer-function model by sampling poles and zeros in the complex plane. We
sample the poles, uniformly in magnitude and phase, in the annulus of radii 0.4
and 0.8. We sample the zeros uniformly in the disk of radius 0.92. We generate
the nonlinear transformation as a finite combination of Legendre polynomials
defined as
\begin{equation}
  \varphi_j(x) = 2^j\cdot \sum_{k=0}^j
  x^k\binom{j}{x}\binom{\tfrac{j+k-1}{2}}{j}.
\end{equation}
We sample the coefficients of the combination independently and uniformly in
the interval $\intcc{-1,1}$.

In each Monte Carlo experiment, we consider Hammerstein systems with
different orders for both the nonlinear system and the polynomial nonlinearity.
In Table~\ref{tab:hammerstein_orders}, we present the orders of the systems
considered in the various experiments.

\begin{table}[ht]
  \centering
  \caption{Orders of the Hammerstein systems used in the
  simulations.}\label{tab:hammerstein_orders}
  \begin{tabular}{ccc}
    \toprule
    Dataset & $S$ & $f(\cdot)$ \\
    \midrule
    LOLO (Low-Low) & $\{3,\ldots, 5\}$ & $\{5,\ldots, 10\}$\\
    HILO (High-Low) & $\{9,\ldots, 20\}$ & $\{5,\ldots, 10\}$\\
    LOHI (Low-High) & $\{3,\ldots, 5\}$ & $\{15,\ldots, 20\}$\\
    HIHI (High-High) & $\{9,\ldots, 20\}$ & $\{15,\ldots, 20\}$\\
    \bottomrule
  \end{tabular}
\end{table}

We simulate the responses of the systems in the datasets to a uniform white
noise input in the interval $\sbr{-1,1}$. We collect $N=200$ samples of the
output, from zero initial conditions, and we estimate the static nonlinearity
and the impulse response.

As described in Section~\ref{sec:hammersteins}, the linear blocks are modeled
as zero-mean Gaussian processes with first order stable-spline kernels. We
consider both a parametric model and a nonparametric model for the static
nonlinearity.  All the methods are initialized with $\rho_1 = 1$, $\rho_2=
0.6$. The noise variances are initialized from the prediction error of an
overparameterized least-squares estimate
(see~\cite{bai1998optimal,risuleo2015new}).

In the simulation, we compare the performance of the following estimators:

\begin{description}
  \item[H-P] A semiparametric model for the Hammerstein
    system. It uses the Legendre polynomial basis to construct a linear
    parameterization (with the correct order) of the input:
    \begin{equation}
      \mu_w(\theta)=\Phi\theta,\qquad \sbr[1]{\Phi}_{i,j} = \varphi_j(u_i).
    \end{equation}
    The dynamical system is modeled as a zero-mean Gaussian process with
    covariance matrix given by the first order stable-spline kernel.
  \item[H-MCEM] A nonparametric model for the Hammerstein system with Gibbs
    sampling from the joint posterior with $B=200$ and $M=500$. It uses the
    radial-basis-function kernel~\eqref{eq:rbf} to model the input
    nonlinearity. Note that, because the Hammerstein system is not
    identifiable, we fix $\theta_1=1$ in the algorithm.
  \item[H-VBEM] A nonparametric model for the Hammerstein system with
    variational-Bayes approximation of the joint posterior. It uses the same
    kernel as H-MCEM to model the input  nonlinearity.
  \item[NLHW] The parametric model in Matlab with the default parameters. It
    corresponds to the maximum-likelihood estimator of the model with the
    correct parameterization.
\end{description}
In all methods, the EM iterations are stopped once the relative change in the
parameter values is below $10^{-2}$.

To evaluate the performance of the methods, we use the standard goodness-of-fit
criterion~\eqref{eq:goodnessoffit} for the impulse response of the linear
system. For the input nonlinearity, we compute the
estimated value $\hat w_j$ on a uniform grid of 300 values between -1 and 1 and
we compare it to the true value $w_j$ according to
\begin{equation}
  \mathrm{Fit}^f_j = 1 - \frac{\enVert{w_j - \hat w_j}}{\enVert{w_j -
  \mathrm{mean}(w_j)}},
\end{equation}
where $w_j$ is the vector of values of the true nonlinearity ot the $j$th Monte
Carlo run, and $\hat w_j$ is an estimate of the same vector of values.

\begin{figure*}[htb]
  \centering
\includegraphics{boxplot_hammerstein.pdf}
  \caption{Boxplot of the estimation result}\label{fig:boxplot_hammerstein}
\end{figure*}

The result of the experiment are presented in
Figure~\ref{fig:boxplot_hammerstein}. The figure shows the boxplots of the fit
of the estimated impulse responses (upper pane) and of the static
nonlinearities (lower pane) over the systems in the datasets.

From this simulation, it appears that the proposed nonparametric models are
capable of recovering the system better than the fully parametric NLHW\@. In
addition, it appears that using the correct parametric model for the
input nonlinearity is beneficial in terms of accuracy. As was the case in the
cascaded-system estimation problem, the two approximation methods have comparable performance.

\section{Conclusions}
In this work, we have proposed a new model structure, which we have called the
\emph{uncertain-input model}. Uncertain-input models describe
linear systems subject to inputs about which we have limited information. To
encode the information we have available about the input and the system, we
have used Gaussian-process models.

We have shown how classical problems in system identification can be seen as
uncertain-input estimation problems. Among these applications we find classical
PEM, errors-in-variables and blind system-identification problems,
identification of cascaded linear systems, and identification of Hammerstein
models.

We have proposed an iterative algorithm to estimate the uncertain-input model.
We estimate the impulse response of the linear system and the input
nonlinearity as the posterior means of the Gaussian-process models given the
data. The hyperparameters of the Gaussian-process models are estimated using
the marginal-likelihood method. To solve the related optimization
problem, we have proposed an iterative method based on the EM method.

In the general formulation, the model depends on the convolution of two
Gaussian processes. Therefore, the joint distribution of the data is not
available in closed form.  To circumvent this issue, we have proposed specialized
models, namely the semiparametric and the parametric models, for which the
integrals defining the posterior distributions are available. In the more
general case, we have proposed two approximation methods for the joint
posterior distribution. In the first method, we have used a particle
approximation of the posterior distribution. The particles are drawn using the
Gibbs sampler from Gaussian full-conditional distributions. In the second
method, we have used the variational-Bayes approach to approximate the
posterior distribution. Using a mean-field approximation, we have found that
the posterior distribution can be approximated as a product of two independent
Gaussian random variables.

We have tested the proposed model on two problems: the estimation of cascaded
linear systems and of Hammerstein models. In both cases, the proposed
uncertain-input formulation is able to capture the systems and to provide good
estimates.

Although hinged on the EM method (which is guaranteed to converge under certain
smoothness assumptions) the approximate methods we have proposed do
not have general convergence guarantees: in the formulation given
by~\eqref{eq:ui_model}, there may instances of uncertain-input models for which
the assumptions required for convergence may not hold.  In future publications,
we plan to analyze whether there exists general conditions on the
uncertain-input model such that the algorithms are guaranteed to converge to
optimal solutions.

In addition, the uncertain-input model can be nonidentifiable in certain
configurations (for instance, consider the general errors-in-variables
problem). We plan to further explore this nonidentifiability. Connections with
other problems sharing the same bilinear
structure~\cite{bai2005least,wang2009revisiting} outside of the system
identification framework are also under investigation.

redpajama_set_name: RedPajamaArXiv


A police investigation is under way after a toddler was allegedly left unattended inside a bus for six hours in southern New South Wales.
The child, who 9News is not naming, was taken to hospital for treatment for dehydration.
The two-year-old boy had been picked up by the Yass Little Learners Centre bus as part of a routine pickup and drop-off service.
But the boy wasn't taken off the bus when it arrived at the daycare centre.
He was discovered by the driver after being left through the day strapped into a child restraint, unattended inside the vehicle.
The bus was parked on hot concrete for some six hours as temperatures in Yass hit 28 degrees.
The toddler was treated for dehydration by paramedics and taken to the nearby Yass Hospital.
The Centre's Director says the staff were left in shock and steps had been taken already to make sure it didn't happen again.
It is understood the toddler was discharged from the hospital a few hours later.
The driver has been suspended and the matter reported to child and community services. A police investigation is now under way.
Investigators are appealing for information and urge anyone with information to come forward.
Meantime the Centre remains open - with other staff from the centre doing the bus run.


Our original event space located in the back of our tasting room, The Krausen Parlour is ideal for smaller parties and get-togethers.
Notice any off flavors? Have a problem with a can, bottle, or packaging? Find outdated beer in the marketplace or anything else unusual? Tell us about it here so we can work to fix it!


text:
Leptasterias arctica is een zeester uit de familie Asteriidae.

De wetenschappelijke naam van de soort werd in 1885 gepubliceerd door Murdoch.

Forcipulatida

redpajama_set_name: RedPajamaWikipedia


Make Money Specialist is a blog site created to filter and present useful information about making money online in one place. This blog site includes reviews, ideas, tips and many other ways to earn money online fast free and easy.
Work in the UK as a respite nurse or home carer. Jobs for people to live in as home nurses or disability carers. Get paid to travel to the UK.
Providing information on The Present and Future Data-Entry Home Jobs, and getting off Unemployment, SSI, Disability, and any 9-5 Job! Bringing Independence, Dignity, and Respect into the lives of people who would not experience these virtues, living hand to foot, or paycheck to paycheck. This is the goal!
Free ad posting site which Upgrade your existing business throughout worldwide, easily.
M.S.I provides various waterproofing solutions like basement waterproofing, waterproofing products. You must try MSI perfect waterproofing products that will be fit for your need and best for basement waterproofing.
Do You Think Attraction Marketing really works in Network Marketing?
Attraction Marketing has become a buzz-word over the last couple of years in the Network Marketing industry, becoming the new standard for online business building to people in all kinds of different MLM companies.
Color Life Painters are local painting contractors with a passion for color. We strive for excellence in all our painting services. As a pre-screened, customer referred paint contractor, we will complete your project for a fair price in a timely manner with the utmost attention to detail.
RealWaysToEarnOnline.com is dedicated to all of those who wish to make money online. Herein you will find hints, tips and suggestions for realistic ways to make money online right from the comfort of your home.
Almost every evening party and awards sceremony, some celebrities will walk the red carpet for their selection, which makes one fashion icon replica Louboutin. Identical copies of replica Louboutin shoes Louboutin shoes look like a real designer shoes.
Original oil paintings for sale online buy landscape portrait impressionist and modern oil paintings.
Offline Arbitrage Review and Bonus. Read Offline Arbitrage right now to get secret insider information on what this course is all about.
Detailed information about the best affiliate programs on the internet. Find the best affiliate niches and the programs that pay the highest commission. A must read for everyone who wants to earn money at home.
Home Based Ceramic Business for Women.
Learn about how to start your own home based business. You can now earn extra money by selling ceramic crafts. Smoothclay.com offers you a great $349 home-based business opportunity by women and for women.
At Niche Profit Classroom, we teach the process of niche marketing, but specifically for the Internet. You'll learn how to find a 'hungry' niche and build profit-producing websites selling your own products or affiliate products.
StudyTube : Online Income | Online Learning is taking off around the World!
We actually pay you to study or teach in our proprietary distribution system, with great online profits everyone participates. We make the business of educating yourself profitable. Learn about our dynamic fractal matrix method to online income and IQ. Become a Student at StudyTube and get paid up to $390.00 per month. An Instructor possibly earns $1,830.00 See how this is possible, StudyTube.biz - Tired of uploading videos to Youtube for free? We pay you to upload or to just watch!
Legit Paying Paid To Click Sites with regular Updates. Earn free money from Internet working from Home. All Sites Are free to join. Best Paid To Click, Paid To Read Sites List.


text:
Q: Can conjugate subgroups be like this? Let G,H,K be three different subgroups of some group and are conjugate with one another. Will the situations below ever happen?
A. The intersection of G and H is a non trivial group but the intersection of H and K is trivial. More generally, can the intersection of each two of them be of different orders?
B. The intersection of G and H, H and K, K and G are non trivial groups of the same order but are different with one another.
Thank you for your time.

redpajama_set_name: RedPajamaStackExchange


Robert (Seymour, CT) is an artist/educator/author and developer of materials for artists. His work can be seen in numerous books, publications, museum and private collections around the world. He has taught workshops around the world and maintains an active teaching studio.


text:
Q: How to count co-occurrence of phenomena from Data with repeated values indicating overlap I am currently looking at the co-occurrence of various phenomena (gestures, intonation in speech) in time. As such, the data appears with each variable as it's own column, and phenomena are shown as repeating values while they are co-occuring, as :
Begin Time        End Time       g-phasing   apex   syllable    words   tones
00:00:02.000    00:00:04.266                         Zia        j'avais
00:00:04.266    00:00:05.390    Preparation          Zia        j'avais
00:00:05.390    00:00:05.519    Preparation           vE        j'avais
00:00:05.519    00:00:05.852    Preparation           vE        j'avais     H*
00:00:05.852    00:00:05.910    Preparation           de        des
00:00:05.910    00:00:05.970    Preparation           de        des
00:00:05.970    00:00:06.236    Preparation           de        des
00:00:06.236    00:00:06.276    Preparation           di        dizaines
00:00:06.276    00:00:06.650    Preparation           di        dizaines
00:00:06.650    00:00:06.795    Preparation          zEn        dizaines
00:00:06.795    00:00:06.835    stroke               zEn        dizaines
00:00:06.835    00:00:07.480    stroke               zEn        dizaines
00:00:07.480    00:00:07.630    stroke        apex   zEn        dizaines
00:00:07.630    00:00:07.857    stroke               zEn        dizaines    H*
00:00:07.857    00:00:08.080    stroke               zEn        dizaines
00:00:08.080    00:00:08.120    stroke             ddeux        de
00:00:08.120    00:00:08.226    Preparation        ddeux        de
00:00:08.226    00:00:08.290    Preparation        ddeux        de
00:00:08.290    00:00:08.900    Preparation           sy        sujets
00:00:08.900    00:00:12.396    Preparation           sy        sujets
00:00:12.396    00:00:12.410    stroke                sy        sujets
00:00:12.410    00:00:12.628    stroke                ZE        sujets
00:00:12.628    00:00:12.776    stroke        apex    ZE        sujets
00:00:12.776    00:00:12.924    stroke                ZE        sujets
00:00:12.924    00:00:12.990    stroke                ZE        sujets      H*
00:00:12.990    00:00:13.400    stroke                ZE        sujets

This dataset shows that there are two strokes (one from 00:00:06.795 to 00:00:08.120, and a second one from 00:00:12.396 to 00:00:13.400)
Ideally I would like to be able to count the number of strokes in the dataset, determine how many overlap with a pitch accented syllable (here, the "H*" value in the "tones" column that correspond to the syllables "zEn" and "ZE"), how many do not co-occur with a pitch-accented syllable, etc.
I'm not sure if I should iterrate over rows and create counters, if I should make use of the begin and end times, or if I should restructure the data.. Any help would be greatly appreciated!

A: Probably helpful: Get statistics for each group (such as count, mean, etc) using pandas GroupBy?
Closely looking at the dataset, sounds like a range aggregation problem is embedded. Once you filter for rows containing stroke data only using something like: How to select rows from a DataFrame based on column values? I recommend extracting out the first two columns into an isolated dataframe/list and work on generating a list of merged intervals.

initialize output_dataframe_list as an empty list
For each (start, end) item in merged_intervals_list:
do
        extract the subset of rows that have column 1>= start and column 2<=end and have the word stroke.
        use a groupby aggregation method on it and obtain your needed dataframe.
        append this dataframe containing the results of the interval into output_dataframe_list
concatenate all the dataframes together to create a new dataframe using pd.concat


What Individual Federal Income Tax Form Should I Use?
What Documents Do You Need to File a Form 1041?
For most taxpayers, federal income tax is filed through the standard 1040 forms. There are four standard 1040 forms that apply to most taxpayers: 1040EZ, 1040A, 1040 and the 1040NR. Depending on your particular financial circumstances one form may be more advantageous than the rest. Speaking with competent tax advisors or doing your own research is important to ensure you receive all the tax credit and deductions to which you are entitled.
While IRS Form 1040EZ is, by the far, the easiest form to use when completing your tax return, the standard 1040 form offers more opportunities for itemizing deductions. Nonresident aliens can use Form 1040NR to complete their tax returns.
The federal income tax form 1040EZ is the simplest form to use to file your federal tax return. According to the Internal Revenue Service, not everyone qualifies to file with the 1040EZ. The 1040 EZ is only available for taxpayers who are either married and filing a joint return or single and under 65 years of age. In addition, filers using this form should not have any dependents and income must be under $100,000. As of 2018, you may be able to claim the Earned Income Credit, though in many cases a 1040 form is required instead.
If you qualify to claim certain credits but don't want to itemize deductions, you may be able to use form 1040A. The IRS criteria for form 1040A includes taxable income under $100,000 that originates from traditional wages, interest, capital gains scholarships, any form of pension distributions, unemployment or distributions from the Alaska Permanent Fund. In addition, form 1040A is ideal if your claim any one of the following income adjustments: IRA contribution deduction, student loan interest, tuition or educator expenses. You may also use form 1040A if you claim certain credits like child tax, education and earned income credit.
The traditional 1040 is the standard federal income tax form that fits most tax situations. A standard 1040 form is the appropriate form to file for a refund of federal tax withholding on your income from a sale of property or other income subject to income tax. The IRS has minimal standards to help taxpayers determine if filing a 1040 is appropriate. If you have certain types of income, such as self-employment income from a farm or business, unreported tips, investment dividends or income from a trust or estate, you must use form 1040. An income greater than $100,000 also requires you to use 1040. Additional reasons for using 1040 include itemizing deductions, claiming certain tax credits or owing household employment taxes.
What Is the Difference Between a 1041, 1040A & 1040EZ?
The IRS requires nonresidents of the U.S. to file a federal income tax return with form 1040NR. This form is required for certain taxpayers under certain conditions. According to the IRS, form 1040NR is required if you are involved in business in the U.S. If sufficient tax is not withheld during the year from employment or business, you must file a 1040NR to settle the outstanding tax. Representatives of deceased persons or estates required to file a 1040NR must do so for the tax year in which the income is generated.
IRS: Which Form – 1040, 1040A, or 1040EZ?


PRICE: 5800 SEK per person, 2 persons sharing one snowmobile. 7600 SEK per person, snowmobile single use.
DEPARTURE: Saturday every week. 10 Jan - 12 April. Other dates possible against request with min 4 participants.
MEALS: Full board from lunch day 1 to lunch day 2. Information on diet (vegetarianism, pescetarianism, allergies etc) needed at time of booking.
INCLUDED: Thermo overall, boots, mittens, hat and helmet. Warm sleeping bag and towel.
NOT INCLUDED: Travel insurance, international and domestic fare, alcoholic beverages.
TOUR CONDITIONS: Min age 18 yrs. General good health and fitness required. Min age to drive own snowmobile 18 yrs and a valid driving licence for car. An excess of a maximum SEK 8,000 per incident/ snowmobile will be charged to the person who was driving the damaged snowmobile.
ACCOMMODATION: Dormitory in wilderness hut 1n.
MEETING POINT: : 7.30 p.m. in Kiruna City/Jukkasjärvi.
GROUP SIZE: Minimum 2 persons, maximum 8 persons per departure.
INCLUDED: Warm thermo overall, warm shoes, gloves and hat, helmet. Transfer from Kiruna/Jukkasjärvi and back again after the tour.
LEADERSHIP: English or Swedish language.
TOUR CONDITIONS: Minimum age to drive your own snowmobile: 18 years and valid driving licence for car. An excess of a maximum SEK 8,000 per incident/snowmobile will be charged to the person who was driving the damaged snowmobile.


text:
Lithuania: Undiscovered Eastern Europe
By A. Gracia, July 19, 2010
Still off the beaten path, Lithuania is one of the three Baltic republics that gained their independence after the fall of the Soviet Union. Its lush forests, beautiful national parks, the medieval Kaunas and Vilnius city centers (both UNESCO World Heritage Sites), and rich historic heritage make this northeastern European country of 3.5 million inhabitants an attractive destination.
Roughly the size of Ireland, Lithuania is bordered by Latvia to the north, Belarus to the east and south, Poland to the south and the Kaliningrad region of Russia to the southwest. The Baltic Sea kisses the surprisingly lovely 99-km stretch of coastline to the west, covered with majestic pines and long beaches with abundant Caribbean-like sand.
The country offers vast expanses of undisturbed forests and a tradition of nature conservancy that can be traced back to pre-Christian times. In this last European country to renounce Paganism in order to embrace Christianity, you'll find odd behavior (such as embracing trees in order to gain strength) that could perhaps be traced back to these beliefs. Moreover, Lithuania's former status as the Soviet Union's breadbasket has evolved into a rich agricultural tradition that ensures fresh produce and tasty local delicacies are available on every corner.
Trakai is a picturesque little town surrounded by several lakes, 30 km west of Vilnius. The castle on Lake Galvė is a beautiful sight, surrounded by water and a source of symbolism – it was once the site of the Lithuanian monarchy and fell only once to enemy attack. It later served as a residence for Lithuania's grand dukes and today is a popular museum of medieval Lithuanian history as well as a stage for concerts, festivals and films. The town also offers wonderful surroundings to explore, such as the Užtrakai and Lentvaris estates with their lovely parks and views. The countryside in and around Trakai provides opportunities for horse riding, swimming, diving, boating and yachting in the summer.
In Lithuania you can have your picture taken in the official center of Europe, determined by the French National Geographic Institute (and recognized by the European Union and the Guinness Book of Records) to be 26 km to the north of Vilnius, close to the village of Purniškės.
Lithuania's capital, Vilnius, offers a well-preserved old town and a Bohemian atmosphere close to its university centers. Kaunas, also the Lithuanian capital at different points in time, is still its second-largest city and another historic, architectural, intellectual and economic center.
Klaipėda, Lithuania's third-largest city, is located along the coast and is noted for being the northernmost ice-free port in the Baltic states. Its architecture departs from traditional Lithuanian styles and seems influenced by the country's northern neighbors. Here you'll find an old town with quaint, reconstructed merchants' houses and converted warehouses with their original exterior timber beams. The coastal towns are cherished by locals for their white sandy beaches reminiscent of Caribbean beaches rather than Nordic countries, as well as for their abundant spas and resorts.
If you visit during the winter, opportunities for cross-country skiing and ice skating abound. (If you're used to warmer climates, you might enjoy stepping indoors more often to enjoy a warm drink or meal.)
The Berneliu Uzeiga restaurants offer not only a meal, but an experience inspired by Lithuanian folk traditions – not just in their extensive menu, but in details such as the pottery in which many dishes are served, the live folk music, and employees in traditional dress. At the time of my visit there were five locations throughout the country, but you can find the latest information at berneliuuzeiga.lt.
corozalwien Oct. 20, 2010 @ 9:49 p.m.
I loved the stories by this writer. Why hasn't he or she published any more? This person must travel a lot, as stories go from the Caribbean and Latin America to everywhere in Europe and back. They make for nice escapist fantasies without exaggerated romanticism
corozalwien July 14, 2011 @ 9:42 p.m.

redpajama_set_name: RedPajamaCommonCrawl


text:
Q: Trouble wrting to Excel I' am new to Python and trying to write into a merged cell within Excel. I can see the data that is already stored within this cell/row, so I know its there. However when I try to overwrite it nothing happens.
I have tried messing with the index and header as well but nothing seems to work.

import pandas as pd
from openpyxl import load_workbook
Read the excel file into a pandas DataFrame
df = pd.read_excel(file here', sheet_name='Sheet1')
print(df.iloc[8, 2])
Make the changes to the DataFrame
df.iloc[8, 2] = "Bob Smith"
Load the workbook
book = load_workbook(file here)
writer = pd.ExcelWriter(file here, engine='openpyxl')
writer.book = book
Write the DataFrame to the first sheet
df.to_excel(writer, index=False)
Save the changes to the Excel file
writer.save()

A: import pandas as pd
from openpyxl import *
file="C:/Users/OneDrive/Bureau/draftExcel.xlsx"
df = pd.read_excel(file,sheet_name='sheet1')
df.iat[5,0]='cell is updated'
print(df) # to check first in the terminal if the content of the cell is updated
book=load_workbook(file)
writer=pd.ExcelWriter(file, engine='openpyxl')
df.to_excel(writer,sheet_name='sheet1',index=False)
writer.close()

I tried to make an example from what you explained because you didn't show your code, so I hope it was helpful.
Instead of using .iloc I used .iat so you can update the data in a specific cell in your DataFrame using column_index instead of column_label.
Remember that the Excel file you are working on must be closed while you are editing data with python, if it is open you will get an error.