Sampling-based Causal Inference in Cue Combination
and its Neural Implementation
Zhaofei Yua , Feng Chena,b , Jianwu Donga , Qionghai Daia

arXiv:1509.00998v1 [cs.NE] 3 Sep 2015

a

Center for Brain-Inspired Computing Research,
Department of Automation, Tsinghua University, Beijing 100084, China
b
Beijing Key Laboratory of Security in Big Data Processing and Application,
Beijing 100084, China

Abstract
Causal inference in cue combination is to decide whether the cues have a
single cause or multiple causes. Although the Bayesian causal inference model
explains the problem of causal inference in cue combination successfully, how
causal inference in cue combination could be implemented by neural circuits,
is unclear. The existing method based on calculating log posterior ratio with
variable elimination has the problem of being unrealistic and task-specific.
In this paper, we take advantages of the special structure of the Bayesian
causal inference model and propose a hierarchical inference algorithm based
on importance sampling. A simple neural circuit is designed to implement
the proposed inference algorithm. Theoretical analyses and experimental
results demonstrate that our algorithm converges to the accurate value as
the sample size goes to infinite. Moreover, the neural circuit we design can
be easily generalized to implement inference for other problems, such as the
multi-stimuli cause inference and the same-different judgment.
Keywords: Causal inference, importance sampling, cue combination,
neural circuit
1. Introduction
Human brain receives cues from multiple sensory modalities and integrates them in an optimal way [1]. The cues from the outside world are noisy
observations of stimuli reflecting uncertainty. It has been demonstrated that,
if all cues have the same cause, the optimal process of cue combination is a

Preprint submitted to arxiv

September 4, 2015

process of Bayesian inference [2, 3, 4, 5]. However, the truth is that, we receive information from various sources simultaneously in our daily life, which
means the cues may come from different causes. How to decide whether a
single cause or multiple causes is responsible for the cues, known as causal
inference in cue combination, is an important problem. This problem is the
precondition of cue combination and is quite common in our daily life [6, 7].
For example, at a cocktail party, we need to decide whether the face and
voice belong to the person who calls our name [8]. Recently, the problem of
causal inference in cue combination is partially answered by Kording et al.[9]
and Sato et al.[10], who propose the Bayesian causal inference model. Their
causal inference model successfully explains the problem of causal inference
in cue combination. Yet, how causal inference in cue combination could be
implemented by neural circuit, is unclear. Solving this problem benefits not
only theoretical researches but also practical applications. On the one hand,
causal inference is the basis for cue combination. On the other hand, if the
causal inference could be implemented by neural circuits, the neural circuits
could be used to perform causal inference in cue combination for robots.
Over the past decade, several methods with different probability codes
have been proposed to perform probability inference with neural circuits.
Rao [11, 12, 13] establishes the relationship between the dynamic equation
of neural circuits and the inference of probabilistic graphical models. He
proves that the process that the firing rate of neurons in the recurrent neural
circuit varies with respect to time is a process of posterior probabilities inference in a hidden Markov model, under the condition that the firing rate is
proportional to the log of posterior probabilities. Ott and Stoop [14] build the
relationship between the dynamical equation of continuous Hopfield network
and belief propagation on a binary Markov random field. Sampling is another
commonly accepted way to perform inference by neural circuits. Based on
Monte Carlo sampling, Huang and Rao [15] build a spiking network model
to perform approximate inference for any hidden Markov model. Maass et
al.[16, 17, 18] propose that stochastic networks of spiking neurons could implement inference for graphical models by Markov chain Monte Carlo. Shi
and Griffiths [19] apply importance sampling to perform inference of chain
Bayesian model and design neural circuits to implement it. Another important framework is Probabilistic population coding (PPC), the core idea of
which is that the neurons are encoders of distributions, instead of the values
of variables [20, 21, 22]. Ma et al.[20] present that the inference of cue integration can be conducted simply by linear combinations of each population
2

activity with PPC. The method is exploited thereafter by Beck et al.[23] to
realize the Bayesian decision making and the inference of marginalization
[24].
To the best of our knowledge, the only work implementing causal inference
in cue combination with neural circuits is proposed by Ma et al. [25] in 2013.
They calculate the ratio of the posterior probabilities of both situations (a
single cause or multiple causes) with variable elimination and then design a
neural circuit to implement it. This method suffers from three shortcomings.
Firstly, the circuits they design are task-specific and only work on two stimuli.
If we want to implement multi-stimuli causal inference [26] with the same
method, the circuit will be completely different. What’s more, the required
number of operations increases faster than linear with respect to the number
of stimuli, which makes the neural circuit unrealistic [25]. Secondly, it is hard
to generalize the circuit to implement a similar task called same-different
judgment [27]. Thirdly, since how to implement logarithmic operations with
neurons remains unknown, approximations are taken in their neural circuit
so that they could only get near-optimal results.
In this paper, different from calculating the posterior ratio with variable
elimination in [25], we propose a hierarchical inference algorithm based on
importance sampling, which takes advantages of the special structure of the
causal inference model. A neural circuit with hierarchical structure is then
designed corresponding to the bottom-up inference process. The proposed
method has three advantages. Firstly, the neural circuit is simple and it is
easy to be realized by PPC and some simple plausible neural operations.
Secondly, it is easy to generalize this neural circuit to implement inference
for other problems, such as the multi-stimuli cause inference and the samedifferent judgment. Thirdly, a theoretical proof is given that the samplingbased method converges to the accurate value with probability one as sample
size tends to infinity.
The rest of this paper is organized as follows. Section 2 briefly reviews
the causal inference in cue combination. In section 3 we present a samplingbased inference algorithm and design the corresponding neural circuit. The
experimental results are shown in section 4. We generalize our method to
solve other two problems in section 5 and make a conclusion in section 6.

3

𝐶𝐶

𝐶𝐶 = 1
𝑋𝑋1

𝑆𝑆

𝑆𝑆1

𝑋𝑋2

𝑋𝑋1

𝐶𝐶 = 2
𝑆𝑆2

𝑋𝑋2

Fig. 1. The causal inference model in cue combination.

𝐶𝐶

2. The Causal Inference Model In Cue Combination
The problem of causal inference
𝑆𝑆1 in cue𝑆𝑆2combination is to infer whether
cues come from a single or multiple causes. Kording et al. [9] and Sato et al.
[10] propose a causal inference model of cue combination respectively, which
could explain physiological and psychological experiments successfully. Here,
𝑋𝑋2 considered here only include
𝑋𝑋1the stimuli
we briefly review this model and
visual and auditory ones. The multi-stimuli problem will be explained in
section 5. In Fig. 1, node C represents the common-cause variable, S, S1 , and
S2 express the stimuli. X1 and X2 are cues received by the sensory system.
The state of cause C is 1 or 2, where C = 1 means the cues have the same
cause and C = 2 means the cues have two different causes. For simplicity, we
assume that P (C = 1) is equal to P (C = 2), both of which have a probability
0.5. When C = 1 , there is a stimulus S with distribution P (S) corresponding
to the common cause, where P (S) is a Gaussian distribution with mean 0
and variance σS2 . Two measurements X1 and X2 are generated from two
Gaussian distributions with different variances σ12 and σ22 , but with the same
mean S. When C = 2, there are two different stimuli S1 and S2 , which are
drawn from the same Gaussian distribution with mean 0 and variance σS2 .
Then two measurements X1 and X2 are drawn from two different Gaussian
distributions with their means being S1 and S2 , and their variances being
σ12 and σ22 respectively. Based on the definitions above, the causal inference
problem is to decide whether C = 1 or C = 2 according to the measurements
4

zation

synaptic
weight: w1

∑

𝐶𝐶 2 ⋯ 𝐶𝐶 𝑀𝑀

1
( C ) 𝐶𝐶

40

𝐶𝐶 = 1

𝑆𝑆21 ⋯ 𝑆𝑆2𝑁𝑁

𝑆𝑆11 ⋯ 𝑆𝑆1𝑁𝑁

∑

2

1

I=
( C j 1)

w1 =

max

∑

∑

𝑋𝑋1

synaptic
weight: w2

(

P S 1i , S 2i | C j

∑ P(S , S
j =1

i

i
2

𝑋𝑋2

)

M

1

𝑆𝑆

|Cj

𝑆𝑆1

)

𝑆𝑆1

𝐶𝐶

𝐶𝐶 = 2

𝑋𝑋1

𝑆𝑆2

𝑋𝑋2

𝑆𝑆2

𝑋𝑋2

𝑋𝑋1

Fig. 2. optimal
The three-layer Bayesian network equivalent to the causal inference
n=100
model n=300
in Fig. 1.
n=1000

X1 and X2 .
3. Sampling-Based Causal Inference

-30

-20

In this section, we first convert the causal inference model to a three-layer
Bayesian network. Then we propose a sampling-based hierarchical inference
method
and30 design
the corresponding neural circuit. We demonstrate that
40
20
10
0
-10
Stimulus disparity
this circuit can be realized by PPC and simple plausible neural operations.
3.1. The three-layer Bayesian network model
In this paper, the problem is to infer the state of node C. In order to
simplify inference, we convert the causal inference model above to a threelayer Bayesian model (Fig. 2) with some appropriate prior probabilities
and conditional probabilities. In the new model, node C is the commoncause variable, which is similar to that in the causal inference model. S1
and S2 refer to two different stimuli, such as visual and auditory stimuli. The conditional probability of S1 and S2 under C is expressed
as


S12
1
√
P (S1 , S2 |C). We define P (S1 , S2 |C = 1) = δ (S1 − S2 ) 2πσ exp − 2σ2
S
S
 2 2
S1 +S2
1
and P (S1 , S2 |C = 2) = 2πσ2 exp − 2σ2 , where δ (S1 − S2 ) is the Dirac
S
S
Delta distribution. X1 and X2 are measurements from S1 and S2 respectively. The
 conditional
 probability of X1 under S1 is defined by P (X1 |S1 ) =
(X1 −S1 )2
1
√
and the conditional probability of X2 under S2 is deexp − 2σ2
2πσ
1

1

5



−S2 )2
1
fined by P (X2 |S2 ) = √2πσ
exp − (X22σ
. It is easy to verify that this
2
2
2
Bayesian network is equivalent to the causal inference model in Fig. 1.
3.2. Sampling-based inference algorithm
Several methods have been developed to perform inference for the Bayesian
model shown in Fig. 2, such as belief propagation (BP) [28] and Markov chain
Monte Carlo (MCMC) [29], all of which are able to be implemented by neural circuits [30, 31, 16, 17]. However, all the circuits have the shortcoming
of being task-specific. Specifically, the neural circuit for belief propagation
[30, 31] requires pools of spiking neurons to represent function nodes of the
factor graph. It is hard to generalize the circuit of the Bayesian model in Fig.
2 to implement multi-stimuli causal inference. Similarly, the neural circuit
[16, 17] based on MCMC should meet the neural computability condition
(NCC) and the circuit will be completely different for multi-stimuli causal
inference. In this paper, we aim to build a general-purpose neural circuit for
causal inference in cue combination. Here we utilize importance sampling to
perform inference. Importance sampling is a kind of Monte Carlo methods
in statistics, which is used to estimate the intractable integrals by random
sampling. Different from other Monte Carlo methods, importance sampling
generates samples from a simple distribution rather than the original distribution [32, 33]. Here we give a simple example.
E(f
R (X))P (X)
= f (X) P (X)dX
X
R
(X)
g (X) dX
= f (X)P
g(X)
X 

(X)
= E f (X)P
g(X)
g(X)
m
P
f
(X
i )P (Xi )
≈ m1
g(Xi )

(1)

i=1:
Xi ∼g(X)

In equation (1), the goal is to calculate the mathematical expectation
of f (X), where X follows the distribution P (X). There are cases where
we can’t sample from the original distribution P (X) of variable X directly.
Instead, we can calculate the expectation of f (X) P (X) /g (X) with X following the simple distribution g (X). Note that for the region with larger
value of g (X), the sampling points should be denser, which means the samples are more important.
6

Using importance sampling to perform inference has its neural basis. The
responses of neurons have been interpreted as Monte Carlo samples by Hoyer
and Hyvarinen [34], which means that the state of each neuron is drawn randomly from a special distribution. Shi and Griffiths [19] have used importance sampling to perform inference of chain Bayesian network. By taking in
the idea of hierarchical inference, we generalize importance sampling to our
model. What’s more, we will prove the convergence of the sampling-based
inference method. We first consider a Bayesian network with only two nodes
A and B, where A is the parent node of B. It is easy to obtain the conditional
expectation of A given B with importance sampling:
P

E (f (A) |B) =
=
≈

P

f (A) P (A|B) =

f (A)P (A|B)P (B)

A

P (B)
PA
f (A)P (B|A)P (A)
E(f (A)P (B|A))
AP
= E(P (B|A)) P (A)
P (B|A)P (A)
P (A)
AP
f (Ai )P (B|Ai )

(2)

Ai :Ai ∼P (A)

P (B|Ai )

P
Ai :Ai ∼P (A)

=

P

f (Ai )

Ai :Ai ∼P (A)

P (B|Ai )
P
P (B|Ai )
Ai :Ai ∼P (A)

In equation (2), the approximation holds as we use importance sampling
to estimate the expectation. Ai ∼ P (A) means that the sample Ai is drawn
from the distribution P (A). Here A is discrete and the sums are replaced by
integrals when A is continuous. It should be noted that the same samples of
Ai are used in both sum. When we calculate the last term in (2), we first
calculate the sum of Ai in the denominator and then calculate the sum of Ai
in the numerator. Equation (2) could be generalized to solve the inference
problem in Fig. 2:

7

P (C = 1|X1 = x1 , X2 = x2 )
= PR (C = 1|x1 , x2 )
=
P (C = 1, S1 , S2 |x1 , x2 ) dS1 , S2
S1R
,S2
=
P (C = 1|S1 , S2 ) P (S1 , S2 |x1 , x2 ) dS1 , S2
S1 ,S
R 2

=
=
≈

P (C=1|S1 ,S2 )P (x1 ,x2 |S1 ,S2 )P (S1 ,S2 )dS1 ,S2

S1 ,S2

R

P (x1 ,x2 |S1 ,S2 )P (S1 ,S2 )dS1 ,S2
S1 ,S2
E(P (C=1|S1 ,S2 )P (x1 ,x2 |S1 ,S2 ))P (S ,S )
1 2
E(P (x1 ,x2 |S1 ,S2 ))P (S ,S )
1 2
N
P

i=1
S1i ,S2i ∼P (S1 ,S2 )

=

N
P
i=1

S1i ,S2i ∼P (S1 ,S2 )

(

P(

N
P

P (x1 ,x2 |S1i ,S2i )

)

P S1i ,S2i |C=1
i
i
S1 ,S2 |C=1 +P S1i ,S2i |C=2

)

(

)

P (x1 ,x2 |S1i ,S2i )
N
P

i=1
S1i ,S2i ∼P (S1 ,S2 )

I S1i = S2i



(3)

P (x1 ,x2 |S1i ,S2i )



i=1
S1i ,S2i ∼P (S1 ,S2 )

N
P
i=1
S1i ,S2i ∼P (S1 ,S2 )

=

P C = 1|S1i , S2i

P (x1 ,x2 |S1i ,S2i )

P (x1 ,x2 |S1i ,S2i )
N
P

i=1
S1i ,S2i ∼P (S1 ,S2 )

P (x1 ,x2 |S1i ,S2i )

In equation (3), the sample S1i , S2i is drawn from P (S1 , S2 ). We abbreviate X1 = x1, X2 = x2 to x1 , x2 and this will hold in the rest of the
paper. I S1i = S2i is an indicator function, it equals to 1 only when S1i =
S2i . The last equality holds due to the definitions of P (S1 , S2 |C = 1) and
P (S1 , S2 |C = 2). Note that equation (3) also holds for P (C = 1) 6= 0.5. It
is easy to find that equation (3) remains the hierarchical structure of Bayesian
model in Fig. 2. Based on this, a neural circuit with a hierarchical structure
could be designed corresponding to the bottom-up process of inference. We
will discuss this part in detail in the following subsection.
An important index for sampling-based algorithm is its accuracy. The
following theorem elucidates that our algorithm converges to the accurate
value with probability one as the sample size goes to infinity. The proof of
theorem 1 is provided in Appendix A.
Theorem 1. The distributions P (C), P (S1 , S2 |C), P (X1 |S1 ) and P (X2 |S1 )
are defined on the Bayesian network in Fig.2. S1i , S2i ∼ P (S1 , S2 ), then for
arbitrary small number ε
8

45
40
35

Firing rate

30
25
20
15
10
5
0
-40

-30

-20

10
-10
0
Stimulus S (orientation angle in degrees)

20

30

40

Fig. 3. The tuning curve of a neuron in the primary visual cortex (V1) and
its neural variability.


lim P 

N →∞

N
P
i=1


P C = 1|S1i , S2i



P (x1 ,x2 |S1i ,S2i )
N
P
i=1

P (x1 ,x2 |S1i ,S2i )

− P (C = 1|x1 , x2 ) < ε = 1
(4)

3.3. Implementation with Neural Circuits
In this section, we design a neural circuit to implement sampling-based
causal inference in cue combination. According to recent studies, one of the
most accepted neural circuits to implement probability inference is based on
PPC and some plausible neural operations. Related researches include that
of Ma et al.[20], presenting that the inference of cue integration can be conducted by PPC and a plausible neural operation-linear combinations. Beck et
al.[24] implement the inference of marginalization with PPC, quadratic nonlinearity and divisive normalization. Our method adopts the same structure
as that of Ma and Beck, and the neural circuit is designed based on PPC and
three types of plausible neural operations, including multiplication, normalization and linear combinations. Here we first give a simple explanation of
PPC. PPC takes advantages of the variability in neuronal responses and considers a population of neurons as the encoders of probability distributions,
9

Mean

0.01
0
10

2

10
Sample size

3

10

4

0.25
Our's
Ma et al.

Error rate

0.2

0.15

0.1

0.05

0
10

2

10
Sample size

3

10

4

=
w2 I ( S1i ≠ S 2i )

=
w1 I=
( S1i S2i )

Coding of the result
max

2

1

⋯

∑
synaptic
weight: w1

∑
synaptic
weight: w2

normalization
1
PPC: r1

r12

S11 , S 21

r1m
S12 , S 22

⋯

S1 , S 2 ~ P ( S1 , S 2 )
i

S1N , S2N

i

Input : X 1 , X 2

Fig. 4. The neural circuit of sampling-based method for causal inference in
cue combination.

10

rather than the values of variables. Specifically, for independent Poisson
spiking neurons, the distribution of the responses r = {r1 , r2 , ..., rN } to the
ri
Q −fi (s)
input stimulus S is P (r|S) = e rif!i (s) , where fi (s) is the tuning curve
i

of the neuron i. The tuning curve is a function of S, which represents the
average firing rate to stimulus S over trials (In theory, an infinite number
of trials). Fig. 3 shows an example, the blue curve is the Gaussian-like
tuning curve of a neuron in the primary visual cortex (V1). This neuron
is sensitive to the moving direction of the stimulus. The red circles are the
firing rates with respect to different moving direction of the stimulus in a
trial. The red circles are not always on the blue curve because the tuning
curve is the average firing rate and neural response has variability. The neural circuits for equation (3) are shown in Fig. 4. We suppose that there
are Poisson spiking neurons S11 S21 , S12 S22 ,..., S1N S2N with their states sampling
from P (S1 , S2 ). The tuning curve of the neuron S1i S2i is proportional to
P (X1 , X2 |S1i , S2i ). These assumptions are reasonable as physiologically studies [35, 36, 37] have demonstrated that the quantity of neurons in human
brain follows some prior distributions. The Poisson spiking neurons S11 S21 ,
S12 S22 ,..., S1N S2N are used to code the input stimuli X1 , X2 and the output firing rates are r11 , r12 ,..., r1N respectively. The firing rates are then normalized.
Note that the normalization operation can be realized by inhibitory
neurons
P
[38]. If we use R to express the total firing rate, where R =
ri , then we
i


P
i
i
i
i
P (X1 , X2 |S1 , S2 ) , which
can get E (ri /R|R = n) = P (X1 , X2 |S1 , S2 ) /
i

is proved in [19]. The equation above means the expectation of the normalized firing ratefor Poisson neurons equals to normalized probability
N


P
P X1 , X2 |S1i , S2i /
P X1 , X2 |S1i , S2i . These neural activities are then
i=1

fed into the third layer with synaptic weights w1 and w2 , where w1 =
I (S1i = S2i ) and w2 = I (S1i 6= S2i ). In the fourth layer, a max operation
is taken to decide whether the cause is 1 or 2. Note that, the precondition
of the inference is that we have known the prior probability and conditional
probability. We suppose that the prior probability is presented by the distribution of Poisson spiking neurons, which means the states of the Poisson
spiking neurons follow the prior distribution. We also assume that the tuning curves are proportional to conditional probability. With the benefit of
sampling-based inference, massive number of neurons could sample in paral-

11

normalized firing rate

Rate

1.5

1

0.5

0
0

100

200

300

400
500
600
Index of the neuron

700

800

900

1000

100

200

300

400
500
600
Index of the neuron

700

800

900

1000

60
50

Firing rate

40
30
20
10
0
0

Fig. 5. The firing rate of the Poisson spiking neurons in a trial.
lel and calculate without iteration. This means that the neural circuit could
trade space for time thus the inference would be quite rapid.
4. Experiments
In this section, we demonstrate the merits of the proposed method with
experiments. Samples of X1 and X2 are generated according to the prior
probabilities P (C) and conditional probabilities P (S1 , S2 |C), P (X1 |S1 ), P (X2 |S2 )
in three steps. First, we generate samples of variables C with equality probabilities 0.5 for each state. Then for each sample C i , if C i = 1, S1i will be
generated from a Gaussian distribution with mean 0 and variances σS2 , and
S2i = S1i . If C i = 2, S1i and S2i will be drawn from the same Gaussian distribution whose mean is 0 and variances is σS2 . At last X1i and X2i are generated
from two different Gaussian distributions, whose means are S1i and S2i and
variances are σ12 and σ22 respectively.
4.1. Experiment 1: Simulating the Poisson spiking neurons and their firing
rates
Here we simulate the behaviors of the Poisson spiking neurons in Fig.
4. We first generate the input X1i and X2i randomly with the method
proposed above. The parameters are specific to σS = 4, σ1 = σ2 = 6.
Then we generate 1000 Poisson spiking neurons S11 S21 , S12 S22 ,..., S11000 S21000
by sampling from P (S1 , S2 ). The tuning curve of the neuron S1i S2i is set to
10000 × P (X1 , X2 |S1i , S2i ). Fig. 5 represents the firing rate r11 , r12 , ..., r11000
of the 1000 Poisson spiking neurons in a trial (the firing rate could vary
12

0
0

200

300

400

500

600

700

800

900

1000

20

Firing rate

100

15

10

5

2

x 10

0

-3

2

normalized probability
normalized firing rate

Rate

1.5

1

0.5

0
0

100

200

300

400
500
600
Index of the neuron

700

800

900

1000

60
Fig. 6. The normalized
firing rate compared with the normalized distribution
in a trial.
50
40

Firing rate

in different trials due to neural variability). Note that here we only show
30
the results of neurons
whose indexes range from 1 to 1000 with a uniform
spacing 30. In Fig.
6,
the circles on the blue curve represents the nor20
1000
P i
i i
i
r1 . Similarly, we only
malized firing rate
10 for neuron S1 S2 , which is r1 /
i=1

show the results of00 neurons
whose
indexes
range
from
1 to10001000 with a uni100
200
300
400
500
600
700
800
900
of the neuron
form spacing 30. The plus on theIndexred
curve is the normalized probability
1000
P
P (X1 , X2 |S1i , S2i ). We can see that the normalized firP (X1 , X2 |S1i , S2i ) /
i=1

ing rate is close to the normalized probability.
4.2. Experiment 2: Testing on the convergence and accuracy of our method
Here we simulate and present the behaviors of the neurons in the last two
layers and show the convergence and accuracy of our method. We first generate 1000 inputs of X1i and X2i randomly with the method proposed above.
For each inputs X1i and X2i , σS , σ1 and σ2 are drawn randomly from a uniform
distribution on [3 7]. Then we calculate P (C = 1|X1i , X2i ) with the samplingbased method and express the result as P sample (C = 1|X1i , X2i ). Meanwhile,
the truth of P (C = 1|X1i , X2i ) is expressed as P truth (C = 1|X1i , X2i ), which is
calculated with the elimination method in [16]. The error of samples X1i and
X2i is defined by δ = P sample (C = 1|X1i , X2i ) − P truth (C = 1|X1i , X2i ) . This
index expresses the gap between the sampling value and optimal value of posterior probability. The mean error is calculated from 1000 different inputs.
The error rate represents the proportion of false results in the 1000 different
13

4

6
Stimulus S

8

Mean error of the posterior probability

0.07
0.06
0.05
0.04
0.03
0.02
0.01

0.07 0

Mean error of the posterior probability

10

2

10
Sample size

3

10

4

0.06

0.25

Fig. 7. Mean error of the posterior probability varies Our's
with
Ma
et al. sample size.
0.05

0.04 0.2

Error rate

0.03

0.15

0.02

0.1

0.01
0 0.05

10

2

10
Sample size

0.25 0
10

2

3

10

10
Sample size

3

Our's
Ma et al.

4

10

4

Error rate

0.2

0.15

=
w1 I=
( S1i S2i )

0.1

=
w2 I ( S1i ≠ S 2i )

Coding of the result

0.05

max

0
10

2

10
Sample size

3

10

4

2
Fig. 8. Error1 rate of causal inference.
i
i
=
w2 I ( S∑
⋯
1 ≠ S2 )

=
w1 I=
( S1i ∑S2i )

Coding of the result synaptic
weight: w2
14max

synaptic
weight: w1

normalization
1
PPC: r1 1

S11 , S 21
i

i

∑

r12

m
2r1

S12 , S 22

⋯⋯

∑

S1N , S2N

inputs when we infer the cause. We calculate the mean error and the error
rate with different sample sizes and repeat the experiments 10 times. The
mean error of the posterior probability varying with sample size is shown in
Fig. 7. We find that the mean error decreases as the sample size increases
and converges to zero when the sample size tends to infinity. This result
demonstrates that if we have enough samples, or we have enough neurons,
our algorithm will get the optimal value. Fig. 8 plots the error rate obtained
from our method and that from the method of Ma et al [16]. The method of
Ma et al. is not related to sample size while ours could get stretchable results
(different accuracies) with different sample sizes. Obviously, our method is
superior to Ma’s when the sample size is larger than 500 (500 neurons for
each variable). We can see that in order to keep the error rate under 0.05
for two-stimuli causal inference, we need at least 1000 neurons to represent
each variable. This means N = 1000 in equation (3).
4.3. Experiment 3: Testing on the applicability of the method with different
parameters.
Experiment 1 and 2 indicate that our inference method could get the
optimal solution given enough samples. However, the parameters σ1 , σ2
and σS are drawn randomly for each trial. In experiment 3, we will make a
concrete analysis of the applicability of our method with different parameters.
In this experiment, σ1 , σ2 and σS can vary from 1 to 8. We test the error
rate for different σ1 , σ2 and σS with the sample size being 1000 and show
the result in Fig. 9. In each sub-figure, σS is set to a fixed value while both
σ2 and σS vary from 1 to 8. We can see that the error rate is less than 0.1
for most of the parameters. However, when σS is in close proximity to zero,
the error rate turns out to be very high. This could be explained as follows.
Since σS is close to zero, the difference between S1 and S2 remains very small
no matter C = 1 or 2. Then the posterior probability of a common cause and
two different causes both will be quite near to 0.5. Due to this, a very small
error could lead to incorrect inference results, making the error rate very
high. Nevertheless, if there are adequate samples, the error rate could be
arbitrarily small, which means our method is robust to different parameters.
4.4. Experiment 4: Testing on the probability of reporting a common cause
with respect to stimulus disparity.
Stimulus disparity refers to the space difference between different stimuli,
which is defined by S2 − S1 , where S1 and S2 are two different stimuli. In15

σS =1

8

σS = 2

8

0.25

0.16
7

7

0.14

0.2

6

6

0.15

4

Sigma 1

Sigma 1

0.12
5

5
0.1
4
0.08

0.1

3

3
0.06
2

2

0.04

0.05
1
1

2

3

4

5

6

7

1
1

8

2

3

4

Sigma 2

7
6

0.14

8

0.12

7

7

8

0.08
4

σS = 4
0.1

6

0.1

5

0.08

Sigma 1

Sigma 1

6

Sigma 2

σS = 3

8

5

5
0.06

4

0.06
3

3
0.04

2
1
1

0.02
2

3

4

5

6

7

1
1

8

Sigma 2

0.02
2

3

4

5

6

7

8

Sigma 2

σS = 5

8

0.04

2

σS = 6

8
0.1

7

0.08

7

0.07

0.08

5

0.07
0.06

4

6

Sigma 1

Sigma 1

0.09
6

0.06

5
0.05
4
0.04

0.05
3

0.03

0.03

2
1
1

3

0.04

2
0.02

0.02
2

3

4

5

6

7

1
1

8

2

3

4

Sigma 2

6

7

8

Sigma 2

σS = 7

8

5

σS = 8

8
0.07

7

7
0.06

6

0.06

6

0.05

4

0.04

3

Sigma 1

Sigma 1

0.05
5

5
0.04

4
3

0.03

0.03
2

2
0.02

0.02
1
1

2

3

4

5
Sigma 2

6

7

8

16

1
1

2

3

4

5

6

7

Sigma 2

Fig. 9. Error rate for different parameters with sample size equals 1000.

8

0.8

optimal
n=100
n=300
n=1000

Proportion reports C=1

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-30

-20

-10
0
10
Stimulus disparity

20

30

Fig. 10. Comparison of proportion reports C = 1 with respect to stimulus
disparity between 8accurate equation and our sampling-based method.0.22
0.2

7

Sigma S

tuitively, it is more likely that there is a common cause if stimulus disparity
0.18
6
is small, while there
are two different causes if stimulus disparity is0.16large.
In this experiment,
0.14
5 200000 samples are generated with parameters σ
1 = 3,
i
i
i
i
σS = σ2 = 10. For each S1 and S2 , stimulus disparity is defined by S0.12
2 − S1 .
4
The state of variable C i is inferred by optimal equation [25] and our method
0.1
respectively. Then3 for the samples with the same stimulus disparity, 0.08
we calculate the proportion
of
reporting
a
common
cause.
Fig.
10
shows
the
result,
0.06
2
the red curve is obtained by optimal equation. The black, green and
0.04 blue
1
curves are calculated
by 2our method
with
sample
size
being
100,8 300 and
1
3
4
5
6
7
Sigma 2
1000 respectively. The result shows that as the sample size becomes larger,
the sampling-based curve tends to be closer to the accuracy curve. When
sample size equals 1000, the sampling-based curve is almost the same as the
accuracy curve. This result indicates the accuracy of our method.
0.2

Error rate for three-stimuli causal inference
Error rate for ten-stimuli causal inference

0.18

5. Generalization
0.16
0.14

Error rate

In this section, we generalize our sampling-based method to implement
0.12
inference for other two important problems: the multi-stimuli causal infer0.1
ence in cue combination
and the same-different judgment.
0.08
0.06
0.04
0.02
0
10

2

17
10
Sample size

3

Error rat

0.2
0.1
0
10
8

10
8

6

6

4

𝑆𝑆1

𝑋𝑋1

4
2

Sigma 2

𝐶𝐶

𝑆𝑆2 ⋯

2

Sigma S

𝑆𝑆𝑛𝑛

𝑋𝑋2 ⋯ 𝑋𝑋𝑛𝑛

𝑆𝑆1

𝐶𝐶

𝑋𝑋1

Fig. 11. The Bayesian model for the multi-stimuli causal inference.
5.1. Multi-stimuli causal inference
In the above experiments, situations where there are only two stimuli
were taken into consideration while in our daily life cues may come from
multiple sensory modalities, such as visual, auditory, and tactile. Despite
the fact that Bayesian model could also explain the causal inference problem
with multi-stimuli [26], how to implement inference with neural circuits is
unclear. With our sampling-based method, it is easy to generalize the neural
circuits to implement inference for multi-stimuli.
The steps are similar to that in section 3. First we convert the problem
to a Bayesian network (shown in Fig. 11). The prior probabilities and
conditional probabilities are defined according to the causal inference model,
where different states of C reflect different situations of the cause. The causal
inference problem is then converted to the inference of posterior probability,
which could be calculated by importance sampling:

18

𝜇𝜇
𝑆𝑆2 ⋯

𝑆𝑆𝑛𝑛

𝑋𝑋2 ⋯ 𝑋𝑋𝑛𝑛

P (C = 1|X1 = x1 , X2 = x2 , ..., Xn = xn )
= P (C
R = 1|x1 , x2 , ..., xn )
=
P (C = 1, S1 , S2 , ..., Sn |x1 , x2 , ..., xn ) dS1 , S2 , ..., Sn
S1 ,S2R,...,Sn
=
P (C = 1|S1 , S2 , ..., Sn ) P (S1 , S2 , ..., Sn |x1 , x2 , ..., xn ) dS1 , S2 , ..., Sn
S1 ,S2R,...,Sn

=
=
≈

P (C=1|S1 ,S2 ,...,Sn )P (x1 ,x2 ,...,xn |S1 ,S2 ,...,Sn )P (S1 ,S2 ,...,Sn )dS1 ,S2 ,...,Sn

S1 ,S2 ,...,Sn

R

P (x1 ,x2 ,...,xn |S1 ,S2 ,...,Sn )P (S1 ,S2 ,...,Sn )dS1 ,S2 ,...,Sn
S1 ,S2 ,...,Sn
E(P (C=1|S1 ,S2 ,...,Sn )P (x1 ,x2 ,...,xn |S1 ,S2 ,...,Sn ))P (S ,S ,...,Sn )
1 2
E(P (x1 ,x2 ,...,xn |S1 ,S2 ,...,Sn ))P (S ,S ,...,Sn )
1 2
N
P

i=1
i ∼P (S ,...,S )
S1i ,...,Sn
n
1

=

N
P
i=1
i ∼P (S ,...,S )
S1i ,...,Sn
n
1

=

N
P
i=1

P C = 1|S1i , S2i , ..., Sni

i
P (x1 ,...,xn |S1i ,S2i ,...,Sn
)



N
P

i=1
i ∼P (S ,...,S )
S1i ,...,Sn
n
1

i
P (x1 ,...,xn |S1i ,S2i ,...,Sn
)

i
P (x1 ,...,xn |S1i ,...,Sn
)

i |C=1
P (S1i ,...,Sn
)
i |C=1 +P S i ,...,S i |C=2
P (S1i ,...,Sn
) ( 1 n
)

N
P

i=1
i ∼P (S ,...,S )
S1i ,...,Sn
n
1

I S1i = S2i = ... = Sni

i ∼P (S ,...,S )
S1i ,...,Sn
n
1



i
P (x1 ,...,xn |S1i ,...,Sn
)

i
P (x1 ,...,xn |S1i ,S2i ,...,Sn
)
N
P

i=1
i ∼P (S ,...,S )
S1i ,...,Sn
n
1

i
P (x1 ,...,xn |S1i ,S2i ,...,Sn
)

(5)
Equation (5) is similar to (3) except that the stimuli here are S1 , S2 , ..., Sn
rather than S1 , S2 . According to (5), the neural circuit of multi-stimuli causal
inference is similar to that of two-stimuli causal inference except three differences. Firstly, the states of the Poission spiking neurons sample from
P (S1 , S2 , ..., Sn ), rather than P (S1 , S2 ). Secondly, the tuning curve
of the

i
i
i
neuron marked as i is proportional to P x1 , x2, ..., xn |S1 , S2 , ..., Sn . Thirdly,

the synaptic weights are I S1i = S2i = ... = Sni instead of I S1i = S2i .
We also test on the convergence and accuracy of our method for multistimuli causal inference to show its feasibility. Three-stimuli casual inference is tested because the situation is quite common in our daily life, such
as the integration of visual, auditory, and tactile input. We also test tenstimuli to show that our method applies to higher dimensions. We assume
that P (C = 1) is equal to P (C = 2), both of which have a probability 0.5.
Note that C = 1 means the three cues have the same cause and C = 2
means the cues have diverse causes. We define P (S1 , S2 , ..., St |C = 1) =

19

0.06

2

0.04
1
1

2

3

4

5

6

7

8

Sigma 2

0.2

Error rate of three-stimuli causal inference
Error rate of ten-stimuli causal inference

Error rate

0.15

0.1

0.05

0
10

2

10

3

Sample size

Fig. 12. Error rate of the multi-stimuli causal inference.
Q

t
S2
exp − 2σ12
δ (S1 − Sj ) and P (S1 , S2 , ..., St |C = 2) = √ 1 2 t
S
2πσS
j=1
 2 2

2
S +S2 +...+St
exp − 1 2σ
. Besides, the conditional probability of Xi under Si
2
S


−Si )2
1
exp − (Xi2σ
. The experi(i = 1, 2, ..., t) is defined by P (Xi |Si ) = √2πσ
2
i
i
mental procedure is similar to that of Experiment 2 in section 3 and the result
is shown in Fig. 12. We can find that the error rate decreases as the sample
size increases and convergences to zero when sample size tends to infinity.
We also find that we don’t need to scale up the samples when dimensions
become higher. 1000 samples (neurons) are required for each variable to keep
the error rate under 0.05. These results are in good agreement with the fact
that importance sampling does not scale up with higher dimensions.
√ 1
2πσS

5.2. Same-different judgment
When faced with multiple objects, probably the first thing our brain needs
to do is to decide whether they are the same or not. Thus the same-different
judgment could be critical in perception and cognition. A straightforward
example is object classification. Human brains are able to recognize the same
object and assign them to the semantic classes. Berg et al. [27] propose
the optimal-observer model and prove that the same-different judgment is
a process of probability inference. As illustrated in Fig. 13, variable C
represents the judgment, C = 1 means the objects are the same while C = 2
means they are different. µ is a single value parameter variable generated
20

Error rate

0.2
0.1
0
10
8

10
8

6

6

4

Sigma 2

𝑆𝑆1

𝑋𝑋1

4
2

𝐶𝐶

𝑆𝑆2 ⋯

2

Sigma S

𝑆𝑆𝑛𝑛

𝑋𝑋2 ⋯ 𝑋𝑋𝑛𝑛

𝑆𝑆1

𝑋𝑋1

𝐶𝐶

𝜇𝜇
𝑆𝑆2 ⋯

𝑆𝑆𝑛𝑛

𝑋𝑋2 ⋯ 𝑋𝑋𝑛𝑛

Fig. 13. Optimal-observer model of the same-different judgment problem.

from a uniform distribution ranging from −L to L. When the objects are
the same, Si equals to µ. When they are different, Si is drawn from a
Gaussian distribution with mean µi and variance σS2 . The distribution of
Xi is a Gaussian distribution with its mean being Si and its variance being
σi2 . Based on these definitions, the same-different judgment problem can be
converted to the posterior probability inference problem of variable C, which
could be calculated by importance sampling:

21

P (C = 1|X1 = x1 , X2 = x2 , ..., Xn = xn )
= P R(C = 1|x1 , x2 , ..., xn )
=
P (C = 1, S1 , ..., Sn |x1 , ..., xn ) dS1 , ..., Sn
S1 ,...,S
n
R
=
P (C = 1|S1 , ..., Sn ) P (S1 , ..., Sn |x1 , ..., xn ) dS1 , ..., Sn
S1 ,...,Sn

≈

N
P
i=1
i ∼P (S ,...,S )
S1i ,...,Sn
n
1

=

(

N
P

P(

)

(

i
P (x1 ,...,xn |S1i ,...,Sn
)

i
P (x1 ,...,xn |S1i ,...,Sn
)

)

i |C=1
P S1i ,...,Sn
i
i
i |C=2
S1 ,...,Sn |C=1 +P S1i ,...,Sn

)

N
P

i=1
i ∼P (S ,...,S )
S1i ,...,Sn
n
1
R

N
P

µ

R
i=1
i ∼P (S ,...,S )
S1i ,...,Sn
n
1

i
P (x1 ,...,xn |S1i ,...,Sn
)



i=1
i ∼P (S ,...,S )
S1i ,...,Sn
n
1

N
P
i=1
i ∼P (S ,...,S )
S1i ,...,Sn
n
1

=

P C = 1|S1i , ..., Sni

µ

i
P (x1 ,...,xn |S1i ,...,Sn
)

i ,µ|C=1 dµ
P (S1i ,...,Sn
)

R
i ,µ|C=1 dµ+ P S i ,...,S i ,µ|C=2 dµ
P (S1i ,...,Sn
( 1 n
)
)
µ


i
P (x1 ,...,xn |S1i ,...,Sn
)
N
P

i=1
i ∼P (S ,...,S )
S1i ,...,Sn
n
1

=




i 
P (x1 ,...,xn |S1i ,...,Sn
)
R

N
P

µ

R
i=1
i ∼P (S ,...,S )
S1i ,...,Sn
n
1

µ

i |C=1,µ dµ
P (S1i ,...,Sn
)

R
i |C=1,µ dµ+ P S i ,...,S i |C=2,µ dµ
P (S1i ,...,Sn
( 1 n
)
)
µ


i
P (x1 ,...,xn |S1i ,...,Sn
)
N
P

=

i=1
i ∼P (S ,...,S )
S1i ,...,Sn
n
1
N
P
i=1
i ∼P (S ,...,S )
S1i ,...,Sn
n
1




i
i
P (x1 ,...,xn |S1 ,...,Sn ) 

I −L ≤ S1i = S2i = ... = Sni ≤ L



i
P (x1 ,...,xn |S1i ,...,Sn
)
N
P

i=1
i ∼P (S ,...,S )
S1i ,...,Sn
n
1

i
P (x1 ,...,xn |S1i ,...,Sn
)

(6)
In equation (6), I −L ≤ S1 = S2 = ... = Sn ≤ L is a indicative function.
It equals to 1 only when all the stimuli are the same and between −L and
L. Equation (6) differs from equation (5) in the indicative function. Due
i

i

i

22



0.4

Error rate of same-different judgment with three objects
Error rate of same-different judgment with ten objects

0.35
0.3

Error rate

0.25
0.2
0.15
0.1
0.05
0
10

2

3

10
Sample size

10

4

Fig. 14. Error rate of the same-different judgment
to this, the neural circuit for the same-different judgment are similar to
that for multi-stimuli causal inference except one difference. That is, the
synaptic weights of the neural circuit for the same-different judgment
are

i
i
i
i
i
i
I −L ≤ S1 = S2 = ... = Sn ≤ L , instead of I S1 = S2 = ... = Sn .
We test the accuracy and convergence of our method for the same-different
judgment problem. In this experiment, we compute the error rate of samedifferent judgment with three objects and ten objects respectively. Samples
of X1 , X2 ...., and Xt are generated according to the similar method in the
section of multi-stimuli causal inference. Note that µi is generated from the
uniform distribution ranging from −10 to 10. σS and σi are drawn randomly
from a uniform distribution on [1 3]. The inference result with our method
is present in Fig. 14, which is similar to that of Fig. 13. We can find that
the error rate of same-different judgment with three objects and ten objects
both decrease as the sample size increases and convergence to zero when
sample size tends to infinity. Besides, the samples needed don’t scale up
when dimensions become higher. 5000 samples (neurons) are required for
each variable to keep the error rate under 0.05.
6. Conclusion
In this paper, we propose an inference algorithm for causal inference in
cue combination based on importance sampling and design a corresponding
neural circuit to implement this inference algorithm. The neural circuit is
23

plausible as it is based on PPC and three types of plausible neural operations.
Theoretical analysis and experimental results show that our algorithm can
converge to the accurate value as sample size goes to infinite. It is worth noting that our method provides a general solution to the other two important
problems, namely the multi-stimuli causal inference and the same-different
judgment.
Different from Markov chain Monte Carlo [16, 17], which represents the
distribution with variability of a neuron over time, our method utilizes the
variability over neurons to represent the distribution. This means that massive number of neurons could sample in parallel and calculate without iteration, thus the inference would be quite rapid.
Despite the plausible neural implementation of inference, the question
of how to learn the prior probabilities and conditional probabilities with
learning rules found in biological studies requires considerable future work.
Besides, learning and inference should be implemented by the same neural circuits. Some recent works have provided reference experiences for implementing learning. For example, Maass et al. prove that Spike-Timing-Dependent
Plasticity (STDP) is able to approximate a parameter estimation algorithm–
expectation maximization (EM) algorithm [39, 40, 41]. This principle may
be used to solve the learning problem in our paper.
Proof of Theorem 1
Lemma 1. Supposing that random variables X 1 , X 2 , ..., X n are pairwise
independent and X i ∼ P (X). Similarly, Y 1 , Y 2 , ..., Y n are pairwise independent and Y j ∼ P (Y ). Besides, E (X) = µ1 , E (Y ) = µ2 , µ1 , µ2 6= 0,
2
V ar (X) = σ12 and
 for arbitrary small number ε, we can
 V ar (Y ) = σ2 . Then
N
P

conclude that P 

i=1
N
P

Xi

−
Yj

µ1
µ2

< ε > 1 −

16σ12
N µ22 ε2

−

16µ21 σ22
.
N µ42 ε2

j=1

Proof: As random variables X 1 , X 2 , ..., X n are pairwise independent and
X i ∼ P (X), Y 1 , Y 2 , ..., Y n are pairwise independent and Y j ∼ P (Y ), We
can get
 N

N
P i
P
1
E N
X = N1
E (X i ) = µ1 ,
i=1

E

1
N

N
P
j=1

i=1

!
Yj

=

1
N

24

N
P
j=1

E (Y j ) = µ2 ,


V ar

1
N

1
N

V ar

N
P

Xi



1
N2

=

i=1
N
P

!
Yj

1
N2

=

j=1

N
P

V ar (X i ) =

σ12
,
N

V ar (Y j ) =

σ22
.
N

i=1
N
P
j=1

For arbitrary small number ε1 and ε2 = µµ12 ε1 , application of the Chebyshev’s
Inequality yields the inequality:


N
P
σ2
1
i
P N
X − µ1 < ε1 ≥ 1 − N ε12 ,
1

i=1

1
N

P

N
P

Because

1
N

µ 2 − ε2 <

N
P

X i and
i=1

1
P  µµ12 −ε
+ε2

1
N

Y j − µ2 < ε2

≥1−

j=1

which is equivalent to

P µ 1 − ε1 <

P

!

1
N

1
N

N
P

σ22
,
N ε22



i

X < µ1 + ε1

≥1−

σ12
,
N ε21

≥1−

σ22
.
N ε22

i=1
N
P

!
Y j < µ2 + ε2

j=1

N
P

Y j are independent, so it is trivial to show that
j=1

N
P



Xi
σ12
σ22
µ1 +ε1 
i=1
< µ2 −ε2 ≥ 1 − N ε2
< P
1 − N ε2 .
N
1

Yj

2

j=1

An application of Taylor’s formula yields:


ε1
ε
2 µ1
µ1 −ε1
µ1 1− µ1
µ1
= µ2 1+ ε1 = µ2 1 − 1+ ε11
µ2 +ε2
µ1
µ1


 2  3
µ1
µ1 2ε1
ε1
ε1
ε1
= µ2 − µ2 µ1 1 − µ1 + µ1 − µ1 + ...
>

µ1
µ2

−

2ε1
µ2

>

µ1
µ2

−

4ε1
,
µ2

25

ε

µ1 +ε1
µ2 −ε2

1
µ1 1+ µ1
ε1
µ2 1− µ

=



ε
2 µ1
1 + 1− ε11
µ1


 2  3
ε1
ε1
ε1
1 + µ1 + µ1 + µ1 + ...

µ1
µ2

=

1

=

µ1
µ2

+

µ1 2ε1
µ2 µ1

<

µ1
µ2

+

µ1 4ε1
µ2 µ1

=

µ1
µ2

+

4ε1
.
µ2

The equations above indicate that


N
P

P  µµ12 −

4ε1
µ2

<



Xi

i=1
N
P

<
Yj

µ1
µ2

+

4ε1 
µ2



≥ 1−

σ12
N ε21


1−

σ22
N ε22



.

j=1

Next, we rewrite the equation above as
 N

Xi

P

P

−

i=1
N
P

Yj

µ1
µ2

<

4ε1 
µ2

≥1−

σ12
N ε21

−

µ21 σ22
N µ22 ε21

+

µ21 σ12 σ22
.
N 2 µ22 ε41

j=1

Then for arbitrary small number ε, we have
 N

P

P

i=1
N
P

Xi

−
Yj

µ1
µ2

< ε ≥ 1 −

16σ12
N µ22 ε2

−

16µ21 σ22
N µ42 ε2

+

256µ21 σ12 σ22
.
N 2 µ62 ε4

j=1

We conclude that for arbitrary small number ε,
 N

P

P

i=1
N
P

Xi

−
Yj

µ1
µ2

< ε > 1 −

16σ12
N µ22 ε2

−

16µ21 σ22
.
N µ42 ε2

j=1

Lemma 2. Supposing that random variables X 1 , X 2 , ..., X n are independent pairwise and X i ∼ P (X).Besides, we also know that E (X) = µ1 ,
µ1 6= 0, V ar (X) = σ12 . Then for arbitrary small number ε,


P

1
N

1
N
P

−
Xi

1
µ1

< ε ≥ 1 −

σ12
.
N µ21 ε2

i=1

Proof: The proof is similar to that of Lemma 1.
Theorem 1. The distributions P (C), P (S1 , S2 |C), P (X1 |S1 ) and P (X2 |S1 )
are defined on the Bayesian network in Fig.2. S1i , S2i ∼ P (S1 , S2 ), then for
26

arbitrary small number ε,

N

P
lim P 
P C = 1|S1i , S2i
N →∞


P(

x1 ,x2 |S1i ,S2i

N
P

i=1

i=1

P(

)

x1 ,x2 |S1i ,S2i

)

− P (C = 1|x1 , x2 ) < ε = 1

Proof: Supposing that
N
P
P (C=1|S1i ,S2i )P (x1 ,x2 |S1i ,S2i )
i=1
f1 (x1 , x2 ) =
, f2 (x1 , x2 ) =
N
P
P (x1 ,x2 |S1i ,S2i )

N
P
i=1

P (C=1|S1i ,S2i )P (x1 ,x2 |S1i ,S2i )

i=1

N
P
j=1

P (x1 ,x2 |S̃1j ,S̃2j )

where S1i , S2i ∼ P (S1 , S2 ) and S̃1j , S̃2j ∼ P (S1 , S2 ), then

 N

P
1
i
i
P x1 , x2 |S1 , S2
E N
i=1

=

1
N

=

1
N

N
P
i=1
N
P

E P x1 , x2 |S1i , S2i



P x1 , x2 |S1i , S2i P S1i , S2i dS1i , S2i

R

i=1 S i ,S i
1



2

= P (x1 , x2 ) ,
E

V ar
=

1
N2

=

1
N2

=

1
N

1
N

1
N
N
P
i=1

N
P



j=1

P x1 , x2 |S̃1j , S̃2j
i

i

P x1 , x2 |S1 , S2



!
= P (x1 , x2 ) ,




N
P


V ar P x1 , x2 |S1i , S2i
i=1
N  
2 

P
i
i 2
i
i
E P x1 , x2 |S1 , S2
− E P x1 , x2 |S1 , S2
i=1
R
P (x1 , x2 |S1 , S2 )2 P (S1 , S2 ) dS1 , S2 − P (x1 , x2 )2

S1 ,S2

= σ12 ,
V ar

1
N

N
P
j=1



P x1 , x2 |S̃1j , S̃2j



!
= σ12 .

Note that the variance is denoted as σ12 . It is easy to use Lemma 2 to show
that for arbitrary small number ε,

27

,


P


1
1
N

N
P

−

P (x1 ,x2 |S̃1j ,S̃2j )

j=1

1
P (x1 ,x2 )

< ε ≥ 1 −


P

σ12
,
N P (x1 ,x2 )2 ε2


1
1
N

N
P
i=1

−

P (x1 ,x2 |S1i ,S2i )

1
P (x1 ,x2 )

< ε ≥ 1 −

σ12
.
N P (x1 ,x2 )2 ε2

The equations above indicate that

P


1
1
N

N
P
j=1

1

−

P (x1 ,x2 |S̃1j ,S̃2j )

1
N

N
P
i=1

As |f1 (x1 , x2 ) − f2 (x1 , x2 )| ≤

P (x1 ,x2 |S1i ,S2i )


< 2ε ≥ 1 −

1
1
N

N
P
j=1

x1 ,x2 |S̃1j ,S̃2j

P(

1

−
)

σ12
N P (x1 ,x2 )2 ε2

1
N

N
P
i=1

P (|f1 (x1 , x2 ) − f2 (x1 , x2 )| < 2ε) ≥ 1 −

σ12
N P (x1 ,x2 )2 ε2

We also can get that

 N


P
1
i
i
i
i
P C = 1|S1 , S2 P x1 , x2 |S1 , S2
E N
i=1

1
N

N
P



E P C = 1|S1i , S2i P x1 , x2 |S1i , S2i
i=1 (
N


R
P
= N1
P C = 1|S1i , S2i P x1 , x2 |S1i , S2i
i=1 S i ,S i
1 2


P S1i , S2i dS1i , S2i
= P (C = 1|x1 , x2 ) P (x1 , x2 ) ,

 N


P
1
i
i
i
i
P C = 1|S1 , S2 P x1 , x2 |S1 , S2
V ar N
i=1
(
R
P (C = 1|S1 , S2 )2 P (x1 , x2 |S1 , S2 )2
= N1
=

S1 ,S2

P (S1 , S2 ) dS1 , S2 ) − P (C = 1|x1 , x2 )2 P (x1 , x2 )2
= σ22 .
28

.

, we can

P (x1 ,x2 |S1i ,S2i )

get for arbitrary small number ε,


2

2

.

1 ,x2 )P (x1 ,x2 )
=
Note that the variance is denoted as σ22 here. Since P (C=1|x
P (x1 ,x2 )
P (C = 1|x1 , x2 ), it is easy to use Lemma 1 to show that for arbitrary small
number ε,

P (|f2 (x1 , x2 ) − P (C = 1|x1 , x2 )| < ε) > 1 −

16σ22
N P (x1 ,x2 )2 ε2

−

16P (C=1|x1 ,x2 )2 σ12
,
N P (x1 ,x2 )2 ε2

We also know that


P (|f1 (x1 , x2 ) − f2 (x1 , x2 )| < 2ε) ≥ 1 −

σ12
N P (x1 ,x2 )2 ε2

2

,

then for arbitrary small number ε,
P (|f1 (x1 , x2 ) − P (C = 1|x1 , x2 )| < 3ε)
2

16P (C=1|x ,x )2 σ 2
16σ 2
σ12
− N P (x ,x2 )2 ε2 − N P (x ,x1 )22ε2 1 .
≥ 1 − N P (x ,x
2 2
) ε
1

2

1

2

1

2

This term holds as P (A ∩ B) ≥ P (A) + P (B) − 1. Then for arbitrary small
number ε, we have
P (|f1 (x1 , x2 ) − P (C = 1|x1 , x2 )| < ε)

2
9σ 2
144P (C=1|x ,x )2 σ 2
144σ 2
≥ 1 − N P (x ,x1 )2 ε2 − N P (x ,x2 )2 ε2 − N P (x ,x 1)2 ε22 1 .
1

2

1

2

1

2

When N goes to infinite, we calculate the limits for both sides and conclude
that, for arbitrary small number ε,

lim P 

N →∞

N
P
i=1



i

P C = 1|S1i , S2

P(

x1 ,x2 |S1i ,S2i

N
P
i=1

)

P (x1 ,x2 |S1i ,S2i )

− P (C = 1|x1 , x2 ) < ε = 1.

Acknowledgments
The authors would like to thank Y. Gao, F. Deng and S. Guo for helpful
discussion. This work is supported in part by the National Natural Science
Foundation of China under Grant 61271388 and Grant 61327902, and in part
by the Research Project of Tsinghua University under Grant 2012Z01011 and
Grant 20141080934.

29

References
References
[1] A. Pouget, J. M. Beck, W. J. Ma, P. E. Latham, Probabilistic brains:
knowns and unknowns, Nature Neuroscience 16 (9) (2013) 1170–1178.
[2] D. Alais, D. Burr, The ventriloquist effect results from near-optimal
bimodal integration, Current biology 14 (3) (2004) 257–262.
[3] M. O. Ernst, M. S. Banks, Humans integrate visual and haptic information in a statistically optimal fashion, Nature 415 (6870) (2002) 429–433.
[4] D. C. Knill, J. A. Saunders, Do humans optimally integrate stereo and
texture information for judgments of surface slant?, Vision research
43 (24) (2003) 2539–2558.
[5] R. J. van Beers, A. C. Sittig, J. J. van der Gon Denier, How humans
combine simultaneous proprioceptive and visual position information,
Experimental Brain Research 111 (2) (1996) 253–261.
[6] R. L. Seilheimer, A. Rosenberg, D. E. Angelaki, Models and processes
of multisensory cue combination, Current opinion in neurobiology 25
(2014) 38–46.
[7] M. Ursino, C. Cuppini, E. Magosso, Neurocomputational approaches
to modelling multisensory integration in the brain: a review, Neural
Networks 60 (2014) 141–165.
[8] C. Kayser, L. Shams, Multisensory causal inference in the brain, PLoS
biology 13 (2) (2015) e1002075.
[9] K. P. Körding, U. Beierholm, W. J. Ma, S. Quartz, J. B. Tenenbaum,
L. Shams, Causal inference in multisensory perception, PLoS one 2 (9)
(2007) e943.
[10] Y. Sato, T. Toyoizumi, K. Aihara, Bayesian inference explains perception of unity and ventriloquism aftereffect: identification of common
sources of audiovisual stimuli, Neural computation 19 (12) (2007) 3335–
3355.

30

[11] R. P. Rao, Bayesian computation in recurrent neural circuits, Neural
computation 16 (1) (2004) 1–38.
[12] R. P. Rao, Hierarchical bayesian inference in networks of spiking neurons, in: Advances in neural information processing systems, 2004, pp.
1113–1120.
[13] K. Doya, Bayesian brain: Probabilistic approaches to neural coding,
MIT press, 2007.
[14] T. Ott, R. Stoop, The neurodynamics of belief propagation on binary
markov random fields, in: Advances in neural information processing
systems, 2006, pp. 1057–1064.
[15] Y. Huang, R. P. Rao, Neurons as monte carlo samplers: Bayesian inference and learning in spiking networks, in: Advances in Neural Information Processing Systems, 2014, pp. 1943–1951.
[16] L. Buesing, J. Bill, B. Nessler, W. Maass, Neural dynamics as sampling:
a model for stochastic computation in recurrent networks of spiking
neurons, PLoS computational biology 7 (11) (2011) e1002211.
[17] D. Pecevski, L. Buesing, W. Maass, Probabilistic inference in general
graphical models through sampling in stochastic networks of spiking
neurons, PLoS computational biology 7 (12) (2011) e1002294.
[18] R. Legenstein, W. Maass, Ensembles of spiking neurons with noise support optimal probabilistic inference in a dynamically changing environment 10 (10) (2014) e1003859.
[19] L. Shi, T. L. Griffiths, Neural implementation of hierarchical bayesian
inference by importance sampling, in: Advances in neural information
processing systems, 2009, pp. 1669–1677.
[20] W. J. Ma, J. M. Beck, P. E. Latham, A. Pouget, Bayesian inference
with probabilistic population codes, Nature neuroscience 9 (11) (2006)
1432–1438.
[21] W. J. Ma, M. Jazayeri, Neural coding of uncertainty and probability,
Annual review of neuroscience 37 (2014) 205–220.

31

[22] A. E. Orhan, W. J. Ma, Neural population coding of multiple stimuli,
The Journal of Neuroscience 35 (9) (2015) 3825–3841.
[23] J. M. Beck, W. J. Ma, R. Kiani, T. Hanks, A. K. Churchland, J. Roitman, M. N. Shadlen, P. E. Latham, A. Pouget, Probabilistic population
codes for bayesian decision making, Neuron 60 (6) (2008) 1142–1152.
[24] J. M. Beck, P. E. Latham, A. Pouget, Marginalization in neural circuits
with divisive normalization, The Journal of Neuroscience 31 (43) (2011)
15310–15319.
[25] W. J. Ma, M. Rahmati, Towards a neural implementation of causal
inference in cue combination, Multisensory Res 26 (2013) 159–176.
[26] D. R. Wozny, U. R. Beierholm, L. Shams, Human trimodal perception
follows optimal statistical inference, Journal of vision 8 (3) (2008) 24.
[27] R. van den Berg, M. Vogel, K. Josić, W. J. Ma, Optimal inference
of sameness, Proceedings of the National Academy of Sciences 109 (8)
(2012) 3178–3183.
[28] D. Koller, N. Friedman, Probabilistic graphical models: principles and
techniques, MIT press, 2009.
[29] C. Andrieu, N. De Freitas, A. Doucet, M. I. Jordan, An introduction to
mcmc for machine learning, Machine learning 50 (1-2) (2003) 5–43.
[30] A. Steimer, W. Maass, R. Douglas, Belief propagation in networks of
spiking neurons, Neural Computation 21 (9) (2009) 2502–2523.
[31] S. Litvak, S. Ullman, Cortical circuitry implementing graphical models,
Neural computation 21 (11) (2009) 3010–3056.
[32] H. Hachiya, M. Sugiyama, N. Ueda, Importance-weighted least-squares
probabilistic classifier for covariate shift adaptation with application to
human activity recognition, Neurocomputing 80 (2012) 93–101.
[33] J. Cheng, M. J. Druzdzel, Ais-bn: An adaptive importance sampling
algorithm for evidential reasoning in large bayesian networks, Journal
of Artificial Intelligence Research 13 (2000) 155–188.

32

[34] P. O. Hoyer, A. Hyvarinen, Interpreting neural response variability as
monte carlo sampling of the posterior, Advances in neural information
processing systems (2003) 293–300.
[35] R. L. De Valois, E. W. Yund, N. Hepler, The orientation and direction
selectivity of cells in macaque visual cortex, Vision research 22 (5) (1982)
531–544.
[36] D. M. Coppola, L. E. White, D. Fitzpatrick, D. Purves, Unequal representation of cardinal and oblique contours in ferret visual cortex, Proceedings of the National Academy of Sciences 95 (5) (1998) 2621–2623.
[37] C. S. Furmanski, S. A. Engel, An oblique effect in human primary visual
cortex, Nature neuroscience 3 (6) (2000) 535–536.
[38] M. Kouh, T. Poggio, A canonical neural circuit for cortical nonlinear
operations, Neural computation 20 (6) (2008) 1427–1451.
[39] B. Nessler, M. Pfeiffer, W. Maass, Stdp enables spiking neurons to detect hidden causes of their inputs, in: Advances in neural information
processing systems, 2009, pp. 1357–1365.
[40] B. Nessler, M. Pfeiffer, L. Buesing, W. Maass, Bayesian computation emerges in generic cortical microcircuits through spike-timingdependent plasticity 9 (4) (2013) e1003037.
[41] D. Kappel, B. Nessler, W. Maass, Stdp installs in winner-take-all circuits an online approximation to hidden markov model learning, PLoS
computational biology 10 (3) (2014) e1003511.

33

