STATISTICS IN MEDICINE
Statist. Med. (in press)
Published online in Wiley InterScience(www.interscience.wiley.com) DOI: 10.1002/sim.2941

A probabilistic algorithm for robust interference suppression
in bioelectromagnetic sensor data
Srikantan S. Nagarajan1, ∗, † , Hagai T. Attias2 , Kenneth E. Hild II1 andKensuke Sekihara3
1 Department of Radiology, University of California at San Francisco, 513 Parnassus Avenue 5362,

San Francisco, CA 94122, U.S.A.
2 Golden Metallic Inc., San Francisco, CA 94147, U.S.A.
3 Department of Electronic Systems and Engineering, Tokyo Metropolitan Institute of Technology,

191-0065, JapanSUMMARY
Magnetoencephalography (MEG) and electroencephalography (EEG) sensor measurements are often contaminated by several interferences such as background activity from outside the regions of interest, by
biological and non-biological artifacts, and by sensor noise. Here, we introduce a probabilistic graphicalmodel and inference algorithm based on variational-Bayes expectation-maximization for estimation ofactivity of interest through interference suppression. The algorithm exploits the fact that electromagnetic
recording data can often be partitioned into baseline periods, when only interferences are present, andactive time periods, when activity of interest is present in addition to interferences. This algorithm isfound to be robust and efficient and significantly superior to many other existing approaches on real andsimulated data. Copyright q 2007 John Wiley & Sons, Ltd.
KEY WORDS:

magnetoencephalography; electroencephalography; graphical models

1. INTRODUCTION
Bioelectromagnetic data are obtained by measuring electric and magnetic fields, which arise inbiological tissues using a sensor array. This paper is focused on electromagnetic fields arisingfrom the brain, but the techniques presented here apply to other biological systems, such as theheart. For brain tissues, electroencephalography (EEG) data are obtained by measuring electricfields using an electrode array placed on the scalp, and magnetoencephalography (MEG) data
∗ Correspondence to: Srikantan S. Nagarajan, Department of Radiology, University of California at San Francisco,†

513 Parnassus Avenue 5362, San Francisco, CA 94143-0628, U.S.A.
E-mail: srikantan.nagarajan@radiology.ucsf.edu, sri@radiology.ucsf.edu

Contract/grant sponsor: NIH; contract/grant numbers: R01DC004855, R01DC006435

Copyright q

2007 John Wiley & Sons, Ltd.

Received 19 April 2007
Accepted 23 April 2007

S. S. NAGARAJAN ET AL.

are obtained by measuring magnetic fields using a SQUID array surrounding the head. Amongexisting techniques for non-invasive mapping of brain functions, both MEG and EEG have thehighest temporal resolution. Both are used by basic neuroscientists in studies of brain functions.
They are also used by clinicians, most commonly in patients suffering from brain tumors andepilepsy. In brain tumor patients, MEG is used to map the cognitive function of the tumor area
and of neighboring areas, in order to guide neurosurgical planning, navigation, and tumor resection. Similarly, in epilepsy patients, MEG and EEG are often used to map where epileptic
activity originates and to map the cognitive function of brain regions surrounding epileptogenic
zones.
However, current techniques for functional brain mapping using MEG and EEG suffer fromimportant shortcomings. The data captured by the sensor array arise not only from signal frombrain sources located in areas of interest, but also from other sources, termed interference sources.
These include sources in other brain areas, such as spontaneous brain activity, biological sources
outside the brain, such as eye blinks, and non-biological sources, such as power lines. Signals frominterference sources overlap with those from the brain sources of interest, making it difficult toaccurately reconstruct the activity of desired brain areas. The task of removing interference signals
from the sensor data is termed interference suppression.
This paper focuses on the stimulus-evoked experimental paradigm, which is extremely popularin EEG and MEG studies. In this paradigm, a stimulus is presented to the subject at a series of
equally spaced time points. Each presentation produces activity in a set of brain sources, which
generates an electromagnetic field captured by the sensor array. Those data constitute the stimulusevoked response, and analyzing them can help to gain some insight into the mechanism used by
the brain to process the stimulus and similar sensory inputs. Perhaps the most important use of
stimulus evoked responses is to identify the brain locations of the sources evoked by the stimulus.
Unfortunately, the presence of interference often results in very inaccurate estimates of thoselocations.
Many approaches to the problem of interference suppression in stimulus-evoked responses havebeen taken, with varying degrees of success. One common method is using a large number ofstimulus presentations (100–200), also termed trials, and averaging the response across trials.
The underlying assumption there is that interference signals in different trials are statisticallyindependent, whereas evoked signals are not. Hence, averaging over sufficiently many trials would
minimize interference and reveal the clean-evoked response. However, the required large numberof trials results in several drawbacks. Since subjects can typically tolerate only 1–2 h of recordingsin the sensor array, the number of stimulus conditions that can be obtained within an experimentis limited. Furthermore, although the evoked response may vary a little across a small number oftrials, it could be non-stationary over a large number of trials. In such cases, averaging would yield
an inaccurate estimate of the evoked response. Moreover, many rapid brain processes that occur
within the course of single trials or a small number of trials cannot be examined by averaging
across many trials.
Data-driven approaches such as principal component analysis (PCA), Wiener filtering, matchedfiltering, and more recently, independent component analysis (ICA), have also been used forinterference suppression [1–3]. Some disadvantages of such approaches include the need to makesubjective choices when running them, such as setting the threshold in PCA and selecting relevant
components in ICA. An important drawback of most of those methods is their inability to exploit
the pre-/post-stimulus partition of the data (see below). In the experiments section, we demonstratethat the new technique presented here significantly outperforms those methods.
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

A PROBABILISTIC ALGORITHM FOR ROBUST INTERFERENCE SUPPRESSION

This paper presents a new technique for interference suppression in stimulus-evoked EEG/MEG
data. Our approach to this problem is formulated in the framework of probabilistic graphical models
with hidden variables, which has been developed considerably during the last decade in the fields
of machine learning and statistics. In this approach, we describe the observed sensor data interms of three types of unobserved signals, arising from evoked sources, interference sources,
and sensor noise. Those signals are described in our model by hidden variables with their own
probability distribution and depend on the sources via an appropriate probability distribution,
derived from the physics of the problem. The model exploits the fact that the data are partitioned
into two periods: pre-stimulus period, where the data include just the response of interference andsensor noise sources, and post-stimulus period, where the data also include the response of evoked
sources. Combining those distributions, we obtain a probabilistic model for the sensor data. We
present a variational Bayesian expectation-maximization (VB-EM) algorithm that infers the model
parameters from data. VB-EM is an extension of standard EM that has two major advantages: (1)it automatically infers the optimal number of interference and evoked sources required to explain
the sensor data and (2) it computes a full posterior distribution over model parameters, rather thana point estimate, which effectively prevents overfitting.
The paper is organized as follows. The probabilistic graphical model, termed partitioned factor
analysis (PFA), is defined in mathematical terms in the next section. Section 3 presents the VB-EM
algorithm for inferring this model from data. Section 4 provides an estimator for the clean-evoked
response, i.e. the contribution of the evoked sources alone to the sensor data, using the model toremove the contribution of the interference sources. This section also presents an automaticallyregularized estimator of the correlation matrix of the clean-evoked response. Section 5 demonstrates,
using real and simulated data, that the algorithm provides interference-robust estimates of the timecourse of the stimulus-evoked response. Section 6 concludes with a discussion of our results andof extensions to PFA.

2. PFA PROBABILISTIC GRAPHICAL MODEL
This section presents the PFA probabilistic graphical model, which is the focus of this paper.
The PFA model describes observed EEG/MEG sensor data in terms of three types of underlying,
unobserved signals: (1) signals arising from stimulus-evoked sources; (2) signals arising frominterference sources; and (3) sensor noise signals. The model is inferred from data by an algorithm
presented in the next section. Following inference, the model is used to separate the evoked source
signals from those of the interference sources and from sensor noise, thus providing a clean versionof the evoked response. In addition, it produces a regularized correlation matrix of the clean-evoked
response, which facilitates localization.
Let yin denote the signal recorded by sensor i = 1 : M y at time n = 1 : N . We assume that thesesignals arise from Mx evoked factors and Mu interference factors that are combined linearly. Letx jn denote the signal of evoked factor j = 1 : Mx , and let u jn denote the signal of interference
factor j = 1 : Mu , both at time n. We use the term factor rather than source for a reason explained
below. Let Ai j denote the evoked mixing matrix, and let Bi j denote the interference mixing matrix.
Those matrices contain the coefficients of the linear combination of the factors that produces thedata. They are analogous to the factor loading matrix in the factor analysis model. Let vin denote
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

S. S. NAGARAJAN ET AL.

the noise signal on sensor i. Mathematically
yin =

Mx


Ai j x jn +

j=1

Mu


Bi j u jn + vin

(1)

j=1

We use an evoked stimulus paradigm, where a stimulus is presented at a specific time, termed
the stimulus onset time. The stimulus onset time is defined as n = N0 + 1. The period preceding theonset n = 1 : N0 is termed pre-stimulus period, and the period following the onset n = N0 + 1 : N
is termed post-stimulus period. We assume that the evoked factors are active only post-stimulus
and satisfy x jn = 0 before its onset. Hence, using vector notations
Bu n + vn ,n = 1 : N0
yn =(2)
Axn + Bu n + vn , n = N0 + 1 : N
To turn (2) into a probabilistic model, each signal must be modelled by a probability distribution.
Here, each evoked factor is modelled by a Gaussian distribution‡ with zero mean and unit precisionp(x jn ) = N(x jn | 0, 1)(3)

We model the factors as mutually statistically independent, hencep(xn ) =M
xp(x jn ) = N(xn | 0, I )

(4)

j=1

For interference signals, we also employ a Gaussian model. Each interference factor is modelled
by a zero-mean Gaussian distribution with unit precision, p(u jn ) = N(u jn | 0, 1). PFA describes
the factors as independent:p(u n ) =Mup(u jn ) = N(u n | 0, I )

(5)

j=1

The sensor noise is modelled by a zero-mean Gaussian distribution with a diagonal precision
matrix ,p(vn ) = N(vn | 0, )(6)From (2) we obtain p(yn | xn , u n ) = p(vn ), where we substitute vn = yn − Axn − Bu n with xn = 0
for n = 1 : N0 . Hence, we obtain the distribution of the sensor signals conditioned on the evoked
and interference factors,
n = 1 : N0
N(yn | Bu n , ),p(yn | xn , u n , A, B) =(7)
N(yn | Axn + Bu n , ), n = N0 + 1 : N‡

A Gaussian distribution over a random vector z with mean  and precision matrix  is defined by



1/2


N(x |,) =  2 1exp[− 2 (x−)T (x−)]

The precision matrix is defined as the inverse of the covariance matrix.
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

A PROBABILISTIC ALGORITHM FOR ROBUST INTERFERENCE SUPPRESSIONPFA also makes an i.i.d. assumption, meaning the signals at different time points are independent.
Hence,
p(y | x, u, A, B) =Np(yn | xn , u n , A, B)

n=1

p(x) =
p(u) =

N
n=N0 +1Np(xn )(8)p(u n )

n=1

where y, x, u denote collectively the signals yn , xn , u n at all time points. The i.i.d. assumption ismade for simplicity, and implies that the algorithm presented below can exploit the spatial statisticsof the data but not their temporal statistics.
To complete the definition of PFA, we must specify prior distributions over the model parameters.
For the noise precision matrix , we choose a flat prior, p() = const. For the mixing matrices
A, B, we use a conjugate prior. A prior distribution is termed conjugate w.r.t. a model when its
functional form is identical to that of the posterior distribution (see the discussion below equation(A15)). We choose a prior where all matrix elements are independent zero-mean Gaussiansp(A) = N(Ai j | 0, i  j )ijp(B) =(9)

N(Bi j | 0, i  j )ij

and the precision of the i jth matrix element is proportional to the noise precision i on sensor i.
It is the  dependence which makes this prior conjugate. (It can be shown that in the limit of zero
sensor noise  → ∞; the impact of the prior on the posterior mean of A, B would vanish in theabsence of this dependence, which would be undesirable.) The proportionality constants  j and j constitute the parameters of the prior, a.k.a. hyperparameters. Equations (8), (9) together withequations (4), (5), (7) fully define the PFA model.

3. INFERRING THE PFA MODEL FROM DATA: A VB-EM ALGORITHM
This section presents an algorithm that infers the PFA model from data. PFA is a probabilistic
model with hidden variables, since the evoked and interference factors are not directly observable.
We use an extended version of the expectation maximization (EM) algorithm to infer the modelfrom data. This version is termed VB-EM.
Standard EM computes the most likely parameter value given the observed data, a.k.a. themaximum a posteriori (MAP) estimate. In contrast, VB-EM considers all possible parameter values,and computes the probability of each value conditioned on the observed data. VB-EM therefore
treats hidden variables and parameters on equal footing by computing posterior distributions forboth quantities. One may, however, choose to compute a posterior only over one set of model
parameters, while computing just a MAP estimate for the other set.
VB-EM is an iterative algorithm, where each iteration consists of an E-step and an M-step. The
E-step computes the sufficient statistics (SS) of the hidden variables, and the M-step computes
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

S. S. NAGARAJAN ET AL.

the SS of the parameters. (SS of an unobserved variable are quantities that define its posterior
distribution.) The algorithm is iterated to convergence, which is guaranteed.
The VB-EM algorithm has several advantages compared with standard EM. It is more robust to
overfitting, which can be a significant problem when working with high-dimensional but relatively
short time series, as we do in this paper. It produces automatically regularized estimators, suchas for the evoked response correlation matrix, whereas standard EM produces under-conditioned
ones. In addition, the variance of the posterior distribution it computes (essentially the estimator’svariance or squared error) provides a measure of the range of parameter values compatible withthe data.
We now describe the VB-EM algorithm for the PFA model. A full derivation is provided inAppendix A.
3.1. E-step
The E-step of VB-EM computes the SS for the hidden variables conditioned on the data. For the
pre-stimulus period n = 1 : N0 , the hidden variables are the interference factors u n . Compute theirposterior mean ū n and covariance  byū n =  B̄ T yn = ( B̄ T  B̄ + I + M y  B B )−1(10)

where B̄ are  B B are computed in the M-step by equations (15)–(17). B̄ is the posterior meanof the interference mixing matrix, and  B B is related to its posterior covariance (specifically, theposterior covariance of the ith row of B is  B B /i ; see Appendix A).
For the post-stimulus period n = N0 + 1 : N , the hidden variables include the evoked andinterference factors xn , u n . To simplify the notation, we combine the evoked and interference
factors into a single vector, and their mixing matrices into a single matrix. Let L  = Mx + Mu
be the combined number of evoked and interference factors. Let A denote the M y × L  matrix
containing A and B, and let xn denote the L  × 1 vector containing xn and u n xn, A = (A B)
xn =
(11)
un
The SS are computed as follows. At time n, compute the posterior means x̄ n and ū n of theevoked and interference factors, and their posterior covariance , byx̄n =  ĀT yn = ( ĀT  Ā + I + M y )−1

(12)Here, as in (11), we have combined the posterior means of the factors into a single vector x̄ n , andthe posterior means of the mixing matrices into a single matrix Ā , x̄nx̄n =, Ā = ( Ā B̄)
(13)
ū n
where Ā, B̄,  are computed in the M-step by equations (15)–(17). As explained in Appendix A,
/i is the posterior covariance of row i of A .
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

A PROBABILISTIC ALGORITHM FOR ROBUST INTERFERENCE SUPPRESSION

The covariances x x and uu of the evoked and interference factors, and their cross-covariancexu , are obtained by appropriately dividing  into quadrantsx x xu
=(14)
Txu uu
where x x is the top left Mx × Mx block of , xu is the top right Mx × Mu block, and uu isthe bottom right Mu × Mu block. These covariances are used in the M-step.
3.2. M-step
The M-step of VB-EM computes the SS for the model parameters conditioned on the data. We
divide the parameters into two sets. The first set includes the mixing matrices A, B, for which
we compute full posterior distributions. The second set includes the noise precision  and thehyperparameters matrices , , for which we compute MAP estimates.
Compute the posterior means of the mixing matrices byĀ = R yx 

(15)B̄ = R yu 
where


=

Rx x + 

Rxu

T
Rxu

Ruu + −1
(16)

The quantities R yx , R yu , Rx x , Rxu , Ruu are posterior correlations between the factors and the data
and among the factors themselves, and are computed below. ,  are diagonal matrices with thehyperparameters  j ,  j on the diagonal.
The covariances  A A and  B B corresponding to the evoked and interference mixing matrix(see Appendix A), and  AB corresponding to their cross-covariance, are obtained by appropriately
dividing  into quadrants A A  AB
=(17)
TAB  B B
where  A A is the top left L × L block of ,  AB is the top right L × M block, and  B B is thebottom right M × M block.
Next, use those covariances to update the hyperparameter matrices ,  by


1 T
−1
 = diagĀ  Ā +  A A
My
(18)


1 T
−1
B̄  B̄ +  B B
 = diagMy
and to update the noise precision matrix  by−1 =Copyright q

1diag(R yy − Ā R Tyx − B̄ R Tyu )N2007 John Wiley & Sons, Ltd.

(19)
Statist. Med. (in press)
DOI: 10.1002/sim

S. S. NAGARAJAN ET AL.

3.2.1. Posterior means and correlations of the factors. Here we compute the posterior correlations,used above, between the factors and the data and among the factors themselves. Let x̄ n = xn  andū n = u n  denote the posterior mean of the evoked and interference factors. During the pre-stimulus
period n = 1 : N0 , x̄n = 0 and ū n is given by (10). During the post-stimulus period n = N0 + 1 : N ,they are givenby (12), (13).

Let R yx = n yn xnT  and R yu = n yn u Tn  denote the data-evoked and data-interference posterior correlations. Then
R yx =
R yu =

N
n=N0 +1N
n=1

yn x̄nT(20)

yn ū Tn
Let Rx x = n xn xnT , Rxu = n xn u Tn , and Ruu = n u n u Tn  denote the evoked–evoked,
evoked–interference, and interference–interference posterior correlations. Then
Rx x =
Rxu =
Ruu =

N


(x̄n x̄nT + x x )n=N0 +1N(x̄n ū Tn + xu )(21)

n=N0 +1N0(ū n ū Tn + ) +

n=1

N


(ū n ū Tn + uu )n=N0 +1

using the factors covariances (14).
Finally, let R yy denote the data–data correlation
R yy =

N
n=1

yn ynT(22)

4. ESTIMATING CLEAN-EVOKED RESPONSE AND ITS CORRELATION MATRIXIn this section, we present two sets of estimators computed by the PFA model after inferring itfrom data. The first estimator computes the clean-evoked response. The second estimator computes
a well-conditioned correlation matrix for the signals obtained by the first estimator.
Let z in denote the combined contribution from all evoked factors to sensor signal i. Thenz in =Mx

Ai j x jn

(23)

j=1

Let z̄ in denote the estimators of z in . This means that z̄ in = z in , where the average is w.r.t. theposterior over A, x. Computing this estimate amounts to obtaining a clean version of the combined
contribution of the evoked factors, removing contributions from interference factors and sensor
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

A PROBABILISTIC ALGORITHM FOR ROBUST INTERFERENCE SUPPRESSION

noise. We obtainz̄ in =MxĀi j x̄ jn

(24)

j=1

Next, consider the correlation matrix of the evoked response, which is a required input for localization algorithms such as beamforming. Let C denote the correlation of the combined contribution
from all evoked factors. ThenC=N

n=N0 +1z n z nT(25)

Let C̄ denote the estimator of C. This means, as above, that C̄ = C. We obtain
C = Ā Rx x ĀT + −1 Tr(Rx x  A A )(26)

We point out an important fact about the estimated correlation matrix C̄. It is always wellconditioned, due to the diagonal  A A term. Hence, the VB-EM approach automatically produces
regularizedcorrelation matrix. Note that the correlation matrix obtained directly from the signal
estimates, n z̄ n z̄ nT , is under-conditioned.
5. MODEL-ORDER SELECTION, INITIALIZATION AND COMPLEXITY
One advantage of the algorithm presented here is that it offers a principled method of modelorder selection. Model-order selection in PFA algorithm refers to the choice of Mx and Mu .
The MAP estimates of the hyperparameters of the mixing matrices can be used to estimate thenumber of factors by thresholding. Alternatively, we can compute the maximum of the posterior
over model structure q(Mx , Mu |y), which is equivalent to maximizing the marginal log likelihood
log p(y|Mx , Mu ). The marginal log likelihood obtained by integrating over all hidden variables
is also referred to as the evidence. The evidence penalizes complexity and corresponds to theBayesian information criterion (BIC) and the minimum description length (MDL) for infinite data[4]. It can be shown that the evidence is lower bounded by a free energy objective function F, as
defined in equations (A5) and (A6). Therefore, after computing F for different model orders Mx
and Mu , we can chooseM̄x , M̄u = argmax F(Mx , Mu )Mx ,MuAlthough, the proposed algorithm is fairly robust to initialization, the specific initializations of
the parameters that we use in the Results section are as follows. We initialize the mixing matrix
B to the dominant eigen-vectors of the data obtained in the pre-stimulus period. The evokedfactor mixing matrix A is initialized as the dominant eigen-vectors of the post-stimulus data after
pre-whitening with the pre-stimulus data covariance.  is initialized to be uniform across sensors
and equal to the inverse of the least-significant eigenvalue of the pre-stimulus data.  and  are
initially assumed to be identity matrices.
For each iteration of the algorithm, the computational complexity of estimation of the PFAgraphical model is O(N ∗ (Mx + Mu ) ∗ M y ). So, PFA is linear in the number of time samples,
sensors and factors.
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

S. S. NAGARAJAN ET AL.

6. RESULTS
6.1. Simulations
Figure 1 shows an example of performance for the proposed interference suppression algorithm
on simulated data. The top row shows simulated noisy MEG data created assuming three brain
sources and 25 interference sources and 275 sensors. The middle row shows the true signal thatis present in the post-stimulus period within the noisy MEG data. The bottom row shows theestimated signal extracted by PFA. When the true signal y ∗ is known, denoising performance canbe quantified using the output signal-to-noise/interference ratio (SNIR)
N
∗2
My
1 
n=N0 y m,nSNIR =
10 log10  N(dB)∗
2
M y m=1
n=N0 (y m,n − ȳm,n )For the example shown, the input SNIR is −13 dB and the output SNIR is −2 dB.
In more extensive simulations, we compare interference suppression performance for the proposed probabilistic algorithm with five other standard methods used in practice—PCA [5], Wiener
Filtering [2], ICA using TDSEP [6] and/or FastICA [7]), and trial averaging. TDSEP and FastICA
were chosen as the representative ICA methods based on their low computational complexity.
Furthermore, when there are more than about 50 sensors, as is typically for high-resolution EEG,
MEG, or magnetocardiography (MCG) systems, TDSEP and FastICA do not require additional
dimensionality reduction. We report the better results between these two ICA algorithms. With
the exception of the trial mean, all the above interference suppression methods are spatial filtering
methods that apply a linear transformation that is applied to the observed data to obtain an estimateof the underlying signal.
The proposed algorithm, and the comparison methods mentioned above, could be applied either
to concatenated single-trial data or to trial-averaged data. For interference suppression performance,

Figure 1. Example of performance of the proposed algorithm.
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

A PROBABILISTIC ALGORITHM FOR ROBUST INTERFERENCE SUPPRESSION

Figure 2. Output SNIR as a function of the input SIR for 10 trials and input SNR = 0 dB.

we first apply each method to the trial-averaged data so that we can directly compare it with thetrial mean. In some cases (as noted), we also apply the interference suppression on single-trial
data and then compute the trial average.
For the simulation results below, there are 1000 data points per trial (the first 63% of which corresponds to the inactive period), M y = 132 sensors, Mx = 2 factors, Ms = 2 sources, and Mu = 1000
interference signals. Results shown represent the mean over 10 Monte Carlo repetitions and theerror bars are used to indicate one standard error of the mean. The input signal-to-interference ratio(SIR) is the ratio of the power of the factors to the power of the interferences (measured in sensor
space). Likewise, the input signal-to-noise ratios (SNR) is the ratio of the power of the factors to
the power of the additive noise. The number of factors, Mx , must be specified for all denoising
methods except the trial mean. The proposed method must also be supplied with a known numberof interference signals, Mu . To simplify the comparisons, it is assumed that the number of factors
is the true number and the number of interference signals Mu = 50.
Figure 2 shows the interference suppression performance as a function of the input SIR. Only
The input SNR is held constant at 0 dB and the number of trials is 10. All of the methods performbetter than the trial mean. PFA performs the best across all input SIR. The performances of both
PCA and Wiener approach that of PFA as the input SIR increases.
Figure 3 shows the interference suppression performance as a function of the number of trials.
The input SIR and input SNR are held constant at −5 and 0 dB, respectively. In this figure thetrial mean outperforms TDSEP and PFA outperforms the other four methods.
6.2. Model-order selection
We demonstrate robustness to model-order selection using the PFA criterion with simulated data.
Figures 4 and 5 show the results of using the PFA criterion, the evidence under the variational
approximation, as a function of model order. For these two figures there are Mx∗ = 2 sources,
Mu∗ = 10 interferences, 1000 data per trial, and 10 trials of data. Figure 4 plots the PFA criterion
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

S. S. NAGARAJAN ET AL.

Figure 3. Output SNIR as a function of the number of trials for input SIR = −5 dB and input SNR = 0 dB.

Figure 4. Performance as a function of Mu , where Mx is assumed to be 2. The true values of Mx
and Mu are 2 and 10, respectively.

as a function of Mu , where it is assumed that Mx = 2. Also shown are the plots of output SNIRand the amplitude of the estimated (inverse) hyperparameters corresponding to the scaled varianceof the columns of the mixing matrix. Note that the output SNIR asymptotes for higher modelorders because the columns of the mixing matrices comprise elements with values near zero. The
PFA criterion matches the output SNIR quite well (note that the precise value of the PFA criterion
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

A PROBABILISTIC ALGORITHM FOR ROBUST INTERFERENCE SUPPRESSION

Figure 5. Performance as a function of Mx , where Mu is assumed to be 15. The true values of Mx
and Mu are 2 and 10, respectively.

can be obtained for real data, whereas the precise value of the output SNIR cannot). Figure 5
plots the PFA criterion as a function of Mx , where it is assumed that Mu = 15. The plot of the
PFA criterion versus model-order peaks at the correct value of M̄x = 2, where the output SNIR
also peaks. Moreover, increasing the specified model order beyond the true model order does notcontribute to significant deterioration in performance, hence our use of the term ‘robust interference
suppression’.
6.3. PFA as preprocessing for ICA
The stimulus-evoked factors in PFA can be subsequently separated using ICA algorithms. Here, we
compare the performance on source separation using ICA after preprocessing with the proposedand comparative algorithms. Source extraction performance is measured using the output sourceto-distortion ratio (SDR), where the distortion for source estimate m includes noise, interference,
and all sources except one. For the case of no permutations, the SDR is defined by
SDR =

Ms
1 
SDRm (dB)
Ms m=1

where⎛
SDRm = 10 log10

Copyright q

Ms ⎜
1 
⎜
Ms m  =1 ⎝2007 John Wiley & Sons, Ltd.

⎞
⎟
1
(dB)

⎟
2
 N
⎠
2−
 n=N0 +1 sm,n s̄m  ,n 
N − N0Statist. Med. (in press)
DOI: 10.1002/sim

S. S. NAGARAJAN ET AL.

sn = W −1 xn is the true source vector at time n, and both sm,n and s̄m,n are normalized to have unit
variance. The definition above is easily extended to account for any possible permutation. This
metric reflects the performance of both the interference suppression/dimension reduction algorithm
and the ICA algorithm. The interference suppression method accounts for all differences in SDR
performance below, since, for each experiment, the same ICA algorithm is used. In general, we
found that TDSEP performed better than FastICA for denoising and FastICA performed better
than TDSEP for source extraction. For TDSEP and FastICA, the source subspace is automaticallydetermined by selecting the components that have the largest ratio of active power to inactive
power. The first component is given byN
2
n=N s̄m,nm̄ 1 = argmax  N −1002mn=1 s̄m,n
and the subsequent Mx − 1 components are found in a similar manner.
Figure 6 shows the source extraction performance as a function of the input SIR. The input SNR
is 0 dB and the number of trials is 10. The non-ICA denoising methods are used to reduce thedimensionality of the data from 132 to 2 prior to applying the ICA algorithm, which in this case is
FastICA. Also shown are the results for FastICA when no dimension reduction method is used. PFA
produces the best overall results and is the least sensitive to input SIR. The results reported here for
FastICA (with no dimension reduction) indicate that denoising/dimension reduction preprocessing
is advantageous when the input SIR is low (<10 dB).
Figure 7 shows the source extraction performance as a function of the number of trials. As
before, the input SIR and input SNR are −5 and 0 dB, respectively, and the results of using
FastICA with no dimension reduction are included. PFA (combined with FastICA) performs thebest and FastICA (with no dimension reduction) performs the worst. These results indicate that 10
trials of 1000 data points per trial is already sufficiently large so that no improvement in separation

Figure 6. Output SDR as a function of the input SIR for 10 trials and input SNR = 0 dB.
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

A PROBABILISTIC ALGORITHM FOR ROBUST INTERFERENCE SUPPRESSION

Figure 7. Output SDR as a function of the number of trials for input SIR = −5 dB and input SNR = 0 dB.

performance is obtained by additional increases in data length. This is not expected to be the caseif the input SIR and/or input SNR are increased.
6.4. Real data
An example of performance of the proposed interference suppression algorithm on auditoryevoked magnetic fields measured across the whole head obtained from a 275-channel sensor arrayin response to a 1 kHz tone pip is shown in Figure 8. For these data, M y = 274, the data lengthis 720 samples per trial (170 pre-stim), and there are 109 trials. Averaged data from 20 trial
averages are noisy as shown in the top left for select channels. The output of PFA is shown inthe middle left. Also shown is the response obtained from averaging 109 trials. It can be seen thatthe response from 20 trials does not resemble the 109 trial average (shown in the bottom left)
suggesting trial-to-trial variability or non-stationarity in the evoked response over 109 trials. The
right column shows waveforms for the noisy input (thin lines) and interference-suppressed output(thick lines) for selected individual channel waveforms.
Figure 9 shows the denoising for a different MEG data set. Here, we examine the stimulusevoked response to a somatosensory stimulus with M y = 274 and the data length is 361 samples
per trial. In this figure, 0 ms corresponds to N0 , which is the onset of the stimulus. We assume
that Mx is 2, Mu is 50. For comparison, PCA denoising on the 10-trial average is also shown,as well as the average across 525 trials. It can be seen that PFA performs adequate interference
suppression of the evoked response.
Quantifying performance of interference suppression with real data is difficult because output
SNIR and SDR can be easily computed only for simulated data since y ∗ , s are not known for realdata. The output SNIR can, however, be used with real data if yn∗ = Ax n can be approximated. In
this latter example, the average response obtained from 525 trials appears to be more similar to the
response to 10 trials, suggesting stationarity in the evoked response. Furthermore, five principal
components explain 97% of the total energy of the trial-averaged data. Therefore, for this real data
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

S. S. NAGARAJAN ET AL.

Figure 8. Example of performance on real auditory-evoked magnetic field data.

we replace yn∗ with the sensor signals due to the five principal components of the trial-averaged
sensor data obtained from 525 trials.
Figure 10 shows denoising performance as a function of the number of trials using the abovementioned procedure. None of the estimated sources produced by the ICA method resembled
the desired signals even when the number of trials was increased to 50. The performance of
ICA denoising depends on being able to correctly select the Mx sources and the results showpoor performance of ICA on this data. Results for PFA, PCA, and Wiener are better than those
produced by the trial mean (when the trial mean uses the same number of trials). PFA performs
the best of these methods, although PCA performs almost identically when the number of trials
equals or exceeds 30.
Figures 11 and 12 show the results of model-order selection for the auditory MEG data set
shown above. For these two figures, only the 20 trials of data are used. Figure 11 plots theevidence as a function of Mu , where it is assumed that Mx = 2. Figure 12 plots the evidence asa function of Mx , where it is assumed that Mu = 25. Also shown are the plots of the amplitudeof the associated inverse hyperparameters. The model orders that maximize the PFA criterion areM̄u = 25 and M̄x = 5. It can be seen that the evidence peaks for small model orders and that theposterior estimates of many of the inverse hyperparameters are zero, thereby demonstrating thebuilt-in model-order robustness of the PFA inference algorithm.
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

A PROBABILISTIC ALGORITHM FOR ROBUST INTERFERENCE SUPPRESSION

Figure 9. Time series after interference suppression of real MEG data from somatosensory cortex.
(A) Trial Mean (10 trials); (B) PCA (10 trials); (C) PFA (10 trials); and (D) Trial mean (525 trials).

Figure 10. Output SNIR as a function of the number of trials for real MEG data.

Figure 13 shows the interference suppression of real EEG data, where M y = 119, the data lengthis 720 samples per trial (170 in the pre-stimulus period), the number of trials is 120, and the data
were the response to an auditory 1 kHz tone. Results for the trial mean are shown for both 10 and
120 trials. Notice that data contain a large 60 Hz contribution (further examination reveals that
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

S. S. NAGARAJAN ET AL.

Figure 11. Model-order selection as a function of Mu , where Mx is assumed to be 2.

Figure 12. Model-order selection as a function of Mx , where Mu is assumed to be 25.

most of the line noise is concentrated in three channels). The proposed algorithm uses 10 trials
and assumes that there are Mx = 5 sources and Mu = 25 interferences. Since the 60 Hz oscillations
occur during both pre-stimulus and post-stimulus periods, our algorithm treats it as an interference
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

A PROBABILISTIC ALGORITHM FOR ROBUST INTERFERENCE SUPPRESSION

Figure 13. Interference suppression of real EEG data. (A) Trial Mean (10 trials); (B) PFA(10 trials); and (C) Trial mean (120 trials).

signal and is therefore able to remove it successfully. Simple temporal filtering, which can also
be used to remove the line noise, will necessarily repress other activity in and near the 60 Hz
frequency, whereas this does not occur with our algorithm. The trial mean, on the other hand,is unable to remove the line noise since the stimulus onset is approximately synchronous withthe line noise. The P1, N1, and P2 responses are clearly visible in the output after interference
suppression (the convention of inverting the polarity, commonly used in EEG analyses, is not usedhere).

7. DISCUSSION
The robustness of the proposed algorithm to the choice of the maximal model orders, Mx , andMu , can be explained using a process known as automatic relevance determination (ARD). The
hyperparameters represent the inverse power of the associated factor/interference signal. When themodel order is chosen larger than necessary, the hyperparameters associated with the redundantsignals approach infinity [8–10]. As a hyperparameter approaches infinity, the observations canbe explained without the associated factor/interference signal. The cost of overestimating themodel order is that the computational complexity increases as either Mu or Mx is increased.
The tendency of the hyperparameters to approach infinity can be used to estimate the two modelorders. The most straightforward way to estimate Mu , for example, is to count the number of
diagonal elements of  that have an inverse value less than a given threshold. In the previous
section, we showed results using the evidence, which does not require an arbitrary selection of athreshold.
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

S. S. NAGARAJAN ET AL.

The proposed algorithm, as given above, has a potential problem of identifiability between theestimation of A and B, especially if the amount of data in the pre-stimulus period is small and both
A and B are primarily estimated from the post-stimulus period (data may have equal likelihood toarise from a source factor and from an interference factor). To avoid this problem, we perform a
two-step procedure for PFA. In the first step, we estimate the interference factor mixing matrix andthe sensor noise precision from data in the pre-stimulus period. In this case, the update equationfor  uses only the pre-stimulus data and is1
diag(R yy − B̄ R yu T )(27)
N
The update rules for B and  are the same as listed in equation (20), with a modified  B B = (Ruu +)−1 . Subsequently, for post-stimulus data, we freeze the above parameters and estimate A usinga modified update rule,−1 =Ā = (R yx − B R yx ) A A(28)

where  A A = (Rx x + )−1 . All other update rules are identical to those listed above.
Furthermore, the proposed model currently assumes that the interferences are statistically stationary between the pre- and post-stimulus periods. However, we can relax this assumption andmodel non-stationary changes in the power of interference if there are no changes in the location of the interferences. We assume that, in the post-stimulus period, the probability distribution
 Muof the interference factors is p(u n ) = m=1p(u m,n ) = N(u n |0, ), where  is a diagonal precision matrix that is equal to the inverse of the power fluctuations of the interference in the
post-stimulus period. In this case, we can learn  from the post-stimulus period using the update
rule −1 = diag((1/N )Ruu ), where Ruu is calculated only for the post-stimulus period.
The algorithm currently assumes that the prior distributions for evoked and interference factors
are i.i.d. and invariant to the time-index permutation. However, this does not appear to impact
performance because in all the simulations presented in the paper both the background sources wereassumed to be sinusoidal (with bimodal distributions) or damped sinusoids (with super-Gaussian
distributions), rather than Gaussians as assumed in the model. Moreover, since the performanceof the algorithm is also good on real bioelectromagnetic data, where the interference factors areindeed oscillatory, the algorithm has some degree of robustness with respect to assumptions about
the prior distribution of interference and evoked factors. Since estimation is data dependent, if
the data suggest that factors have temporal continuity, then the estimated factors will have somesmoothness. Nevertheless, an algorithm that is able to exploit temporal correlation in factors could
potentially be more powerful. We are currently pursuing such an extension, using several different
models that incorporate temporal statistics of the evoked and interference factors, whose parameters
are inferred from data. Algorithms derived from such models perform interference suppression
using not just spatial but also spatio-temporal filtering. On a separate note, since bioelectromagnetic
data are often non-Gaussian, we are currently extending the model to incorporate non-Gaussian
factor models.
APPENDIX A: THE VB-EM ALGORITHM
This section outlines the derivation of the VB-EM algorithm that infers the PFA model fromdata.
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

A PROBABILISTIC ALGORITHM FOR ROBUST INTERFERENCE SUPPRESSION

A.1. Model
The full joint distribution of the PFA model is given byp(y, x, u, A, B) = p(y | x, u, A, B) p(x) p(u) p(A) p(B)(A1)

together with equations (5), (7), (8).
A.2. Variational Bayesian inference
The Bayesian approach, as discussed above, treats hidden variables and parameters on equal
footing: both are unobserved quantities for which posterior distributions must be computed. A
direct application of Bayes rule to the PFA model would compute the joint posterior over thehidden variables x, u and parameters A, B
p(x, u, A, B | y) =1p(y, x, u, A, B)
p(y)

(A2)

where the normalization constant p(y), termed the marginal likelihood, is obtained by integratingover all other variables
p(y) = dx du dA dB p(y, x, u, A, B)(A3)
However, this exact posterior is computationally intractable, because the integral above cannot be
obtained in closed form.
The VB approach approximates this posterior using a variational technique. The idea is torequire the approximate posterior to have a particular factorized form, and then optimize it byminimizing
the Kullback–Leibler (KL) distance from the factorized form to the exact posterior
q log( p/q) [11].
Here, we choose a form which factorizes the hidden variables from the parameters given thedatap(x, u, A, B | y) ≈ q(x, u, A, B | y) = q(x, u | y)q(A, B | y)(A4)

It is worth emphasizing that (1) beyond the factorization assumption, we make no further approximation when computing q, and (2) the factorized form still allows correlations among x, u, aswell as among the matrix elements of A, B, conditioned on the data.
Rather than minimize the KL distance directly, it is convenient to start from an objective function
defined byF[q] = dx du dA dBq(x, u, A, B | y)[log p(y, x, u, A, B) − log q(x, u, A, B | y)](A5)
It can be shown that
F[q] = log p(y) − K L[q(x, u, A, B | y) || p(x, u, A, B | y)](A6)and, since the marginal likelihood p(y) is independent of q, maximizing F w.r.t. q is equivalent
to minimizing the KL distance. Furthermore, F is upper bounded by log p(y) because the KL
distance is always non-negative. Hence, any algorithm that successively maximizes F, such as
VB-EM, is guaranteed to converge.
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

S. S. NAGARAJAN ET AL.

A.3. Derivation of VB-EM
VB-EM is derived by alternately maximizing F w.r.t. the two components of the posterior q. In
the E-step one maximizes w.r.t. the posterior over hidden variables q(x, u | y), keeping the second
posterior fixed. In the M-step one maximizes w.r.t. the posterior over parameters q(A, B | y),keeping the first posterior fixed. When performing maximization, normalization of q must be
enforced by adding two Lagrange multiplier terms to F in (A5).
Maximization is performed by setting the gradients to zero:*F= log p(y, x, u, A, B)2 − log q(x, u | y) + C1 = 0*q(x, u | y)*F= log p(y, x, u, A, B)1 − log q(A, B | y) + C2 = 0*q(A, B | y)(A7)

where C1 , C2 are constants depending only on the data y. ·1 denotes averaging only w.r.t.
q(x, u | y), and ·2 denotes averaging only w.r.t. q(A, B | y). Hence, the posteriors are given by1
exp[log p(y, x, u, A, B)2 ]Z1
1exp[log p(y, x, u, A, B)1 ]q(A, B | y) =Z2q(x, u | y) =(A8)where Z 1 , Z 2 are normalization constants.
A.4. E-step
It follows from (A8) that the posterior over u, x factorizes over time, and has different pre- andpost-stimulus forms,
q(u, x | y) =

N0

n=1

q(u n | yn ) ·N
n=N0 +1

q(u n , xn | yn )(A9)It also follows that in the pre-stimulus period q(u n | yn ) is Gaussian in u n , and in the poststimulus period q(u n , xn | yn ) is Gaussian in u n , xn . To see this, consider log q(x, u | y) in(A8) and observe that it is a sum over n, where the nth element depends only on x n , u n and thedependence is quadratic.
For the pre-stimulus period we obtainq(u n | yn ) = N(u n | ū n , −1 )(A10)with mean ū n and covariance matrix  given by (10). (One first obtains  = (B T B + I )−1 ,and then performs the average using (A18).) For the post-stimulus period, the posterior is also
Gaussianq(xn , u n | yn ) = q(xn | yn ) = N(xn | x̄n , −1 )(A11)with mean x̄n and covariance matrix −1 given by (12) (as for  above, one first obtains
 = (AT A  + I )−1 , then applies (A18)).
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

A PROBABILISTIC ALGORITHM FOR ROBUST INTERFERENCE SUPPRESSION

It is useful to make explicit the correlations among the factors implied by their posteriors (A10),(A11). For the pre-stimulus period, we obtainu n u Tn  = ū n ū Tn + (A12)For the post-stimulus period, we obtain x n  = x̄n and xn xnT  = x̄n x̄nT + . In terms of xn , u nxn  = x̄n
u n  = ū n
xn xnT  = x̄n x̄nT + x x

(A13)

u n u Tn  = ū n ū Tn + uu
xn u Tn  = x̄n ū Tn + xu
where we have used (13), (14).
A.5. M-step
It follows from (A8) that the parameter posterior factorizes over the rows of the mixing matrices,and correlates their columns. Let wi denote a column vector containing the ith row of the combined
mixing matrix A = (A, B)⎛⎞· · · w1 · · ·⎜⎟⎜ · · · w2 · · · ⎟
⎜
⎟
⎟
(A14)
A = ⎜⎜⎟⎜· · · · · · · · ·⎟⎝⎠· · · w My···

so wij = Ai j . Then, the posterior over each row is Gaussianq(A, B | y) = q(A | y) =MyN(wi | w̄i , i −1 )

(A15)

i=1

with mean w̄ij = Āi j computed by (15). The precision matrix i −1 is computed using (16). To
see this, consider log q(A, B | y) in (A8) and observe that it is a sum over i, where the ith element
depends only on the ith rows of A, B and the dependence is quadratic.
It is now evident that p(A, B) of equation (9) is indeed a conjugate prior. Rewriting it in theformp(A, B) = p(A ) =MyN(wi | 0, i  )

(A16)

i=1

where  is a diagonal matrix with the hyperparameter matrices ,  on its diagonal, shows thatits functional form is identical to that of the posterior (A15), with −1 replacing  .
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

S. S. NAGARAJAN ET AL.

It is useful to make explicit the correlations among the elements of the mixing matrices impliedby their posterior (A15). They are Ai j Akl  = Āi j Ākl + ik  jl /i , or, in terms of A, B,
1( A A ) jl
i1
Bi j Bkl  = B̄i j B̄kl + ik ( B B ) jl
i1
Ai j Bkl  = Āi j B̄kl + ik ( AB ) jl
i

Ai j Akl  = Āi j Ākl + ik(A17)

where we used (17). It follows thatAT A = ĀT  Ā + M y  A A
B T B = B̄ T  B̄ + M y  B B(A18)

AT A  = ĀT  Ā + M y which are needed for (10), (12).
To obtain the update rules for the hyperparameters (18), observe that the part of the objectivefunction F (A5) that depends on ,  is
log p(A) + log p(B)

(A19)

where the averaging is w.r.t. the posterior q. Next, compute the derivative of this expression w.r.t.
,  and set it to zero. The solution of the resulting equation is (18). It is easier to first compute thederivative and then apply the average. Similarly, to obtain the update rule for the noise precision(19), observe that the part of F that depends on  is
log p(y | x, u, A, B) + log p(A) + log p(B)

(A20)

and set its derivative w.r.t.  to zero.
ACKNOWLEDGEMENT

This work was supported by NIH grants R01DC004855 and R01DC006435 to S. S. N.
REFERENCES
1. Ossadtchi A, Baillet S, Mosher JC, Thyerlei D, Sutherling W, Leahy RM. Automated interictal spike detection
and source localization in magnetoencephalography using independent components analysis and spatio-temporalclustering. Clinical Neurophysiology 2004; 115(3):508–522.
2. Ungan P, Basar E. Comparison of Wiener filtering and selective averaging of evoked potentials.
Electroencephalography and Clinical Neurophysiology 1976; 40(5):516–520.
3. Nagarajan S, Attius HT, Sekihara K, Hild II KE. Partitioned factor analysis for interference suppression andsource extraction. International Workshop on Independent Component Analysis and Signal Separation (ICA ’06),
Charleston, SC. Lecture Notes in Computer Science, vol. 3889. Springer: Berlin, 2006; 189–197.
4. Attias H. A variational Bayesian framework for graphical models. Advances in Neural Information Processing
Systems (NIPS ’99). MIT Press: Cambridge, MA, 2000; 209–215.
5. Jackson JE. A User’s Guide to Principal Components. Wiley: Hoboken, NJ, 2003.
Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim

A PROBABILISTIC ALGORITHM FOR ROBUST INTERFERENCE SUPPRESSION

6. Ziehe A, Muller KR. TDSEP—an efficient algorithm for blind separation using time structure. International
Conference on Artificial Neural Networks (ICANN ’98), Skovde, Sweden, September 1998; 675–680.
7. Hyvarinen A. Fast and robust fixed-point algorithms for independent component analysis. IEEE Transactions on
Neural Networks 1999; 10(3):626–634.
8. MacKay DJC. Bayesian non-linear modeling for the energy prediction competition. ASHRAE Transactions 1994;
100(2):1053–1062.
9. Sahani M, Linden J. Evidence optimization techniques for estimating stimulus–response functions. Advancesin Neural Information Processing Systems (NIPS ’02), vol. 15. MIT Press: Cambridge, MA, December 2002;
301–308.
10. MacKay DJC. Bayesian interpolation. Neural Computation 1992; 4(3):415–447.
11. Cover TM, Thomas JA. Elements of Information Theory. Wiley: New York, 1991.

Copyright q

2007 John Wiley & Sons, Ltd.

Statist. Med. (in press)
DOI: 10.1002/sim