Multiscale quantile segmentation

TL;DR

This paper presents a multiscale quantile segmentation method for serial data that accurately detects change points in quantile functions with strong statistical guarantees, without relying on distributional assumptions beyond independence.

Contribution

It introduces a novel multiscale statistic for quantile segmentation that controls the probability of correct segment detection and achieves minimax optimal estimation rates.

Findings

Method reliably detects quantile changes in simulations.

Provides asymptotic confidence bands for quantile functions.

Effective in genetic sequencing and ion channel data analysis.

Abstract

We introduce a new methodology for analyzing serial data by quantile regression assuming that the underlying quantile function consists of constant segments. The procedure does not rely on any distributional assumption besides serial independence. It is based on a multiscale statistic, which allows to control the (finite sample) probability for selecting the correct number of segments S at a given error level, which serves as a tuning parameter. For a proper choice of this parameter, this tends exponentially fast to the true S, as sample size increases. We further show that the location and size of segments are estimated at minimax optimal rate (compared to a Gaussian setting) up to a log-factor. Thereby, our approach leads to (asymptotically) uniform confidence bands for the entire quantile regression function in a fully nonparametric setup. The procedure is efficiently implemented using dynamic programming techniques with double heap structures, and software is provided. Simulations and data examples from genetic sequencing and ion channel recordings confirm the robustness of the proposed procedure, which at the same hand reliably detects changes in quantiles from arbitrary distributions with precise statistical guarantees.

Figures

Tables

Table 1. Table C.1: MISE and V-measure of the two proposed estimators in equation ( 21 ) in the main text (Koenker) and Definition C.1 (Runs), respectively, for a simple bump function as described in the text, for additive Gaussian and Cauchy noise and quantiles β = 0.25 , 0.5 , 0.75 𝛽 0.25 0.5 0.75 \beta=0.25,0.5,0.75 . The outperforming method is marked bold in each case.

Setting	$\beta$	Estimator	MISE	V-m.
Normal	0.5	Runs	0.09	0.74
		Koenker	0.09	0.72
	0.25	Runs	0.43	0.69
		Koenker	0.42	0.68
	0.75	Runs	0.43	0.72
		Koenker	0.49	0.66
Cauchy	0.5	Runs	0.03	0.90
		Koenker	0.01	0.98
	0.25	Runs	0.05	0.85
		Koenker	0.02	0.95
	0.75	Runs	0.05	0.84
		Koenker	0.08	0.78

Table 2. Table E.1: Probability of estimating correctly the number of segments, frequency of confidence interval coverage of the true c.p’s (CI cov.), and frecuency of confidence band coverage of the true function (CB cov.), for data as in Figure G.2 with different levels α 𝛼 \alpha and n = 2485 𝑛 2485 n=2485 .

	$\beta$	$1-\alpha$	$\mathbb{P}\left(\hat{S}\le S\right)$	CI cov.	CB cov.		$\mathbb{P}\left(\hat{S}\le S\right)$	CI cov.	CB cov.
Normal ${\sigma }^2=0.04$	0.5	0.99	1.000	1.000	0.998	t (3 d.f.) ${\sigma }^2=0.04$	1.000	1.000	0.996
		0.95	1.000	1.000	0.977		1.000	1.000	0.974
		0.90	0.997	1.000	0.936		0.998	1.000	0.947
	0.25	0.99	1.000	1.000	1.000		1.000	1.000	1.000
		0.95	1.000	1.000	0.994		1.000	1.000	0.996
		0.90	0.999	1.000	0.951		0.998	1.000	0.948
	0.75	0.99	1.000	1.000	0.998		1.000	1.000	0.998
		0.95	0.999	1.000	0.962		0.999	1.000	0.972
		0.90	0.998	1.000	0.953		0.998	0.999	0.954
Cauchy	0.5	0.99	1.000	1.000	0.998	${\chi }_3^2$ ${\sigma }^2=0.04$	1.000	1.000	0.993
		0.95	0.999	1.000	0.965		0.999	1.000	0.966
		0.90	0.995	0.998	0.928		1.000	0.998	0.936
	0.25	0.99	1.000	1.000	1.000		1.000	1.000	1.000
		0.95	1.000	1.000	0.997		1.000	1.000	0.994
		0.90	1.000	1.000	0.943		0.998	1.000	0.955
	0.75	0.99	1.000	1.000	0.995		1.000	1.000	0.992
		0.95	1.000	1.000	0.971		1.000	1.000	0.977
		0.90	0.999	0.999	0.932		0.998	0.999	0.952

Table 3. Table G.1: Summary table of the characteristics of MQS versus SMUCE (Frick et al., 2014 ) , HSMUCE (Pein et al., 2017 ) , WBS (Fryzlewicz, 2014 ) , R-FPOP (Fearnhead and Rigaill, 2017 ) , QS (Eilers et al., 2005 ) , NOT (HT contrast) (Baranowski et al., 2019 ) , and NWBS (Padilla et al., 2019 ) . The computation time comparison is based on normal noise ( n = 1988 𝑛 1988 n=1988 ), mean as in Figure G.2 and variance 0.04. With ”robust to heterogeneity” we mean that the median/mean estimator is robust again changes in variance.

	MQS	SMUCE	HSMUCE	WBS	R-FPOP	QS	NOT	NWBS
Confidence statements	✓	✓	✓	$\times$	$\times$	$\times$	$\times$	$\times$
Consistency results	✓	✓	✓	✓	✓	$\times$	✓	✓
No distributional assumptions	✓	$\times$	$\times$	$\times$	✓	✓	$\times$	✓
Target a specific quantile	✓	$\times$	$\times$	$\times$	$(\times )$	✓	$\times$	($\times$)
Robust to outliers	✓	$\times$	✓	$\times$	✓	✓	✓	✓
Robust to heterogeneity	✓	$\times$	✓	✓	$\times$	✓	✓	$\times$
Computation time [s] (see caption)	0.65	0.01	0.02	0.07	0.01	38	0.35	63.8

Table 4. Table G.2: Frequencies of estimated number of segments in [ % ] [\%] , MIAE ( × 100 ) (\times 100) , and V-measure ( × 10 ) (\times 10) for data as in Figure G.2 . Here, the true number of segments equals 7. The proposed MQS estimator is compared with SMUCE (Frick et al., 2014 ) , HSMUCE (Pein et al., 2017 ) , WBS (Fryzlewicz, 2014 ) , R-FPOP (Fearnhead and Rigaill, 2017 ) , QS (Eilers et al., 2005 ) , NOT (Baranowski et al., 2019 ) , and NWBS (Padilla et al., 2019 ) .

	Method	$\le 5$	6	$7$	$8$	$\ge 9$	MIAE	V-m.		$\le 5$	6	$7$	$8$	$\ge 9$	MIAE	V-m.
Normal ${\sigma }^2=0.04$	MQS(0.5)	0.00	0.70	99.20	0.10	0.00	3.83	9.29	t (3 d.f.) ${\sigma }^2=0.04$	0.00	10.30	89.50	0.20	0.00	3.38	9.33
	SMUCE	0.00	0.00	99.70	0.30	0.00	0.90	9.96		0.00	0.00	0.00	0.00	100.00	2.41	7.72
	HSMUCE	0.00	0.00	98.40	1.60	0.00	0.91	9.95		0.00	0.00	99.50	0.50	0.00	0.91	9.96
	R-FPOP	0.00	0.00	96.90	3.10	0.00	0.92	9.95		0.00	0.00	95.30	4.20	0.50	0.72	9.97
	WBS	0.00	0.00	98.00	1.90	0.10	0.97	9.95		0.00	0.00	1.00	0.00	99.00	1.56	8.40
	QS(0.5)	0.00	0.00	0.10	0.10	99.80	7.50	3.54		0.00	0.00	0.10	0.10	99.80	7.49	3.54
	NOT(HT)	0.00	0.00	98.80	1.20	0.00	0.91	9.95		0.00	0.00	99.10	0.80	0.10	0.89	9.97
	NOT(VAR)	0.00	0.00	99.60	0.40	0.00	20.96	9.97		0.00	0.00	7.90	8.70	83.40	20.50	9.64
	NWBS	0.00	0.00	96.60	1.30	2.10	1.26	9.86		0.00	0.00	98.70	0.50	0.80	0.80	9.93
Cauchy	MQS(0.5)	0.00	17.50	82.30	0.20	0.00	2.73	9.34	${\chi }_3^2$ ${\sigma }^2=0.04$	0.00	0.00	99.70	0.30	0.00	3.53	9.34
	SMUCE	0.00	0.00	0.00	0.00	100.00	58.52	5.38		0.00	0.00	0.10	0.70	99.20	30.03	8.32
	HSMUCE	0.00	0.20	99.60	0.20	0.00	1.93	9.81		0.00	0.00	86.20	13.30	0.50	30.49	9.89
	R-FPOP	0.00	0.00	91.00	7.00	2.00	0.22	9.98		0.00	0.00	84.40	7.80	7.80	32.53	9.92
	WBS	0.00	0.00	0.00	0.00	100.00	81.71	6.94		0.00	0.00	26.20	0.30	73.50	30.30	9.42
	QS(0.5)	0.10	0.00	0.00	0.00	99.90	5.07	3.54		0.00	0.10	0.20	0.40	99.30	7.40	3.54
	NOT (HT)	56.20	1.80	41.60	0.20	0.20	56.25	6.62		0.00	0.00	99.40	0.30	0.30	30.29	9.96
	NOT (VAR)	0.00	0.00	0.00	0.10	99.90	121.06	8.95		0.00	0.00	61.20	15.20	23.60	29.26	9.89
	NWBS	0.00	0.10	99.10	0.40	0.40	0.19	9.97		0.00	0.00	98.20	0.90	0.90	1.03	9.93
Normal ${\sigma }^2=0.04$	MQS(0.25)	81.40	18.60	0.00	0.00	0.00	8.19	8.64	t (3 d.f.) ${\sigma }^2=0.04$	82.80	17.20	0.00	0.00	0.00	7.76	8.64
	MQS(0.75)	0.00	0.00	99.90	0.10	0.00	6.06	9.02		0.00	0.00	100.00	0.00	0.00	5.68	9.06
	QS(0.25)	0.30	0.10	1.80	3.60	94.20	8.75	3.55		0.00	0.10	1.00	2.30	96.60	7.58	3.55
	QS(0.75)	0.00	0.00	0.10	0.30	99.60	6.83	3.54		0.00	0.00	0.10	0.10	99.80	6.02	3.54
Cauchy	MQS(0.25)	82.80	17.20	0.00	0.00	0.00	7.76	8.64	${\chi }_3^2$ ${\sigma }^2=0.04$	82.50	17.50	0.00	0.00	0.00	7.36	8.74
	MQS(0.75)	0.00	0.00	100.00	0.00	0.00	5.68	9.06		0.00	0.50	99.50	0.00	0.00	6.84	8.94
	QS(0.25)	0.00	0.10	1.00	2.30	96.60	7.58	3.55		0.00	0.10	1.00	1.80	97.10	7.07	3.54
	QS(0.75)	0.00	0.00	0.10	0.10	99.80	6.02	3.54		0.00	0.00	0.00	0.30	99.70	8.23	3.54

Table 5. Table G.3: Frequencies of estimated number of segments (in percentage), MIAE ( × 100 ) (\times 100) and V-measure ( × 10 ) (\times 10) for simultaneous changes in mean and variance, for data as in Figure G.3 . The true number of segments for the median is 4 and for the other quantiles 6.

	Method	$\le 3$	4	5	6	$\ge 7$	MIAE	V-m.		$\le 3$	4	5	6	$\ge 7$	MIAE	V-m.
Normal	MQS(0.5)	16.08	83.11	0.81	0.00	0.00	6.13	8.67	t (3 d.f.)	0.00	99.39	0.61	0.00	0.00	3.92	8.88
	SMUCE	0.30	33.06	37.41	19.82	9.40	1.97	9.77		0.00	0.00	0.00	0.00	100.00	13.38	5.41
	HSMUCE	24.57	73.61	1.82	0.00	0.00	1.96	9.97		13.09	85.99	0.92	0.00	0.00	2.19	9.95
	R-FPOP	0.10	65.52	8.70	14.36	11.32	1.31	9.77		0.00	80.16	7.67	8.79	3.37	1.80	8.45
	WBS	0.61	90.60	1.62	6.17	1.01	2.02	9.77		0.00	0.00	0.00	0.72	99.28	2.72	7.75
	QS(0.5)	0.00	0.00	0.10	0.51	99.39	9.93	2.55		0.00	0.00	0.00	0.51	99.49	6.49	2.55
	NOT (HT)	0.61	97.37	0.81	1.11	0.10	2.05	9.63		1.53	97.34	0.10	0.61	0.41	2.19	9.95
	NOT (VAR)	0.00	0.00	0.00	99.70	0.30	28.40	9.31		0.00	0.00	0.00	8.18	91.82	29.20	9.11
	NWBS	2.83	29.52	36.50	20.32	10.82	3.04	7.58		0.20	64.93	16.36	10.12	8.38	1.76	7.83
Normal	MQS(0.25)	0.51	65.82	33.47	0.20	0.00	11.98	7.33	t (3 d.f.)	2.04	57.98	35.79	4.19	0.00	12.78	6.81
	MQS(0.75)	0.00	0.30	57.53	42.16	0.00	19.19	6.95		0.00	2.76	70.25	26.99	0.00	11.84	7.33
	QS(0.25)	0.20	0.71	2.63	7.28	89.18	13.82	3.64		0.31	3.07	6.54	9.10	80.98	11.42	3.64
	QS(0.75)	0.00	0.00	0.00	0.00	100.00	10.98	3.66		0.00	0.00	0.00	0.00	100.00	9.73	3.64

Equations

$$\begin{gathered} \Sigma =\{\vartheta =\sum _{s=1}^S{\theta }_s\text{mathds}{1}_{[{\tau }_{s-1},{\tau }_s)}:{\theta }_s\ne {\theta }_{s+1}, \\ 0={\tau }_0<{\tau }_1<...<{\tau }_S=1,S<\infty \}. \end{gathered}$$

$${\vartheta }_{\beta }(x_i)=\inf \{\theta :\text{P}(Z_i\le \theta )=\beta \}\text{mbox}fori=1,\dots ,n.$$

$$W_i=W_i(Z_i,\vartheta (x_i)):=\{\begin{array}{ll}0& \text{mbox}if\vartheta (x_i)-Z_i<0,\\ 1& \text{mbox}if\vartheta (x_i)-Z_i\ge 0.\end{array}$$

$$H_{ij}:W_i,\dots ,W_j\stackrel{i.i.d.}{\sim }B(\beta )\text{vs.}{K}_{ij}:W_i,\dots ,W_j\stackrel{i.i.d.}{\sim }B({\beta }^{\prime })\text{ with }{\beta }^{\prime }\ne \beta ,$$

$$\hat{S}=\hat{S}(q):=\underset{\vartheta \in \Sigma }{\inf }\#J(\vartheta )\text{ s.t }{T}_n(W(Z,\vartheta ),\vartheta )\le q,$$

$$\mathcal{H}(q):=\{\vartheta \in \Sigma :\#J(\vartheta )=\hat{S}(q)\text{ and }{T}_n(W(Z,\vartheta ),\vartheta )\le q\}.$$

$$M_n:=\underset{1\le i\le j\le n}{\max }\sqrt{2T_i^j(X,\beta )}-P_{\ell ,n},$$

$$\mathbb{P}\left({\hat{S}}_{n,\alpha }>S\right)\le \alpha ,$$

$$\mathbb{P}\left({\hat{S}}_{n,\alpha }>S+s\right)\le {\alpha }^{\lfloor s/2\rfloor +1},$$

$${\xi }_{F,\beta }(\delta )=\left|F({\theta }_{\beta }+\delta )-\beta \right|.$$

$$\Xi :=\underset{s=1,\dots ,S}{\min }\min \left\{{\xi }_{F_{\min },\beta }({\theta }_s-{\theta }_{s-1}),{\xi }_{F_{\max },\beta }({\theta }_{s-1}-{\theta }_s)\right\}\text{mbox}and\Lambda :=\underset{s=1,\dots ,S}{\min }\left|{\tau }_s-{\tau }_{s-1}\right|.$$

$$\mathbb{P}\left({\hat{S}}_{n,\alpha }<S\right)\le 4(S-1){\text{e}}^{-n\Lambda {\Xi }^2}\left[{\text{e}}^{2\sqrt{n\Lambda }\Xi \left(q_n(\alpha )/\sqrt{2}+\sqrt{\log (\text{nicefrac}2\text{e}\Lambda )}\right)}+1\right].$$

$$\mathbb{P}\left({\hat{S}}_{n,\alpha }=S\right)\ge 1-\alpha (q_n)-4(S-1){\text{e}}^{-n\Lambda {\Xi }^2}\left[{\text{e}}^{2\sqrt{n\Lambda }\Xi \left(q_n/\sqrt{2}+\sqrt{\log (\text{nicefrac}2\text{e}\Lambda )}\right)}+1\right],$$

$$\underset{\begin{array}{c}\vartheta \in \Sigma \\ \Lambda >{\Lambda }_n,\Xi >{\Xi }_0\end{array}}{\sup }\mathbb{P}\left(\underset{\tau \in J(\vartheta )}{\max }\underset{\hat{\tau }\in J(\hat{\vartheta }(q))}{\min }|\hat{\tau }-\tau |>\frac{\log (n)}{n}\frac{1}{{\Xi }^2}\right)\to 0,$$

$${\vartheta }_{\beta }(x_i)=\inf \{\theta :\text{P}(Z_i\le \theta )\ge \beta \}.$$

$$\begin{array}{r}{H}_{ij}:\mathbb{P}(Z_k\le {\theta }_s)=\beta \forall i\le k\le j\text{vs.}{K}_{ij}:\mathbb{P}(Z_k\le {\theta }_s)\ne \beta \forall i\le k\le j.\end{array}$$

$$\begin{array}{rl}{T}_i^j(W)& =T_i^j\left(W\left(Z,{\theta }_s\right)\right)=\log \left(\frac{{\sup }_{{\beta }^{\prime }}{\prod }_{l=i}^j{f}_{{\beta }^{\prime }}(W_l)}{{\prod }_{l=i}^j{f}_{\beta }(W_l)}\right)\\ & =(j-i+1)\left({\overline{W}}_i^j\log \left(\frac{{\overline{W}}_i^j}{\beta }\right)+\left(1-{\overline{W}}_i^j\right)\log \left(\frac{1-{\overline{W}}_i^j}{1-\beta }\right)\right),\end{array}$$

$${\Phi }_{ij}\left(Z\right)=\{\begin{array}{ll}0& \text{mbox}if\sqrt{2T_i^j(W(Z,{\theta }_s))}\le q_{i,j},\\ 1& \text{mbox}otherwise,\end{array}$$

$$T_n(Z,\bullet ):=\underset{\begin{array}{c}1\le i\le j\le n\\ \bullet \text{mbox}isconstanton[x_i,x_j]\end{array}}{\max }\sqrt{2T_i^j(W(Z,\bullet ),\beta )}-P_{\ell ,n}.$$

$$q_{ij}(\alpha )=q_n(\alpha )+P_{\ell ,n}.$$

$$\underset{\begin{array}{c}1\le i\le j\le n\\ {\vartheta }_{\beta }\text{mbox}isconstanton[x_i,x_j]\end{array}}{\max }\sqrt{2T_i^j(X,\beta )}-P_{\ell ,n},$$

$$q_n(\alpha ):=\inf \{q:\mathbb{P}(M_n\le q)\ge 1-\alpha \}.$$

$$\hat{\vartheta }:={\mathrm{argmin}}_{\hat{\vartheta }\in \mathcal{H}(q)}\sum _{i=1}^n\left(Z_i-\hat{\vartheta }(x_i)\right)\left(\beta -\text{mathbbm}{1}_{\{Z_i<\hat{\vartheta }(x_i)\}}\right).$$

$${\xi }_s:=\min \left\{{\xi }_{F_{\min },\beta }(-{\delta }_s),{\xi }_{F_{\max },\beta }({\delta }_s)\right\}s=1,\dots ,S,$$

$${\gamma }_{n,s}(q)={\left(1-2\exp \left(-\frac{{\left(\sqrt{2n{\lambda }_s}{\xi }_s-q-\sqrt{2\log (2\text{e}/{\lambda }_s)}\right)}_{+}^2}{2}\right)-2\exp \left(-n{\lambda }_s{\xi }_s^2\right)\right)}^2.$$

$$\mathbb{P}\left(\hat{S}(q)\ge S\right)\ge \prod _{s=1}^S{\gamma }_{n,s}(q).$$

$$\underset{\begin{array}{c}\vartheta \in \Sigma :\\ \Lambda >{\Lambda }_0,\Xi >{\Xi }_0\end{array}}{\inf }\mathbb{P}\left({\vartheta }_{\beta }\in \mathcal{H}(q_n(\alpha ))\right)$$

$$\ge 1-\alpha -\mathbb{P}\left(\hat{S}(q)<S\right)\ge 1-\alpha +o(1)$$

$$\underset{n\to \infty }{\lim }\underset{\begin{array}{c}\vartheta \in \Sigma :\\ \Lambda >{\Lambda }_0,\Xi >{\Xi }_0\end{array}}{\inf }\mathbb{P}(\vartheta \in \mathcal{H}(q_n(\alpha )))\ge 1-\alpha .$$

$$\sqrt{n{\Lambda }_n}{\Xi }_n>2\sqrt{-\log ({\Lambda }_n)}.$$

$$\mathbb{P}\left(\underset{\tau \in J(\vartheta )}{\max }\underset{\hat{\tau }\in J(\hat{\vartheta }(q))}{\min }|\hat{\tau }-\tau |>{ϵ}_n\right)\le 4(S-1){\text{e}}^{-n{ϵ}_n{\Xi }^2}\left[{\text{e}}^{2\sqrt{n{ϵ}_n}\Xi \left(\text{nicefrac}q\sqrt{2}+\sqrt{\log (2\text{e}/{ϵ}_n)}\right)}+1\right].$$

Full text

Keywords: Change-points, Double heap, Dynamic programming, Multiscale methods, Quantile regression, Robust segmentation.

2010 Mathematics Subject Classification: 62G08, 62G15, 62G30, 62G35, 90C39.

1 Introduction

The analysis of serial data with presumably abrupt underlying distributional changes, for example, in its mean, median, or variance, is a long-standing issue and relevant to a magnitude of applications, e.g. to econometrics and empirical finance (Preuss et al., 2015; Shen, 2016; Russell and Rambaccussing, 2019), evolutionary and cancer genetics (Liu et al., 2013; Zhang and Siegmund, 2007; Jónás et al., 2016) or neuroscience (Cribben and Yu, 2017), to mention a few (see Section 1.4 for a more comprehensive discussion). The present work proposes a new methodology for this task, denoted as Multiscale Quantile Segmentation (MQS) which on the one hand, is extremely robust as it does not rely on any distributional assumption (apart from independence) and on the other hand, still has high detection power with (even non-asymptotic) statistical guarantees. More precisely, our methodology is based on quantile segments, i.e., quantile regression functions which are modeled as right-continuous piecewise constant functions $\vartheta$ with finitely but arbitrary many (unknown) segments $S$. We stress that even in such situations where the quantiles are not piecewise constant, this may serve as a reasonable proxy to cartoonize the quantile function in a simple but meaningful way, resulting in a “quantilogram” similar in spirit to Tukey’s regressogram (Tukey, 1961). To fix our setting, we assume that the underlying regressor (e.g., time) is in the interval $[0,1)$ and is sampled equidistantly at $n$ sampling points $x_{i}\mathrel{\mathop{:}\!\!=}(i-1)/n$. As our main results are nonasymptotic, extensions to general (ordered) sampling domains and non-equidistant sampling points are immediate. All such segment functions $\vartheta:[0,1)\rightarrow\mathds{R}$ are then comprised in the space

[TABLE]

Hence, our quantile function consist of $S$ (unknown) distinct segments with unknown segment values $\theta_{s}\in\mathds{R}$ and segment lengths $\tau_{s}-\tau_{s-1}$ (see Figure 1.1 for illustration).

QSR-model (Quantile Segment Regression model) Let $\beta\in(0,1)$ and $\vartheta_{\beta}\in\Sigma$ a segment function. In the QSR-model one observes $n$ independent random variables $Z_{i}$ at equidistant sampling points $x_{i,n}=x_{i}=(i-1)/n$ for $i=1,\ldots n$ such that its $\beta$-quantiles are given as

[TABLE]

Note that for a particular observed signal $Z_{i}$, $i=1,\ldots,n$, different quantiles $\beta\neq\beta^{\prime}$ will, in general, correspond to different segment functions $\vartheta_{\beta}\neq\vartheta_{\beta^{\prime}}\in\Sigma$ in (1). Here we focus on finding the corresponding segment function for a particular, but arbitrary, given quantile $\beta$, e.g., the median ($\beta=0.5$). Screening for several $\beta$-values then will allow to specify in which part of the distribution significant changes occur. This is in contrast to non-parametric distributional segmentation (see Section 1.4 for references), where one would aim to find changes in the distribution per se, i.e. a segment change is detected when (at least) some of the quantile functions changes. Therefore, we do not consider $\beta$ as a tuning parameter, but rather as a user-specific input, that depends on the particular task of interest.

Example 1.1**.**

Figure 1.1 shows a data example from an ion channel recording experiment performed by the Steinem lab (Insitute of Organic and Biomolecular Chemistry, University of Göttingen) (see more details in Section 5.2). Using the patch-clamp technique (Sakmann and Neher, 1995), one can measure the current flow of ions transported through an individual channel in a cell’s membrane. Figure 1.1 shows data from the bacterial porin PorB, an outer membrane porin of Neisseria meningitidis – a pathogenic bacteria well known for being the agent of epidemic meningitis (Virji, 2009). Over time, the ion channel changes its gating behavior by closing and reopening its pore. This leads to a piecewise constant current flow structure, shown as black dots in Figure 1.1. The opening and closing behavior are of immediate physiological interest, and we can consider this within the QSR-model, where in this data example $n=2,801$. By targeting different quantiles, one can recover different properties of the channel. For example, the median $\beta=0.5$ provides a measure of average current flow and can thus quantify the channel’s overall configuration. The second row in Figure 1.1 shows in red our MQS segmentation for $\beta=0.5$ and $\alpha=0.01$, with fully automatic estimated number of segments $\hat{S}=10$ (recall (1)).

In addition, upper and lower tail quantiles provide more information on the variability of the channel behavior. The third and fourth rows in Figure 1.1 show in light red our MQS segmentation for $\beta=0.25$ and $\beta=0.75$, respectively. Moreover, the first row in Figure 1.1 shows together the MQS’ estimates for the 0.25, 0.5, and 0.75-quantiles. We refer to this visualization of the three quantiles as the multiscale segment boxplot (MSB). This shows that the variance on the lower state is lower than on the higher state, a common characteristic of many channels, known as open channel noise, see Sakmann and Neher (1995), which is well known to aggravate data analysis. By targeting different quantiles, MQS provides a direct approach for its quantification. Besides estimates for the respective quantile segmentation, our procedure also comes with a precise uncertainty quantification. This is illustrated with confidence intervals for the segment locations (blue lines) and confidence bands for the segment function (grey areas). We stress that these confidence statements do not rely on any distributional assumptions for the data generating process (besides independence), which is in contrast to other available approaches for this purpose, see Section 1.4. This makes the MQS procedure particularly useful in practice, where parametric distributional assumptions are often problematic. For example, although a Gaussian assumption is widely used for ion channel data, see e.g., Pein et al. (2017), performing a Shapiro-Wilk’s (Shapiro and Wilk, 1965) test on the first segment of the present ion channel data rejects with p-value = 0.0016.

As illustrated in Example 1.1, the aim of this work is to provide a statistical methodology for multiscale quantile segmentation (MQS) in the general QSR-model. In particular, we stress that MQS is robust to arbitrary distributional changes that keep the respective quantile under consideration unchanged, including changes of variance, as demonstrated in Example 1.1. See also Appendix B, which further illustrates this robustness with some synthetic simulated data example. For each $\beta\in(0,1)$, MQS provides estimates and confidence statements for:

the number $S$ of segments $S(\vartheta_{\beta})$, 2. 2.

the segment locations $\tau_{1},\ldots,\tau_{S-1}\in[0,1)$, 3. 3.

and, based on these, the segment values, i.e., the quantiles $\theta_{1},\ldots,\theta_{S}\in\mathds{R}$.

Our approach is based on a simple transformation: Given a candidate segment function $\vartheta\in\Sigma$ (which will depend on data $Z_{i}$), we consider the transformed binary (pseudo) data

[TABLE]

Note that if and only if the candidate function $\vartheta$ equals the true underlying regression function $\vartheta_{\beta}$ in the QSR-model, then $W_{1},\ldots,W_{n}$ are i.i.d. Bernoullis with success probability $\beta$. Based on this, our methodology “tests” in a multiscale fashion any possible candidate function in $\Sigma$ to be valid: it selects a segment quantile function which does not contradict the i.i.d. Bernoulli assumption and, among those, has the smallest number of segments.

To this end, the unknown number and locations of segments of $\vartheta_{\beta}$ will be detected by a certain multiscale statistic $T_{n}(Z,\vartheta)=T_{n}(W(Z,\vartheta),\vartheta)$ (see (18)), which combines on intervals $[x_{i},x_{j}]$ where $\vartheta|_{[x_{i},x_{j}]}$ is constant, with $1\leq i\leq j\leq n$, the corresponding log-likelihood ratio tests for the hypothesis testing problems

[TABLE]

where $B(\beta)$ denotes a Bernoulli distribution with success probability $\beta$ (see (18) for the precise definition of $T_{n}$). In a first step MQS determines the number of segments from the data given such $\beta\in(0,1)$. For a given threshold $q\in\mathds{R}$ only depending on $\beta$, $n$, and the nominal level $\alpha$ (to be specified later), the estimated number of segments $\hat{S}$ is then the solution of a (non-convex) optimization problem with convex constrains given by the multiscale statistic $T_{n}$, namely,

[TABLE]

where $J(\vartheta)$ denotes the set of segments and $\#J(\vartheta)$ the number of segments of $\vartheta\in\Sigma$. That is, we choose the smallest number of segments such that the multiscale test $T_{n}(W(Z,\vartheta),\vartheta)\leq q$ still accepts. As, in a certain sense, this is a model selection step, we call $\hat{S}$ the MQS selector of S. In a second step, we then constrain all candidate segment quantile functions to those with $\hat{S}$ segments, that is, to the set

[TABLE]

The MQS estimate (MQSE) $\hat{\vartheta}$ is then a particular simple segment function in $\mathcal{H}(q)$ based on the Wald-Wolfowitz runs tests (see Section 2 for details and Figure 1.1 for illustration).

1.1 Main results: Theory

For the MQSE $\hat{\vartheta}$, with $\vartheta_{\beta}$ being the true underlying regression function in the QSR-model, we immediately get from (5) that $\mathbb{P}\left(\hat{S}(q)>S\right)\leq\mathbb{P}\left(T_{n}\left(Z,\vartheta_{\beta}\right)\leq q\right).$ In Section 2.1 we show that $T_{n}(Z,\vartheta_{\beta})$ can be bounded in distribution with a random variable $M_{n}=M_{n,\beta}$ that only depends on $\beta$ and $n$ and hence does not depend on $\vartheta_{\beta}$ or any other characteristics of the underlying distribution of the observations $Z_{i}$, namely

[TABLE]

where $X=(X_{1},\ldots,X_{n})$ are i.i.d. Bernoulli distributed with mean $\beta$, and a penalization $P_{\ell,n}$ depending only on $n$ and $\ell=j-i+1$, see (18). Therefore, for given $n$ and $\beta$ the finite sample-quantiles of $M_{n}$ can be computed by Monte-Carlo simulations in a universal manner. Hence, choosing $q_{n}=q_{n}(\alpha)$ as the $(1-\alpha)$-quantile of $M_{n}$ gives us control for the overestimation error of the number of segments included in our final estimator by a desired level $\alpha$, namely

[TABLE]

where $\hat{S}_{n,\alpha}\mathrel{\mathop{:}\!\!=}\hat{S}\left(q_{n}(\alpha)\right)$. In Theorem 1.1 we present a refinement of this fact. Note that $q_{n}(\alpha)$ depends on $\beta$, however, asymptotically this dependency vanishes, see Remark 2.2.

Theorem 1.1** (Overestimation error).**

Consider the QSR-model. For $q=q_{n}(\alpha)$ as in (5) and (20) the MQS-selector $\hat{S}_{n,\alpha}:=\hat{S}(q)$ satisfies uniformly over all segment functions $\vartheta\in\Sigma$

[TABLE]

where $\lfloor\cdot\rfloor:=\max\{m\in\mathbb{Z}:m\leq x\}$, for all $x\in\mathbb{R}$, i.e. the floor function.

This means, in addition to controlling the overall error to include at least one extra segment ($s=0$), the error of overestimating the number of segments by more than $s$ decays exponentially fast. This reveals that segments detected by MQS are, with very high probability, indeed present in the signal. Therefore, (9) guarantees that $\alpha$ controls the false positives in a strong family-wise error sense. The overestimation bound in Theorem 1.1 is complemented by an explicit bound for the underestimation error, i.e. to miss a segment, see Theorem 1.2. Together, this allows for precise fine-tuning of the error of a wrong number of detected segments via the choice of the error level $\alpha$. Clearly, such an underestimation bound for $\mathbb{P}(\hat{S}_{n,\alpha}<S)$ has to depend on some characteristics of the function $\vartheta_{\beta}$, specifically on the length and height of the jumps, as no method can detect arbitrary small changes for a fixed number of data. Moreover, note that even a large jump of $\vartheta_{\beta}$ is not identifiable from the QSR-model, if it does not induce a sufficiently large jump in the respective distribution functions, as the following example shows.

Example 1.2**.**

For sufficiently large $L>0$ and small $\epsilon>0$ consider random variables $X,Y$ such that $\textbf{P}(X=L)=0.5$, $\textbf{P}(X=-L)=0.5-\epsilon$, $\textbf{P}(X=-L+\epsilon)=\epsilon$ and $\textbf{P}(Y=-L)=0.5-\epsilon$, $\textbf{P}(Y=L)=0.5$, $\textbf{P}(Y=L-\epsilon)=\epsilon$. Then, if half of the sample comes from $X$, i.e., $Z_{1},\ldots,Z_{n/2}\;\sim\;X$ and the other half from $Y$, i.e., $Z_{n/2+1},\ldots,Z_{n}\;\sim\;Y$ as in the QSR-model, the underlying median regression function $\vartheta_{0.5}$ has a large jump of size $2L-2\epsilon$ at $\tau=0.5$. However, for sufficiently small $\epsilon$ (depending on $n$), because the probability of sampling $-L+\epsilon$ or $L-\epsilon$ is very small, this jump is not detectable from the observations $Z_{1},\ldots,Z_{n}$. See Figure 1.2 for illustration. The following definition specifies situations where a quantile segment is detectable.

Definition 1.1**.**

For a distribution function $F$ and $\beta\in(0,1)$, let $\theta_{\beta}\mathrel{\mathop{:}\!\!=}\inf\{\theta:F(\theta)\geq\beta\}$ be the $\beta$-quantile. Then the quantile jump function $\xi_{F,\beta}\colon\mathbb{R}\rightarrow[0,1]$ is defined as

[TABLE]

The quantile jump function quantifies how much a quantile jump in the QSR-model influences data around a jump. The bound for the underestimation error for the QSR-Model naturally depends on the minimal length of constant segments and on the minimal quantile jump for the distributions $F_{1},\ldots,F_{n}$ of $Z_{1},\ldots,Z_{n}$, that is,

[TABLE]

where $F_{\text{min}}\mathrel{\mathop{:}\!\!=}\min\{F_{1},\ldots,F_{n}\}$ and $F_{\text{max}}\mathrel{\mathop{:}\!\!=}\max\{F_{1},\ldots,F_{n}\}$ (pointwise). In Theorem 2.1 we show a general exponential bound for underestimating the number of segments by the MQSE, which provides the following result.

Theorem 1.2**.**

Consider the QSR-model, $\hat{S}_{n,\alpha}$ as in (8), and $\Lambda,\Xi$ as in (10). Then, for $q_{n}(\alpha)>0$ the $1-\alpha$ quantile of $M_{n}$ in (7), it holds that

[TABLE]

Note that whenever $\beta\to 0$ or $\beta\to 1$ then $\Xi\to 0$ and hence, the bound on the r.h.s. in (11) becomes trivial. This reflects the fact that for very high or very small quantiles it is arbitrarily difficult to capture changes from finitely many samples. Combining (8) and (11), we obtain an explicit bound for estimating the number of segments correctly, depending on $\alpha(q_{n})=\mathbb{P}(M_{n}>q_{n})$, $n$, $\Lambda$ and $\Xi$ in (10), namely,

[TABLE]

uniformly over all possible segment functions in $\Sigma$ satisfying (10). As $\alpha(q_{n})$ converges to zero as $q_{n}\to\infty$ (see Section 2), the MQS selector is consistent, that is $\mathbb{P}\left(\hat{S}_{n,\alpha}=S\right)\rightarrow 1$ as $n\rightarrow\infty$, whenever the threshold parameters are chosen such that $q_{n}\to\infty$ and $q_{n}/\sqrt{n}\rightarrow 0$. Moreover, $q_{n}$ can be chosen such that the r.h.s of (12) is maximized leading to exponentially fast selection consistency. In Theorem 2.3 we refine these consistency results to the situation where $\Xi=\Xi_{n}$ and $\Lambda=\Lambda_{n}$ vanish as $n\to\infty$. We give a sharp condition on $\Xi_{n},\Lambda_{n},q_{n}$ under which MQS consistently estimates the number of segments $S$ and show that (up to constants) these conditions cannot be improved, in general.

With the choice of $q_{n}=q_{n}(\alpha)$ as in (8), the multiscale approach of MQS directly yields (asymptotically) honest simultaneous confidence bands (i.e. uniformly over all possible segment functions in $\Sigma$ with minimal scale $\Lambda$ and minimal distribution jump $\Xi$) for the underlying regression function $\vartheta_{\beta}$ and simultaneous confidence intervals for the change-points $\tau$ via the set $\mathcal{H}(q_{n})$ in (6), see last three plots of Figure 1.1 and Theorem 2.2 for further details. Besides model selection consistency and confidence statements for all quantities, the MQS procedure also attains minimax optimal estimation rates (up to a log-factor) for the change-point locations. In Theorem 2.4 we show that whenever $q_{n}=o(\sqrt{\log(n)})$ and $\Lambda_{n}^{-1}=o(n/\log(n))$ then for any $\Xi_{0}>0$

[TABLE]

where the minimax rate is lower bounded by the sampling rate $1/n$ and hence, the rate in (13) is optimal (up to the log-factor).

1.2 Implementation

It has been exploited for a long time that global optimization procedures to detect segment changes in parametric models (mainly Gaussian), can be computed efficiently using dynamic programming, see e.g., (Bellman, 1954; Friedrich et al., 2008; Boysen et al., 2009; Killick et al., 2012; Davies et al., 2012; Frick et al., 2014; Zou et al., 2014; Pein et al., 2017; Haynes et al., 2017; Celisse et al., 2018; Wang et al., 2018). However, the exact computation of the nonparametric quantile segments of the MQSE in (5) leads to an additional computational burden compared to the case where the underlying jump signal corresponds to a parameter of a specific data distribution. Whereas the parametric case typically involves updating a running mean, which can be done in $\mathcal{O}(1)$ time, the quantile case involves updating the empirical quantile, which depends on the ordering of the data and hence, is computationally more involved. We incorporate an efficient running quantile computation as described in (Astola and Campbell, 1989), which is based on double heap structures. Thereby, the worst case computation time of MQSE only increases by a log-factor compared to the parametric case, as e.g., discussed in (Frick et al., 2014), being of the order $\mathcal{O}(n^{2}\log n)$. However, in many situations, depending on the underlying signal, the actual complexity will be almost linear. We give more details on the implementation of MQSE in Section 3. An R package mqs is available at https://github.com/ljvanegas/mqs.

1.3 Simulation results and data examples

In Section 4 we explore MQS in a comprehensive simulation study, including a comparison with several other state of the art segmentation methods, which have been designed to be robust, namely R-FPOP (Fearnhead and Rigaill, 2017), WBS (Fryzlewicz, 2014), HSMUCE (Pein et al., 2017), QS (Eilers et al., 2005), NWBS (Padilla et al., 2019), and NOT (Baranowski et al., 2019). NOT has two contrast functions for different types of robustness, namely ”pcwsConstMeanHT” (abbreviated HT) and ”pcwsConstMeanVar” (abbreviated VAR) for heavy-tailed noise and changes in variance, respectively. In order to compare the detection power of MQSE in a benchmark scenario, we also compare with SMUCE (Frick et al., 2014) which is tailored to i.i.d. normal error. A major finding is that while other methods only work well in some specific cases, MQS reliably detects changes in the quantiles for arbitrary distributions.

Figure 1.3 illustrates a benchmark scenario which shows data from different distributional regimes together with the MQS, SMUCE, R-FPOP, WBS, HSMUCE, QS, NOT, and NWBS estimates (red lines) for the median. In order to show robustness of the MQSE to changes in distribution we choose data from 3 different distributional regimes. The first 350 observations are drawn from a normal distribution with variance 2, the next 1190 from a t distribution with 1 d.f. and variance 0.1 and the last 945 from a $\chi^{2}$ distribution with 1 d.f. and variance 0.1. MQS turns out to have good detection power (see second regime), to be robust to heavy tails and change in variance (see first and second regime), and more generally to arbitrary distributional variations and skewness (see third regime). In contrast, SMUCE (specifically designed for Gaussian, homogeneous noise) is very sensitive to changes in variance and heavy-tailed distributions, adding a lot of artificial changes. WBS and NOT(VAR) are more robust to high variances, but also add artificial changes in the presence of outliers. R-FPOP is robust to heavy-tailed distributions, but is susceptible to changes in variance, adding artificial changes in the first 350 observations. On the contrary, HSMUCE, QS, and NOT(HT) suffer from oversmoothing and miss important data features. In this case, MQSE is capable of exploring the concentration of points around the median from the heavy-tailed $t$ (1 d.f.) distribution, despite low signal to noise ratio. NWBS is very robust to different noise families, but adds change-points in every change of distribution function, in this case in the first change from the normal to the t distribution. Moreover, we found NWBS (and QS) to be several magnitudes slower than the other methods, for instance, for this data example NWBS took $64s$ and MQS $0.64s$.

Recall that MQS does not just provide some reconstruction of the underlying $\beta$-quantile, but it also comes with precise statistical guarantees, such as honest confidence statements. This is particularly valuable for many real data examples, where there is uncertainty about the precise observational distribution. We demonstrate this in Figure 1.4 (see also Figure G.1) which shows an example from cancer genetics Behr et al. (2018), where one aims to detect copy-number aberrations (changes in the number of copies of certain regions of the genome) in tumor DNA (see Section 5 for details and a further data set on ion channel recordings). While Figure G.1 shows that most of the other methods include artificial jumps in their reconstruction (which in this data example are known to be not present in the underlying signal and result from sequencing artifacts, see Section 5.1 for details), MQS reliably recovers most of the copy-number aberrations correctly. Moreover, copy-number aberrations are known to often omit heterogeneity on different segments, and MQS adds a powerful visualization tool for this via the multiscale box plot, see Figure 1.4.

1.4 Related work

This work extends previous methodology developed in (Frick et al., 2014) on change-point regression under a parametric model (e.g., for a Gaussian mean signal) to a semi-parametric model. A major and novel feature of our approach is that it comes with statistical guarantees without making any such parametric assumptions. In particular, our confidence statements for all quantities and minimax optimal estimation rates, expand those offered in Frick et al. (2014) for parametric models to quantile segmentation. The technical challenges in such a semi-parametric model are tackled by a novel exponential inequality (see (11)) based on the quantile jump function (see Definition 1.1), which is fundamental to our theory. Algorithmically, we are able to exploit a double heap structures for the optimization problem in ((5)), which overall only leads to a loss of a $\log(n)$ factor in runtime compared to the parametric situation in (Frick et al., 2014). The transformation (3) extends the use of residual signs of Dümbgen (1998) and Davies and Kovac (2001) and was already mentioned in (Frick et al., 2014), however, without any theoretical analysis or efficient implementation.

The estimation of step functions with unknown number and location of segments are widely discussed problems; we mention (Fearnhead, 2006; Spokoiny, 2009; Boysen et al., 2009; Harchaoui and Lévy-Leduc, 2010; Jeng et al., 2010; Killick et al., 2012; Niu and Zhang, 2012; Siegmund, 2013; Matteson and James, 2014; Du et al., 2016; Gao et al., 2019; Fryzlewicz, 2018) for a variety of methods and related statistical theory. Methods for segment regression problems that offer statistical guarantees in general restrict those to (minimax optimal) estimation rates which typically concern the mean function and assume normality, see e.g., (Harchaoui and Lévy-Leduc, 2010; Cai and Xiong, 2012; Fryzlewicz, 2014; Li et al., 2016; Pein et al., 2017; Baranowski et al., 2019). In contrast, we target the quantile function (see e.g., Koenker (2005) for a survey) and do not assume a specific distributional model. One scenario related to our work is nonparametric distribution segmentation, see e.g. (Zou et al., 2014; Chu and Chen, 2019; Padilla et al., 2019). In this case, the aim is to pick up any distributional changes, not only in a particular quantile. Conceptually, mostly related to our work are methods which explicitly aim to provide robust methodology for segmentation. From a general perspective, segment regression may be considered as a particular case of a high dimensional linear model, with a very specific design matrix. Belloni and Chernozhukov (2011) consider sparse high dimensional quantile regression and show that minimizing the asymmetric absolute deviation loss together with $L_{1}$ penalization yields almost optimal estimation rates (in $L_{1}$ loss). However, their results require certain regularity conditions on the design matrix (similar to restricted eigenvalue conditions) which are not fulfilled for the specific design which corresponds to a multiscale segmentation setting. Other methods are concerned with particular types of robustness. For example Fearnhead and Rigaill (2017); Baranowski et al. (2019) consider the case of heavy-tailed symmetric distributions. Eilers et al. (2005); Li and Zhu (2007) consider minimizing absolute loss with an $L_{1}$-penalty for quantile regression for the segment regression setting, while not providing any theoretical results. More generally, Dümbgen and Kovac (2009) consider quantile regression via minimizing a convex loss with a total variation penalty and provide certain consistency results. Finally, Aue et al. (2014, 2017) consider quantile segmentation in a time series context based on a minimum description length criterion to select the number of segments and show consistency of their method, whereas in our (simpler) setting we obtain exponentially fast selection consistency, see (12).

2 Theory

In Section 2.1 we introduce the MQS methodology in more detail and in Section 2.2 we present our results on model selection consistency, confidence statements, and minimax optimal estimation rates for MQS. Proofs are postponed to Appendix A. We start with a discussion of the model assumptions.

At first glance, it might seem restrictive that the QSR-model requires for the underlying $\beta$-quantile regression function $\vartheta_{\beta}$ that $\textbf{P}(Z_{i}\leq\vartheta_{\beta}(x_{i}))=\beta$ (and not, more generally, that $\vartheta_{\beta}(x_{i})=\inf\{\theta:\textbf{P}(Z_{i}\leq\theta)\geq\beta\}$). However, such assumption is unavoidable for efficient (i.e., with non-trivial power) control of overestimation for the number of segments, as it is achieved by the proposed method, without assuming a specific observational distribution. To see this, note that without making any distributional assumptions, for given independent observations $Z_{1},\ldots,Z_{n}$ and candidate quantiles $\vartheta(x_{1}),\ldots,\vartheta(x_{n})$, all information from $Z_{i}$ about the candidate $\vartheta$ is captured in the transformation $W_{i}=W_{i}(Z_{i},\vartheta(x_{i}))$, $i=1,\ldots,n$ in (3). Now assume that in the QSR-model the condition in (2) is replaced by

[TABLE]

Consider some (discontinuous) distribution function $F$ for which (14) and (2) differ, such that there exists $\theta_{0}\in\mathds{R}$ with $F(\theta_{0})=1$ and $\lim_{x\nearrow\theta_{0}}F(x)<\beta$. Consider observations in the QSR-model with $Z_{1},\ldots,Z_{n}\overset{i.i.d.}{\sim}F$. Then, the true $\beta$-quantile $\vartheta_{\beta}$ is constant with $\vartheta_{\beta}\equiv\theta_{0}$ and for any data $Z_{i}$ the transformation in (3) with truth $\vartheta_{\beta}$ yields $W_{i}=W_{i}(Z_{i},\theta_{0})=1$ for $i=1,\ldots,n$. However, for any other data $Z_{i}$ the constant candidate $\vartheta\equiv\theta_{0}$ with $\theta_{0}\geq\max(Z_{i})$, results in exactly the same transformation $W_{i}=W_{i}(Z_{i},\theta_{0})=1$ for $i=1,\ldots,n$. Consequently, if one allows in the QSR-model for generalized quantiles as in (14) controlling segment overestimation as in Theorem 1.1 appears too ambitious, as it rules out any reasonable estimator, i.e., this can only be achieved when setting $S=0$ always. In fact, quantile estimation of non continuous distributions is well known to require more specific model assumptions in general, see e.g. (Machado and Silva, 2005). The MQS procedure allows quantile regression for discrete distributions, as long as the assumptions of the QSR-model hold.

2.1 Multiscale procedure

Given $\beta$, observations $Z=(Z_{1},\ldots,Z_{n})$ from the QSR-model and a function $\vartheta\in\Sigma$ with $S$ segments as in (1) (which may depend on $n$), the MQS methodology is based in a first step on a multiscale test to decide whether or not $\vartheta$ is a good candidate for the underlying unknown quantile function $\vartheta_{\beta}$. To this end, fix an interval $[x_{i},x_{j}]\subseteq[\tau_{s-1},\tau_{s})$ for some $1\leq s\leq S$, that is, an interval where the candidate $\vartheta$ is constant with value $\theta_{s}$. When the true quantile function $\vartheta_{\beta}$ is constant on $[x_{i},x_{j}]$ as well, the decision problem of whether or not $\vartheta$ coincides with the truth on the interval $[x_{i},x_{j}]$ translates to the testing problem

[TABLE]

Equivalently, one can write (15) using the transformed data in (3) as the hypothesis testing problem in (4). The log-likelihood-ratio test for the hypothesis testing problem (4) and (15), respectively, is given by the test statistic

[TABLE]

where $f_{\beta}(x)=\beta^{x}(1-\beta)^{(1-x)}$ denotes the probability mass function of the Bernoulli distribution and $\overline{W}_{i}^{j}\mathrel{\mathop{:}\!\!=}(j-i+1)^{-1}\sum_{l=i}^{j}W_{l}$. That is, the test is then of the form

[TABLE]

where the threshold $q_{i,j}:=q_{i,j}(\alpha)$ determines the level $\alpha$ of the test (note that, under $H_{ij}$ the distribution of $T_{i}^{j}(W)$ is independent of $\vartheta_{\beta}$). However, we do not know where and on which scales the true $\vartheta_{\beta}$ is constant and thus, we have to consider all intervals on all different scales simultaneously. This means that a candidate function $\vartheta$ is accepted if and only if it gets accepted by all local tests $\Phi_{ij}$ in (17), for appropriately chosen thresholds $q_{ij}$. More precisely, define the penalized multiscale statistic (as a functional on $\Sigma$) as

[TABLE]

with $P_{\ell,n}\mathrel{\mathop{:}\!\!=}\sqrt{2\log\left({\textup{e}\,n}/{\ell}\right)}$, where $\ell=j-i+1$ denotes the scale and $\textup{e}=\exp(1)$. For some given threshold $q=q_{n}(\alpha)$, a candidate function $\vartheta$ is accepted if and only if $T_{n}(Z,\vartheta)\leq q_{n}(\alpha)$. This means that we choose local thresholds $q_{ij}(\alpha)$ in (17) as in (Dümbgen and Spokoiny, 2001; Dümbgen and Walther, 2008; Frick et al., 2014) of the form

[TABLE]

Remark 2.1**.**

A heuristic reasoning for the particular penalization $P_{\ell,n}$ in (18) is as follows. In order to put different scales on equal footing, one has to chose larger thresholds for smaller scales. On the other hand, it follows from Wilk’s theorem (see (Frick et al., 2014) for a more precise argument) that for sufficiently large intervals the local statistics $\sqrt{2T_{i}^{j}}$ are well approximated by the absolute value of a standard normal. As the maximum of $\ell/n$ independent standard normals grows as $\sqrt{2\log(\ell/n)}$, the particular choice in (19) ensures that $T_{n}(Z,\bullet)$ in (18) remains finite even when $n\to\infty$ for any fixed $\alpha\in(0,1)$.

Note that for $X=(X_{1},\ldots,X_{n})$ i.i.d. Bernoulli distributed with mean $\beta$ and $\vartheta=\vartheta_{\beta}$, the statistic $T_{n}(Z,\vartheta_{\beta})$ in (18) follows the same distribution as

[TABLE]

which is bounded by $M_{n}$ in (7) in probability. The distribution of $M_{n}$ does not depend on any characteristics of the unknown $\vartheta_{\beta}$ and hence, we can determine its quantiles via Monte-Carlo simulations. Thus, we can define $q_{n}(\alpha)$ as the $(1-\alpha)$-quantile of $M_{n}$, i.e.,

[TABLE]

This choice guarantees that the true quantile function $\vartheta_{\beta}$ gets accepted by the above testing procedure with probability at least $1-\alpha$, which leads to confidence statement for all quantities of $\vartheta_{\beta}$, see Theorem 2.2.

Remark 2.2**.**

It is shown in (Frick et al., 2014; Dümbgen and Spokoiny, 2001) that $M_{n}$ converges in distribution to an almost surely finite random variable $M$ (which is independent of $\beta$) and hence, $\limsup_{n}q_{n}(\alpha)<\infty$ for any $\alpha\in(0,1)$ . In particular, for large $n$, one may choose the parameter $q(\alpha)$ as a quantile of $M$, which can be simulated once and then stored, lowering the overall computation time.

Finally, we define the MQS estimator (MQSE) $\hat{\vartheta}$ as an element in $\mathcal{H}(q_{n}(\alpha))$ from (6) that minimizes a particular cost function. One cost function that has also been considered in several other works (Eilers et al., 2005; Li and Zhu, 2007; Dümbgen and Kovac, 2009; Belloni and Chernozhukov, 2011; Lee et al., 2018) is the asymmetric absolute deviation loss. That is,

[TABLE]

In Appendix C, we give details on an alternative choice which is based on the Wald-Wolfowitz runs statistic (Wald and Wolfowitz, 1940). We stress that all our theoretical results, detailed in the following, hold for any estimator in $\mathcal{H}(q_{n}(\alpha))$ and are thus, completely independent of this choice. Both cost functions, the one in (21) and in Appendix C, are available in our R package implementation. In simulations in Appendix C we find that they typically perform comparably, with the latter slightly outperforming (21) for higher quantiles.

2.2 Consistency results and confidence statements

From the construction of the MQSE $\hat{\vartheta}$ it follows with $q_{n}(\alpha)$ as in (20), that the corresponding number of segments does not exceed the true number of segments at given error level $\alpha$, as described in equation (8) and Theorem 1.1. As already argued in Section 1, any bound on the underestimation error must depend on further characteristics of the underlying quantile function $\vartheta_{\beta}$ and on the respective quantile jumps. To this end, let

[TABLE]

with $\delta_{s}=\theta_{s}-\theta_{s-1}$, $F_{\min}\mathrel{\mathop{:}\!\!=}\min\{F_{1},\ldots,F_{n}\}$, and $F_{\max}\mathrel{\mathop{:}\!\!=}\max\{F_{1},\ldots,F_{n}\}$ and

[TABLE]

Theorem 2.1** (Underestimation error).**

Consider the QSR-model and $\hat{S}(q)$ as in (5). Then, for $q>0$

[TABLE]

From Theorem 2.1 we obtain that for any fixed $q>0$ and $\Lambda,\Xi>0$ in (10) the probability of underestimating the number of segments vanishes exponentially fast, see Theorem 1.2. Additionally, we obtain from Remark 2.2 and Theorem 2.1 for $\Lambda_{0}>0$ and $\Xi_{0}>0$, that for any $\alpha\in(0,1)$ as $n\rightarrow\infty$

[TABLE]

and, therefore, we can state the following theorem.

Theorem 2.2** (Confidence statements).**

Consider the QSR-model and $\hat{S}(q)$ as in (5), $q_{n}(\alpha)$ as in (20), $\Lambda$ and $\Xi$ as in (10). Then, the set $\mathcal{H}(q_{n}(\alpha))$ in (6) constitutes an asymptotically honest $(1-\alpha)$-confidence band for $\vartheta_{\beta}$ in the QSR-model uniformly over segment functions with minimal $\Lambda_{0},\Xi_{0}>0$, i.e.,

[TABLE]

As detailed in Appendix D, $\mathcal{H}(q)$ can be computed easily simultaneously with the MQS estimator $\hat{\vartheta}$ and from $\mathcal{H}(q_{n}(\alpha))$ confidence intervals for the segment locations $\tau$ and for the quantile values $\theta$ can be constructed.

From Theorems 1.2 and the fact that if $q_{n}\to\infty$, $\alpha_{n}\mathrel{\mathop{:}\!\!=}\textbf{P}(M_{n}>q_{n})\to 0$, it follows directly that for any fixed $\Lambda,\Xi>0$ and some sequence $q_{n}\to\infty$ such that $\nicefrac{{q_{n}}}{{\sqrt{n}}}\rightarrow 0$ MQS performs consistent model selection for the number of c.p.’s. The following result goes beyond this and considers the situation of a sequence of regression function $\vartheta_{\beta}(n)$ in the QSR-model, where the minimal scale $\Lambda=\Lambda_{n}$ and the minimal quantile jump $\Xi=\Xi_{n}$ can vanish as $n\to\infty$.

Theorem 2.3** (Model selection consistency).**

For a sequence $\vartheta_{\beta,n}\in\Sigma$ with $\Lambda_{n}$ and $\Xi_{n}$ as in (10), consider the QSR-model and $\hat{S}(q)$ as in (5). For some sequence $q_{n}\to\infty$, assume the following.

For signals with $\liminf_{n}\Lambda_{n}>0$ and $\liminf_{n}\Xi_{n}>0$, that $\nicefrac{{\sqrt{n}}}{{q_{n}}}\rightarrow\infty$. 2. 2.

For signals with $\liminf_{n}\Lambda_{n}>0$ and $\Xi_{n}\rightarrow 0$, that $\nicefrac{{\sqrt{n}\Xi_{n}}}{{q_{n}}}\rightarrow\infty$. 3. 3.

For signals with $\Lambda_{n}\rightarrow 0$, that $\sqrt{n\Lambda_{n}}\Xi_{n}\geq(2+\epsilon_{n})\sqrt{-\log(\Lambda_{n})}$, for some $\epsilon_{n}$ with $\epsilon_{n}{{\sqrt{-\log\Lambda_{n}}}/{q_{n}}\rightarrow\infty}$.

Then, for such sequences $\vartheta_{\beta,n}$, the MQS selector is consistent, that is, $\mathbb{P}\left(\hat{S}(q_{n})=S\right)\rightarrow 1$.

Theorem 2.3 shows that for a sequence of signals $\vartheta_{\beta,n}$ the number of segments is estimated consistently as long as for the respective minimal scale $\Lambda_{n}$ and minimal quantile jump $\Xi_{n}$ it holds that

[TABLE]

Frick et al. (2014) showed that in the case of Gaussian observations with piecewise constant mean, no method can consistently estimate the number of segments for a sequence of signals with minimal scale $\Lambda_{n}$ and minimal jump height $\Delta_{n}$ whenever $\sqrt{n\Lambda_{n}}\Delta_{n}<\sqrt{-2\log(\Lambda_{n})}+o(1).$ Furthermore, for any continuous distribution $F$ with density $f$ we find that $\lim_{\delta\to 0}\frac{\xi_{F,\beta}(\delta)}{\delta}=f(\theta_{\beta})$. For Gaussian observations, the mean and the median coincide, hence, one obtains that $\Xi_{n}\in\mathcal{O}(\Delta_{n})$. Consequently, possibly up to the constants, (24) cannot be improved in general. Besides model selection consistency and confidence statements for all quantities, MQS also yields minimax optimal estimation rates for the location of segments (up to a log-factor) as the following theorem shows.

Theorem 2.4** (Estimation rates).**

Consider the QSR-model, $\hat{\vartheta}\in\mathcal{H}(q)$ as in (6), and $\Lambda,\Xi$ as in (10). Then, for any $q>0$ and sequence $\epsilon_{n}\searrow 0$

[TABLE]

Note that for $\Lambda^{-1}=o\left(n/\log(n)\right)$ and $q=o\left(\sqrt{\log(n)}\right)$ a sufficient condition for the right hand side to vanish as $n\rightarrow\infty$ is $\epsilon_{n}\geq\log n/(\Xi^{2}n)$, which, up to the log-term equals the minimax optimal sampling rate $1/n$ under a normal error assumption.

Remark 2.3** (Consistency for increasing number of change-points).**

We stress that for all our results, the number of change-points $S$ can depend on $n$, that is, we can consider a sequence $(S_{n})_{n\in\mathbb{N}}$. Note that, the number of change-points $S_{n}$ is upper bounded by the minimal scale $\Lambda_{n}$ via $\Lambda_{n}^{-1}\geq S_{n}-1$ for all $n\in\mathbb{N}$. As argued above, our MQSE remains consistent when $\Lambda_{n}$ converges to [math] at rate $\Lambda^{-1}=o\left(n/\log(n)\right)$, i.e., even when the number of change-points $S_{n}$ goes to infinity as $n$ increases, as long as $S_{n}\in o\left(n/\log(n)\right)$.

3 Implementation

The MQSE and its associated confidence bands can be computed with dynamic programming, employing the double heap structure underlying the computation of quantiles as in (Astola and Campbell, 1989). The major idea of the dynamic programming scheme is analog to the one in (Frick et al., 2014): One successively computes the segmentation (including confidence statements) for the sub-problems consisting of the first $j$ data points, with $j=1,\ldots,n$. Thereby, one makes use of overlapping structure, such that the $j$th solution is updated efficiently from the solutions of the first $j-1$ sub-problems. The basic Bellman equation in this dynamic programming scheme is as follows: For the $j$th sub-problem with observations $Z_{1},\ldots,Z_{j}$, if one knew that the last change-point of the MQSE is at location $i$, then the MQSE for the $j$th sub-problem equals the MQSE for the $i$th sub-problem on the first $i$ data points. There are $\mathcal{O}(n)$ possible locations for the last change-point and there are $n$ successive sub-problems. Thus, overall this dynamic programming scheme has a worst case complexity of $\mathcal{O}(n^{2})$ times the cost to compute the constant solution between the last change-point $i$ and the last data point $j$ on interval $[x_{i},x_{j}]$. Note that, in order to compute a constant MQS solution on some interval $[x_{i},x_{j}]$, one has to intersect all confidence intervals associated with the level $\alpha$-tests $\phi_{k,l}$ in (17) for $i\leq k\leq l\leq j$ . When the intersection is empty, this means that there is no constant function for which all local tests in (17) accept and hence, a new change-point has to be added. For sake of brevity, here we only outline the major differences to the algorithm in (Frick et al., 2014). More details on the implementation, including confidence statements and the specific cost function for the MQSE, are given in Appendix D. For the MSB, which plots estimates for the $0.25$, $0.5$ and, $0.75$ quantiles simultaneously, in order to avoid crossing of different quantiles, we make some further minor modifications detailed in Appendix D.1.

The local level $\alpha$-tests in (17) can be inverted into $(1-\alpha)$-confidence statements for the underlying parameter $\theta_{s}$ in (15) as follows. For $q>0$ define $l(q)$ and $u(q)$ as the two unique solutions of

[TABLE]

such that $0\leq l(q)<u(q)\leq 1$. Then, some straight forward calculations show that $T_{n}(Z,\vartheta)\leq q$ if and only if

[TABLE]

with $\underline{m}_{i,j}=\max\left(1,\lceil(i-j+1)l(\tilde{q})\rceil\right)$, $\overline{m}_{i,j}=\min\left(n,\lfloor(i-j+1)u(\tilde{q})\rfloor\right)$ and $\tilde{q}=(q+P_{\ell,n}^{2}/(2(j-i+1))$. Note that $\underline{m}_{i,j}=\underline{m}_{1,j-i+1}\mathrel{=\!\!\mathop{:}}\underline{m}_{j-i+1}$ and $\overline{m}_{i,j}=\overline{m}_{1,j-i+1}\mathrel{=\!\!\mathop{:}}\overline{m}_{j-i+1}$ only depend on the length of interval $[x_{i},x_{j}]$. Just as in (Frick et al., 2014), the computation of MQS is based on these confidence boxes $\left\{\left[\underline{b}_{ij},\overline{b}_{ij}\right)\;:\;1\leq i\leq j\leq n\right\}$ in (26). In (Frick et al., 2014) the boxes $\left[\underline{b}_{ij},\overline{b}_{ij}\right)$ depend on the local sums of observations $\sum_{l=i}^{j}Z_{l}$ and its dynamic program explores that these local sums and hence, boxes of an interval $[i,j]$, can be updated in $\mathcal{O}(1)$ time from the boxes of intervals $[i+1,j]$ and $[i,j-1]$, respectively. In contrast, for MQS the boxes depend on the particular quantiles of the observations $\{Z_{i},\ldots,Z_{j}\}$. Thus, in order to adapt the programming scheme of Frick et al. (2014), one has to update the running $\overline{m}_{1},\underline{m}_{1},\ldots,\overline{m}_{n},\underline{m}_{n}$-quantiles efficiently. Here, we use double heap structures as in (Astola and Campbell, 1989) to update the running quantile in $\mathcal{O}(\log(n))$ time. In total, this increases the overall computation time by a log-factor, with worst case complexity of order $\mathcal{O}(n^{2}\log(n))$. However, depending on the reconstructed signal, pruning steps and the use of smaller interval systems often lead to a computation time which is almost linear in $n$, see Frick et al. (2014) for details.

4 Simulations

In the following, we explore MQS in a simulation study. Thereby, the choice of the threshold parameter $q$ is essential as it balances detection and overestimation of the number of segments and hence false positives. Thus, $q$ can be seen as a tuning parameter of MQS. Via the one-to-one correspondence of a confidence level $\alpha$ and $q_{n}(\alpha)$ in (20), in the following, we choose $\alpha=0.1$ and hence $q=q_{n}(0.1)$. In this way, the probability that MQS overestimates the number of segments does not exceed $10\%$. Such a choice depends on the application. However, we see a great advantage of our methodology, to rely only on this parameter, which has an immediate statistical meaning. For a more refined discussion on the choice of threshold parameter $q$ and possible data driven model selection procedures, we refer to (Frick et al., 2014). For example, another possible parameter choice for $q$ is via minimizing the right hand side of (12) together with Monte-Carlo simulations of $M_{n}$.

In the following, we consider five different competitors for MQS for which software is available online: SMUCE from (Frick et al., 2014), HSMUCE from (Pein et al., 2017), wild binary segmentation (WBS) from (Fryzlewicz, 2014), R-FPOP from (Fearnhead and Rigaill, 2017), quantsmooth (QS) from (Eilers et al., 2005), NOT from (Baranowski et al., 2019), and NWBS (Padilla et al., 2019). SMUCE provides a multiscale methodology for normal observations with homogeneous variance. HSMUCE is also designed for normally distributed data, but is robust against changes in variance. WBS can be seen as a “greedy” procedure which successively adds changes based on a localized CUMSUM statistic. The theoretical results implicitly assume normally distributed observations with change in mean. R-FPOP is a penalized cost approach that uses the biweight loss, which is designed to be particularly robust to extreme outliers. It considers arbitrary changes in the underlying distribution, but cannot be tuned to search for changes in some specific quantiles. In particular, it considers that observations in segments are i.i.d.. QS is a smoothing method which performs minimization of the asymmetric absolute deviation loss (recall (21)) together with $L_{1}$ penalization and hence, can compute arbitrary quantile curves. However, QS is not designed for change-point detection, but for data visualization through smoothing. This will be clear in simulations, where the number of segments is very often very large. NOT choose random subsamples and then uses a tailor-made contrast function to find possible change-points. The contrast functions more related to our work are the (HT) for heavy-tailed noised and (VAR) for changes in variance. NWBS is a non parametric method based on a CUSUM-like statistic for changes in the distribution function. It then uses wild binary segmentation ideas, as in WBS. For all competitors, we choose tuning parameters as default in the available software. We summarize some features of these methods in Table G.1.

As measures of evaluation for the simulation study we use the number of estimated segments, the mean absolut squared error (MIAE) $\sum_{i=1}^{n}|\hat{\vartheta}(x_{i})-\vartheta_{\beta}(x_{i})|/n$, and the entropy-based V-measure introduced in (Rosenberg and Hirschberg, 2007). The latter takes values in $[0,1]$ and measures whether given clusters include the correct data points of the corresponding class. Larger values indicate higher accuracy with $1$ corresponding to a perfect segmentation. All results were obtained from $1,000$ Monte Carlo runs. We also performed simulations on the coverage properties for the confidence intervals and bands. Details can be found in Appendix E. Similar to other multiscale approaches (Frick et al., 2014; Behr et al., 2018), we find that the level is typically exceeded, meaning that, in general, these confidence statements are conservative.

4.1 Additive error

First, we consider an additive model with i.i.d. error terms, that is,

[TABLE]

with $\vartheta\in\Sigma$ and $\varepsilon_{1},\ldots,\varepsilon_{n}$ i.i.d. according to some distribution. For observations as in (27) in the QSR-model the quantile functions $\vartheta_{\beta}$ are shifted versions of $\vartheta$, namely, $\vartheta_{\beta}=\vartheta+\theta_{\beta}$, where $\theta_{\beta}$ is the $\beta$-quantile of $\epsilon_{1}$. In particular, for any quantile $\beta\in(0,1)$, $\vartheta_{\beta}$ has the same number and locations of segments. Here, we consider $\vartheta$ as in Figure G.2 (top row), which has $7$ segments and $n=1,988$. For the error terms $\varepsilon_{i}$ we consider normal distribution $\varepsilon_{i}\sim\mathcal{N}(0,\sigma^{2})$, $t$-distribution with $3$ degrees of freedom $\epsilon_{i}\sim t_{3}\sigma/\sqrt{3}$ and variance $\sigma^{2}$, rescaled Cauchy distribution (heavy tails) $\varepsilon_{i}\sim 0.02\,\text{Cauchy}(0,1)$, and rescaled chi-square distribution (skewed) with $3$ degrees of freedom and median [math], that is, $\varepsilon_{i}\sim(\chi_{3}^{2}-\beta_{0.5}(\chi_{3}^{2}))\sigma/\sqrt{6}$ with variance $\sigma^{2}$, with $\sigma^{2}=0.04$ and $\beta_{0.5}(\chi_{3}^{2})$ the median of the distribution $\chi_{3}^{2}$. (see Figure G.2).

The results for MQS$(\beta)$, $\beta=0.25,0.5,0.75$ are shown in Table G.2. Note that, in general, the MQS selector of $S$ for $\beta=0.5,0.75$ seems to have a higher detection power as for $\beta=0.25$ in this example. This is due to the fact that the test signal $\vartheta$ in Figure G.2 (top row) has 4 jumps upwards but just 2 jumps downwards. It is easy to check that jumps upwards have a stronger influence on higher (overall) empirical quantiles. MQS is a reasonable estimator in the four scenarios presented in this section. It is robust against outliers, as well as to skewness of the distributions. For normally distributed data is not surprising that methods tailored for this scenario outperform MQS, such as SMUCE, HSMUCE, WBS, and NOT(VAR), in particular in terms of MIAE and segmentation accuracy. However, MQS(0.5) is comparable in terms of calculating the number of segments. On the other hand, for heavy tailed and skewed distributions SMUCE, WBS, and NOT(VAR) fail completely, while MQS(0.5) retains a very high segmentation accuracy, also outperforming HSUMCE for skewed data. The nonparametric methods tailored for heavy-tails, R-FPOP, NOT(HT), and NWBS, perform comparably to MQS in this additive i.i.d. setting. Although, NOT(HT) highly overestimated the number of segments for Cauchy noise. However, as shown in Section 4.2, these methods are highly sensitive to any distributional changes (not present in this set-up), in contrast to MQS which robustly estimates quantile segments among all scenarios. QS performs comparably to MQS in terms of MIAE, however, it highly overestimates the number of segments and performs worse in terms of estimation accuracy and computational speed. This is partly due to the fact that it is designed as a smoothing method and not a segmentation method.

4.2 Changes in variance

The QSR-model further allows for changes in variance or other characteristics which are independent of changes in the respective $\beta$ quantile $\vartheta_{\beta}$. Here, we consider changes in variance for normally and $t$ distributed (with $3$ d.f.) observations, $n=2,000$. In Figure G.3 (top row) the underlying median (mean) (solid line) and variance (dashed line) functions are displayed. The second and third row show the true $0.25$, $0.5$, and $0.75$ quantile functions (black lines) and the MQS box plot (red lines). Note that in this example, the $0.5$ quantile has $4$ and the $0.25$,$0.75$ quantiles have $6$ segments.

Simulation results are shown in Table G.3. MQSE appears very robust against changes in variance within a segment even for heavy tailed distributions. It estimates the correct number of change-points with high probability. At the same time, the MQSE for the $0.25$ and $0.75$ quantiles depict changes in variance. Similar as before, methods designed for Gaussian (homoscedastic) distributions, such as SMUCE and WBS fail completely for the heavy tailed setting. Note that, although, HSMUCE was particularly designed for Gaussian distributions with change in variance, MQS(0.5) outperforms it in terms of estimated number of change-points. As expected, nonparametric methods, such as NWBS and R-FPOP, are sensitive to these changes in variance and perform poorly compared to MQS in terms of estimated number of segments. Just as in the previous section, QS highly overestimates the number of segments. NOT(HT) performs comparably to MQS(0.5) in both settings, being better in the normal case, but worse for t distribution. NOT(VAR) performs bad in both scenarios. In summary, we find that, whereas other methods only work well in specific situations, MQS is very robust to all scenarios, estimating quantile segments with high accuracy independent of any distributional assumptions.

5 Real data examples

5.1 Copy Number Aberrations

Copy Number Aberrations (CNA’s) are sections of DNA in the genome of cancer cells that are either multiplied or deleted, relative to the state present in normal tissue. CNA’s are important factors of tumor progression, through the deletion of tumor suppressing genes and the multiplication of genes involved for example in cell division. The number of copies of DNA sections, depending on chromosomal loci, corresponds to a segment function, where a segment corresponds to a different copy number. A common measurement technique is via whole genome sequencing (WGS). Thereby, the tumor DNA is fragmented into several pieces. Then the single pieces are sequenced using short “reads”, and finally these reads are aligned to a reference genome by a computer. Statistical modeling of WGS data is particularly difficult as random variations and systematic biases, such as mappability and CG bias, lead to violations of parametric model assumptions, such as normal or Poisson, see e.g. (Liu et al., 2013). Quantile segmentation with MQS does not require any such specific model assumptions and hence, is particularly suited for this setting.

Figure 1.4 shows (pre-processed111Sequencing produces spatial artifacts in the data and waviness, which can be pre-process using standard procedures of smoothing filter, baseline correction and binning, see, (Behr et al., 2018) for details.) WGS data of cell line LS411 from colorectal cancer. Sequencing was performed by Complete Genomics in collaboration with the Wellcome Trust Centre for Human Genetics at the University of Oxford. For this particular data set, it is known that the underlying CNA’s only take values in the natural numbers. This is because it was collected under special conditions, where cells come from a single homogeneous tumor-clone, see (Behr et al., 2018) for more details. This allows certain validation of the estimated segments, something which is not feasible for most real patient tumors.

Figure 1.4 shows the MSB (MQS’ estimated $0.25$, $0.5$, and $0.75$ quantiles) at confidence level $1-\alpha=0.99$. MQSE recovers most of the signal structure correctly. In particular, MQSE is way more robust than SMUCE (Frick et al., 2014), WBS (Fryzlewicz, 2014), R-FPOP (Fearnhead and Rigaill, 2017), NOT(HT), and NOT(VAR) (Baranowski et al., 2019), see Figure G.1. SMUCE, WBS, R-FPOP, NOT introduce many artificial changes, which cannot be present in the signal as, in this particular example, it is known to only take integer values. HSMUCE (Pein et al., 2017) and NWBS (Baranowski et al., 2019) (fifth and tenth rows in Figure G.1) are more robust, but still adds artificial changes. The sixth row of Figure G.1 shows the estimated $0.25$, $0.5$, and $0.75$ quantiles of QS (Eilers et al., 2005). Similar to MQS, it correctly recovers most of the signal structure, but it misses some underlying changes, see, in particular, the last change at data point $7148$. Moreover, QS has a much higher running time compared to MQS: while MQS took 31 seconds to run each quantile for this data set, QS took 54 minutes.

5.2 Ion channel data

Ion channels are pore-forming proteins that allow ions to pass through a cell membrane. They are vital for several processes like excitation of neurons and muscle cells. The pores of an ion channel can open and close, a process called gating, often as a result of external stimuli. Therefore, the amount of ions that can pass through a channel is not constant in time (Chung et al., 2007). A major tool for a quantitative analysis of the gating dynamics is the patch clamp technique, which allows to measure the conductance of a single ion channel in time (Sakmann and Neher, 1995). Roughly speaking, this kind of data is obtained by inserting a single ion channel in an (often artificial) membrane surrounded by an electrolyte with an electrode to measure the current while constant voltage is applied. These recordings can be modeled as a segment function disturbed by an error, see e.g. (Pein et al., 2017; Gnanasambandam et al., 2017). The particular data set considered in Figure 1.1 comes from a single channel of the bacterial porin PorB from the Steinam lab (Institute of Organic and Biomolecular Chemistry, University of Göttingen). The measurement protocol incorporates a lowpass filter which leads to local dependencies of the error terms, see (Pein et al., 2017). To remove these dependencies, which violate modeling assumptions of the QSR-model, we subsampled every 11th observation. The MQS multiscale box plot is shown in the top row of Figure 1.1. A common feature of ion channel data, called open channel noise, is that the noise variance in open states is often much higher than in closed states (Sakmann and Neher, 1995, Section 3.4.4). MQS is very robust against this heterogeneity while at the same time reliably detects most of gating events. Figure G.4 shows that SMUCE (Frick et al., 2014), R-FPOP (Fearnhead and Rigaill, 2017), and WBS (Fryzlewicz, 2014) introduce a lot of artificial changes, because of increased variance for open channel noise. QS (Eilers et al., 2005) (fifth row in Figure G.4) misses most of the structural segment changes. Although particularly tailored to this application, HSMUCE (Pein et al., 2017) (and NOT(VAR) (Baranowski et al., 2019)), which assume heteroscedastic normal observations, do not seem to be superior to MQS here. NOT(HT) shows a similar reconstruction to MQS(0.5). NWBS misses considerable spikes that other methods can estimate. Moreover, in contrast to the other methods, MQS explicitly quantifies the change in variance via the inter-quantile distance in the multiscale segment boxplot.

6 Discussion

In this work we proposed a new approach for quantile segmentation, that does not just provide an estimate with optimal detection rates, but also comes with honest confidence statements for the number of change-points, their locations, as well as confidence bands for the full underlying quantile function. These results hold under the minimal assumption of independence between different observations and do not make any other distributional assumptions. In simulations and real data examples we also observe empirically that the MQSE is very robust and consistently estimates changes in quantiles in various settings. While other methods usually only work well in a particular setting (for which they are designed for), MQSE can be applied universally for any quantile segmentation task (under the independence assumption). In the following we discuss a few possible extensions.

Multivariate change-point estimation

One may wonder whether the MQS procedure can also be applied in a multivariate setting, i.e., for $Z_{i}\in\mathbb{R}^{d}$ for $d>1$. However, it is well known that the univariate concept of quantiles as points of division of the mass of a probability distribution is not immediate to generalize for higher dimension. No single point in a multidimensional ($d>1$) space can divide the space in this way, what makes it difficult to extend the QSR-model (1) in a straightforward manner to the multivariate setting (for a survey on multidimensional medians see (Small, 1990)). Some attempts have been made in terms of depth curves (Hallin et al., 2010), and in this sense one could attempt to look for changes in these curves. We believe this is an interesting venue to take, however outside of the scope of this paper. Note, that even for the more simple concept of component-wise quantiles, there is a severe computational challenge, since the inversion of the local tests leads to ellipsoids (instead of intervals as in the one-dimensional case), which are hard to intersect repetitively.

Change-point estimation for dependence structure

The only assumption that we make for the theoretical analysis of the MQSE is independence of the sequential observations. One might wonder whether this assumption can be further weakened. Clearly, without any specific assumptions on the dependence structure there is no hope for consistent change-point estimation. For the Gaussian case Tecuapetla-Gómez and Munk (2017) provide some strategy to estimate the underlying covariance structure for piece-wise constant signals with m-dependent errors. Taking such an estimated covariance structure, one might correct the local likelihood ratio tests for dependence (see (Frick et al., 2014) for the Gaussian mean case) and thus, preserve the confidence statements of MQS. However, any consistent estimation of such covariance structure will unavoidably come at the cost of making additional assumptions on the underlying distribution. The major strength of our approach is that such assumptions are not required (and instead we only rely on the independence assumption). In practice, we still found that as long as the dependence structure is sufficiently weak, MQSE is robust to this. Interestingly, a slightly negative autocorrelation can even be beneficial for the detection power of MQS, which is consistent with results on minimax detection boundaries in the autocorrelated case Enikeeva et al. (2020). We provide some simulation results in Appendix F.

Acknowledgments

The authors acknowledge support of DFG-RTG 2088, DFG-SFB 803 Z02, and DFG Cluster of Excellence 2067 MBExC. MB was supported by DFG postdoctoral fellowship BE 6805/1-1. Helpful comments from Chris Holmes, Housen Li, and Florian Pein are gratefully acknowledged.

Appendix A Proofs of Section 2

The proofs of this section are similar in spirit to those in (Frick et al., 2014). However, there are important diferences due to the discrete nature of the transformation (3) and the corresponding convergence rates of the empirical quantiles. Before we prove the main results of Section 2, we require a couple of auxiliary results.

To this end, let $W_{1},W_{2},\ldots,W_{n}$ be i.i.d. Bernoulli distributed random variables with mean $\beta$ and let for $x\geq k\geq 0$

[TABLE]

For fixed $n$ define the random variables

[TABLE]

Theorem A.1**.**

Let $k\in\mathbb{N}$ with $k\geq 1$ and $q_{n}(\alpha)$ as in (20). Then

[TABLE]

Proof.

Define the stopping times

[TABLE]

Note that

[TABLE]

Consider the Markov process $\left(\sum_{i=1}^{n}W_{i}\right)_{n\in\mathbb{N}}$. By the strong Markov property, for any stopping time $\tau$ the process $(\sum_{i=1}^{n}W_{i+\tau})$ is independent of $W_{1},\ldots,W_{\tau}$ conditioned on $\tau<\infty$. Note also that $(\sum_{i=1}^{n}W_{i+\tau})$ is identically distributed for all stoping times $\tau$, because it is the sum of $n$ i.i.d. random variables. This implies that

[TABLE]

are independent and identically distributed. Therefore, for any $k\geq 1$ and $x>0$ it follows that

[TABLE]

By definition $\zeta_{1}(q)\leq n$ implies that $M_{n}>q$, for $M_{n}$ as in (7). Therefore, it follows from (20) that

[TABLE]

∎

For the following two theorems let $Z_{1},\ldots,Z_{n}$ be independent random variables with distribution functions $F_{1},\ldots,F_{n}$ and $\beta\in(0,1)$ be given. Assume that the beta quantile is identical for all random variables, i.e.,

[TABLE]

Define the empirical quantile and the minimal distribution function, the maximal distribution function, and the mean distribution function as:

[TABLE]

Note that $F_{\min}$, $F_{\max}$, and $\overline{F}$ are again distribution functions with $\beta$-quantile $\theta_{\beta}$.

Lemma A.1**.**

Let $X_{1},\ldots,X_{n}$ and $Y_{1},\ldots,Y_{n}$ be independent Bernoulli random variables s.t. $X_{i}\sim B(p_{i})$, $Y_{i}\sim B(p_{\min})$, and $Z_{i}\sim B(p_{\max})$ for all $i=1,\ldots,n$ and $p_{\min}=\min\{p_{1},\ldots,p_{n}\}$, $p_{\max}=\max\{p_{1},\ldots,p_{n}\}$. Then

[TABLE]

Proof.

Choose $W_{i}\sim B(p_{\min}/p_{i})$ independently for all $i=1,\ldots,n$ and independently of $X_{1},\ldots,X_{n}$. Let

[TABLE]

Then $\tilde{Y}_{i}\sim B(p_{\min})$ and $\sum_{i=1}^{n}\tilde{Y}_{i}\leq\sum_{i=1}^{n}X_{i}$ a.s. Therefore,

[TABLE]

The other inequality follows analogously. ∎

Theorem A.2**.**

Let $\xi_{F,\beta}$ be as in Definition 1.1. Then for any $\delta>0$

[TABLE]

Proof.

[TABLE]

Now, let $X_{1},\ldots,X_{n}$ be i.i.d. random variables with distribution function $F_{\min}$. Then it holds

[TABLE]

From Lemma A.1 it follows that

[TABLE]

where the last inequality follows from the Dvoretzky-Kiefer-Wolfowitz inequality, see e.g. (Massart, 1990). The other inequality follows analog. ∎

Theorem A.3**.**

Let $T_{1}^{n}$ be as in (16) and $W$ as in (3). Then, for any $\delta,\epsilon>0$

[TABLE]

Proof.

For any $p,q\in(0,1)$, Pinsker’s inequality, see e.g. (Tsybakov, 2009, Lemma 2.5), implies

[TABLE]

$W_{1},\ldots,W_{n}$ are independent, but not identically, Bernoulli distributed random variables with mean $F_{i}(\theta_{\beta}+\delta)$. Define $\overline{W}=n^{-1}\sum_{i=1}^{n}W_{i}$, $\overline{F}(\theta_{\beta}+\delta)=n^{-1}\sum_{i=1}^{n}F_{i}(\theta_{\beta}+\delta)$, and assume that ${\xi_{\overline{F},\beta}(\delta)>\sqrt{\nicefrac{{\epsilon}}{{2n}}}}$ (otherwise the assertion follows trivially). Then (29) implies

[TABLE]

where the last inequality follows from Hoeffding’s inequality. ∎

Let $\delta>0$ and note that $\xi_{\overline{F},\beta}(\delta)\geq\xi_{F_{\min},\beta}(\delta)$, since $F_{\min}\leq\overline{F}$ implies

[TABLE]

where the last inequality follows from the monotonictiy of $F_{\min}$ and $\delta>0$. Similarly, $\xi_{\overline{F},\beta}(-\delta)\geq\xi_{F_{\max},\beta}(-\delta)$.

Now we are ready to prove the main results of Section 2.

Proof of Theorem 1.1.

Let $J(\vartheta)$ denote the set of change-points of $\vartheta$. First, note that

[TABLE]

with $\tilde{Z}=Z-\vartheta$ and $\tilde{S}_{n,\alpha}$ as in (5) with $Z$ replaced by $\tilde{Z}$. Thus, we can assume w.l.o.g. that $\vartheta\equiv\theta_{0}$ and $S=1$.

Observe that $\hat{S}_{n,\alpha}\geq 2s+2$ implies that the multiscale constraint for the true regression function $\vartheta$ is violated on at least $s+1$ disjoint intervals $[\nicefrac{{i_{1}}}{{n}},\nicefrac{{j_{1}}}{{n}}],\ldots,$ $[\nicefrac{{i_{s+1}}}{{n}},\nicefrac{{j_{s+1}}}{{n}}]$ $\subseteq[0,1]$, that is

[TABLE]

and it follows from Theorem A.1 that

[TABLE]

∎

Proof of Theorem 2.1.

Define for $s=1,\ldots,S-1$ the intervals

[TABLE]

Note that these intervals are pairwise disjoint because $\lambda_{s}\leq\min\{\tau_{s}-\tau_{s-1},\tau_{s+1}-\tau_{s}\}$.

Let $\theta_{s}^{+}=\max\{\theta_{s},\theta_{s+1}\}$ and $\theta_{s}^{-}=\min\{\theta_{s},\theta_{s+1}\}$ and split the interval $I_{s}$ accordingly, i.e.

[TABLE]

We are interested in the event that a function exists which is constant on $I_{s}$ and fulfills the multiscale constraints in both $I_{s}^{-}$ and $I_{s}^{+}$, i.e.

[TABLE]

Observe that either $\hat{\theta}\leq\theta_{s}^{+}-\nicefrac{{\delta_{s}}}{{2}}$ or $\hat{\theta}\geq\theta_{s}^{-}+\nicefrac{{\delta_{s}}}{{2}}$ and define

[TABLE]

[TABLE]

Due to the independence of $\Omega_{s}^{+}$ and $\Omega_{s}^{-}$ and the fact that $\Omega_{s}\subseteq\Omega_{s}^{+}\cup\,\Omega_{s}^{-}$, we get that

[TABLE]

Next, we prove an upper bound for $\mathbb{P}(\Omega_{s}^{-})$; a bound for $\mathbb{P}(\Omega_{s}^{+})$ follows from symmetry.

Let $F_{1},\ldots,F_{\left|I_{s}^{-}\right|}$ be the distribution functions of the random variables in $I_{s}^{-}$, then $F_{i}^{-1}(\beta)=\theta_{s}^{-}$ for all $i=1,\ldots,\left|I_{s}^{-}\right|$ . Moreover, let $F_{min}$ and $\overline{F}$ denote the minimum and mean distribution function of the random variables in $I_{s}^{-}$. Let

[TABLE]

be the empirical quantile of the observations on the interval $I_{s}^{-}$. Thus, for all $\hat{\theta}\geq\theta_{s}^{-}+\nicefrac{{\delta_{s}}}{{2}}$, if $\zeta_{I_{s}^{-}}\leq\theta_{s}^{-}+\nicefrac{{\delta_{s}}}{{2}}$, we get that

[TABLE]

Moreover, the function

[TABLE]

is strictly convex with minimum $\beta$ and hence, f\big{|}_{[\beta,1]} is strictly increasing. Thus, (32) implies that for all $\hat{\theta}\geq\theta_{s}^{-}+\nicefrac{{\delta_{s}}}{{2}}$, if $\zeta_{I_{s}^{-}}\leq\theta_{s}^{-}+\nicefrac{{\delta_{s}}}{{2}}$

[TABLE]

and hence,

[TABLE]

where the last inequality follows from Theorems A.2 and A.3. Moreover, note that $\xi_{\overline{F},\beta}(\nicefrac{{\delta_{s}}}{{2}})\geq\xi_{s}(\nicefrac{{\delta_{s}}}{{2}})$. Hence,

[TABLE]

with $\gamma_{n,s}$ as in (23). For $s=1,\ldots,S-1$ define the random variables

[TABLE]

Observe that $X_{s}=1$ implies that any function $\hat{\vartheta}\in\Sigma$ with $T_{n}(W(Z,\vartheta),\beta)\leq q$ has at least one change-point on $I_{s}$. Since $I_{1},\ldots,I_{S-1}$ are pairwise disjoint this implies $\hat{S}(q)-1\geq\sum_{s=1}^{S-1}X_{s}$ . Therefore

[TABLE]

∎

Proof of Corollary 1.2.

Let $\gamma_{n,s}(q)$ be as in (23), $\Xi,\Lambda$ as in (10), and $\gamma_{n}(q)$ be defined as $\gamma_{n}(q)=\min_{1\leq s\leq S-1}\gamma_{n,s}(q)$. Then

[TABLE]

From Theorem 2.1 it follows that

[TABLE]

Using the inequality $(1-x)^{m}\geq 1-mx$ for all $x\in(0,1)$ and $m\in\mathbb{N}$ and (33) it follows that

[TABLE]

where the last inequality comes from the fact that $a^{2}-(a-b)^{2}\leq 2ab$, with $a=\sqrt{n\Lambda}\,\Xi$ and $b=\nicefrac{{q}}{{\sqrt{2}}}+\sqrt{\log(\nicefrac{{2\textup{e}}}{{\Lambda}})}$. The result follows from the fact that $S\leq\Lambda^{-1}$ ∎

Proof of Theorem 2.3.

From the proof of Corollary 1.2 we get the following bound

[TABLE]

with $\Gamma_{1,n}=\left(\sqrt{n\Lambda_{n}}\,\Xi_{n}-\nicefrac{{q_{n}}}{{\sqrt{2}}}-\sqrt{\log(\nicefrac{{2\textup{e}}}{{\Lambda_{n}}})}\right)_{+}^{2}+\log\Lambda_{n}\,$ and $\,\Gamma_{2,n}=n\Lambda_{n}\,\Xi_{n}^{2}+\log\Lambda_{n}$. Then a sufficient condition for $\mathbb{P}\left(\hat{S}(q)<S\right)\rightarrow 0$ is that $\Gamma_{1,n}\rightarrow\infty$ and $\Gamma_{2,n}\rightarrow\infty$ as $n\rightarrow\infty$.

Cases 1 and 2: If $\liminf\Lambda_{n}>0$, then a sufficient condition for $\Gamma_{1,n}\rightarrow\infty$ is that

[TABLE]

This holds if $\nicefrac{{\sqrt{n}\Xi_{n}}}{{q_{n}}}\rightarrow\infty$. Moreover, if this is the case then $\Gamma_{2,n}\rightarrow\infty$ as $n\rightarrow\infty$ and the proof is finished.

Case 3: If $\Lambda_{n}\rightarrow 0$, assume that $\sqrt{n\Lambda_{n}}\Xi_{n}\geq(2+\epsilon_{n})\sqrt{-\log\Lambda_{n})}$ for a sequence $\epsilon_{n}$ such that $\epsilon_{n}\sqrt{-\log\Lambda_{n}}/q_{n}\rightarrow\infty$. Using the inequality $\sqrt{x+y}-\sqrt{x}\leq\nicefrac{{y}}{{(2\sqrt{x})}}$ for $x,y\geq 0$ we obtain

[TABLE]

where the last inequality comes from the fact that $(a+b)^{2}-a^{2}\geq b^{2}$, for $a,b\geq 0$. A sufficient condition for $\Gamma_{1,n}\rightarrow\infty$ is then that $\epsilon_{n}\sqrt{-\log\Lambda_{n}}/q_{n}\rightarrow\infty$, as it was assumed. Note that $\sqrt{n\Lambda_{n}}\Xi_{n}\geq(2+\epsilon_{n})\sqrt{-\log\Lambda_{n}}$ implies

[TABLE]

which goes to infinity as by definition $\liminf_{n}\epsilon_{n}\geq 0$. ∎

Proof of Theorem 2.4.

As in the proof of Theorem 2.1, define $S$ disjoint intervals

[TABLE]

and define $I_{s}^{+}$, $I_{s}^{-}$, $\theta_{s}^{+}$ and $\theta_{s}^{-}$ accordingly. Now assume an estimator $\hat{\vartheta}$ of $\vartheta$, with estimated number of segments $\hat{S}$, such that $T_{n}(Z,\hat{\vartheta})\leq q$ and

[TABLE]

In other words, there exists $0\leq s\leq S-1$ such that $|\hat{\tau}_{l}-\tau_{s}|>\epsilon_{n}$ for all $0\leq l\leq\hat{S}-1$, i.e. $\hat{\vartheta}$ does not have a change-point in the interval $I_{s}$. Then, as in the proof of Theorem 2.1,

[TABLE]

Finally, the assertion follows by replacing $\lambda_{s}$ by $\epsilon_{n}$ in the proof of Theorem 2.1. ∎

Appendix B Illustrative example with simulated data

In the following we provide another example with synthetic simulated data which highlights MQS robustness to arbitrary distributional changes. The data, shown as black dots in the top row of Figure B.1, comes from 3 different distributional regimes. The first 350 observations are drawn from a normal distribution with variance 1, the next 1190 observations from a t distribution with 3 d.f. and variance 0.05 and the last 945 from a $\chi^{2}$ distribution with 1 d.f. and variance 0.05. This data setting is covered by the QRS-model, which only assumes independence of the observations, but does not make any other distributional assumptions. The black line in the top row of Figure B.1 shows the true median of the underlying signal and the gray lines show the $0.25$ and $0.75$ quantiles. For example, the median is a segment function $\vartheta_{0.5}\in\Sigma$ (black line) with $S=7$ segments. Note that the QSR-model allows that different quantile curves may have entirely different segments, for example, the $0.25$- and $0.75$-quantile functions $\vartheta_{0.25},\vartheta_{0.75}$ (gray lines) have an additional segment in the first regime ($S=8$). The second row in Figure B.1 shows the MQS estimates (MQSE) for the median (red line) and the 0.25- and 0.75-quantiles (light red lines) for the data as in the first row. The following three rows show the three estimates together with the corresponding confidence bands. Note that, MQS accurately estimates change-point locations even when changes in the distribution occur, either at a change-point location or within a constant segment.

Appendix C MQSE with Wald-Wolfowitz runs statistic

In the following, we present an alternative approach to (21) for MQSE as a particular element in (6), which uses the transformed data in (3) and is based on the Wald-Wolfowitz runs statistic (Wald and Wolfowitz, 1940). Recall that, for $\vartheta=\vartheta_{\beta}$ the true underlying regression function in (3), the transformed observations $W_{1},\ldots,W_{n}$ are i.i.d. Bernoulli distributed with mean $\beta$. Let $R$ be the number of runs of the sequence $W_{1},\ldots,W_{n}$, i.e.,

[TABLE]

and $N_{0}(W)\mathrel{\mathop{:}\!\!=}\sum_{i=1}^{n}W_{i}$. Then, conditioned on $N_{0}=n_{0}$ and $n-N_{0}=n_{1}$, it holds that $R$ is asymptotically normally distributed with mean $\mu={2n_{1}n_{0}}/{n}-1$ and variance $\sigma^{2}={2n_{1}n_{0}(2n_{1}n_{0}-n)}/({n^{2}(n-1)})$ see (Wald and Wolfowitz, 1940). Hence, for $\vartheta=\vartheta_{\beta}$ in (3) it follows that $\textbf{P}(R(W)=r,\;N_{0}(W)=k)=\textbf{P}(R(W)=r,\;|\;N_{0}(W)=k)\textbf{P}(N_{0}(W)=k)$ is asymptotically well approximated by

[TABLE]

Definition C.1**.**

For a threshold $q>0$, $\mathcal{H}(q)$ as in (6), the MQS estimator (MQSE) is defined as

[TABLE]

Remark C.1**.**

Note that the estimator MQSE may not be unique, as for two estimators $\hat{\vartheta},\tilde{\vartheta}\in\mathcal{H}(q)$ it is possible that $R(\hat{W})=R(\tilde{W})$ and $N_{0}(\hat{W})=N_{0}(\tilde{W})$, for $\hat{W}\mathrel{\mathop{:}\!\!=}W(Z,\hat{\vartheta})$ and $\tilde{W}\mathrel{\mathop{:}\!\!=}W(Z,\tilde{\vartheta})$ as in (3). For $n$ sufficiently large, we found it to be unique in most of the cases. When it is not unique we choose the first feasible change-point location.

In the following, we provide a small simulation study to compare the segmentation performance of MQSE for the Koenker loss from equation (21) in the main text and the runs based loss from equation (35). To this end, we considered a simple bump function with $n=700$ equally spaced observations in an additive regression model for Gaussian and Cauchy noise, respectively. The underlying median signal is [math] for the first $300$ and last $300$ observations and $1$ elsewhere. Table C.1 shows the MISE and the V-measure for $\beta=0.25,0.5,0.75$. For Gaussian noise, the estimate from Definition C.1 sightly outperforms the Koenker based estimate from equation (21) in the main text, in particular for the higher quantile $\beta=0.75$. For the heavy tailed Cauchy noise, it also outperforms the Koenker based estimate for the $0.75$ quantile, but not for the other two quantiles. In summary, we find that both estimates perform comparable, with the estimate from Definition C.1 slightly outperforming the Koenker based estimate for higher quantiles. Therefore, we conclude that this choice is not of major concern in practice. Both estimates are available in the online implementation at https://github.com/ljvanegas/mqs.

Appendix D Details of the Implementation

Consider some given observations $Z_{1},\ldots,Z_{n}\in\mathds{R}$, quantile $\beta\in(0,1)$, and threshold $q\in\mathbb{R}$ (corresponding to some confidence level $\alpha$ as in (7)). In this section we give more details on the implementation of MQS. Pseudocode for MQS is given in Algorithm 1, where the general dynamic programming structure is similar to the algorithm for SMUCE in (Frick et al., 2014).

Bellman equation:

Following a dynamic programming approach, for $j=1,\ldots,n$ we successively solve the estimation problem for $Z_{1},\ldots,Z_{j}$ adding additional change-points when necessary. When solving sub-problem $j$, that is, calculating the respective MQSE $\hat{\vartheta}_{1,j}$ for observations $Z_{1},\ldots,Z_{j}$ we have already solved sub-problems $l$ with estimates $\hat{\vartheta}_{1,l}$ for $l=1,\ldots,j-1$, where $\hat{\vartheta}_{1,j-1}$ had $k$ change-points. If we knew, that the last change-point of $\hat{\vartheta}_{1,j}$ was at location $i<j$, then

[TABLE]

Note that (36) corresponds to a dynamic programming Bellman equation, which connects the solution of a previous sub-problem $i$ to the solution of the current sub-problem $j$. For an illustration see Figure D.1, where the red line represents $\hat{\vartheta}_{1,i}$ and the blue line $\hat{\theta}_{i+1,j}$. Thus, in order to find $\hat{\vartheta}_{1,j}$ we need to find its last change-point $i$.

Adding a change-point:

Recall that MQSE has a minimal number of change-points among all valid solutions in the set $\mathcal{H}$ from equation (6) (meaning that all local confidence statements are satisfied). Therefore, when solving sub-problem $j$, in a first step, the algorithm checks whether there is a valid solution for $\hat{\vartheta}_{1,j}$ with $k$ change-points. Note that, as $\hat{\vartheta}_{1,j-1}$ has $k$ change-points, this implies that any valid solution for $\hat{\vartheta}_{1,j}$ must have at least $k$ change-points, too. Further, note that there is a valid solution for $\hat{\vartheta}_{1,j}$ with $k$ change-points if and only if there exists some $i<j$ with valid constant estimator on the interval $[i+1,j]$ such that $\hat{\vartheta}_{1,i}$ has $k-1$ change-points. When no such $i$ exists, this implies that $\hat{\vartheta}_{1,j}$ has exactly $k+1$ change-points. Searching through all indexes $i$ where $\hat{\vartheta}_{1,i}$ has $k$ change-points, the algorithm finds the optimal location, in terms of the statistic $D$ in (34), for $\hat{\vartheta}_{1,j}$’s last change-point. The indexes $R_{k}$ (see pseudocode) keep track of those candidate locations that need to be considered when there are $k+1$ change-points, see next paragraph.

Confidence intervals for change-points:

While iterating over $j$ (and thereby dynamically increasing $k$), we keep track of indices $R_{k}$ and $L_{k}$. $R_{k}\in\{1,\ldots,n\}$ is such that $\hat{\vartheta}_{1,R_{k}}$ has $k$ change-points and $\hat{\vartheta}_{1,R_{k}+1}$ has $k+1$ change-points, i.e., $R_{k}$ is the greatest index $j$ such that $\hat{\vartheta}_{1,j}$ has $k$ changes. Similarly, $L_{k}\in\{1,\ldots,n\}$ is such that $\hat{\vartheta}_{L_{k},R_{k}}$ has no change-point and $\hat{\vartheta}_{L_{k}+1,R_{k}}$ has one change-point, i.e., $L_{k}$ is the smallest index such that $\hat{\vartheta}_{L_{k},R_{k}}$ is constant. Note that this implies that for any $\hat{\vartheta}\in\mathcal{H}$ the location $\hat{\tau}_{k}$ of its $k$-th change-point is contained in the interval $[L_{k},R_{k}]$. As the true underlying change-point function is contained in $\mathcal{H}(q)$ with probability at least $1-\alpha(q)+{o}(1)$ (recall Theorem 2.2), this implies that $[L_{k},R_{k}]$ is a valid confidence interval for its $k$th change-point location.

Double-heap structures:

In order to check whether a valid constant solution exists on a given interval, one has to intersection the confidence boxes from equation (26) of the main text. This can be done efficiently by updating the $\overline{m}_{1,j-i+1}$ and $\underline{m}_{1,j-i+1}$ quantiles for observations $Z_{i},\ldots,Z_{j}$ from the respective quantiles for $Z_{i-1},\ldots,Z_{j-1}$. To implement this dynamic quantile update, we use double-heaps ($DH$) as described in (Astola and Campbell, 1989). Double-heaps are graphs that consist of two ordered binary trees, one decreasingly and one increasingly, with a common root (see Figure D.2 for an example). For a given quantile $\beta\in(0,1)$ and observations $Z_{1},\ldots,Z_{n}$, whenever we have the observations stored in a double-heap with the right proportion of nodes stored below and above the root, we immediately get that the root is a $\beta$-quantile. For example, if we have the same number of data points stored above and below the root (as in Figure D.2), the root is always the median. The benefit of working with double-heaps, is that they can be updated in $\log(n)$ time, since to maintain the ordering structure by adding (or replacing) one data point, one only has to go through a single branch.

Intersecting confidence boxes:

As detailed in Section 3 of the main text, the local tests can be inverted into intervals as in (26). Recall that their bounds correspond to the $\overline{m}_{j-i+1}$ and $\underline{m}_{j-i+1}$ quantiles of $Z_{i},\ldots,Z_{j}$. Therefore, we can store $2n$ double heaps (denotes as $\overline{DH}_{l}$ and $\underline{DH}_{l}$, $l=1,\ldots,n$) for the upper and lower bounds of the intervals, updating them in each step by replacing one element. In other words, to find the bounds for the interval $[i,j]$, we use the $(j-i+1)$ double-heap ($DH_{j-i+1}$), that has previously stored the bounds for $[i-1,j-1]$ and update it accordingly. Then, we can intersect this bounds recursively to check whether a constant valid solution exists.

Updating the candidates and cost:

An additional double-heap structure $DH_{\beta}$ allows us to efficiently update the $\beta$-quantile for each interval $[i,j]$, and therefore, calculate the candidate function with $k$-th change-point $i$ (see Figure D.1 blue line). If the $\beta$-quantile is inside the multiscale bound (i.e., if it is in the intersection of all local intervals) we choose it as $\hat{\theta}_{i+1,j}$, otherwise we choose either the upper or lower bound of the intersection. The cost of (36) for $i$ can also be calculated dynamically. The two quantities needed to calculate the statistic $D$ in (34) are the sum of transformed data and the number of runs. Both can be updated efficiently in this dynamic programming scheme.

Confidence bands:

For the constant segments $[R_{k}+1,L_{k+1}-1]$, the confidence bands are given by the intersection of all its local confidence boxes in (26). For the confidence intervals $[R_{k},L_{k}]$, the band is calculated as the union over all possible change-points $i$ of the corresponding intersected local boxes in (26).

Computation time and space:

The worst case computation time of this algorithm is $\mathcal{O}(n^{2}\log n)$, and this time is achieved with a constant estimator. In many cases the computation time is smaller, in particular when there are $O(n)$ jumps, in which case the computation time is of order $\mathcal{O}(n\log n)$. Compared to the SMUCE algorithm in Frick et al. (2014) this is an increase of order $\log(n)$, which is due to the additional $\log(n)$ cost of quantile updates. We also note that the required space of MQS is bigger than for SUMCE, since here we need to store double heap structures for each of the $2n+1$ quantiles that need to be calculated. SMUCE needs $\mathcal{O}(n)$ space, while MQS requires $\mathcal{O}(n^{2})$.

D.1 The Multiscale Segment Boxplot

When calculating the Multiscale Segment Boxplot (MSB) (0.25-, 0.5-, and 0.75-quantiles) and two or more of the quantiles have an estimated change-point at a similar location, we might prefer an estimator that chooses the same location, instead of slightly different locations. This will also avoid, with high probability, intersections in the three different quantiles which eases interpretation (see e.g. (He, 1997; Chernozhukov et al., 2010)). Therefore, for the MSB we will modify the estimated MQSE change-point locations whenever the changes in different quantiles are sufficiently close. Thereby, we ensure that the resulting estimator is still an element of the set $\mathcal{H}$ in (6) for each quantile, thus all theoretical guarantees given in Section 2 in the main text still hold. Essentially, whenever the confidence intervals for change-points of two quantiles overlap, we select a common change-point for the MSB.

To do this, we first run the MQS algorithm to obtain the confidence statements, recall Section D. To ensure locally simultaneous changes, if the confidence intervals of two or three quantiles intersect, we choose a common change-point inside the intersection. The particular choice we make is the mean of the corresponding MQS change-point estimators if it is inside the intersection, otherwise we choose the upper or lower limit of the interval. Then, after choosing the change-point, the estimated signal value is the empirical quantile of the delimited segment whenever this is contained in the confidence band. Otherwise we choose the upper or lower limit of the band.

By restricting to locally simultaneous changes, one typically also avoids the problem of quantile-crossing, that is, the estimated quantile function for some quantile $\beta$ exceeds in some region the estimate for some quantile $\beta^{\prime}>\beta$. More generally, when estimating $k$ quantile functions simultaneously, using a Bonferroni correction, we obtain simultaneous coverage of the underlying quantile functions by our confidence bands with probability at least $1-k\alpha$. As the underlying quantile functions do not cross, by making an appropriate choice within the confidence band, quantile crossing happens with probability at most $k\alpha$. In practice, we found that by making these adjustment quantile crossing for the MSB has a very low probability of happening. In 1,000 Monte Carlo simulations of MSB with $\alpha=0.05$ for the Gaussian additive model as in Figure G.2, we did not find any quantile crossing.

Appendix E Simulation results on confidence statements

MQS also provides confidence statements for the number of segments, the change-point locations, and finally the underlying segment function itself (see Theorem 2.2). In this section, we investigate the behavior of such confidence statements in the additive error models as considered in Section 4.1 in the main text, see Figure G.2, with $n=2485$.

Simulation results are shown in Table E.1. As demonstrated in Section 4.1 in the main text, the number of segments is estimated correctly with high probability for $\beta=0.5,0.75$ and slightly underestimated for $\beta=0.25$. Our method gives the theoretical guarantee that $P(\hat{S}\leq S)\geq 1-\alpha$ (see (8)). In Table E.1, we give the frequency of $\hat{S}\leq S$. For all considered examples, (8) is fulfilled for a much larger level than required. Moreover, given that the number of segments is estimated correctly, we construct simultaneous confidence intervals for the true change-point locations. In Table E.1 column ”CI cov.” (short for Confidence Interval coverage), we give the frequency that all the change-points $\tau_{s}$ were inside the respective confidence intervals (given that $\hat{S}=S$). We find that with very high probability, all the change-points are covered by the confidence intervals. We also provide confidence bands that cover the true function with probability at least $1-\alpha$ (see Theorem 2.2), given that the number of change-points is estimated correctly. In Table E.1 column ”CB cov.” (short for Confidence Band coverage), we give the frequency that the true function $\vartheta$ was fully inside the confidence band (given that $\hat{S}=S$). Here, we see that the nominal level is kept, but the confidence bands are, in general, too conservative.

Appendix F MQS under serial dependence

The only assumption that is required for our theoretical results to be valid is independence of the observations as specified in the QSR model (2). In this section we explore empirically the MQS performance when this assumption is violated. To this end, we consider an additive regression model with autoregressive noise. More specifically, for $n=500$ we consider observations of the form $Y_{i}=\mu_{i}+\sqrt{1-\theta^{2}}X_{i}$, where $\mu$ is a simple bump function with its three segments at $[1,125]$, $[126,375]$ and $[376,500]$ taking values [math], $1$, [math], respectively, see Figure F.1. For $X_{i}$ we consider an AR(1) process, that is $X_{i+1}=\theta X_{i}+\epsilon_{t}$ with $\epsilon_{t}\sim\mathcal{N}(0,1)$, for different values of $\theta$ with $|\theta|<1$. The autocorrelation function of the process $X_{i}$ is then given by $\rho_{h}=\theta^{h}$, where $h$ represents the lag. Thereby, positive values of $\theta$ correspond to a positive autocorrelation, negative values of $\theta$ to negative autocorrelation, and $\theta=0$ corresponds to the independent case. Note that since $X_{i}$ has variance $1/(1-\theta^{2})$ we multiply the error by the variance stabilizing scaling $\sqrt{1-\theta^{2}}$.

Figure F.1 shows examples for three different values of $\theta$. Figure F.2 illustrates estimation performance as the parameter $\theta$ deviates from 0 for the example described above. For small $\theta$ values (i.e., small auto-correlation), MQS performs comparably to the independent case (red line) in terms of number of segments, MISE, and V-measure. We highlight that for negative $\theta$ the probability of estimating the correct number of segments can be even higher than in the independent case (see upper-left corner of Figure F.2). This is consistent with recent general results on minimax detection boundaries for bump detection with stationary Gaussian error Enikeeva et al. (2020). When $\theta$ gets very large the strong postive autocorrelation results in overestimation of the number of segments for the MQS (see upper-right corner of Figure F.2).

Appendix G Additional Tables and Figures