Bayesian Function Learning Using MCMC Methods.
AbstractThe paper deals with the problem of reconstructing a continuous one-dimensional function from discrete noisy samples. The measurements may also be indirect in the sense that the samples may be the output of a linear operator applied to the function (linear inverse problem, deconvolution). In some cases, the linear operator could even contain unknown parameters that are estimated from a second experiment (joint identification-deconvolution problem). Bayesian estimation provides a unified treatment of this class of problems, but the practical calculation of posterior densities leads to analytically intractable integrals. In the paper it is shown that a rigourous Bayesian solution can be efficiently implemented by resorting to a MCMC (Markov chain Monte Carlo) simulation scheme. In particular, it is discussed how the structure of the problem can be exploited in order to improve computational and convergence performances. The effectiveness of the proposed scheme is demonstrated on two classical benchmark problems as well as on the analysis of IVGTT (IntraVenous Glucose Tolerance Test) data, a complex identification-deconvolution problem concerning the estimation of the insulin secretion rate following the administration of an intravenous glucose injection.
Introduction
The problem of reconstructing (learning) an unknown function from a set of experimental
data plays a fundamental role in engineering and science. In the present paper,
the attention is restricted to scalar functions of a scalar variable
the main ideas may apply to more general maps as well. In the most favourable case, it is
possible to directly sample the function. Rather frequently, however, only indirect measurements
are available which are obtained by sampling the output of a linear operator
applied to the function. For instance, in deconvolution problems what is sampled is the
convolution of the unknown function with a known kernel, see e.g. [1], [2], [3], [4], [5].
The approaches used to solve the function learning problem can be classified according
to three major strategies. The parametric methods assume that the unknown function
belongs to a set of functions which are parameterized by a finite-dimensional parameters
vector. For instance, the function can be modelled as the output of a multi-layer
perceptron which, for a given topology, is completely characterized by the values of its
weights [6]. Another possibility is to use a polynomial spline with fixed knots. In both
cases, function learning reduces to the problem of estimating the model parameters, a task
that can be performed by solving a (possibly nonlinear) least squares problem. If it were
true that the unknown function belongs to the given function space, statistical estimation
theory could be invoked in order to find minimum variance estimators and compute confidence
intervals [7]. Moreover, it would also be possible to compare parametric models
of increasing complexity using statistical tests (F-test) or complexity criteria (Akaike's
criterion).
The second strategy, namely regularization [8], [9], [10], [11], avoids introducing heavy
assumptions on the nature of F(\Delta) but rather classifies the potential solutions according
to their "regularity" (typically by using an index of smoothness such as the integral
of the squared k-th derivative of the function). The relative importance of the sum of
squared residuals against the regularity index is controlled by the so-called regularization
parameter. The key problem is finding an optimal criterion for the selection of the
regularization parameter, although empirical criteria such as ordinary cross validation
and generalized cross validation [12], [13] perform satisfactorily in many practical cases.
Moreover, regularization does not provide confidence intervals so that it is not possible
to assess the reliability of the reconstructed function.
The present paper deals with the third strategy which is based on Bayesian estimation.
The unknown function is seen as an element of a probability space whose probability
distribution reflects the prior knowledge. For instance, the prior knowledge that F(\Delta)
is smooth is translated in a probability distribution that assigns higher probabilities to
functions whose derivatives have "small" absolute values. A practical way to do that is
to describe F(\Delta) as a Gaussian stochastic process whose k-th derivative is a white noise
process with intensity - 2 . Provided that both - 2 and the variance oe 2 of the (Gaussian)
measurement noise are known, the Bayes formula can be used to work out the posterior
distribution of F(\Delta) given the data [2], [14], [15], [16]. The posterior provides a complete
description of our state of knowledge. In particular, the mean of the posterior can be used
as a point estimate (Bayes estimate) whereas the variance helps assessing the accuracy.
It is notable that, if the regularization parameter is taken equal to the ratio oe 2 =- 2 the
regularized estimate coincides with the Bayes one.
The main advantage of the Bayesian approach is the possibility to address the selection
of the regularization parameter in a rigourous probabilistic framework. In fact, when - 2
is not known, it can be modelled as a random variable and two different approaches are
possible. The simpler one is based on the following observation: if the prior distribution
of - 2 is very flat, the maximum of its posterior given the data is close to the maximum
likelihood estimate - 2
ML . Then, if the posterior of - 2 is very peaked around its max-
imum, it is reasonable to estimate F(\Delta) as if - 2
ML were the true value of - 2 [14], [17], [4].
MAGNI ET AL.: BAYESIAN FUNCTION LEARNING USING MCMC METHODS 3
However, if the posterior of - 2 is not "very peaked" (which is likely to happen especially
for medium and small data sets), neglecting the uncertainty on - 2 would lead to under-estimated
confidence margins for F(\Delta). The truly Bayesian approach, conversely, calls for
the computation of the posterior of F(\Delta) taking into account also the random nature of
. Since the involved integrals are analytically intractable, one has to resort to Monte
Carlo methods.
A first purpose of the present paper is to show how a truly Bayesian solution of the
function learning problem can be efficiently worked out. We discuss the various stages of
the procedure starting from the discretization of F(\Delta) to arrive at the practical implementation
of Markov chain Monte Carlo (MCMC) methods [18]. Our approach is similar to the
one proposed in [15] where, however, the case of indirect measurements (deconvolution
problem) is not treated.
A further issue addressed in the paper is the joint identification-deconvolution problem,
which arises when the convolution kernel is not a priori known but is to be identified by a
separate experiment. The standard (suboptimal) approach is to identify the convolution
kernel and then use it as if it were perfectly known in order to learn the unknown function
F(\Delta). As shown in the paper, the use of MCMC methods allows to learn both functions
jointly.
The paper is organized as follows. Section II contains the statement of the problem. In
Section III, after a concise review of MCMC methods, a numerical procedure for solving
the Bayesian function learning problem is worked out. In section Section IV the proposed
method is illustrated by means of simulated as well as real-world data coming from the
analysis of metabolic systems. Some conclusions (Section V) end the paper.
II. Problem statement
In this paper we consider the problem of reconstructing a function ~
discrete and noisy samples y k such that
~
denotes the measurement error and ~
L k is a linear functional. In increasing order
of generality, we have:
~
~
~
~
denote the sampling instants and t is the initial time.
The first definition of ~
corresponds to the function approximation problem based on
samples of the function itself. Using the second and third definitions, the problem (1)
becomes a deconvolution problem, or an integral equation of the first kind (also called
Fredholm equation), respectively. As for the noise v k , letting v := [v 1
assumed that is a known matrix and oe 2 is a (possibly
unknown) scalar. In the following, y := [y 1
~
f ~
~
non parametric estimator based on Tychonov regularization is given by
~
~
f
P is a suitable operator and k \Delta k is a norm in a
suitable function space. A typical choice for ~
is:
~
Then, if in addition ~
L k is as in (2) and "L 2 norm" is used, ~
turns out to be a smoothing
spline [13]. The basic idea behind (5) is to find a balance between data fit and smooth-
MAGNI ET AL.: BAYESIAN FUNCTION LEARNING USING MCMC METHODS 5
ness of the solution, the relative weight being controlled by the so-called regularization
parameter fl. Various criteria have been proposed to tune fl. Among them, we may mention
ordinary cross-validation [1], generalized cross-validation [12], and the L-curve [19].
It is worth noting that without the aid of the smoothness penalty flk ~
would be
impossible to learn the function ~
f (belonging to an infinite dimensional space) from the
finite data set fy k g.
Rather interestingly, (5) can also be interpreted as a Bayesian estimator. To this pur-
pose, assume that ~
f is a stochastic process such that ~
f is a white noise with intensity
~
w(t)]=0, 8t and E[ ~
w(t) ~
assuming that both v k and ~
f
are normally distributed and letting
, the regularized estimate ~
coincides with the
conditional expectation of ~
f given the observations y k , i.e. E[ ~
f fl . According to
the Bayesian paradigm, the probabilistic assumptions on the unknown function ~
f should
reflect our prior knowledge. For instance, if the operator ~
P is as in (6), then ~
f is the
double integral of a white noise: as such it is a relatively smooth signal (the smaller ~
the smoother ~
f will be).
The above considerations suggest that ~
f fl is the "optimal" estimator provided that
. Unfortunately, ~ - 2 (and sometimes also oe 2 ) is very unlikely to be known a priori.
In the literature [20], [14], [4] it has been proposed to regard ~
as an unknown parameter
and compute its maximum likelihood estimate:
~
In this way ~
f fl with
becomes a sub-optimal estimator.
In this paper, conversely, we pursue a rigourous Bayesian approach. The unknown
parameter ~ - 2 is modelled as a hyper-parameter having a suitable prior distribution, which
is taken into account in the computation of the posterior density
Z
~
f)p( ~
6 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE
Since the analytic evaluation of the posterior is made intractable by the presence of
the hyper-parameter, the calculations are carried out by means of Monte Carlo sampling
algorithms. In this context it is convenient to discretize the original problem. Without
loss of generality, it is assumed that t 0 =0. Then, given a sufficiently small discretization
interval T , consider f
. For simplicity, it is
assumed that the sampling instants t k are multiples of T , i.e. t
the operator ~
P is suitably discretized. For instance, with approximated by:
Moreover, ~
f is approximated by Lf , where L is a suitable n \Theta N matrix:
if ~
L k is as in (2) then L kj
if ~
L k is as in (3) then L
R T~
if ~
L k is as in (4) then L
R T~
These approximations are obtained under the assumption that the unknown function ~
f
is constant in between the sampling instants. After the discretization, (5) is approximated
by:
f
denotes the usual norm in ! n ), whose closed form solution is:
MAGNI ET AL.: BAYESIAN FUNCTION LEARNING USING MCMC METHODS 7
In the Bayesian setting it is assumed that the entries of the vector are taken
from the realization of a (discrete-time) white noise with variance - 2 I, that is E[w]= 0,
A first purpose of the present paper is to investigate the applicability of Markov chain
Monte Carlo [18] sampling methods to estimate p(f j y) assuming that - 2 is a parameter
with suitable prior distribution. Another issue concerns the Bayesian solution of the
joint identification-deconvolution problem. In fact, in many deconvolution problems the
impulse response ~ h(-) that enters the convolution integral (3), depends on one or more
parameters -, that are estimated by a separate experiment: ~ putting
together the two experiments, we obtain:
where M(-) is a known function of -, z := [z is the data vector used to identify
- and ffl := is the corresponding measurement error.
The standard approach is to estimate -
using (10) and then estimate f from using -
as if it were the true value of -. On the other hand, a truly Bayesian approach describes
- as a random variable. Then, p(f j y) should be evaluated by considering (9) and (10)
simultaneously. As a particular case, it is possible to consider (9) alone with - modelled
as a random variable to allow for its uncertainty. Again, the standard "suboptimal"
approach is to compute f fl using the nominal value of - [4] and then assess the sensitivity
of the estimate with respect to parameters uncertainty.
III. Markov chain Monte Carlo Methods in Bayesian estimation
problems
Probabilistic inference involves the integration over possibly high-dimensional probability
distributions. Since this operation is often analytically intractable, it is common to
8 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE
resort to Monte Carlo techniques, that requires sampling from the probability distribution
to be integrated. Unfortunately, sometimes it is impossible to extract samples directly
from that distribution. Markov chain Monte Carlo (MCMC) methods [18] provide a
unified framework to solve this problem.
MCMC methods are based on two steps: a Markov chain and a Monte Carlo integra-
tion. By sampling from suitable probability distributions, it is generated a Markov chain
that converges (in distribution) to the target distribution, i.e. the distribution to be in-
tegrated. Then, the expectation value is calculated through Monte Carlo integration over
the obtained samples.
The MCMC methods differ from each other in the way the Markov chain is generated.
However, all the different strategies proposed in the literature, are special cases of the
Metropolis-Hastings [21], [22] framework. Also the well-known Gibbs sampler [23] fits in
the Metropolis-Hastings scheme.
In the following sub-sections, we will describe the application of the Metropolis-Hastings
algorithm to Bayesian function learning as well as discuss about its possible variants.
A. The method
In order to describe the Metropolis-Hastings algorithm we will use the following notation

the vector of the model parameters
the vector of the data (i.e. the observations)
ffl \Theta i - the i th samples drawn
the target distribution (proportional to the posterior distribution)
where p(') is the prior distribution of the model parameters and p y
the likelihood.
MAGNI ET AL.: BAYESIAN FUNCTION LEARNING USING MCMC METHODS 9
The Markov chain derived by the Metropolis-Hastings method is obtained through the
following steps:
1. at each time t, a candidate sample \Theta is drawn from a proposal distribution
2. the candidate point \Theta is accepted with probability:
3. if the candidate point \Theta is accepted, the next sample of the Markov chain is \Theta
else the chain does not move and \Theta
It is important to remark that the stationary distribution of the chain (i.e. the distribution
to which the chain converges) is independent of the proposal distribution [18], and
coincides with the target distribution p ';y (').
Although any proposal distribution, on long-run, will deliver samples from the target
distribution, the rate of convergence to the stationary distribution of the generated
Markov chain crucially depends on the relationships between the proposal and the target
distributions; moreover, the number of samples necessary to perform the Monte Carlo
steps depends on the speed with which the algorithm "mixes" (i.e. spans the support of
the target distribution).
When the vector of the model parameters is large, it is often convenient to divide '
into K components and update the samples \Theta of these components one by one [24]. This
scheme is called Single-Component Metropolis-Hastings
Let \Theta (i)
t the i th component of \Theta at time t and \Theta (\Gammai)
\Theta (K)
t g, the Metropolis-Hastings scheme turns to:
1. at each time t, the next sample \Theta (i)
t+1 is derived by sampling a candidate point \Theta (i)
from a proposal distribution q (i)
2. the candidate point \Theta (i)
is accepted with probability:
ff(\Theta (i)
3. if the candidate point \Theta (i)
is accepted, the next sample of the Markov chain is \Theta (i)
\Theta (i)
, else the chain does not move and \Theta (i)
t .
The Gibbs sampler (GS) is just a special case of the single-component Metropolis-
Hastings. The GS scheme exploits the full conditional (the product of the prior distribution
and the likelihood) as the proposal distribution. In this case, it is easy to verify
that the candidate point is always accepted, so that the Markov chain moves at every
step. When the full conditionals are standard distributions (easy to sample from), the
GS represents a suitable choice. On the contrary, when it is not possible to draw samples
directly from the full-conditional distributions, it is convenient to resort to mixed schemes
Metropolis-Hastings). In this setting, a portion of the model parameters
is estimated using the Gibbs Sampler, while the other ones are treated using "ad-hoc"
proposal distributions.
These algorithms have been extensively used in the field of probabilistic graphical modelling
[25]. Using this kind of models a suitable partition (blocking) of the vector of model
parameters is naturally obtained and it is also easy to derive the full conditional distribu-
tion. The convergence rate and the strategies for choosing the proposal distribution are
described in [26], [27].
B. MCMC in Function Reconstruction
In this sub-section we will describe how the problems defined in Section II can be tackled
using MCMC methods. In order to explain the probabilistic models used in the different
sampling schemes, we will resort to a Bayesian Network (BN) representation. BNs are
Directed Acyclic Graphs (DAGs) in which nodes represent variables, while arcs express
direct dependencies between variables. These models are quantified by specifying the
conditional probability distribution of each node given its parents. They will help us in
MAGNI ET AL.: BAYESIAN FUNCTION LEARNING USING MCMC METHODS 11
expressing the conditional independence assumptions underlying the different function
reconstruction problems. For further details see [28], [29].
B.1 Function approximation based on direct sampling
Consider the function approximation problem based on samples of the function itself
(smoothing problem). Its formal definition is given in equations (1) and (2). Our goal is
to provide a Bayesian estimate of the vector f (i.e. the discretized unknown function).
As shown in Section II, the discretized form of the problem may be written as:
Referring to Section II, f :=P \Gamma1 w with E[w]=0, E[ww T
To apply the MCMC strategy described in Section III-A, we must assign a suitable
probabilistic model to the parameter set g. By exploiting a set of standard
choices in Bayesian estimation problems, it is assumed that w is normally distributed
given - 2 and that the precision parameter 1
has a Gamma distribution. More formally,
and
. Moreover, we suppose that the noise v has a Normal distribution,
covariance matrix oe 2 \Psi. This implies that the data model can
be written as: p(y
In this setting the target distribution becomes:
The model is described by the simple BN of Fig. 1
It is easy to see that, in order to apply MCMC integration, it is useful to adopt the
partition =wg. In this context, it is convenient to adopt the Gibbs
Sampler, since the full conditional distributions assume the following standard form:
y
I
Y
Fig. 1. Function approximation: probabilistic model of the smoothing problem, deconvolution problem
and Fredholm equation problem.
where
From the point estimate -
derived by the MCMC algorithm, it is trivial to
reconstruct the unknown function f as -
Moreover, having samples from the
joint posterior distribution, it is possible to derive any statistics of interest, including
confidence intervals (or more appropriately: Bayes intervals).
B.2 Function approximation in inverse problems
In deconvolution problems the unknown function has to be reconstructed on the basis
of indirect measurements: a convolution integral expresses the relationships between the
samples and the unknown function, see (1) and (3). The structure of the problem is
analogous when considering integral equations of the first kind (Fredholm equations), see
(1) and (4)
Again, our goal is to provide a Bayesian estimate of the (discretized) unknown f func-
MAGNI ET AL.: BAYESIAN FUNCTION LEARNING USING MCMC METHODS 13
tion. The discretized form of the problem still has the form (11). Since the functions ~ h(-)
and ~ h(t; -) that enter in the definitions (3) and (4) are assumed to be known, the matrix
L is completely specified. Also in this case, f := P \Gamma1 w with
so that the parameter set to be estimated is the same as in the smoothing problem:
g. Thus, the probabilistic model for the inverse problems (1),(3) and (1),(4) is
again described by the BN of Fig. 1. In fact, in our setting, the smoothing problem, the
deconvolution problem and the Fredholm equation problem differ only in the computation
of the matrix L.
B.3 Deconvolution problems with uncertain impulse response
An interesting extension of the problem described in the previous section is to relax the
assumption of complete knowledge of the impulse response ~ h(-) appearing in equation
(3). As anticipated in Section II, we suppose that ~ is a function of a set
of unknown parameters -, which have to be estimated from experimental data. In this
way, the problem becomes the simultaneous estimation of the unknown function f and
the parameter set -, given the model described by (9) and (10).
Once again, f := P \Gamma1 w with I. In this case, however, the
probabilistic model to be specified becomes more complex. The parameter set to be
estimated is now -g. The corresponding probabilistic model is assumed as:
- 0 and \Sigma - . Moreover, we suppose that the noise signals v and ffl are independent and
normally distributed with known covariance matrices, namely v - N(0; oe 2 \Psi) and ffl -
Then, the data model can be written as: p(y j w;
In this setting the target distribution becomes:
The resulting model is described by the BN of Fig. 2.
14 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE
y
I
x
x
x
z ~ N(M(
1/ l
YFig. 2. Function approximation in deconvolution problems in presence of uncertainty on the impulse
response: probabilistic model.
In order to devise an MCMC scheme for estimating ', it is convenient to partition ' as
=-g. The full conditional distributions are as follows:
where
and
The full conditional (16) is clearly a non-standard distribution, so that it cannot be
sampled directly. In order to use the GS strategy it is necessary to apply sampling
algorithms for general distributions, like rejection sampling or adaptive rejection sampling
MAGNI ET AL.: BAYESIAN FUNCTION LEARNING USING MCMC METHODS 15
[30]. Unfortunately, such algorithms may impair the overall efficiency of the stochastic
simulation machinery.
A valuable alternative is to resort to a mixed MCMC scheme in which different proposal
distributions are used to extract samples from the different partitions of '.
ffl The proposal for ' (1) ; ' (2) are the full conditional distributions, as in the GS. In this
way the candidate point for these partitions is always accepted.
ffl The proposal distribution for ' (3) can be chosen as the prior distribution
In this case, the proposal distribution is independent from the past sample
drawn by the Markov chain, so that This scheme is also known as
independence sampler [31]. The acceptance probability for the candidate sample \Theta (3)
simplifies as:
ff(\Theta (3)
so that the acceptance rate depends only on the ratio of the likelihood in the candidate
point to the likelihood for the current one. On the basis of the specific problem, other
proposal distributions can be chosen in order to obtain the best performance in the
computational speed. For example, when possible, a good choice can be the use of a
standard distribution that is a good approximation of the full conditional distribution.
This is a way to preserve the advantages of the single-component Metropolis-Hastings,
avoiding the additional computational burden that would be entailed by the GS in
presence of non-standard full conditional distributions.
IV. Bayesian function learning at work
In this section we will show how the above presented methodology is able to cope with
three different benchmark problems taken from the literature.
A. Function approximation based on direct sampling
To test the performance of the MCMC function approximator in the smoothing problem,
we consider an example proposed by Wahba [13]. The function to be approximated is:
~
and the noisy samples y k are:
Since we are interested in reconstructing the function only in correspondence of the
measurement times, in equation (11) L is the identity matrix. Following Section III-B.1
we take the following prior
where I is the identity matrix. The choice of the prior distribution parameters for the- 2
reflects the absence of reliable prior information on the regularization parameter. In
fact with this choice - with probability of 0.9.
As in [13], we assume that the second derivative of the function is regular; taking into
account the discretization, the operator P is chosen as:
MAGNI ET AL.: BAYESIAN FUNCTION LEARNING USING MCMC METHODS 17
-0.50.5t (arbitrary units)
f
(arbitrary
units)
Fig. 3. Function approximation using the MCMC-smoother.: noisy data (stars), true function (dash-dot
line) and reconstructed function (continuous line).
The starting point of the Markov chain was extracted from the prior distribution
of the parameters fw; - 2 g. After a 1000 samples run of the MCMC scheme, convergence
of the estimates was verified by the method described in [26]. In particular, after
choosing the quantiles 0:975g to be estimated with precision
respectively with probability burn-in (N) of
samples and a number of required samples (M) of 870 were calculated. The results
are shown in Fig. 3.
Although the samples are rather noisy, the smoothed signal is close to the true function.
The RMSE (root mean square error) obtained is 0.073. Moreover, the performances of the
MCMC smoother are similar to the ones obtained in [13], where a cubic smoothing spline
was used and the regularization parameter was tuned through ordinary cross-validation
(OCV). Hence, the MCMC smoother is as good as OCV-tuned smoothing splines in
avoiding under- and over-smoothing problems. The main advantage of MCMC smoother is
that it provides also the a-posteriori sampling distribution of the regularization parameter
and the confidence intervals for the reconstructed function in a rigourous Bayesian setting.
B. Function approximation in deconvolution problems
In order to test the performances of the MCMC deconvolution scheme, we consider a
well-known benchmark problem [32], [2], [4]. The input signal given by:
~
has been convoluted with the impulse response:
Then, by adding measurement errors v k simulated as a zero-mean white Gaussian noise
sequence with variance equal to 9, 52 noisy samples are generated at time t
Fig. 4(a) shows the true function ~
while Fig. 4(b) depicts the convoluted
function together with the 52 noisy samples.
Our goal is to reconstruct the unknown function with a sufficiently fine resolution. This
means that we are interested in the function estimates not only in correspondence of the
measurement times, but also in other "out of samples" time points. In particular, we
consider a 208 points on an evenly-spaced time grid in the interval [0 1035]. Then, the
entries of L are
R T~
with
We take the following prior distributions, similar to the ones used in the previous section:- 2 - \Gamma(0:25; 5e \Gamma 7)
MAGNI ET AL.: BAYESIAN FUNCTION LEARNING USING MCMC METHODS 19
(a)
(arbitrary units)
f
(arbitrary
units)
(arbitrary units)
y
(arbitrary
units)
Fig. 4. Simulated deconvolution problem. Panel (a): the function to be reconstructed. Panel (b): the
convoluted function (dash-dot line) and the noisy samples (stars).
Moreover the operator P is selected as:
This choice corresponds to a penalty on the squared norm of the first derivative, (which
is approximated by the first difference of the discretized signal).
The bottleneck of the MCMC scheme is the computation of B \Gamma1 in equation (15)(this
matrix inverse has to be performed at each step of the MCMC scheme). However, by a
proper change of coordinates, it is possible reduce the size of the matrix has to be inverted
from N \Theta N to n \Theta n (i.e. from 208 \Theta 208 to 52 \Theta 52 in our problem). This goal can be
achieved through the following steps:
1. Let hence an n \Theta N matrix), and compute the
SVD (singular value n) and V (N \Theta N)
are orthogonal matrices (UU and D is an n \Theta N diagonal matrix
(D
2. In view of (11)
It is easy to verify that v is
distributed as N(0; I n ), and w as N(0; - 2 I N ).
3. We apply the same MCMC scheme described in the previous section to the reformulated
problem (17). In the new coordinates B is a block diagonal matrix, in which
the first one is an N \Theta N block and the other ones are 1 \Theta 1 blocks.
4. The final estimate is obtained by re-transforming the variables in the original co-
ordinates; in particular we need to compute:
The starting point of the Markov chain was extracted from the prior distribution of
the parameters fw; - 2 g. After 5000 steps of the MCMC scheme, the convergence of
the estimates was verified by using the method described in [26]. In particular after
choosing the quantiles 0:975g to be estimated with precision
respectively with probability burn-in (N) of
132 samples and a number of required samples (M) of 3756 were calculated. The results
are shown in Fig. 5.
The performance of our approach is comparable with the one proposed in [4], where the
regularization parameter is estimated according to a maximum likelihood criterion (see
Fig. 6, Fig. 7).
MAGNI ET AL.: BAYESIAN FUNCTION LEARNING USING MCMC METHODS 21
(arbitrary units)
f
(arbitrary
units)
(arbitrary units)
y
(arbitrary
units)
Fig. 5. Simulated deconvolution problem. Panel (a): the true function (dash-dot line) and the re-constructed
one (continuous line), with its 95% confidence interval (dashed line), obtained with the
MCMC scheme. Panel (b): the true noiseless output (dash-dot line), the estimated output (continu-
ous line) and the noisy samples (*).
The RMSE obtained with our approach is 0.065 while the one obtained by the method
of [4] is 0.059. Again, the advantage of the MCMC scheme is its ability to provide
the a-posteriori sampling distribution of the regularization parameter and the confidence
intervals for the reconstructed function in a rigourous Bayesian setting.
C. Deconvolution with uncertain impulse response
In this subsection, the MCMC scheme is applied to a real-world problem taken from
[33]; in particular, we demonstrate that deconvolution and impulse response identification
can be addressed jointly, as described in Section III-B.3.
The goal is to quantify the Insulin Secretion Rate (ISR) in humans after a glucose stim-
ulus; the experimental setting is related to the so-called IntraVenous Glucose Tolerance
22 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE
(arbitrary units)
f
(arbitrary
units)
Fig. 6. Simulated deconvolution problem. Comparison of the MCMC scheme with the maximum
likelihood one (see text). True function (dashed line), MCMC estimate (continuous line), maximum
likelihood estimate (thick line).
Test (IVGTT), where an impulse dose of glucose is administered in order to assess the
subject capability of bringing Blood Glucose Levels within normal ranges through the
endogenous release of Insulin, the main glucoregolatory hormone.
The ISR cannot be directly measured since insulin is secreted by the pancreas into the
portal vein which is not accessible in vivo. It is possible measure only the effect of the
secretion in the circulation (the plasma concentration of insulin). But, because of the large
liver extraction, the plasma insulin concentration reflects only the post-hepatic delivery
rate into the circulation. This problem can be circumvented by measuring C-peptide (CP)
concentration in plasma. The CP is co-secreted with insulin on an equimolar basis, but is
not extracted by liver, so that it directly reflects the pancreatic ISR. Thus, the problem
turns into the estimation of the ISR on the basis of the (noisy) measurements of CP in
plasma. Since the CP kinetics can be described by a linear model, we obtain the following
MAGNI ET AL.: BAYESIAN FUNCTION LEARNING USING MCMC METHODS 23
samples
lambda^2
samples
of
lambda^2
(b)
Fig. 7. Simulated deconvolution problem. Panel (a): Markov chain of the - 2 parameter. Panel (b):
frequency histogram proportional to the posterior distribution of the parameter - 2 estimated. The
derived with the maximum likelihood approach (see text) was 0.00665.
~
where y k are the CP plasma measurements (pmol/ml), ~
f(\Delta) is the ISR (pmol/min), ~ h(\Delta)
is the CP impulse response (ml \Gamma1 ) and v k is the measurement error (pmol/ml). The CP
impulse response is [33]:
~ h(-) =X
where the parameters A i s and ff i s have to be estimated through an ad-hoc experiment; in
particular an intravenous bolus of biosynthetic CP is delivered to the patient and a number
of plasma measurements of CP are collected (a somatostatin infusion is administered in
order to avoid the endogenous pancreatic secretion). The impulse response parameters
depend on the single patient but are considered constant over time for a specific patient.
We can apply the MCMC scheme of Section III-B.3 in order to jointly perform the CP
impulse response identification and the ISR reconstruction.
Following Section II the problem can be written as:
is the matrix obtained from the discretization of
the deconvolution integral, z are the noisy measurements of the CP plasma concentration
during the impulse response identification experiment, v and ffl are the measurements
errors and are the sampling instants of the
identification experiment. L(-) has elements:
In this case, some prior knowledge on the signal is available [33]: the ISR is known
to exhibit a spike just after the external glucose stimulus and a more regular profile
thereafter. This knowledge can be modelled by considering two different regularization
parameters for the two phases.
For computational reasons, a nonuniform discretization for f has been adopted. Ac-
cordingly, the regularization operator P has been chosen so as to satisfy:
-T
Then:
MAGNI ET AL.: BAYESIAN FUNCTION LEARNING USING MCMC METHODS 25
y
~ N(0, L)
x
x
~ N( ,S x
x
z ~ N(M(
1/ l
e
s
Fig. 8. Probabilistic model for the reconstruction of Insulin Secretion Rate
Moreover, the variances oe 2
v and oe 2
ffl of the measurements errors v and ffl are only imprecisely
known and must be estimated as well. The complete model is shown in Fig. 8.
The derived sampling distributions for the MCMC scheme described in Section III-B.3
are:
26 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE
where
In the above formulas N 1 is the number of point used in the regularization of the first
region (the spike) while N 2 is the number of point used for the second region
the vector w is such that
the number of measurements in the identification experiment.   is a diagonal matrix with
the first N 1 elements equal to - 2
1 and with the other N 2 elements equal to - 2
2 . Finally,
\Psi and \Psi ffl are diagonal matrices with elements y 2 and z 2 respectively; in this way, the
parameters oe v and oe ffl represent the CV (coefficient of variation) of the measurements
errors in the two experiments.
To completely specify the MCMC scheme the following prior distributions must be
3 0:25 0:125 0:05] T and \Sigma - is a diagonal matrix with elements
such that \Sigma - (i;
MAGNI ET AL.: BAYESIAN FUNCTION LEARNING USING MCMC METHODS 27
The value of - 0 was derived on the basis of the prior knowledge on the dynamics of the
impulse response, while the value of the hyper-parameters for oe v and oe ffl are assessed by
knowing that the measurements error has a CV that ranges from 4% to 6%. The priors
for the -s reflect the knowledge on the signal shape in the two response phases.
We perform our test on the data set described in [34], [33]. The data set used for the
impulse response identification are collected after a CP bolus of 49650 pmol, while the
data set used for the ISR reconstruction are taken after an intravenous glucose bolus of
0:5 g=kg. The basal value of the ISR is estimated on the basis of the five CP measurements
taken before the glucose bolus. Our goal is to reconstruct the ISR in correspondence of
the sampling instants in which CP measurements are taken (in this case N=n).
As in [33], we take N corresponding to a first phase ranging
16. The data of the two experiments are shown in Fig. 9.
After a 2500 samples run of the MCMC scheme (convergence was verified by using the
method described in [26] as the same assumption on q, r, s previous reported), the results
shown in Fig. 10 were obtained.
Fig. 10(a) shows the estimate of the CP plasma levels in the IVGTT experiment: the
estimated curve is slightly smoother than the measurements. Fig. 10(b) depicts the ISR
curve as estimated after deconvolution by the MCMC scheme: the reconstructed ISR
reproduces the expected physiological shape, characterized by two regions with different
regularities. The results obtained are comparable with those obtained in [33], where
deconvolution and impulse response identification are treated separately. Fig. 10(c) shows
the estimated CP impulse response identification. It is easy to notice the good quality
of the fit. In Fig. 11 the frequency histograms of the samples generated by the MCMC
estimator for the six impulse response parameters are reported.
The proposed MCMC scheme is able to jointly perform the identification of the impulse
response and the deconvolution of the ISR. In the classical approach [33], the two experi-
28 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE
(a)
time (min)
in
plasma
(pmol/ml)
time (min)
in
plasma
(pmol/ml)
Fig. 9. The data set for the ISR (Insulin Secretion Rate) reconstruction. Panel (a): impulse response
identification experiment, consisting of 32 noisy samples of the CP (C-Peptide concentration) in
plasma, collected after a CP intravenous bolus at time 0 min. Panel (b): Intravenous Glucose
Tolerance Test, consisting of 27 noisy samples of CP in plasma, collected around an intravenous
glucose bolus at time 0 min.
ments are treated in a separate fashion: in the first step, the impulse response is identified,
using the measurements of the "identification set", and only in the second step, the data
of the "deconvolution set" is used to reconstruct the unknown function. The uncertainty
on the impulse response is possibly taken into account only after the deconvolution step.
On the contrary, our scheme combines together the information coming from the two ex-
periments, and uses it in order to provide "optimal" point estimates as well as posterior
moments and confidence intervals.
V. Conclusions
MCMC methods constitute a set of emerging techniques, that have been successfully
applied in a variety of contexts, from statistical physics to image analysis [23], and from
MAGNI ET AL.: BAYESIAN FUNCTION LEARNING USING MCMC METHODS 29
time (min)
in
plasma
(pmol/ml)
(b)
time (min)
ISR
(c)
time (min)
in
plasma
(pmol/ml)
Fig. 10. Solution of the joint deconvolution and impulse response identification problem. Panel (a):
estimated level of CP in plasma during the IVGTT (continuous line) and the collected samples
(stars). Panel (b): the ISR reconstructed (with 95% confidence interval) by the MCMC-estimator.
Panel (c): the estimated impulse response ~ h(\Delta) (continuous line) and the samples (stars) collected
during the identification experiment.
medical monitoring [25] to genetics [35] and archaeology[36]. The powerfulness of MCMC
methodologies lies in their inherent generality that enables the analyst to deal with the
full complexity of real world problems.
In this paper, we have exploited such generality to propose a unified Bayesian framework
for the reconstruction of functions from direct or indirect measurements. In particular,
by using the same conceptual scheme we easily coped with problems that had been previously
solved with ad-hoc methods. The obtained results are, in all cases, at least as
samples
of
(a)
samples
of
alpha2
(b)
samples
of
alpha3
(c)
value (1/ml)
samples
of
(d)
samples
of
samples
of
A3
Fig. 11. Frequency histograms proportional to the posterior distribution of estimated impulse response
parameters.
good as the previously proposed solutions. In addition, since our approach is able to
soundly estimate the posterior probability distribution of the reconstructed function, the
information provided at the end of the estimation procedure is richer than in all other
methods: first and second moments, confidence intervals and posterior distributions are
obtained as a by-product. Finally, our framework has been exploited to implement a new
strategy for the joint estimation of a deconvoluted signal and its impulse response. The
previous approaches were based on a two-step procedure, which is not able to optimally
combine all the information available in the data.
The main limits of MCMC approach is the time required to converge to the posterior
distribution and the difficulty to choose the best sampling scheme. These limitations force
to use MCMC methods only for off-line reconstructon.
In summary, MCMC methods have been shown to play a crucial role in the off-line
MAGNI ET AL.: BAYESIAN FUNCTION LEARNING USING MCMC METHODS 31
function learning problem, since they provide a flexible and relatively simple strategy,
able to provide optimal results in a Bayesian sense.

Acknowledgements

The authors would like to thank Antonietta Mira for her methodological support in
designing the MCMC scheme and Claudio Cobelli and Giovanni Sparacino for having
provided the experimental data for the Insulin Secretion Rate reconstruction problem.
They thank also the anonymous reviewers for their useful suggestions.



--R

"Practical approximate solutions to linear operator equations when the data are noisy,"
"The deconvolution problem: Fast algorithms including the preconditioned conjugate-gradient to compute a map estimator,"
"Linear inverse problems and ill-posed problems,"
"Nonparametric input estimation in physiological systems: Problems, methods, case studies,"
"Blind deconvolution via sequential imputations,"
Bayesian Learning for Neural Networks
Parameter Estimation in Engineering and Science
"A technique for the numerical solution of certain integral equations of the first kind,"
"On the numerical solution of Fredholm integral equations of the first kind by the inversion of the linear system produced by quadrature,"
Solutions of Ill-Posed Problems
"Networks for approximation and learning,"
"Smoothing noisy data with spline functions: Estimating the correct degree of smoothing by the method of generalized cross-validation,"
Spline Models for Observational Data
"Bayesian interpolation,"
"Gaussian processes for regression,"
"Automatic bayesian curve fitting,"
"Nonparametric spline regression with prior information,"
Markov Chain Monte Carlo in Practice
"Numerical tools for analysis and solution of Fredholm integral equation of the first kind,"
"A time series approach to numerical differenti- ation,"
"Equations of state calculations by fast computing machine,"
"Monte Carlo sampling methods using Markov Chain and their applications,"
"Stochastic relaxation, Gibbs distributions, and the bayesian restoration of images,"
"Likelihood analysis of non-gaussian measurment time series,"
"A unified approach for modeling longitudinal and failure time data, with application in medical monitoring,"
"Implementing MCMC,"
"Inference and monitoring convergence,"
Probabilistic Reasoning in Intelligent Systems
"Dynamic probabilistic networks for modelling and identifying dynamic systems: a MCMC approach,"
"Adaptive rejection sampling for Gibbs sampling,"
"Markov Chains for posterior distributions (with discussion),"
"The inverse problem of radiography,"
"A stochastic deconvolution method to reconstruct insulin secretion rate after a glucose stimulus,"
"Peripheral insulin parallels changes in insulin secretion more closely than C-peptide after bolus intravenous glucose administration,"
"Censored survival models for genetic epidemi- ology: a Gibbs sampling approach,"
"An archaeological example: radiocarbon dating,"
--TR

--CTR
Gianluigi Pillonetto , Claudio Cobelli, Brief paper: Identifiability of the stochastic semi-blind deconvolution problem for a class of time-invariant linear systems, Automatica (Journal of IFAC), v.43 n.4, p.647-654, April, 2007
Xudong Jiang , Wee Ser, Online Fingerprint Template Improvement, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.24 n.8, p.1121-1126, August 2002
Ranjan K. Dash , Erkki Somersalo , Marco E. Cabrera , Daniela Calvetti, An efficient deconvolution algorithm for estimating oxygen consumption during muscle activities, Computer Methods and Programs in Biomedicine, v.85 n.3, p.247-256, March, 2007
system identification;inverse problems;Markov chain Monte Carlo methods;dynamic systems;smoothing;bayesian estimation