Modeling time-dependent system

STATE SPACE MODEL

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

1 	Demo	Hippocampal place cells recorded in the Wilson lab	
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter2 	Two-stages analysis of neural signal	

Brown et al. Journal of Neuroscience 1998; 18: 7411-7425.
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

3 	

Discrete-time Poisson-GLM
Δ kth-bin

Covariate signal	xk

Point process	110011000yk

1001

10(01P ( yk | x k , β ) = exp yk log {λk ( x k , β ) Δ} − λk ( x k , β ) Δ

λk ( x k , β ) Δ = e
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	00)

β0 +xTkβState-space model and recursive Bayesian filter

4 	

State-space model
Hidden statex k−1xk

Observationy k−1ykState model: GaussianHyper-parameters:

x1 ~ N (µ, Σ)
x k = x k−1 + ξ kw = [Q, µ, Σ]ξ k ~ N ( 0, Q)

Observation model: Point process (discretized)(P ( yk | x k , β ) = exp yk log {λk ( x k , β ) Δ} − λk ( x k , β ) Δ

λk ( x k , β ) Δ = e
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	)β0 +xTkβState-space model and recursive Bayesian filter

5 	Posterior density of the state
Bayes’ rule:
Posterior

Likelihood

Prior (State model)

P ( y1:K | x1:K , β ) p ( x1:K | w )p ( x1:K | y1:K , β, w ) =P ( y1:K | β, w )Evidenceβ, w parameters

Once we construct the posterior, we can obtain	MAP estimate: the most likely path of the hidden state.
Credible interval: an analog of confidence interval in Bayesian estimation.

We obtain the joint posterior density => Bayesian recursive filter.	
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

6 	

Recursive Bayesian filter
One-step prediction(

p x k | y1:k−1

x k−1Filter(

p x1:k | y1:k)y k−1

Filter at k-th step
Bayes’ rulexk)

ykLikelihood

One-step prediction

P ( yk | x k , β ) p ( x k | y1:k−1, β, w )p ( x k | y1:k , β, w ) =P ( yk | x k−1, β, w )Chapman-Kolmogorov equationp ( x k | y1:k−1, β, w ) = ∫ p ( x k | x k−1, w ) ⋅ p ( x k−1 | y1:k−1, β, w ) dx k−1
Filter density at (k-1)-th step
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

7 	

Recursive Bayesian filter	
xk−1

xkFilter distribution
Filter distribution

observation

Prediction

0t-1k −1

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	

One-step prediction distribution

ktState-space model and recursive Bayesian filter

8 	Methods to obtain the posterior
Methods for obtaining the posterior density to name a few…
Analytical methods
Gaussian approximation (Laplace’s method).
Conjugate prior.
Expectation propagation.
Monte Carlo methods
Sequential importance resampling (Particle filter).
Markov chain Monte Carlo (MCMC).

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

9 	One-step prediction density
Filter distribution
Filter distribution

observation

Prediction

0t-1

One-step prediction distributiontChapman-Kolmogorov equationp ( x k | y1:k−1, w ) = ∫ p ( x k | x k−1, w ) ⋅ p ( x k−1 | y1:k−1, β, w ) dx k−1

N ( x k−1,Q )N ( x k−1|k−1,Wk−1|k−1 )Assumption

One-step prediction is a normal density with	mean	x k|k−1 = x k−1|k−1

covariance	 Wk|k−1 = Wk−1|k−1 + Q
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter10 	Filter density
Filter distribution
Filter distribution

observation

Prediction

0t-1

One-step prediction distributiontFilter density

Likelihood

One-step prediction densityp ( x k | y1:k , β, w ) =P ( yk | x k , β ) p ( x k | y1:k−1, β, w )p ( yk | x k−1, β, w )∝ exp[yk log {λk ( x k , β ) Δ} − λk ( x k , β ) Δ
1
T−1− ( x k − x k|k−1 ) Wk|k−1( x k − x k|k−1 )]2
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter11 	Log-concave functions
Rationale
Log-concave:the function is log-concave if the logarithm of the function is concave.

Concave function.
Exponential family of distributions are in general log-concavewith respect to parameters in ‘natural form’.
A function obtained by multiplication of two log-concave functions
is log-concave.
Log-concave

Log-concave Log-concavep ( x k | y1:k , β, w ) ∝ P ( yk | x k , β ) p ( x k | y1:k−1, β, w )Paninski, L. (2005). Advances in NIPS 17 1025–1032
Boyd S. (2004). Convex Optimization
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

12 	Newton-Raphson method
Newton-Raphson method to find a posterior mode
Posterior mode (MAP estimate):

Posterior modex k|k ≡ argmax f ( x k )xkLog-posterior:f ( x k ) = yk log ( λk Δ ) − λk Δ −1
T−1( x k − x k|k−1 ) Wk|k−1( x k − x k|k−1 )2

Newton-Raphson method:−1

oldx newk = x k − ( ∇∇f ) ∇fold

Hessian Gradient evaluated at x k

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

13 	Gradient and Hessian∂log λk−1− Wk|k−1( x k − x k|k−1 )∂x kgradient∇f ( x k ) = ( yk − λk Δ )Hessian∂f ( x k )∂log λk∂log λk∂log λk−1H≡=−λΔ+y−λΔ−W()
kk
kk|k−1∂x k∂x Tk∂x T∂x k∂x k∂x Tk

If we use a canonical link function

λk = e β 0 + β x kgradient−1∇f ( x k ) = ( yk − λk Δ ) β − Wk|k−1( x k − x k|k−1 )Hessian−1
H = −β T λk Δβ − Wk|k−1

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

14 	Gaussian approximation

Laplace’s method for a Gaussian approximationx k|k∇f ( x k ) = 0

Filter covariance: Wk|kWk|k = −H −1

Filter mean:

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	

x k|k

State-space model and recursive Bayesian filter

15 	

Recursive nonlinear filtering
Filter distribution
Filter distribution

observation

Prediction

0One-step prediction distribution

t-1t

Recursive nonlinear filtering
One-step prediction
Filter mean:x k|k = x k|k−1 + Wk|k−1β ( yk − λk ( x k|k ) Δ )Solved by
Newton-Raphson Method

Observation error (innovation)−1−1T

Filter covariance: W k|k = Wk|k−1 + β λk ( x k|k ) Δβ
One-step prediction
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	

Information given by data at k-th step.
State-space model and recursive Bayesian filter16 	Filtering/Smoothing
Forward recursion

k = 1,2,, K

One-step predictionMeanx t|t−1 = x t−1|t−1Covariance

Wt|t−1 = Wt−1|t−1 + Q

Recursive nonlinear filtering
Meanx k|k = x k|k−1 + Wk|k−1β ( yk − λk ( x k|k ) Δ )Covariance−1TW −1=W+βλk ( x k|k ) Δβk|kk|k−1

Solved by
Newton-Raphson Method

Backward recursion k = K −1, K − 2,,2,1
Fixed interval smoothing
Meanx k|K = x k|k + Ak "#x k+1|K − x k+1|k $%Covariance

Wk|K = Wk|k + Ak "#Wk+1|K −Wk+1|k $% ATkHideaki Shimazaki, Ph.D. http://goo.gl/viSNG	

Ak = Wk|kW −1
k+1|kState-space model and recursive Bayesian filter

17 	Joint estimation of posterior and parameters	Expectation-Maximization algorithm for joint estimation of posterior andparameters.	
(ref. Smith & Brown, 2003)
E-step: Given the parameters, construct the posterior
Posterior

LikelihoodPrior

P ( y1:K | x1:K , β ) p ( x1:K | w )p ( x1:K | y1:K , β, w ) =P ( y1:K | β, w ) Evidence(marginal likelihood)M-step: Given the posterior, optimize the parametersOptimization of the parameters bymaximizing the (marginal) likelihood.	

w MLE = argmax log P ( y1:K | β, w )w

For this goal, we maximize a lower bound of the marginal likelihood (Q-function)	

Q ( w new | w ) ≡ Ex1:K |y1:K , β ,w "#log P ( y1:K , x1:K | β, w new )$%Expected complete data log-likelihood.

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

18 	Derivation of the lower bound	
The marginal log-likelihood function can be bounded as follows. 	
log P ( y1:K | β, w new )= log ∫ P ( y1:K , x1:K | β, w new ) dx1:K= log ∫ p ( x1:K | y1:K , β, w )P ( y1:K , x1:K | β, w new )p ( x1:K | y1:K , β, w )dx1:K" P ( y , x | β, w new ) %
1:K
1:K
'
= log Ex1:K |y1:K , β ,w $$# p ( x1:K | y1:K , β, w ) '&"P ( y1:K , x1:K | β, w new ) %'≥ Ex1:K |y1:K , β ,w $log
p ( x1:K | y1:K , β, w ) '&$#Entropy of the posterior	= Ex1:K |y1:K , β ,w "#log P ( y1:K , x1:K | β, w new )%& − Ex1:K |y1:K , β ,w "#log p ( x1:K | y1:K , β, w )%&Q-function	
Hence we obtain 	Irrelevant to 	w

newlog P ( y1:K | β, w new ) ≥ Q ( w new | w )Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

19 	Model selection	( )()p y1:T = ∫ p y1:T ,x1:T d x1:T
Data

Data Hidden state

Bayesian model selection

Bayes Factor()( ) p y |M()

B12 y1:T =p y1:T | M 1
1:T

2(Jeffrey 61, Kass&Raftery 95)

Model Selection penalized by model dimension()

ABIC = −2log ∫ p y1:T ,x1:T | w d x1:T + 2 × model dimension
Log-quadratic approximation

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

20 	

APPLICATIONS OF THE
STATE-SPACE MODEL
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

21 	

1998

Two stages analysisHideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

22 	Encoding model	

1998

Likelihood function

Conditional intensity function

Hippocampal place cell
receptive field model

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	

Theta phase modulation of spike rate

State-space model and recursive Bayesian filter

23 	

1998

Goodness-of-fit of the encoding model
p-values of observed spike countagainst a null hypothesis of aPoisson model

Cells are more variable thana Poisson

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	

Residual analysis of the state-model

Test rejected a bivariate Gaussianmodel.
Assumption of independence wasalso rejected.
State-space model and recursive Bayesian filter

24 	

Recursive point-process filter	

1998

Posterior (Filter)
Data up to time t_k
Data up to time t_k-1

Observation

One-step prediction

Filter at t_k-1
State equation
Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter25 	

Algorithm for a point process filter	

Filter density
Filter density

Observation
Prediction

0t-1

One-step prediction densityt

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

26 	

1998

DECODING

Abrupt change
was not captured.

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	

TIMEState-space model and recursive Bayesian filter

27 	

1998

DECODING

No difference.

ML was the second best,
however with a lot of outliers.

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

28 	

2001

1D animal position

Change in the receptive fieldåspikes

Conditional intensity model

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	

Update ruleState-space model and recursive Bayesian filter

29 	

Demo	Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	

2001

State-space model and recursive Bayesian filter

30 	

2004

•
•
•

A full mathematical formulation of an adaptive point process filter.
Fast approximation of the adaptive point process filter.
Simulation study on (1) tracking place filed dynamics, (2) simultaneous estimationof receptive field dynamics and arm trajectory (decoding).

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter31 	Selected references of state-space analyses	
•Brown EN, Frank LM, Tang D, Quirk MC, Wilson MA. A statistical paradigm for neural spike train decoding applied to position
prediction from ensemble firing patterns of rat hippocampal place cells, Journal of Neuroscience 1998; 18: 7411-7425.

•Receptive field plasticity–––•

Brown EN, Nguyen DP, Frank LM, Wilson MA, Solo V. An analysis of neural receptive field plasticity by point process adaptive filtering.
Proceedings of the National Academy of Sciences 2001; 98: 12261-12266. PMID: 11593043
Frank LM, Eden UT, Solo V, Wilson MA, Brown EN. Contrasting patterns of receptive field plasticity in the hippocampus and the entorhinal
cortex: an adaptive filtering approach. Journal of Neuroscience 2002; 22: 3817-30. PMID: 11978857
Eden UT, Frank LM, Barbieri R, Solo V, Brown EN, Dynamic analyses of neural encoding by point process adaptive filtering, NeuralComputation, 2004, 16(5): 971-998. PMID: 15070506

Multiple neuron GLM-point process–––

Truccolo W, Eden U, Fellow M, Donoghue JD, Brown EN. A point process framework for relating neural spiking activity to spiking history, neuralensemble and covariate effects. Journal of Neurophysiology, (published online Sept. 8, 2004), 2005, 93: 1074-1089.
Okatan M, Wilson MA, Brown EN. Analyzing functional connectivity using a network likelihood model of ensemble neural spiking activity. NeuralComputation, 2005, 17(9): 1927-1961.
Czanner G, Eden UT, Wirth S, Yanike M, Suzuki WA, Brown EN. Analysis of between-trial and within-trial neural spiking dynamics. Journal ofNeurophysiology, 2008, 99: 2672-2693. PMID: 18216233•EM-algorithm (Joint state-space and parameter optimization)•Behavioral analysis––––•Smith AC, Brown EN. Estimating a state-space model from point process observations. Neural Computation. 2003; 15: 965-91. PMID: 12803953
Smith AC, Frank LM, Wirth S, Yanike M, Hu D, Kubota Y, Graybiel AM, Suzuki W, Brown EN. Dynamic analysis of learning in behavioral
experiments, Journal of Neuroscience, 2004, 15: 965-91. PMID: 14724243
Smith AC, Stefani MR, Moghaddam B, Brown EN. Analysis and design of behavioral experiments to characterize population learning. Journal ofNeurophysiology (published on line Sept. 29, 2004), 2005, 93: 1776-1792.
Smith AC, Wirth A, Suzuki W, Brown EN. Bayesian analysis of interleaved learning and response bias in behavioral experiments. Journal ofNeurophysiology, 2007, Mar; 97(3):2516-24. PMID: 17182907

Motor prosthetics––––Brockwell, AE, Rojas, A L, Kass, RE, Recursive Bayesian decoding of motor cortical signals by particle filtering. Journal of Neurophysiology,
2004, 91(4) 1899-1907
Srinivasan L, Eden UT, Willsky AS, Brown EN. A state-space analysis for reconstruction of goal-directed movements using neural signals.
Neural Computation, 2006, 18(10): 2465-2494. PMID: 16907633
Srinivasan L, Brown EN. A state-space framework for movement control to dynamic goals through brain-driven interfaces. IEEE Transactions on
Biomedical Engineering, 2007, 54(3):526-535.
Srinivasan L, Eden UT, Mitter SK, Brown EN. General purpose filter design for neural prosthetic devices. Journal of Neurophysiology, 2007,
98(4): 2456-2475. PMID: 17522167

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

32 	What we learned	12345• Framework of a state-space model.
• Recursive Bayesian filter (Laplace’s approximation).
• Simultaneous estimation of posterior and parameters (EM-algorithm).
• Model validation in a Bayesian framework (Bayes factor, ABIC),• Applications to neural decoding and plasticity.

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	State-space model and recursive Bayesian filter

33 	

Relations between GLM and Max-ent

Bernoulli-GLM model

Maximum entropy model

Neuron 1
Neuron 2
Neuron 3
Neuron 4

𝑘−1 𝑘
Conditional independence

Hideaki Shimazaki, Ph.D. http://goo.gl/viSNG	

𝑘−1𝑘Joint probability mass function

State-space model and recursive Bayesian filter

34