Detec%ng Simultaneous
Changepoints in Mul%ple Sequences
Joint with Lester Mackey
April 25, 2017

Introduc%on

Three data examples:
Copy number varia%on in cancer

Three data examples:
Historical vola%lity of U.S. stocks

Three data examples:
Conforma%on of a protein molecule

Characteris%cs of data
• Each sequence has “piecewise homogeneous”
behavior separated by changepoints• Changepoints represent boundaries of copy
number aberra%ons, economic events,
conforma%onal changes in protein structure,etc.
• Goal: Iden%fy the changepoints(Shewhart 1925, ES Page 1954)

Characteris%cs of data

In several modern applica%ons:– We have mul%ple (oXen%mes many) sequences– Changepoints tend to occur at the same sequen%alposi%ons– Obtain greater detec%on power and beZer summary ofchangepoints by pooling across sequences

Approaches to inference• Recursive segmenta%on with scan sta%s%cs
Zhang, Siegmund, Ji, Li 2010
Siegmund, Yakir, Zhang 2011
Zhang, Siegmund 2012
Jeng, Cai, Li 2013

• Penalized maximum likelihood
Fan, Dror, Mildorf, Piana, Shaw 2015

• Total-varia%on/rank penaliza%on
Nowak, Has%e, Pollack, Tibshirani 2011
Zhou, Yang, Wan, Zhao, Yu 2013

Our approach
Bayesian Analysis of SImultaneous Changepoints(BASIC):– A full probabilis%c, genera%ve model– Changepoints modeled as unobserved random
events– Prior changepoint distribu%on encodes frequencyof occurrence and simultaneity across sequences– Infer changepoints by es%ma%ng/compu%ng theirposterior probabili%es, given the observed dataModelBASIC ModelObserve J sequences of length T(Xj,t : 1  j  J, 1  t  T ) 2 RJ⇥

T 1J⇥T(q:2tT)
2[0,1]t
Step 1: Generate“changeprobabili%es”
X 2J⇥TRX2R
iid
qt 2 [0, 1]
qt ⇠ ⇡ Q
qt 2 [0, 1]⇡Q = 0.9 0 + 0.1 0.2
according to qt ⇠qt⇡iid
Q⇠ , ⇡e.g.
Q⇡Q =+ 0.1+ 0.2⇡Q0.9= 00.90 0.1(Z: 10.2 j  J, 2  t  T ) 2 {0, 1}
j,tZj,tZ2 {0,21}j,t {0, 1}
indZj,t ⇠ Bernoulli(q)t Zj,t ⇠ Bernoulli(qt)indZj,tt)✓j,t⇠2 Bernoulli(q⇥✓j,td(✓2j,t⇥ : 1  j  J, 1  t  T ) 2 ⇥J⇥TX 2 RJ⇥TJ⇥TX
2R
qt 2 [0, 1]
qtiid
2 [0, 1]qtq ⇠⇠ ⇡⇡QtQ
Step 2: Generate⇡⇡Qchangepoints
==0.9
++0.1
00.9
0 0.1 0.20.2
Q
ZZj,t 22{0,1}{0,1}j,tZj,tind
Bernoulli(qt))
according toZj,t⇠⇠ Bernoulli(qtj,t22⇥⇥✓✓j,tdd⇥⇢R⇥⇢R⇢⇢indind⇠ ⇡⇡⇥✓⇠j,t = ✓ ⇥✓j,tj,t 11= ✓j,ttt==11j,t==11ZZj,tZZj,t==00

BASIC Model⇡Q = 0.9 0 + 0.1⇡Q0.2= 0.9 0 + 0.1 0.2⇡Q = 0.9 0 + 0.1
0.2
Zj,t 2 {0, 1} Zj,t 2 {0, 1} Zj,t 2 {0, 1}
Zj,tZj,t )⇠ Bernoulli(q) ⇠ Bernoulli(qt)tZj,t ⇠ Bernoulli(qtd ✓ 2 ⇥ ⇢ Rdd
Step 3:✓ Generateparameters✓
2⇥⇢
Rj,tj,t 2 ⇥ ⇢ Rnj,tn⇠ ⇡n⇠ ⇡Zj,t = 1⇠⇡1⇥Zj,t = 1 ⇥Z=✓j,t = ✓according
to ✓j,tj,t = ✓Zj,t = 0✓j,t = ✓ ⇥Z=0j,t1j,t1j,tj,t 1 Zj,t = 0⇡5)⇥ = Normal(0, 5)⇡=Normal(0,⇥e.g. ⇡⇥ = Normal(0, 5)
Xj,t 2 RX
2Rj,tXj,t 2 R
indindindX⇠ p(·|✓j,t)X⇠p(·|✓)
j,tj,tXj,t ⇠ p(·|✓j,t)j,tp(x|✓)= Normal(✓, 1)p(x|✓)=Normal(✓,1)
p(x|✓) = Normal(✓,1)
ZZZk (1 q)J k ⇡ (dq)kJkq
⇡Q(dq)Qq k (1 q)J k ⇡qQ(1(dq) q)

BASIC Model⇡=X

Q q=⇡ ⇡(q)X

⇡Qq =
⇡Q(q)

X

⇡Q(q) qj,tn⇠ ⇡dd✓
2⇥⇢
RZj,t = 1
✓j,t 2j,t⇥⇢R⇥✓j,t = ✓n⇠ n
Zj,t = 0⇠ ⇡⇥ Z Z=1j,t1⇡=1j,tj,t⇥✓j,t ✓j,t= ✓j,tZ 1 Z= ✓j,t0 = 0⇡⇥ = Normal(0, 5)
1 j,t =j,t= Normal(0,Step 4: Generatedata⇡⇥ =⇡⇥Normal(0,5) 5)
Xj,t 2 R
Xj,t ⇠ p(·|✓j,t)
Xj,tX2j,tR2 Rindj,t), e.g. p(x|✓) = Normal(✓, 1)
according to Xj,tX⇠ p(·|✓j,t ⇠ p(·|✓j,t) Zp(x|✓)= Normal(✓,1) 1) kp(x|✓)=Normal(✓,Z Zq (1 q)J k ⇡Q(dq)
Xq k (1q k (1q)J q)kJ⇡Qk(dq)⇡Q(dq)⇡Q =⇡Q(q) q
XX⇡Q =⇡ = ⇡Q(q)q2Sq(q) q⇡!
QQ Znq2S q2SY! !
Z ZYp(xi|✓) ⇡⇥(d✓)
n nYp(xi|✓)⇡⇥(d✓)
p(xi|✓)⇡⇥(d✓) i=1BASIC ModelBASIC ModelqZΘX• Observed data: X• Unobserved (latent) variables: q, Z, Θ• Prior distribu%ons: πQ , πΘPrimary inference task: Infer Z, upon observing XThree inference procedures
We develop three inference procedures:
1. Compute changepoint probabili%es
Sample from Pr(Z|X) using Markov Chain Monte Carlo(MCMC).

2. Compute most likely changepoint loca%ons
Maximize Pr(Z|X) over Z

3. Fit the genera%ng model
Es%mate priors πQ , πΘ from the data using Monte Carlo
Expecta%on-Maximiza%on (MCEM)

Simple example

Pr(Zj,t=1|X) computed using 100 MCMC
samples under true πQ , πΘ

Simple example

Pr(Zj,t=1|X) computed using 100 MCMC
samples under (very) wrong πQ , πΘSimple example

Pr(Zj,t=1|X) computed using 100 MCMC samples
under πQ , πΘ es%mated from the dataSimple example

True changepoints in black
Changepoints maximizing Pr(Z|X) in red

Comparison with Hidden Markov
Models (HMMs)
Changepoint (one sequence):• “Inﬁnite” number of states
which never recur• Compute changepointprobabili%es• Compute most likely
changepoint loca%ons• Fit the genera%ng model(πQ , πΘ)

HMM:• Finite number of recurringstates• Compute state probabili%es(forward-backward algo.)• Compute most likely states(Viterbi algo.)• Es%mate Markov transi%onmatrix and data distribu%on
parameters(Baum-Welch/EM)

Inference(Xj,t : 1  j  J, 1  t  T ) 2 RJ⇥TIdea #1: Marginaliza%on
indXj,t ⇠ p(·|✓j,t)p(x|✓)= Normal(✓,Model has latentvariables{qt},1){Zj,t}, {Θj,t}, but…
Z
• We can analy%callymarginalizeout {qt}:kJ
kq (1 q)⇡Q(dq)
Use a discrete mixtureX⇡Q =⇡Q(q) qq2S

• We can analy%cally marginalize out {Θj,t}:
Use a conjugate prior πΘ.
Take-away: Sample from and maximize overPr(Z|X) directly, rather than Pr(q,Z,Θ|X).
2

Idea #2: Row/column-wise opera%ons
Naïve approach:• Sample from Pr(Z|X) (Gibbs):– Arbitrarily ini%alize all Zj,t– For each (j,t): Sample Zj,t from {0,1}, given remaining
variables Z• Maximize Pr(Z|X) (coordinate descent):– Arbitrarily ini%alize all Zj,t– For each (j,t): Maximize Zj,t over {0,1}, given remaining
variables Z

Will get stuck in local maxima of Pr(Z|X), because ofstrong correla%ons between diﬀerent variables Zj,t

Idea #2: Row/column-wise opera%ons
Row/column-wise sampling:
1. For each j, sample Zj, from Pr(Zj,|X,Z(-j),)
2. For each t, sample Z,t from Pr(Z,t|X,Z,(-t))

This exactly samples all Zj,t in a single row orcolumn from their joint condi%onal distribu%on.

Row-wise sampling

Sample Zj,1, marginalizing over Zj,2 to Zj,T

Row-wise sampling

Sample Zj,2, marginalizing over Zj,3 to Zj,T

Row-wise sampling

Normaliza%on constants computed recursively
for t=T, T-1, … (Fearnhead 2006)

Column-wise sampling
To sample Z,t from Pr(Z,t|X, Z,(-t)):
1. Sample Z1,t, marginalizing over Z2,t to ZJ,t.
2. Sample Z2,t, marginalizing over Z3,t to ZJ,t.
3. Etc.
Normaliza%on constants computed recursively,
for j=J, J-1, …

Comparison of sampling approaches

Idea #2: Row/column-wise opera%ons
Row/column-wise maximiza%on:
1. For each j, maximize Pr(Z|X) over Zj,2. For each t, maximize Pr(Z|X) over Z,tDynamic programming recursion used in step 1.
(Jackson et al 2005)

Z

q2S!

nY
Fiong the model

Maximum likelihood: i=1p(xi|✓)⇡⇥(d✓)⇡Q, ⇡⇥ = argmax⇡ ,⇡ log p(X; ⇡Q, ⇡⇥)Q ⇥⇤(i+1) Hard(i+1)because of marginaliza%on over⇥Z⇡,⇡= argmax⇡ ,⇡ Elog p(X, Z; ⇡Q, ⇡⇥)⇥
QQ ⇥ Z|X;⇡ (i),⇡ (i)Q ⇥1

ZZp(x|✓) = Normal(✓, 1)

q2S!

Y
Fiong the model
nqk (1⇡Q

Maximum likelihood: i=1q)J k ⇡Q(dq)X|✓) ⇡ (d✓)
p(x= i ⇡ (q) q⇥
Qq2S!
Z Y
⇡Q, ⇡⇥ = argmaxn⇡ ,⇡ log p(X; ⇡Q, ⇡⇥)Q ⇥p(xi|✓) ⇡⇥(d✓) ⇥⇤(i+1) Hard(i+1)because of marginaliza%on over Z⇡,⇡= argmax⇡ ,⇡ i=1E(i) (i) log p(X, Z; ⇡Q, ⇡⇥)⇥Q⇥
Qlog p(X; ⇡ , ⇡ )EM algorithm: ⇡Q, ⇡⇥ = argmax⇡Z|X;⇡
Q,⇡⇥ Q ,⇡⇥Q ⇥⇥⇤(i+1) (i+1)⇡,⇡= argmax⇡ ,⇡ E1(i) (i) log p(X, Z; ⇡Q, ⇡⇥)⇥QQ ⇥Z|X;⇡

Q,⇡

S%ll hard to compute the E-step
1⇥ZZp(x|✓) = Normal(✓, 1)

q2S!

YFiong the model
nqk (1⇡Q

Maximum likelihood: i=1q)J k ⇡Q(dq)X|✓) ⇡ (d✓)
p(x= i ⇡ (q) q⇥
Qq2S!
Z Y
⇡Q, ⇡⇥ = argmaxn⇡ ,⇡ log p(X; ⇡Q, ⇡⇥)Q ⇥p(xi|✓) ⇡⇥(d✓) ⇥⇤(i+1) Hard(i+1)because of marginaliza%on over Z⇡,⇡= argmax⇡ ,⇡ i=1E(i) (i) log p(X, Z; ⇡Q, ⇡⇥)⇥Q⇥
Qlog p(X; ⇡ , ⇡ )EM algorithm: ⇡Q, ⇡⇥ = argmax⇡Z|X;⇡Q,⇡⇥ Q ,⇡⇥Q ⇥⇥⇤(i+1) (i+1)⇡,⇡= argmax⇡ ,⇡ E1(i) (i) log p(X, Z; ⇡Q, ⇡⇥)⇥QQ ⇥Z|X;⇡

Q

,⇡

⇥

S%ll hard to compute the E-stepBX⇥ MCEM: Approximate⇤ 1 the E-step by (b)B
1log p(X, Z; ⇡Q, ⇡⇥) ⇥ ⇡E

⇤ 1ZX ; ⇡Q, ⇡⇥(b)log p(X,)

log p(X, Z
B Z; ⇡Q, ⇡⇥) ⇡ B(i) (i) log p(X,Z|X;⇡ ,⇡b=1Q ⇥b=1
(1)(B) Z (1), . . . , Z (B)where Z,...,Z; ⇡ Q , ⇡⇥ )are B Monte Carlo samples

CNV Simula%on

CNVs in 50 “cancer” samples
Simulator from Louhimo et al, Nature methods
2012:• J = 50 samples (sequences)• T = 1000 probe loca%ons• 6 types of CNVs:– Dele%ons of both or 1 copy
– Ampliﬁca%ons of 1 to 8 copies– Sample heterogeneity from 30% to 70%– Shared between 5% to 50% of all sequences

CNVs in 50 “cancer” samples

Comparison of methods

Applied BASIC using normal data model with changing mean.
Computed most-likely changepoint loca%ons, aXer ﬁong
model priors using MCEM.

BASIC vs. SIMPLE• BASIC: Bayesian inference with “simultaneouschangepoint” prior• SIMPLE (Fan et al 2015): Penalized maximum likelihood
with “simultaneity-inducing” penal%es• Key diﬀerences:- Tuning parameters in SIMPLE are user-speciﬁed; priors
in BASIC are es%mated from the data
- SIMPLE returns only the most likely changepoints;
BASIC also es%mates changepoint occurrence
probabili%es- MCMC sampling/MCEM prior es%ma%on using BASIC isslower

CNVs in NCI-60 cancer cell lines

CNVs in NCI-60 cell lines• 60 cell lines derived from cancers of various%ssues and organs• Copy numbers measured using Agilent 44K
oligonucleo%de microarray• 2 technical replicates for most cell lines• T = 40,217 probe loca%ons (sex chromosomesremoved)

4 Melanoma samples, Chromosome 11

ME:SK MEL 20−1101

ME:SK MEL 28

0−1101

ME:MDA MB 435

0−1101

ME:MDA N

0−110

E[qt | X]

Chromosome 1 position

MbM
bM
b240230M
b220210M
bM
b200M
b190

180M
b170M
b160entromere150M
bCM
bM
10 b0M
11 b0M
b90

MbM
b80

7060

M
b50Mb40Mb30Mb20M
b100Mb0

Cross-replicate agreement

Focal CNVs

CFHR3, CFHR1

FAF1, CDKN2C

Focal CNVs

FBXW7

Focal CNVs

MYC

CDKN2A, CDKN2B

Focal CNVs

CYP2E1

Focal CNVs

A2BP1

URI1

APOBEC3A, APOBEC3B

Vola%lity changes in S&P 500 stocks

S&P 500 stock returns• Daily returns of U.S. stocks in the S&P 500
from 2000 to 2009
• J = 401 total stocks (sequences)• T = 2514 trading days
• Applied BASIC with Laplace data model, ﬁxed
mean = 0, changing scale17,200
6JSe ul 7
p ,21 0
D 5, 2 08
ec 0
9, 08
20
M08
ay
19,2009Nov
10,2009Ju
l24N , 20
ov
0
1, 7
2007

Aug

BKO
ct
14,Fe 200
4bM 9, 2
ay 0
05, 5
2005

BAC

Ju
l1,200N2
ov
29,200Ap2
r29,Jul 2 200
8, 32003

AON

Ap
r20,20Se01
p17,200
1AIG

Ap
r27,2000

4 U.S. stocks, 2000-2010

0.50.0−0.5
100.2

0.0−0.2

100.2−0.20.010.200.0

−0.2
10

E[qt | X]
01JSe ul 7p ,21 0
D 5, 2 08
ec 09, 0820
M08
ay19,200
9N
ov
10,200
9

BK
0.0

Ju
l24N , 20
ov
0
1, 7
20
07

−0.2

17,200
6BAC

Au
g

0.2

O
ct
14,Fe 200
4bM 9, 2
ay 005, 5
2005

0

Ju
l1,200N2ov
29,200Ap2
r29,Jul 2 200
8, 32003

AON

Ap
r20,20Se01
p17,200
1AIG

Ap
r27,2000

4 U.S. stocks, 2000-2010

0.50.0−0.5
100.2

0.0−0.2

1

Sep 17, 2001: First trading dayaXer 9/11 terrorist aZacks

0.2
00.0

−0.2
10E[qt | X]
01JSe ul 7
p ,21 0
D 5, 2 08
ec 0
9, 08
20
M08
ay
19,2009N
ov
10,200
9

BK0.0

Ju
l2
4
N , 20
ov
01, 7
20
07

−0.2

17,200
6

BAC

Aug

0.2

O
ct
14,Fe 200
4bM 9, 2
ay 0
05, 5
2005

0Jul1,200
N2
ov
29,200
Ap2r29,Jul 2 200
8, 32003AONApr20,20
Se01
p17,200
1AIG

Apr27,2000

4 U.S. stocks, 2000-2010

0.5

0.0

−0.5
100.2

0.0

−0.2

1Jul 1, 2002: Collapse ofWorldCom (93% price drop)

0.2
00.0

−0.2
10

E[qt | X]
0

−0.2
1

J
Se ul 7
p ,21 0
D 5, 2 08
ec 09, 0820
M08
ay
19,2009N
ov
10,20090.0Jul24N , 20
ov
01, 7
20
07

BK17,200
60Aug0.2O
ct14,Fe 200
4bM 9, 2
ay 005, 52005

BACJu
l1,200N2ov
29,200
Ap2r29,Jul 2 200
8, 3
2003

AONApr20,20
Se01
p17,200
1AIG

Apr27,2000

4 U.S. stocks, 2000-2010

0.5

0.0

−0.5
10

0.2

0.0

−0.2

10

0.2

−0.2
0.0

1

Oct 14, 2004: AIG / Marsh &McLennan insurance scandal

0E[qt | X]
0BKJSe ul 7
p ,21 0
D 5, 2 08
ec 0
9, 08
20
M08
ay
19,2009Nov
10,20090Jul24N , 20
ov
01, 7
2007117,200
6−0.2Aug

BAC0.0O
ct14,Fe 200
4bM 9, 2ay 005, 5
2005

0.2

Jul1,200
N2ov
29,200
Ap2r29,Jul 2 2008, 32003

AONAp
r20,20Se01
p17,200
1AIG

Ap
r27,2000

4 U.S. stocks, 2000-2010

0.5

0.0

−0.5
10Sep 15, 2008: Lehmann Bros.
ﬁles Chapter 11 bankruptcy

0.2

−0.2
0.0

1

0.2
00.0

−0.2
10E[qt | X]
0Summary• In many applica%ons, changepoints are
simultaneous across mul%ple data sequences• BASIC: A probabilis%c mul%-sequencechangepoint model• Infer changepoint probabili%es & most likely
changepoint loca%ons from model posterior• Fit model priors using the data
• Applicable to con%nuous data (changes in mean,variance), counts, binary observa%ons, etc. byusing diﬀerent data distribu%on models

Thanks!
Preprint: arxiv.org/abs/1508.01280
SoXware: github.com/zfzf21/BASIC_changepoint