You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


---
abstract: 'In this paper we introduce a novel method for linear system identification with quantized output data. We model the impulse response as a zero-mean Gaussian process whose covariance (kernel) is given by the recently proposed stable spline kernel, which encodes information on regularity and exponential stability. This serves as a starting point to cast our system identification problem into a Bayesian framework. We employ Markov Chain Monte Carlo methods to provide an estimate of the system. In particular, we design two methods based on the so-called Gibbs sampler that allow also to estimate the kernel hyperparameters by marginal likelihood maximization via the expectation-maximization method. Numerical simulations show the effectiveness of the proposed scheme, as compared to the state-of-the-art kernel-based methods when these are employed in system identification with quantized data.'
address:
- |
    Department of Electrical Engineering, Eindhoven University of Technology, Eindhoven, The Netherlands\
    (e-mail: <PRESIDIO_ANONYMIZED_EMAIL_ADDRESS>)
- 'Automatic Control Lab and ACCESS Linnaeus Centre, School of Electrical Engineering, KTH Royal Institute of Technology, Stockholm, Sweden (e-mail: <PRESIDIO_ANONYMIZED_EMAIL_ADDRESS>)'
- 'Department of Information Engineering, University of Padova, Padova, Italy (e-mail: <PRESIDIO_ANONYMIZED_EMAIL_ADDRESS>)'
author:
- Giulio Bottegal
- Hakan Hjalmarsson
- Gianluigi Pillonetto
bibliography:
- 'biblio.bib'
title: |
    A new kernel-based approach to system identification\
    with quantized output data
---

,

,

System identification; kernel-based methods; quantized data; expectation-maximization; Gibbs sampler

Introduction
============

Many applications in communications, control systems, bioinformatics, require modeling and prediction of dynamic systems with quantized output data (see e.g. [@bae2004gene], [@carli2010quantized] and [@wang2010system]). In particular, in this paper we assume that an unknown input-output map is defined by the composition of a linear time-invariant dynamic system and a quantizer, as depicted in Figure \[fig:block\_scheme\].

0.2in

![Block scheme of a system with quantized output data.[]{data-label="fig:block_scheme"}](block_scheme_journal){width="\columnwidth"}

-0.2in

Estimation of this type of structure is a challenging task. Standard system identification techniques, such as least-squares or the Prediction Error Method (PEM) [@Ljung], [@Soderstrom], may give poor performances, because the presence of the quantizer can determine a significant loss of information on the system dynamics. Recently, there has been an increasing interest on developing new methods to solve this problem. Particular attention has been devoted to the case of binary measurements [@wang2003system], [@wang2006joint], also studying on-line recursive identification [@jafari2012convergence], [@guo2013recursive]. In other works, e.g. [@colinet2010weighted], the knowledge of a dithering signal is exploited to improve the identification performance. The problem of experiment design is analyzed in [@casini2011input], [@casini2012input] and [@godoy2014novel]. More recently, system identification with non-binary quantizers has been also addressed in [@marelli2013identification], [@moschitta2015parametric], [@guo2015asymptotically], [@wang2015identification]. In [@godoy2011identification], [@godoy2014identification], [@chen2012impulse], the system identification problem is posed as a maximum likelihood/a-posteriori problem. In particular, the framework proposed in [@chen2012impulse] uses the nonparametric kernel-based identification approach proposed in [@SS2010], see also [@pillonetto2014kernel] for a survey.

Similarly to [@chen2012impulse], in this paper we cast the problem of identifying a linear dynamic systems with quantized data in a Bayesian framework. To this end, we model the impulse response of the unknown system as a realization of a zero-mean Gaussian random process. The covariance matrix (in this context also called a *kernel*) corresponds to the recently introduced *stable spline kernel* (see [@SS2010], [@SS2011], [@bottegal2013regularized]). The structure of this type of kernel depends on two *hyperparameters*, that need to be estimated from data. Following an empirical Bayes approach [@Maritz:1989], the hyperparameters, together with the noise variance, are estimated via marginal likelihood maximization. These quantities are then used to compute the minimum mean-square estimate (MMSE) estimate of the impulse response.

The whole procedure, which has been proven effective and relatively simple in the standard (non-quantized) setting, becomes more involved in the scenario under study. The MMSE system estimate requires the computation of an analytically intractable integral. To accomplish this task, in this paper we propose sampling methods based on Markov Chain Monte Carlo (MCMC) techniques. In particular, we design two different sampling techniques using the so called Gibbs sampler [@geman1984stochastic], which often enjoys faster convergence properties compared to standard Metropolis-Hastings sampling techniques [@Gilks].

Another contribution of the paper is to show that the task of estimating the kernel hyperparameters (and the noise variance) can be accomplished using the proposed sampling techniques. Using the Expectation-Maximization (EM) method [@dempster1977maximum], we design an iterative scheme for marginal likelihood maximization, where the E-step characterizing the EM method makes use of the Gibbs sampler (see also [@casella2001empirical]), and the M-step results in a sequence of straightforward optimization problems. Interestingly, the resulting estimation scheme can be seen as a generalization of the method proposed in [@godoy2014identification] to the Bayesian nonparametric framework and, differently from [@bottegal2015bayesian] and [@chen2012impulse], it allows tuning all the kernel hyperparameters.

The paper is organized as follows. In the next section, we introduce the problem of identifying dynamic systems from quantized data. In Section \[sec:Bayesian\_model\], we formulate the proposed Bayesian model and discuss its inference. In Section \[sec:sysid\_method\], we describe the proposed identification method. Section \[sec:experiments\] shows the results of several numerical experiments to assess the performance of the proposed method. Some conclusions end the paper.

Problem formulation {#sec:problem_formulation}
===================

We consider a linear time-invariant discrete-time dynamic system of the form $$\label{eq:sys1}
z_t = \sum_{i=1}^{+\infty} g_i u_{t-i} + e_t \,,$$ where $\{g_t\}_{t=1}^{+\infty}$ is a bounded-input-bounded-output (BIBO) stable impulse response representing the dynamics of the system. We approximate the impulse response by considering the first $m$ samples only, namely $\{g_t\}_{t=1}^{m}$, where $m$ is assumed large enough to capture the system dynamics[^1]. The input $u_t$ is a measurable signal; for the sake of simplicity, we assume that the system is at rest until $t=0$, meaning $u_t = 0,\,t<0$. We shall not make any specific requirement on the input sequence (i.e., we do not assume any condition on persistent excitation in the input [@Ljung]), requiring only $u_t \neq 0$ for some $t$. Notice, however, that even though the algorithm do not require any specifics of the input, the resulting estimate is of course highly dependent on the properties of the input [@Ljung]. The output $z_t$ is corrupted by the additive zero-mean Gaussian white noise $e_t$ with variance $\sigma^2$, assumed unknown. Introducing the vectors $$g := \begin{bmatrix} g_1 \\ \vdots \\ g_{m} \end{bmatrix} ,\,  \varphi_t^T := [ u_{t-1} \, \ldots \, u_{t-m} ] \,,$$ an approximation of can be written in linear regression form, namely $$\label{eq:sys11}
z_t = \varphi_t^T g + e_t \,.$$ The output $z_t$ is not directly measurable, only a quantized version being available, namely $$y_t = \mathcal Q [z_t] \,,$$ where $\mathcal Q$ is a known map of the type $$\mathcal Q[x] = s_k \qquad \mbox{if } x \in (q_{k-1},\,q_k] \,,$$ with $s_k \in \{s_1,\,\ldots,\,s_Q\}$ and $q_k \in  \{q_0,\,\ldots,\,q_Q\}$ (and typically $q_0 = -\infty$ and $q_Q = +\infty$).

A particular and well-studied case is the binary quantizer, defined as $$\mathcal Q[x] = \left\{ \begin{array}{ll} -1 & \mbox{if } x < C \\
                                  1 & \mbox{if } x \geq C
        \end{array}  \right. \,.$$ It is well-known that a condition on the threshold to guarantee identifiability of the system is $C \neq 0$. In fact, when $C=0$, the system can be determined up to a scaling factor [@godoy2011identification].

We assume that $N$ input-output data samples $y_1,\,\ldots,\,y_{N}$, $u_0,\,\ldots,\,u_{N-1}$ are collected during an experiment. The problem under study is to estimate $\{g_t\}_{t=1}^{m}$ using the collected measurements.

It is also useful to write the dynamics in the following vector notation $$\label{eq:sys2}
z = Ug + e\,,$$ where $z$ and $e$ are $N$-dimensional vectors collecting the samples of $z_t$ and $e_t$, respectively, and $U \in \mathbb{R}^{N \times m}$ is a matrix whose $t$-th row corresponds to $\varphi_t^T$. Similarly, we denote by $y:= [ y_1 \,\, \ldots \,\, y_N ]^T$ the vector collecting the available quantized measurements.

Bayesian modeling and inference {#sec:Bayesian_model}
===============================

Impulse response prior
----------------------

In this paper we cast the system identification problem into a Bayesian framework, setting an appropriate prior on $g$. Following a Gaussian regression approach [@Rasmussen], we model the impulse response as a zero-mean Gaussian random vector, i.e. we assume the following prior distribution $$\begin{aligned}
\label{eq:model_gh}
p(g;\,\lambda,\,\beta) \sim \mathcal N (0,\,\lambda K_\beta) \,. 
$$ Here, $K_\beta$ is a covariance matrix whose structure depends on the *shaping hyperparameter* $\beta$, and $\lambda \geq 0$ is a scaling factor. In this context, $K_\beta$ is usually called a [*kernel*]{} and determines properties of the realizations of $g$. In particular, we choose $K_\beta$ from the family of *stable spline kernels* [@SS2010], [@SS2011]. Such kernels are specifically designed for system identification purposes and give clear advantages compared to other standard kernels [@bottegal2013regularized], [@SS2010] (e.g. the Gaussian kernel or the Laplacian kernel, see [@Scholkopf01b]). Motivated by its maximum entropy properties [@chen2016maximum], in this paper we make use of the *first-order stable spline kernel* (or *TC kernel* in [@ChenOL12]). It is defined as $$\label{eq:ssk1}
\{K_\beta\}_{i,j} := \beta^{ \max(i,j)} \quad,\, 0 \leq \beta  <1 \,,$$ with $\beta$ regulating the impulse response decay rate.

Bayesian inference of systems with quantized output data {#sec:Gibbs}
--------------------------------------------------------

In this section we describe the complete Bayesian model that stems from the prior adopted for the impulse response. Our main aim is to devise an MCMC-based sampling mechanism for inferring the model. We introduce the vector of parameters $\eta = [\lambda \,\, \beta \,\, \sigma^2]$, seen as a deterministic quantity.

[![*Bayesian network describing the system model. Non-dashed nodes denote deterministic or given quantities; dashed nodes denote random variables.*[]{data-label="fig:bayesian_network"}](bayesian_network_FIR.pdf "fig:"){width="6cm"}]{}

Figure \[fig:bayesian\_network\] depicts the Bayesian network describing the system identification problem with quantized data. The kernel hyperparameters $\lambda$ and $\beta$ determine the stochastic model of $g$, which, together with $\sigma^2$, in turn determines the distribution of the $z_t$. The output measurements $y_t$ are determined directly by the non-quantized outputs. Any kind of inference on this network crucially depends on the capability of performing Bayesian inference on the hidden nodes $g,\,\{z_t\}_{t=1}^N$, by computing functionals of the posterior density $$\label{eq:target_density}
 p(g,\,z|y;\,\eta)$$ of the form $$\label{eq:functional_generic}
\int f(g,z) p(g,\,z|y;\,\eta) dg\,dz \,,$$ where $f(g,z)$ is a general function of $g$ and $z$. Special cases of that will be used later comprise the MMSE of $g$ given $y$ and the computation of the marginal likelihood of $\eta$.

In the remainder of the section, we assume that $\eta$ is fixed.

Method 1: sampling from the joint posterior {#sec:sampling1}
-------------------------------------------

Quite unfortunately, analytical computation of is intractable, due to the involved structure of the posterior density. Instead, we can resort to Markov Chain Monte Carlo methods, i.e. use the approximation $$\label{eq:MCMC_approx}
\int f(g,z) p(g,\,z|y;\,\eta) dg\,dz \simeq \frac{1}{M} \sum_{k=1}^M f(g^{(k)},\,z^{(k)})\,,$$ where $g^{(k)}$, $z^{(k)}$ are random samples drawn from and $M$ is a large integer. The problem is how to design an effective technique to sample from the distribution . If all the conditional probability densities of such a distribution are available in closed-form, the problem of sampling can be solved efficiently using a special case of the Metropolis-Hastings method, namely the Gibbs sampler (see e.g. [@Gilks]). The basic idea is that each conditional random variable is the state of a Markov chain; then, by drawing samples from the conditional probability densities iteratively, we converge to the stationary state of this Markov chain and generate samples of the desired distribution. The conditionals of are as follows.

1.  $p(z_i|g,\,\{z_j\}_{j \neq i},\,y;\,\eta),\,i=1,\,\ldots,\,N$. First note that from , given $g$, $z_i$ is independent of $\{z_j\}_{j \neq i}$ and $\{y_j\}_{j \neq i}$. Therefore we have $$p(z_i|g,\,\{z_j\}_{j \neq i},\,y;\,\eta)  = p(z_i|g,\,y_i;\,\eta) \,.$$ Now, using again we observe that, if $y_i$ were not given, then $$\label{eq:pz_gaussian}
    p(z_i|g;\,\eta) \sim \mathcal{N}( \varphi_i^T g ,\,\sigma^2) \,.$$ Knowing the quantized output $y_i$ permits to narrow the range of possible values of $z_i$. In particular, if $y_i = s_k$, for ant $k = 1,\,\ldots,\,Q$, we have $$\label{eq:truncated}
    p(z_i|g,\,y_i = s_k;\,\eta)  \sim \mathcal{N}_{q_{k-1}}^{q_k} ( \varphi_i^T g,\,\sigma^2)\,,$$ where $\mathcal N_{a}^{b}(\mu,\,\sigma^2)$ denotes a Gaussian distribution truncated below $a$ and above $b$, whose original mean and variance are $\mu$ and $\sigma^2$ respectively.

2.  $p(g|z,\,y;\,\eta)$. Given $z$, information carried by $y$ becomes redundant and can be discarded. Due to the assumption on the distribution of the noise $e$, the vectors $g$ and $z$ are jointly Gaussian, so that the posterior density of $g$ given $z$ is also Gaussian. Combining the linear model with the prior it is straightforward to obtain (see e.g. [@Anderson:1979]) $$\label{eq:post_g}
    p(g|z;\,\eta) \sim \mathcal{N}(m_g,\,P_g) \,,$$ with $$\begin{aligned}
    P_g & = \left( \frac{1}{\sigma^2}  U^T U + \left( \lambda_g K_\beta \right)^{-1}\right)^{-1} \label{eq:cov_g} \\
    m_g & = \frac{1}{\sigma^2} P_g  U^T z = H z \label{eq:mean_g} \,,\end{aligned}$$ with obvious definition of $H$.

Algorithm 1 summarizes the computation of using the sampling mechanism described in this subsection.

**Algorithm 1**: Method 1 for inference\
Input: $\{y_t\}_{t=1}^N,\,\{u_t\}_{t=0}^{N-1},\,\eta$\
Output: ${{\mathbb E}\,}[f(g,z)|y]$

Initialization: Compute initial value $g^{(0)}$\
For $k=1$ to $M+M_0$:

1.  Sequentially draw the samples $z_i^{(k)}$, $i=1,\,\ldots,\,N$, from $p(z_i|g^{(k-1)},\,y_i;\,\eta)$

2.  Draw the sample $g^{(k)}$ from $p(g|z^{(k)};\,\eta)$

Compute $\frac{1}{M} \sum_{k=M_0}^{M+M_0} f(g^{(k)},z^{(k)})$

In Algorithm 1, the parameter $M_0$ is introduced. It represents the number of initial samples to be discarded and is also known as *burn-in* [@meyn2009markov]. In fact, the conditional densities from which those samples are drawn are to be considered as non-stationary, because the associated Markov Chain takes a certain number of iterations to converge to a stationary distribution.

Method 2: sampling from the marginalized posterior of $z$ {#sec:sampling2}
---------------------------------------------------------

Method 1 for sampling allows for computing expected values (with respect to the posterior $p(g,\,z|y;\,\eta)$) of any kind of function $f(g,\,z)$. For the problem under study however, we are particularly interested in computing the expected value (conditional on $y$) of $g$, and of quadratic functions in $g$ and $z$, i.e. $$\label{eq:quadr_f}
f(g,z) = [g^T\,\,z^T] \begin{bmatrix}
A & B^T \\ B & C
\end{bmatrix} \begin{bmatrix} g \\ z \end{bmatrix} \,,$$ for given $A \in \mathbb R^{N \times N}$, $B \in \mathbb R^{m \times N}$, $C \in \mathbb R^{m \times m}$. These type of functions will play a central role in estimating the vector of parameters $\eta$ (see Section \[sec:sysid\_method\]).

The sampling method proposed in this subsection relies on the following result.

\[th:method2\] Let $f(g,\,z)$ be a quadratic function of the type . Then $$\begin{aligned}
 \label{eq:th_result}
{{\mathbb E}\,}[f(g,\,z)|y]& = {\mbox{\rm tr}}\{A P_g\}  \\ & + \int z^T(H^TAH + 2BH + C)z \, p(z|y;\eta)  dz \,, \nonumber\end{aligned}$$ where $P_g$ and $H$ are defined in and .

#### Proof: {#proof .unnumbered}

We have $$\begin{aligned}
& \int f(g,\,z) p(g,\,z|y;\,\eta) \,dg\,dz = \nonumber \\
& \int f(g,\,z) p(g|z;\,\eta) p(z|y;\,\eta) \,dg\,dz = \nonumber \\
& \int \left( \int f(g,\,z) p(g|z;\,\eta)\,dg \right) p(z|y;\,\eta) \,dz \,.  \end{aligned}$$ We now focus on the internal integral in the above expression. Developing $f(g,\,z)$ and recalling the first and second moments of $p(g|z;\,\eta)$, given in , we obtain the following terms: $$\begin{aligned}
\int g^TAg \, p(g|z;\,\eta)dg &= {\mbox{\rm tr}}\{AP_g\} +  m_g^TAm_g \\ & = {\mbox{\rm tr}}\{AP_g\} + z^TH^TAHz \,, \\
\int z^TBg \, p(g|z;\,\eta)dg& = z^T B m_g = z^T B H z \,,  \\
\int z^TCz \, p(g|z;\,\eta)dg & = z^T C z \,.\end{aligned}$$ Combining the three terms yields . $\blacksquare$

A analogous result holds for the computation of ${{\mathbb E}\,}[g|y]$.

\[th:method2\_2\] It holds that $$\begin{aligned}
 \label{eq:th_result_2}
{{\mathbb E}\,}[g|y]& = H \int z \, p(z|y;\eta) \, dz \,, \end{aligned}$$ where $H$ is defined in .

Lemmas and state that the posterior expectations of $g$ and of quadratic functions in $g$ and $z$, are equivalent to the expectation of specific functions with respect to the marginalized posterior of $z$ given $y$. This distribution corresponds to a multivariate truncated Gaussian; so, analytical computation of the relative integrals is quite challenging. Computing the integrals by sampling from $p(z|y;\,\eta)$ is instead a viable approach.

An effective technique for sampling from $p(z|y;\eta)$ is again based on the Gibbs sampler. The idea is to iteratively draw samples from the conditional densities $p(z_i|\{z_j\}_{j \neq i},\,y_i;\,\eta)$. Generating samples from these distributions is relatively simple, because they are scalar truncated normal distributions. To see this, first consider the marginal distribution $p(z;\,\eta)$, which is a multivariate normal random vector with zero-mean and covariance matrix $\Sigma_z = \lambda U K_\beta U^T + \sigma^2 I$. Then $$p(z_i|\{z_j\}_{j \neq i};\,\eta) \sim N(m_{z_i},\,P_{z_i})\,,$$ where $$\begin{aligned}
m_{z_i} & =  {\mbox{\rm Cov}\,}[z_i,\,\{z_j\}_{j \neq i}] {\mbox{\rm Var}\,}[\{z_j\}_{j \neq i}]^{-1} \label{eq:post_mean_z}\\
& \qquad\qquad \times [z_1\,\ldots z_{i-1}\,z_{i+1}\,\ldots\,z_N]^T \,, \nonumber \\
P_{z_i}& = {\mbox{\rm var}}[z_i] - {\mbox{\rm Cov}\,}[z_i,\,\{z_j\}_{j \neq i}] {\mbox{\rm Var}\,}[\{z_j\}_{j \neq i}]^{-1} \label{eq:post_cov_z} \\
& \qquad \qquad  \qquad \times {\mbox{\rm Cov}\,}[z_i,\,\{z_j\}_{j \neq i}]^T \nonumber \,.\end{aligned}$$ Therefore, $p(z_i|\{z_j\}_{j \neq i},\,y_i;\,\eta) \sim N_{q_{k-1}}^{q_k}(m_{z_i},\,P_{z_i})$.

Algorithm 2 summarizes the sampling procedure described in this subsection, for quadratic functions. The same algorithm can be used for computing the posterior expectation of $g$.

**Algorithm 2**: Method 2 for inference\
Input: $\{y_t\}_{t=1}^N,\,\{u_t\}_{t=0}^{N-1},\,\eta$\
Output: ${{\mathbb E}\,}[f(g,z)|y]$

Initialization: Compute initial value $g^{(0)}$\
For $k=1$ to $M+M_0$:

1.  Sequentially draw the samples $z_i^{(k)}$, $i=1,\,\ldots,\,N$, from $p(z_i|\{z_j\}_{j \neq i},\,y_i;\,\eta)$

Compute ${\mbox{\rm tr}}\{AP_g\} + \frac{1}{M} \sum_{k=M_0}^{M+M_0} z^{(k)T} (H^TAH + 2BH + C)z^{(k)}$

\[rem:sampling\] The quantities required to compute and can be retrieved from $\Sigma_z^{-1}$. Assume we are interested in generating a sample of $z_1$. To this end, partition $\Sigma_z^{-1}$ as follows $$\Sigma_z^{-1} = \begin{bmatrix}
s_1 & s_{j1}^T \\ s_{j1} & S_j
\end{bmatrix}\,,$$ where $s_1 \in \mathbb{R},\,s_{j1} \in \mathbb{R}^{N-1},\,S_{j} \in \mathbb{R}^{N-1 \times N-1}$. Then $P_{z_1} = 1/s_1$ and $m_{z_1} = -\frac{s_{j1}^T}{s_1}[z_2\,\ldots\, z_{N}]^T$ (see e.g. [@noble1988applied]). This procedure turns out computationally convenient, because it is required to invert $\Sigma_z$ only one time, instead of computing the inversion of $N$ matrices of size $N-1 \times N-1$ (which would be necessary for computing ${\mbox{\rm Var}\,}[\{z_j\}_{j \neq i}]^{-1}$ for any $i = 1,\,\ldots,\,N$).

System identification from quantized output data {#sec:sysid_method}
================================================

In the previous section we have defined an adequate Bayesian framework to describe the problem of identification of systems from quantized output data. We have also presented two methods for inference of functions of $z$ and $g$. In this section, we describe the proposed system identification procedure. It relies on the so called empirical Bayes approach [@Maritz:1989] and consists of the following two steps:

1.  Estimate the parameter vector $\eta$ via marginal likelihood maximization;

2.  Compute the MMSE estimate of the impulse response, using the estimated parameter vector $\hat \eta$.

In the following we analyze in detail the two steps composing the proposed system identification scheme.

Marginal likelihood estimation of the parameter vector $\eta$
-------------------------------------------------------------

In the empirical Bayes approach, the estimation of $\eta$ is tackled via marginal likelihood maximization, i.e., by computing $$\begin{aligned}
\hat \eta & = \arg \max_{\eta}  p(y;\,\eta) \label{eq:marg_lik} \\
& = \arg \max_{\eta} \int p(y,\,g,\,z;\,\eta)\, dg\,dz \,. \nonumber\end{aligned}$$ In the standard (non-quantized) scenario, solving is straightforward because the marginal likelihood admits a closed-form expression [@pillonetto2014kernel]. This does not hold in the problem under study. Therefore, to solve we propose an iterative solution scheme based on the EM method, where we treat $g$ and $z$ as *latent variables* that are iteratively estimated together with the parameter vector $\eta$. Let $\eta^{(n)}$ be the parameter estimate obtained at the $n$-th iteration of the EM method. Then, $n+1$-th update is obtained with the following steps:

(E-step)

:   Compute $$\label{eq:estep}
        Q(\eta,\,\hat \eta^{(n)}) := {{\mathbb E}\,}\left[\log p(y,\,g,\,z;\,\eta) \right] \,,$$ where the expectation is taken with respect to the posterior density $p(g,\,z|y;\,\hat \eta^{(n)})$, with $\eta$ fixed at the value $\eta^{(n)}$;

(M-step)

:   Compute $$\hat \eta^{(n+1)} = \arg \max_{\eta}  Q(\eta,\,\hat \eta^{(n)}) \,.$$

We first focus on performing the E-step. It requires the computation of , which corresponds to the integral $$\label{eq:compl_lik_integral}
\int \log p(y,\,g,\, z;\,\eta) p(g,\,z|y;\,\hat \eta^{(n)}) \,dg\,dz \,.$$

\[th:complete\_likelihood\] The complete likelihood admits the decomposition $$\begin{aligned}
-2 \log p(y,\,g,\,z;\,\eta) & = \frac{1}{\sigma^2} f_1(g,\,z) + \frac{1}{\lambda} f_2(g,\,z,\,\beta) \\ 
& \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!  + N \log \sigma^2 + \log \det \lambda K_\beta + \log p(y|z,\,g;\,\eta) \nonumber \,,\end{aligned}$$ where $f_1(g,\,z)$ is a quadratic function of the type , with $A = U^TU,\,B = U,\,C = I$, and $f_2(g,\,z,\,\beta)$ is also a quadratic function in $g$ and $z$, with $A = K_\beta^{-1},\,B = 0,\,C = 0$.

#### Proof {#proof-1 .unnumbered}

Using Bayes’ rule we can decompose the complete likelihood as follows $$\begin{aligned}
 \log p(y,\,g,\,z;\,\eta)&= \log p(y|z,\,g;\,\eta)  \\ & \quad+  \log p(z|g;\,\eta)  + \log p(g;\,\eta)  \,.\nonumber\end{aligned}$$ The term $p(z|g;\,\eta)$ is a vectorized version of , so that $$\begin{aligned}
 \label{eq:log_pz}
\log p(z|g;\,\eta)  =  -\frac{N}{2} \log \sigma^2 - \frac{1}{2\sigma^2} \|z - U g\|_2^2 \,.\end{aligned}$$ The last term on the right hand side gives $f_1(z,\,g)$. Similarly, $$\log p(g;\,\eta) = \log \det \lambda K_\beta + g^T(\lambda K_\beta)^{-1}g \,,$$ where the last term corresponds to $f_2(z,\,g,\,\beta)$. $\blacksquare$

The proposition reveals the nature of the complete likelihood, which is the summation of terms that are either constant or quadratic functions of $z$ and $g$, plus the term $\log p(y|z,\,g;\,\eta)$. As for this term, we note that $p(y|z,\,g;\,\eta)$ factorizes, each factor being of the type $$\label{eq:py_factor}
p(y_i|z,\,g;\,\eta)\! =\! \left\{ \begin{array}{ll}
    \! 1 & \mbox{if } y_i \! = \! s_k  \mbox{ and } z_i \! \in \! (q_{k-1},\,q_k] \\
    \! 0 & \mbox{otherwise} \end{array} \right.\! .$$ When computing the integral of this term using the sampling mechanisms introduced in Section \[sec:Bayesian\_model\], it is ensured that all the generated samples $z_i^{(k)}$ belong to the interval corresponding to the observed quantized value $y_i$. Hence, when we compute the expectation of $p(y|g,\,z;\,\eta)$ techniques of Section \[sec:Bayesian\_model\], it is ensured that each factor is always equal to 1 and thus $\log p(y|g,\,z;\,\eta) = 0$. Therefore, computing reduces to computing the expectation of two quadratic functions, and this can be done using both the sampling mechanisms introduced in Section \[sec:Bayesian\_model\]. We denote by $\hat f^{(n)}_1$ and $\hat f^{(n)}_2(\beta)$, respectively, the expected value of $f_1(g,\,z)$ and $f_2(g,\,z,\,\beta)$ at the $n$-th iteration of the EM method (i.e., when the $\eta^{(n)}$ of the parameter vector is available). Neglecting constant terms, we have $$\begin{aligned}
-2Q(\eta,\,\hat \eta^{(n)}) & = \frac{1}{\sigma^2} \hat f^{(n)}_1 + \frac{1}{\lambda}\hat f^{(n)}_2(\beta) \\
& \quad + N \log \sigma^2 + \log \det \lambda K_\beta \,. \nonumber\end{aligned}$$

\[th:m-step\] Let $$h(\beta) := m \log \hat f^{(n)}_2(\beta) + \log \det K_\beta\,.$$ Then the EM update $\eta^{(n+1)} = [\lambda^{(n+1)}\, \beta^{(n+1)}\, \sigma^{2(n+1)}]$ is obtained as follows: $$\begin{aligned}
\beta^{(n+1)} & = \arg \min_{\beta \in (0,\,1]} h(\beta)\,, \label{eq:betag_update}\\
\lambda^{(n+1)}& = m^{-1} \hat f^{(n)}_2(\beta^{(n+1)})   \label{eq:lambdag_update}\,, \\
\sigma^{2(n+1)}& = N^{-1} \hat f^{(n)}_1 \,. \label{eq:sigma_update}\end{aligned}$$

#### Proof {#proof-2 .unnumbered}

Differentiating $-2 Q(\eta,\,\eta^{(n)})$ with respect to $\lambda$, one obtains that the minimum, as a function of $\beta$, is achieved at $\lambda^* =  m^{-1} \hat f^{(n)}_2(\beta)$. Inserting this value into $-2 Q(\eta,\,\eta^{(n)})$ yields $h(\beta)$, so that and follow. Finally, it is straightforward to see that is the minimizer of $-2 Q(\eta,\,\eta^{(n)})$. $\blacksquare$

The M-step results in a series of computationally attractive operations. The kernel scaling hyperparameter $\lambda$ and the noise variance $\sigma^2$ admit a solution in closed form. The shaping hyperparameter $\beta$ is updated solving a simple scalar optimization problem. This problem can be further simplified by using a factorization of the first-order stable spline kernel, see [@bottegal2016robust] for details.

MMSE estimate of the impulse response
-------------------------------------

Having an estimate $\hat \eta$ available, the MMSE estimate of $g$ corresponds to $$\begin{aligned}
\hat g & := {{\mathbb E}\,}[g|y]  = \int g\, p(g|y;\,\hat \eta)\, dg \label{eq:int_1}\\
& = \int g\, p(g,\,z|y;\,\hat \eta) \,dg\,dz \,. \label{eq:int_2}\end{aligned}$$ Again, this quantity can be computed using both the sampling methods described in the previous section.

We summarize the overall procedure for system identification from quantized data in the following algorithm.

**Algorithm 3**: System identification with quantized output measurements\
Input: $\{y_t\}_{t=1}^N,\,\{u_t\}_{t=0}^{N-1}$\
Output: $\{\hat{g_t}\}_{t=1}^m$

Initialization: Set an initial value of $\hat \eta^{(0)}$\
Repeat until convergence:

1.  E-step: Compute $Q(\eta,\,\eta^{(n)})$ using Algorithm 1 or Algorithm 2

2.  M-step: update $\hat \eta^{(n+1)}$ according to Proposition \[th:m-step\]

Compute $\hat g$ using Algorithm 1 or Algorithm 2

The initial value $\hat \eta^{(0)}$ can be set either randomly or by using values of the kernel hyperparameters and noise variance returned by the standard nonparametric method for non-quantized data (see [@SS2010]).

Which sampling method?
----------------------

Algorithm 3 for quantized system identification works with both the sampling methods presented in Section \[sec:Bayesian\_model\], because it requires the computation of the posterior expectation of quadratic functions of the type and of $g$. Choosing the sampling method depends on the user requirements. Method 1 allows for any type of inference and therefore it can be used to compute useful statistics such as confidence intervals on the estimate of $g$ and quantile estimation. On the other hand, Method 2 requires sampling only from $p(z|y,\,\eta)$; if implemented properly (see Remark \[rem:sampling\]), this method is expected to have a faster convergence to the required expectations. The experiments presented in the next section show that the two sampling methods substantially give the same performance in terms of accuracy in reconstructing the true impulse response.

Numerical experiments {#sec:experiments}
=====================

We test the proposed technique by means of 2 Monte Carlo experiments of 100 runs each. For each Monte Carlo run, we generate an impulse response by picking 10 pairs of complex conjugate zeros with magnitude randomly chosen in $[0,\,0.99]$ and random phase and, similarly, 10 pairs of complex conjugate poles with magnitude randomly chosen in $[0,\,0.92]$ and random phase. The obtained impulse response is rescaled in order to have a random $\ell_2$-gain in the interval $[2,\,4]$. The goal is to estimate $m=50$ samples of this impulse response from $N$ input-output data. The input sequences are realizations of white noise with unit variance. We compare the following estimators.

-   **KB-GS-1**: This is the method described in this paper, making use of the sampling technique of Section \[sec:sampling1\] (Algorithm 1). The parameter $M$, denoting the number of samples generated by the sampler, is set to $100$ for each iteration of the EM method, and to $500$ for the final computation of $g$ (last step of Algorithm 3). The burn-in phase consists of $M_0 = 100$ samples. Convergence of the EM method is established by a threshold on the relative difference of the current and previous parameter estimates, i.e. the method stops if the condition $$\frac{\|\eta^{(n+1)} -\eta^{(n)}\|_2^2}{\|\eta^{(n)}\|_2^2} \leq 10^{-3}$$ is satisfied. In our simulations, we have verified that our choice of $M$ was adequate to guarantee a good degree of approximation (the reader interested in convergence diagnostics is referred to, e.g., [@gelman2014bayesian Ch. 11.4]).

-   **KB-GS-2**: This is the method described in this paper, making use of the sampling technique of Section \[sec:sampling2\] (Algorithm 2). The algorithm settings are the same as the previous method.

-   **KB-St.**: This is the standard nonparametric kernel-based method proposed in [@SS2010] and revisited in [@ChenOL12]. It is not designed for handling quantized data. The kernel adopted for identification is . The kernel hyperparameters are estimated by marginal likelihood maximization, while $\sigma^2$ is estimated via least-squares residuals.

-   **KB-Or.**: Same as KB, with the difference that in this case the vector $z$ is available to the estimator. Hence, this estimator is regarded as an Oracle.

-   **ML-GS**: This is the method proposed in [@godoy2011identification], based on maximum likelihood identification of the impulse response, namely $$\label{eq:MLFIR}
    \hat g = \arg \max_{g,\,\sigma^2}  p(y;\,g,\,\sigma^2) \,.$$ To solve the likelihood problem, a scheme based on the EM method is proposed. In our implementation, the E-step of the EM iterations is computed using the Gibbs sampler (in contrast, [@godoy2011identification] proposes a scenario-based approach). The length of $g$ is also set to 50, while the convergence is established by using the same conditions as the estimator KB-GS-1.

-   **MAP-GS**: This is the method proposed in [@chen2012impulse], based on maximum-a-posteriori (MAP) identification of the impulse response, namely $$\label{eq:MAPFIR}
    \hat g = \arg \max_{g}  p(y|g;\,\sigma^2) p(g;\,\lambda,\,\beta) \,,$$ where the prior $p(g;\,\lambda,\,\beta)$ corresponds to . Again, an EM solution scheme based on the Gibbs sampler is employed to solve this problem, setting $m=50$, while the convergence is established by using the same conditions as the estimator KB-GS-1. To facilitate hyperparameter and noise variance tuning, we plug those that are estimated by the method KB-Or., which has access to the non-quantized output $z$.

The performance of the estimators is evaluated by means of the fitting score, computed as $$FIT_i  = 1-\frac{\|g^i - \hat g^i \|_2}{\|g^i\|_2} \,,$$ where $g^i$ is the impulse response generated at the $i$-th Monte Carlo run and $\hat g^i$ the estimate computed by the tested methods.

Binary quantizer
----------------

The first experiment considers the following binary quantizer $$\mathcal Q_b [x] := \left\{ \begin{array}{ll} 1 & \mbox{ if } x \geq 1 \\ -1 & \mbox{ if } x < 1 \end{array} \right.\,.$$ For each Monte Carlo run, the noise variance is such that $\frac{{\mbox{\rm var}}(Ug) }{\sigma^2} = 10$, i.e. the ratio between the variance of the noiseless (non-quantized) output and the noise is equal to 10. If $\sigma^2 > 2$, the system is discarded. We generate $N = 500$ data samples.

[![*Box plots of the fitting scores for the binary quantizer experiment.*[]{data-label="fig:bin_exp"}](fig_ceil_500_2GS_rev1_swap-eps-converted-to "fig:"){width="\columnwidth"}]{}

Figure \[fig:bin\_exp\] shows the results of the Monte Carlo runs. The advantage of using the proposed identification scheme, compared to the method that does not account for the quantizer, is evident. Despite the large loss of information caused by the quantizer, the proposed method gives a fit which is quite comparable to the oracle method (KB-Or.). The proposed sampling mechanisms yield nearly equivalent performance. Furthermore, there seems to be a substantial advantage in using a Bayesian approach compared to the non-regularized estimator ML-GS, which is tailored for short FIR systems rather than IIR systems. The high number of coefficients to be estimated inevitably leads to high variance in the estimates. If this effect is not suitably alleviated by regularization (i.e., introducing a “good” bias), the performance of the estimator is doomed to be poor. Finally, the proposed estimators, based on computing the MMSE estimate of the impulse response, outperform the estimator MAP-GS, which is based on computing the MAP estimate of the impulse response.

Ceil-type quantizer
-------------------

In the second experiment we test the performance of the proposed method on systems followed by a ceil-type quantizer, which is defined as $$\mathcal Q_c [x] := \lceil x \rceil \,.$$ Again, for each Monte Carlo run, the noise variance is such that $\frac{{\mbox{\rm var}}(Ug) }{\sigma^2} = 10$. We generate $N = 200$ data samples.

As shown in Figure \[fig:ceil\_exp\], in this case, if one compares the oracle-type method (i.e. KB-Or.) with the same method using quantized data (KB-St.), the loss of accuracy is relatively low. This is because this type of quantizer has a mild effect on the measurements. It can be seen, however, that the proposed methods are able to give a fit that is comparable to the oracle KB-Or., regardless of the employed sampling method. We notice also that the Bayesian approach has a major impact on the performance, as seen by the gap in the performance between ML-GS and the other estimators.

[![*Box plots of the fitting scores for the ceil-type quantizer experiment.*[]{data-label="fig:ceil_exp"}](fig_bin_500_2GS_rev1_swap-eps-converted-to "fig:"){width="\columnwidth"}]{}

Conclusions
===========

In this paper, we have introduced a novel method for system identification when the output is subject to quantization. We have proposed a MCMC scheme that exploits the Bayesian description of the unknown system. In particular, we have shown how to design two integration schemes based on the Gibbs sampler by exploiting the knowledge of the conditional probability density functions of the variables entering the system. The two sampling techniques can be used in combination with the EM method to perform empirical Bayes estimation of the kernel hyperparameters and the noise variance. We have highlighted, through some numerical experiments, the advantages of employing our method when quantizers affect the accuracy of measurements.

As a final remark, we note that the cascaded composition of a linear LTI dynamic system followed by a static nonlinear function is known in system identification as a *Wiener system* [@giri2010block], [@hagenblad2008maximum]. However, in Wiener systems the nonlinear function is generally assumed unknown, so a direct extension of the method proposed in this paper to general Wiener systems does not appear immediate, and would require the use of more involved and MCMC techniques (see, e.g., [@lindsten2013bayesian]).

[^1]: All the results obtained in this paper can be also extended to the case $m=\infty$ but we prefer to consider a finite value for $m$ to simplify exposition. In addition, the assumption $m<N$, with $N$ the data set size, is typical of the system identification scenario.