A Nonparametric Bayesian Approach for Sparse Sequence Estimation

TL;DR

This paper introduces a nonparametric Bayesian method using Gaussian mixtures for sparse sequence estimation, achieving optimal asymptotic behavior and superior empirical performance with a variational algorithm.

Contribution

It proposes a novel Gaussian mixture prior for sparse sequence estimation and a deterministic variational algorithm for efficient posterior computation.

Findings

Posterior concentrates at minimax rate

Gaussian mixtures outperform pure Gaussian priors

Algorithm shows superior results on benchmark data

Abstract

A nonparametric Bayes approach is proposed for the problem of estimating a sparse sequence based on Gaussian random variables. We adopt the popular two-group prior with one component being a point mass at zero, and the other component being a mixture of Gaussian distributions. Although the Gaussian family has been shown to be suboptimal for this problem, we find that Gaussian mixtures, with a proper choice on the means and mixing weights, have the desired asymptotic behavior, e.g., the corresponding posterior concentrates on balls with the desired minimax rate. To achieve computation efficiency, we propose to obtain the posterior distribution using a deterministic variational algorithm. Empirical studies on several benchmark data sets demonstrate the superior performance of the proposed algorithm compared to other alternatives.

Tables

Table 1. Table 1: MSE of Simulation Study 1, error of the best method marked as bold

$s_n$	10				20				40				80
${\mu }_0$	1	3	5	7	1	3	5	7	1	3	5	7	1	3	5	7
EBMW	10	54	21	13	20	96	35	25	40	152	61	49	79	234	108	96
EBKM	12	35	15	6	21	53	19	7	32	74	25	7	44	93	30	8
EBmed	10	43	22	14	20	69	36	28	37	103	66	54	64	157	127	107
SURE	14	42	45	43	23	68	69	69	42	104	105	105	74	152	153	153
Soft Thresholding	10	76	116	116	20	153	228	233	40	304	457	464	80	609	916	929
Hard Thresholding	13	61	21	12	24	122	39	23	45	237	75	42	87	473	149	82
FDR $q=0.01$	10	75	25	11	20	143	40	22	41	253	67	44	81	434	112	85
FDR $q=0.1$	11	55	24	20	23	93	39	37	44	141	67	64	88	208	113	111
FDR $q=0.4$	22	63	53	51	36	94	82	80	64	134	119	117	119	175	161	161
DP	11	37	11	3	19	50	17	4	33	71	22	4	46	92	26	6

Table 2. Table 2: MAE of Simulation Study 1, error of the best method marked as bold

$s_n$	10				20				40				80
${\mu }_0$	1	3	5	7	1	3	5	7	1	3	5	7	1	3	5	7
EBMW	10	23	14	13	20	43	26	24	40	74	47	45	80	125	86	83
EBKM	27	35	22	19	40	48	25	20	60	62	29	21	82	76	32	22
EBmed	12	21	12	10	23	38	23	20	44	72	45	39	80	133	96	80
SURE	15	30	32	31	24	49	51	51	44	78	79	79	78	117	119	119
Soft Thresholding	10	27	33	33	20	54	65	65	40	109	130	130	80	217	260	261
Hard Thresholding	11	22	10	9	21	45	19	17	41	88	37	32	82	175	74	64
FDR $q=0.01$	10	27	10	8	20	51	19	17	40	93	36	33	80	163	69	66
FDR $q=0.1$	10	21	12	11	21	36	22	22	41	59	41	41	82	98	77	76
FDR $q=0.4$	14	26	25	25	25	44	43	42	48	72	70	70	94	110	108	108
DP	23	31	18	14	36	42	20	16	61	57	25	17	87	72	27	19

Table 3. Table 3: MSE of Simulation Study 2, error of the best method marked as bold

$s_n$	25			50			100
${\mu }_0$	3	4	5	3	4	5	3	4	5
EBMW	138	99	52	233	156	90	385	248	153
EBKM	81	57	28	120	80	41	174	114	52
EBMed	107	78	50	165	124	91	254	208	163
SURE	102	106	105	163	168	170	257	260	259
ST	201	289	327	403	577	658	806	1155	1308
HT	172	143	64	341	282	129	680	563	251
FDR $q=0.01$	192	151	62	356	241	101	639	385	164
FDR $q=0.1$	139	87	55	224	135	99	355	213	167
FDR $q=0.4$	150	133	126	229	206	201	332	302	294
DP	80	55	25	119	79	35	171	109	49

Table 4. Table 4: MAE of Simulation Study 2, error of the best method marked as bold

$s_n$	25			50			100
${\mu }_0$	3	4	5	3	4	5	3	4	5
EBMW	58	47	36	105	82	66	186	143	118
EBKM	75	54	38	102	67	45	139	84	50
EBMed	49	37	29	88	69	58	176	135	113
SURE	73	76	75	120	124	124	194	196	196
ST	70	83	87	141	166	175	281	332	349
HT	62	44	26	123	86	52	245	172	102
FDR $q=0.01$	67	46	26	127	77	48	233	134	90
FDR $q=0.1$	52	34	29	87	61	56	148	111	103
FDR $q=0.4$	63	61	60	107	106	106	179	175	174
DP	60	42	29	93	58	34	128	74	43

Table 5. Table 5: MSE of Simulation Study 3, error of the best method marked as bold

$s_n$	25			50			100
${\mu }_0$	3	4	5	3	4	5	3	4	5
EBMW	110	91	62	186	150	103	308	245	176
EBKM	77	68	49	121	105	77	185	163	124
EBMed	93	79	59	151	130	100	236	212	176
SURE	96	104	106	155	165	168	243	257	260
ST	196	269	313	389	534	626	780	1075	1251
HT	135	123	79	263	237	155	527	479	311
FDR $q=0.01$	151	129	79	268	216	133	479	368	221
FDR $q=0.1$	112	90	66	187	148	110	300	235	187
FDR $q=0.4$	138	137	129	216	213	201	320	307	301
DP	76	68	53	120	110	84	189	164	123

Table 6. Table 6: MAE of Simulation Study 3, error of the best method marked as bold

$s_n$	25			50			100
${\mu }_0$	3	4	5	3	4	5	3	4	5
EBMW	51	46	38	92	82	69	164	145	126
EBKM	75	63	52	108	91	74	160	135	111
EBMed	44	38	31	79	70	59	147	132	115
SURE	69	74	76	114	122	124	187	195	197
ST	67	79	86	133	158	171	266	317	342
HT	51	43	31	100	84	61	200	169	121
FDR $q=0.01$	55	44	31	101	80	57	188	144	104
FDR $q=0.1$	45	37	31	81	67	58	143	120	108
FDR $q=0.4$	60	62	60	105	107	105	176	175	176
DP	59	54	44	97	84	69	156	127	102

Table 7. Table 7: MSE of Simulation 4, error of the best method is marked as bold

$A$	2	3	4	5	6	7
EBMW	320	421	291	175	134	125
EBKM	206	223	150	79	46	35
EBMed	337	350	240	173	148	138
SURE	278	327	333	337	336	334
ST	504	894	1259	1436	1475	1483
HT	375	671	610	301	134	102
FDR $q=0.01$	375	633	441	199	121	111
FDR $q=0.1$	373	421	259	194	185	181
FDR $q=0.4$	440	451	404	403	399	396
DP	204	220	161	205	151	85

Table 8. Table 8: MAE of Simulation 4, error of the best method marked as bold

$A$	2	3	4	5	6	7
EBMW	179	199	159	130	122	120
EBKM	225	182	115	79	61	61
EBMed	176	161	127	110	101	97
SURE	209	241	245	248	246	247
ST	216	294	347	367	371	372
HT	189	242	184	110	85	80
FDR $q=0.01$	189	232	147	96	85	83
FDR $q=0.1$	188	168	120	110	109	108
FDR $q=0.4$	217	214	208	212	210	209
DP	215	161	91	89	75	57

Equations

$$X_i={\theta }_i+e_i,i=1,2,\dots ,n,$$

$$\Theta_{0}(s_{n})=\Big{\{}\boldsymbol{\theta}\in\mathbb{R}^{n}:\#\{1\leq i\leq n:\theta_{i}\neq 0\}=s_{n}\Big{\}}.$$

$$\pi ({\theta }_i)=w{\delta }_0({\theta }_i)+(1-w)g({\theta }_i).$$

$$\pi ({\theta }_i)=w{\delta }_0({\theta }_i)+(1-w)\sum _{t=1}^T{w}_t\text{N}(m_t,{\sigma }^2),$$

$$w$$

$${\theta }_i\mid w$$

$$Q_n(A)=\frac{{\int }_A\{p_{\mathbf{\theta }}(X^n)/p_{{\mathbf{\theta }}^{\ast }}(X^n){\}}^{\kappa }{\Pi }_{X^n}(d\mathbf{\theta })}{{\int }_{{\mathbb{R}}^n}\{p_{\mathbf{\theta }}(X^n)/p_{{\mathbf{\theta }}^{\ast }}(X^n){\}}^{\kappa }{\Pi }_{X^n}(d\mathbf{\theta })}.$$

$$A_{M{\epsilon }_n}=\{\mathbf{\theta }\in {\mathbb{R}}^n:\Vert \mathbf{\theta }-{\mathbf{\theta }}^{\ast }{\Vert }^2>M{\epsilon }_n\}.$$

$${\theta }_i$$

$$X_i$$

$$G_0=w_0{\delta }_0+(1-w_0)\text{N}(0,{\sigma }_0^2);$$

$$G=\sum _{t=1}^{\infty }{\pi }_t({\xi }_t{\delta }_0(\cdot )+(1-{\xi }_t){\delta }_{{\eta }_t^{\ast }}(\cdot ))=\sum _{t=1}^{\infty }{\pi }_t{\xi }_t{\delta }_0(\cdot )+\sum _{t=1}^{\infty }{\pi }_t(1-{\xi }_t){\delta }_{{\eta }_t^{\ast }}(\cdot )$$

$$\equiv w{\delta }_0(\cdot )+(1-w)\sum _{t=1}^{\infty }{w}_t{\delta }_{{\eta }_t^{\ast }}(\cdot ),$$

$$V_t|{\alpha }_0\sim \text{Beta}(1,{\alpha }_0),t=1,2,\cdots ,T-1,V_T=1;$$

$${\xi }_t\sim \text{Ber}(w_0),t=1,2,\cdots ,T;$$

$${\eta }_t^{\ast }|{\xi }_t\sim \{\begin{array}{rlr}& {\delta }_0& {\xi }_t=1\\ & \text{N}(0,{\sigma }_0^2)& {\xi }_t=0\end{array};t=1,2,\cdots ,T;$$

$${\pi }_t=V_t\prod _{j=1}^{t-1}(1-V_j),t=1,2,\cdots ,T-1,{\pi }_T=\prod _{j=1}^T(1-V_j);$$

$$Z_k|\{V_1,V_2,\cdots ,V_{T-1}\}\sim \text{Multinomial}(\mathbf{\pi });$$

$$X_i|Z_i\sim \text{N}({\eta }_{Z_i}^{\ast },1),i=1,2,\cdots ,n.$$

$$q(\mathbf{Z},\mathbf{V},\mathbf{\eta },\mathbf{\xi })=q_{\mathbf{p},\mathbf{m},\mathbf{\tau }}(\mathbf{\eta },\mathbf{\xi })q_{{\mathbf{\gamma }}_1,{\mathbf{\gamma }}_2}(\mathbf{V})q_{\mathbf{\Phi }}(\mathbf{Z});$$

$$(\sum _{t=1}^T{\hat{\varphi }}_{it}{\hat{p}}_t){\delta }_0(\cdot )+\sum _{t=1}^T{\hat{\varphi }}_{it}(1-{\hat{p}}_t){\delta }_{m_t}(\dot{)}.$$

$${\hat{G}}_{\text{MAP}}={\tilde{w}}_0\cdot {\delta }_0(\cdot )+\sum _{t=1}^T{\tilde{w}}_t{\delta }_{{\hat{m}}_t}(\cdot ).$$

$${\theta }_i$$

$$X_i$$

$$\frac{{\tilde{w}}_0\exp (-\frac{\kappa X_i^2}{2}){\delta }_0+{\sum }_{t=1}^T{\tilde{w}}_t\exp (-\frac{\kappa (X_i-{\hat{m}}_t{)}^2}{2}){\delta }_{{\hat{m}}_t}}{{\tilde{w}}_0\exp (-\frac{\kappa X_i^2}{2})+{\sum }_{t=1}^T{\tilde{w}}_t\exp (-\frac{\kappa (X_i-{\hat{m}}_t{)}^2}{2})}\equiv {\hat{w}}_0\cdot {\delta }_0(\cdot )+\sum _{t=1}^T{\hat{w}}_t{\delta }_{{\hat{m}}_t}(\cdot );$$

$${\hat{\theta }}_i=\sum _{t=1}^T{\hat{w}}_t{\hat{m}}_t.$$

$$D_n={\int }_0^1\prod _{i=1}^n{\int }_{\mathbb{R}}\{\frac{p_{{\theta }_i}(X_i)}{p_{{\theta }_i^{\ast }}(X_i)}{\}}^{\kappa }{\Pi }_{w,\mathbf{w},\mathbf{m}}(d{\theta }_i)\Pi (dw);$$

$$D_n>{\int }_0^1{w}^{n-s_n}(1-w{)}^{s_n}\Pi (dw)\prod _{i\in S^{\ast }}{\int }_{\mathbb{R}}{e}^{\frac{\kappa }{2}\{(X_i-{\theta }_i^{\ast }{)}^2-(X_i-{\theta }_i{)}^2\}}\sum _{t=1}^T\frac{w_t}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{1}{2{\sigma }^2}(m_t-{\theta }_i{)}^2}d{\theta }_i.$$

$${\int }_{\mathbb{R}}{e}^{\frac{\kappa }{2}\{(X_i-{\theta }_i^{\ast }{)}^2-(X_i-{\theta }_i{)}^2\}}\frac{w_t}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{1}{2{\sigma }^2}(m_t-{\theta }_i{)}^2}d{\theta }_i$$

$$\prod _{i\in S^{\ast }}{\int }_{\mathbb{R}}{e}^{\frac{\kappa }{2}\{(X_i-{\theta }_i^{\ast }{)}^2-(X_i-{\theta }_i{)}^2\}}\sum _{t=1}^T\frac{w_t}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{1}{2{\sigma }^2}(m_t-{\theta }_i{)}^2}d{\theta }_i$$

$$\prod _{i\in S^{\ast }}\{\sum _{t=1}^T{w}_t\exp (-\frac{\kappa (X_i-m_t{)}^2}{\kappa {\sigma }^2+1})\}\ge \prod _{i\in S^{\ast }}\{w_{t_i}\exp (-\sum _{t=1}^T\frac{\kappa (X_i-m_{t_i}{)}^2}{\kappa {\sigma }^2+1})\}$$

$$P(|\sum _{i\in S^{\ast }}(X_i-{\theta }_i^{\ast }{)}^2/s_n-1|\ge \log (n/s_n)-1)\le 2\exp (-s_n(\log (n/s_n)-1{)}^2/8).$$

$$\prod _{i\in S^{\ast }}\{\sum _{t=1}^T{w}_t\exp (-\frac{\kappa (X_i-m_t{)}^2}{\kappa {\sigma }^2+1})\}\ge \exp (\log C\cdot s_n-{\epsilon }_n-\frac{\kappa }{\kappa {\sigma }^2+1}{o}_p({\epsilon }_n)).$$

$$\prod _{i\in S^{\ast }}{e}^{\frac{\kappa }{2}\{(X_i-{\theta }_i^{\ast }{)}^2\}}=e^{\frac{\kappa }{2}{\prod }_{i\in S^{\ast }}\{(X_i-{\theta }_i^{\ast }{)}^2\}}=e^{\frac{\kappa s_n}{2}+o_p(s_n)}.$$

$$\prod _{i\in S^{\ast }}{\int }_{\mathbb{R}}{e}^{\frac{\kappa }{2}\{(X_i-{\theta }_i^{\ast }{)}^2-(X_i-{\theta }_i{)}^2\}}\sum _{t=1}^T\frac{w_t}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{1}{2{\sigma }^2}(m_t-{\theta }_i{)}^2}d{\theta }_i\ge (1+\kappa {\sigma }^2{)}^{-\frac{s_n}{2}}\exp (-{\epsilon }_n-o_p({\epsilon }_n)).$$

Full text

Keywords: Empirical Bayes; Variational; Minimax; Posterior Consistency.

1 Introduction

Consider a Gaussian sequence model

[TABLE]

where $\boldsymbol{\theta}=(\theta_{1},\dots,\theta_{n})$ is the unknown mean parameter, often assumed to be sparse, and $e_{i}$’s are independent errors following a standard normal distribution. The primary interest is to reconstruct the sparse vector $\boldsymbol{\theta}$ based on the data $X^{n}=(X_{1},X_{2},\cdots,X_{n})$. Such a simple model arises in many applications, such as astronomy, signal processing and bioinformatics. It is also a canonical model for many modern statistical methods, such as nonparametric function estimation, large-scale variable selection and hypothesis testing.

There are different ways to define the sparsity of a vector. In this paper, we focus on the common definition of sparsity, i.e., many elements of $\theta_{i}$’s are zero. In particular, we assume $\boldsymbol{\theta}$ is from the class of nearly black vectors,

[TABLE]

The parameter $s_{n}$ measures the sparsity of $\boldsymbol{\theta}$, which is usually assumed to be $o(n)$, but unknown. So a desired feature of any reconstruction procedure is to adapt to the unknown sparsity level.

Since $\boldsymbol{\theta}$ is known to be sparse, a natural approach is to threshold $X^{n}$. The thresholding rule can be chosen by the principle of minimizing the empirical fitting error with penalization (Golubev, 2002) or by minimizing the False Discover Rate (Abramovich et al., 2006). In addition, a plethora of research on variable selection and prediction in the context of high dimensional linear regression models can also be applied on this problem(Fan and Lv, 2010; Bühlmann and Van De Geer, 2011; Zhang and Huang, 2008).

Some thresholding procedures are motivated from a Bayesian aspect, such as George and Foster (2000), Johnstone and Silverman (2004), and Castillo et al. (2012). The idea is to model $\theta_{i}$’s with a two-group structure prior, where one component is a point mass at zero due to the sparsity assumption and the other component is a continuous distribution with a smooth symmetric density function $g$, namely,

[TABLE]

Theoretical studies have indicated that $g(\cdot)$, the prior density function on the non-zero component, should have a heavy tail for the posterior distribution to have the desired asymptotic behavior. Therefore, the double exponential distribution and its scaled variants are recommended for the choice of $g$, but not Gaussian distributions. The problem with distributions like the Gaussians which have lighter tails is that they tend to over-shrink $X_{i}$’s for large $\theta_{i}$’s therefore attain a lower posterior contraction rate (Johnstone and Silverman, 2004; Castillo et al., 2012).

In this paper, we propose an Empirical Bayes approach to this problem. Despite the warning message on Gaussian distributions, we still use Gaussians in our prior specification, which takes the following form

[TABLE]

then use data to learn weights $w,w_{1},\cdots,w_{T}$ and centers $m_{1},\cdots,m_{T}$. Our prior choice can be viewed as a special case of the two-group prior (2) with $g$ being a mixture of Gaussian density functions. The Gaussian prior is appealing due to its conjugacy, which simplifies our analysis and also enables tractable computation. As revealed by our asymptotic analysis, the suboptimal behavior of Gaussian distributions as mentioned in Johnstone and Silverman (2004) and Castillo et al. (2012) is on Gaussian distributions with mean zero, which can be avoided by a proper choice on weights $w_{1},\cdots,w_{T}$ and centers $m_{1},\cdots,m_{T}$.

The remainder of this paper is organized as follows. We present our main results in Section 2. We show that with some mild conditions on $m_{t}$’s and $w_{t}$’s, the corresponding posterior will have desired asymptotic behavior: it concentrates around the true parameter $\boldsymbol{\theta}^{*}$ at the minimax rate, its effective dimension adapts to the unknown sparsity $s_{n}$, and the corresponding the posterior mean is asymptotic minimax. A variational implementation of our approach, as well as empirical studies, is presented in Section 3. All the proofs are given at the end after the conclusion.

We will start the remaining of this paper by briefly discussing some related work.

• •

This paper is motivated by a recent work by Martin and Walker (2014), where $g(\theta)=g_{i}(\theta)$ in (2) is set to be a Gaussian distribution centered exactly at the data point $X_{i}$. In Martin and Walker (2014), the resulting estimate of $\theta_{i}$ is still a shrinkage estimate toward zero, while in our approach, the estimate of $\theta_{i}$ is adaptively shrunk toward the mean of nearby data points. This explains why the empirical performance of our approach is better than the one from Martin and Walker (2014).

• •

Estimating $\boldsymbol{\theta}$ for the Gaussian sequence model (1) can be also viewed as a compound decision problem (Robbins, 1951), where $\theta_{i}$’s are assumed to be i.i.d. random variables from a common but unknown distribution $G$. The aforementioned two-group prior (2) can be viewed as a special form of $G$, but in general $G$ can take any form. A number of nonparametric approaches have been proposed to estimate $G$ and then in return provide an estimate for $\boldsymbol{\theta}$, such as the Maximum likelihood approach by Jiang and Zhang (2009), the nonparametric empirical Bayes approach by Brown and Greenshtein (2009), and a convex optimization based approach by Koenker (2014). In particular, Gaussian mixtures are used to estimate $G$ by Jiang and Zhang (2009) and Brown and Greenshtein (2009). However, their asymptotic results do not cover the nearly black class $\Theta_{0}(s_{n})$. In their simulation study, Jiang and Zhang (2009) indeed consider the sparse situation and suggest to add a Gaussian component with mean zero. In our asymptotic study, we have found that to achieve the desired asymptotic properties, a point mass zero, instead of a Gaussian with mean zero, seems necessary for the nearly black class.

2 Main Results

First we introduce some key notations. Data vector is denoted as $X^{n}=(X_{1},X_{2},\cdots,X_{n})$; $\Pi_{X^{n}}(\cdot)$ is used to denote the prior indicating dependence with data; $\boldsymbol{\theta}^{\ast}=(\theta^{\ast}_{1},\cdots,\theta^{\ast}_{n})$ is the unknown true mean of $X^{n}$; $p_{\boldsymbol{\theta}}(X^{n})$ is the likelihood function of $\boldsymbol{\theta}$ given $X^{n}$. $a_{n}\succeq b_{n}$ is used to represent $\lim_{n}b_{n}/a_{n}=0$.

Consider the following hierarchical prior $\Pi_{X^{n}}(\cdot\mid w_{1:T},m_{1:T},\sigma^{2},\alpha)$ on $\boldsymbol{\theta}$ as follows

[TABLE]

where $\sigma^{2}$, $\alpha$, $w_{1:T}$ and $m_{1:T}$ are fixed parameters. $T$ is independent of the dimension $n$.

The posterior distribution is proportional to $p_{\boldsymbol{\theta}}(X^{n})\Pi_{X^{n}}(\boldsymbol{\theta})$. Strong assumptions are usually needed to prove Bayesian consistency. However, following Walker and Hjort (2001) and Martin and Walker (2014), we can obtain our posterior distribution with fractional likelihood, which can weaken the conditions required for Bayesian consistency. Given a fractional parameter $\kappa$, the corresponding posterior distribution is then proportional to $p_{\boldsymbol{\theta}}^{\kappa}(X^{n})\Pi_{X^{n}}(\boldsymbol{\theta})$. In implementation, $\kappa$ is set to be a large number close to 1 to capture most information from data. In all of the simulation studies in Section 3, we set $\kappa=0.99$. In particular, the posterior measure of a Borel set $A\subset\mathbb{R}^{n}$ involving $\kappa$, denoted by $Q_{n}(A)$, can be expressed as

[TABLE]

$\boldsymbol{\theta}^{\ast}$ is assumed to be sparse. Denote $S^{\ast}$ as the support of $\boldsymbol{\theta}^{\ast}$ and $s_{n}=|S^{\ast}|$ is the cardinality. Minimax rate established in (Donoho et al., 1992) for estimating a vector in $\Theta_{0}(s_{n})$ is $\varepsilon_{n}=s_{n}\log(n/s_{n})$. The aim is to prove posterior mean estimator based on the above prior attains asymptotically minimax $L^{2}$ error rate in $\Theta_{0}(s_{n})$. To accomplish this, first we prove our posterior measure $Q_{n}$ has a desired concentration rate and then prove MSE is bounded.

The key conditions are imposed on $w_{1},\cdots,w_{T}$ and $m_{1},\cdots,m_{T}$. To find suitable $m_{1},\cdots,m_{T}$, we apply suitable clustering algorithms on nonzero mean set $\boldsymbol{\theta}_{S^{\ast}}=\{\theta_{i}^{\ast}|i\in S^{\ast}\}$ and retrieve cluster centers as $m_{1},\cdots,m_{T}$. Condition 1 essentially puts an upper bound on within-cluster sum of squares of $\boldsymbol{\theta}_{S^{\ast}}$ if we apply some clustering algorithm. Denote $m_{t_{i}}$ as the corresponding cluster center of $\theta_{i}^{\ast}\in\boldsymbol{\theta}_{S^{\ast}}$, we assume

Condition 1**.**

$\sum_{i\in S^{\ast}}(\theta_{i}^{\ast}-m_{t_{i}})^{2}=o(\varepsilon_{n})$.

Since the summation in the left hand side has $s_{n}$ terms, if $\max_{i\in S^{\ast}}|\theta_{i}^{\ast}-m_{t_{i}}|=o(\log(n/s_{n}))$, that is, the distance between $\theta_{i}^{\ast}$ and the corresponding cluster centers grows slower than $\log(n/s_{n})$, our condition could be satisfied. Therefore Condition 1 is not strict since $s_{n}=o(n)$. In implementation to specify $m_{1},\cdots,m_{T}$, since we do not know $\boldsymbol{\theta}_{S^{\ast}}$, we use a novel clustering algorithm to identify the cluster centered at 0 and estimate $T$ nonzero cluster centers. The details will be provided in Section 3.

The second condition is on the weight $w_{1},w_{2},\cdots,w_{T}$:

Condition 2**.**

$\min_{1\leq t\leq T}w_{t}\geq Cs_{n}/n>0$.

Since $T$ is independent of $n$. This condition could be easily satisfied if we simply set $w_{t}=1/T$ for each $t$. If each cluster size is bounded below, we could also plug in cluster weights.

Given the above 2 conditions, first we introduce Lemma 1, an analogue of Lemma 1 in Martin and Walker (2014).

Lemma 1**.**

Let $D_{n}$ be the denominator of posterior measure $Q_{n}$. If Condition 1 and Condition 2 hold, then $D_{n}>\frac{\alpha}{1+\alpha}\exp\{-2\varepsilon_{n}-o(\varepsilon_{n})\}$ with $P_{\boldsymbol{\theta}^{\ast}}-$probability 1.

After establishing a lower bound for the denominator of $Q_{n}$, in order to prove posterior concentration, we’ll establish an upper bound for that numerator of $Q_{n}$ measuring the complement of a ball centered at $\boldsymbol{\theta}^{\ast}$. Specifically, we show that posterior probability measure $Q_{n}$ concentrates asymptotically on a ball centered at the truth $\boldsymbol{\theta}^{\ast}$ with square radius proportional to $\varepsilon_{n}$, namely,

[TABLE]

Theorem 1**.**

If Condition 1 and Condition 2 hold, then there exists $M>0$ such that $Q_{n}(A_{M\varepsilon_{n}})\rightarrow 0$ with $P_{\boldsymbol{\theta}^{\ast}}$-probability 1 as $n,s_{n}\to\infty$ with $s_{n}=o(n)$.

Theorem 1 implies that our posterior distribution concentrates around the right place at the right rate, so it ought to produce an estimator of $\boldsymbol{\theta}$ with good properties. Next we show that the posterior mean $\hat{\boldsymbol{\theta}}$ is a minimax estimator.

Theorem 2**.**

If Condition 1 and Condition 2 hold, there exists a universal constant $M^{\prime}>0$, such that $\mathbb{E}_{\boldsymbol{\theta}^{\ast}}\|\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{\ast}\|^{2}\leq M^{\prime}\varepsilon_{n}$ for all large $n$.

The proof is the same as the one by Martin and Walker (2014), provided that we have proved Lemma 1 and Theorem 1.

As posterior concentrates around $\boldsymbol{\theta}^{\ast}$ at the minimax rate, we could conclude the majority of the posterior mass concentrates on $s_{n}$-dimensional subspaces of $\mathbb{R}^{n}$. The sparsity level of posterior is measured by the posterior distribution of $w$. Theorem 3 shows posterior distribution of $w$ put large probability mass above $1-s_{n}n^{-1}$.

Theorem 3**.**

If Condition 1 and Condition 2 hold, let $\delta_{n}=K\varepsilon_{n}n^{-1}$, and $K>0$ is a suitably large constant. Then $E_{\boldsymbol{\theta}^{\ast}}\{P(1-w>\delta_{n}|X^{n})\}\rightarrow 0$ as $n\rightarrow\infty$.

Notice that $\delta_{n}=Ks_{n}/n\log(n/s_{n})\approx Ks_{n}/n$, therefore the posterior distribution of $w$ concentrates around $1-s_{n}n^{-1}$. That is, the effective dimension of our posterior distribution adapts to the unknown sparsity level $s_{n}$. Our posterior mean estimator could detect true sparsity pattern with large probability. The proof is the same as the one by Martin and Walker (2014) as long as we have proved Lemma 1 and Theorem 1.

3 Implementation

In Section 2 we have shown using good prior with parameters $(w_{1:T},m_{1:T},\sigma^{2},\alpha)$ could result in posterior concentration and minimax rate of Bayes posterior mean estimator. In implementation, we first specify $\sigma^{2}$ and $\alpha$, then estimate $(w_{1:T},m_{1:T})$ using $X^{n}$. In Condition 1 $\sum_{i\in S^{\ast}}(\theta_{i}^{\ast}-m_{t_{i}})^{2}$ needs to be bounded, however, in practice $\boldsymbol{\theta}_{S^{\ast}}$ is unknown so that we cannot apply clustering algorithms on it. One remedy is to apply the clustering algorithm on $X^{n}$ since $\sum_{i\in S^{\ast}}(\theta_{i}^{\ast}-m_{t_{i}})^{2}\leq 2\sum_{i\in S^{\ast}}(X_{i}-m_{t_{i}})^{2}+2\sum_{i\in S^{\ast}}(X_{i}-\theta_{i}^{\ast})^{2}$. $\sum_{i\in S^{\ast}}(X_{i}-\theta_{i}^{\ast})^{2}=o_{p}(\varepsilon_{n})$, therefore $\sum_{i\in S^{\ast}}(\theta_{i}^{\ast}-m_{t_{i}})^{2}=o(\varepsilon_{n})$ if $\sum_{i\in S^{\ast}}(X_{i}-\theta_{i}^{\ast})^{2}=o_{p}(\varepsilon_{n})$.

we apply some proper clustering algorithm on $X^{n}$ to assign clustering centers to $m_{1:T}$ and cluster weights to $w_{1:T}$. We need to pre-specify 0 as one clustering center so that the major task is to estimate the nonzero cluster cneters. In Bayesian Framework, Dirichlet process mixture model is used to do clustering. We build a Dirichlet process (DP) mixture model on $X^{n}$ and estimate $m_{t}$ by plugging in corresponding cluster centers. The general framework is summarized as follows:

[TABLE]

where $G_{0}$ is the base measure and $\alpha_{0}$ is the concentration measure. Given $G_{0}$ and $\alpha_{0}$, we have stick breaking representation of $G$ as $\sum^{\infty}_{t=1}\pi_{t}\delta_{\eta_{t}}(\cdot)$, where $\eta_{t}$ is drawn independently and identically distributed (i.i.d.) from base measure $G_{0}$, while $\pi_{t}=V_{t}\prod^{t-1}_{l=1}(1-V_{l})$. $V_{l}$ is drawn i.i.d. from $\text{Beta}(1,\alpha)$. From this formulation we could see that random distribution function $G$ is almost surely discrete. Since $\boldsymbol{\theta}^{\ast}$ is a sparse vector, in order to get a random distribution $G$ with a positive point mass at 0, $G_{0}$ ought to have a positive mass at 0. Therefore we model $G_{0}$ as a normal component with a point mass at 0:

[TABLE]

where $w_{0}$ and $\sigma_{0}^{2}$ are 2 pre-specified parameters. Then $G$ could be written as

[TABLE]

where $w=\sum^{\infty}_{t=1}\pi_{t}\xi_{t}$ and $w_{t}=\frac{\pi_{t}(1-\xi_{t})}{1-w}$. In this formulation $\xi_{t}$ is drawn i.i.d. from $\textsf{Ber}(w)$. If $\xi_{t}=0$, then $\eta_{t}^{\ast}$ is drawn from $\textsf{N}(0,\sigma_{0}^{2})$, otherwise $\eta_{t}^{\ast}=0$. Since $w>0$, $G$ always has a positive probability mass at 0 which could induce certain level of sparsity. In the following subsection, we’ll develop a variational Algorithm to estimate $G$.

3.1 Variational Algorithm for Estimating Prior

Once we have specified Dirichlet Process as the prior, we need to compute the posterior distribution of $G$ in order to calculate Maximum A Posteriori (MAP) estimator of $G$. The major challenge here is that we have the infinite sum in the stick breaking form of $G$, which makes it impossible to sample from $G$. One remedy here is to fix $T$ as the upper bound of the number of clusters (Blei et al., 2006). Then we have the following truncated version of stick breaking process using $G_{0}=w_{0}\delta_{0}+(1-w_{0})\textsf{N}(0,\sigma_{0}^{2})$ as the base measure.

[TABLE]

The observed data are $\mathbf{y}$ and the parameters are $\mathbf{Z}_{1\times n},\mathbf{V}_{1\times(T-1)},\boldsymbol{\eta}^{\ast}_{1\times T},\boldsymbol{\xi}_{1\times T}$. $\boldsymbol{\eta}^{\ast}=(\eta_{1}^{\ast},\cdots,\eta_{T}^{\ast})$ contains all unique values of $\boldsymbol{\eta}=(\eta^{\ast}_{Z_{i}})^{n}_{i=1}$.

We could use Markov Chain Monte Carlo (MCMC) (Escobar and West, 1995) to compute the posterior distribution of these parameters. However, due to large $n$, the number of parameters is huge, making MCMC converging very slowly. In this paper, we use a variational algorithm to get posterior distribution. Variational Bayes methods are deterministic. The essence of variational inference is to regard the computation of posterior distribution as an optimization problem. Solving this optimization problem gives an approximation to the posterior distribution. Blei et al. (2006) proposed variational algorithms for Dirichlet process mixture models for exponential family. Although normal distribution with a positive mass at 0 does not belong to exponential family, we could follow the same philosophy to deduce the corresponding algorithm.

We consider mean field variational inference and assume the following fully factorized variational distribution:

[TABLE]

Through the calculation shown in the Appendix, we could prove the optimal $q$ must be further factorized as follows:

• •

$q_{\mathbf{p},\mathbf{m},\boldsymbol{\tau}}(\boldsymbol{\eta}^{\ast},\boldsymbol{\xi})=\prod^{T}_{t=1}q_{p_{t},m_{t},\tau_{t}}(\eta^{\ast}_{t},\xi_{t})$, where $\mathbf{p}=(p_{1},p_{2},\cdots,p_{T})$, $\mathbf{m}=(m_{1},m_{2},\cdots,m_{T})$, $\boldsymbol{\tau}=(\tau_{1},\tau_{2},\cdots,\tau_{T})$, and $q_{p_{t},m_{t},\tau_{t}}(\eta^{\ast}_{t},\xi_{t})=p_{t}1_{\xi_{t}=1}\delta_{0}+(1-p_{t})1_{\xi_{t}=0}q_{m_{t},\tau_{t}}(\eta^{\ast}_{t})$, where $q_{m_{t},\tau_{t}}(\eta^{\ast}_{t})$ is Normal density with mean $m_{t}$ and variance $\tau_{t}^{2}$.

• •

$q_{\boldsymbol{\gamma}_{1},\boldsymbol{\gamma}_{2}}(\mathbf{V})=\prod^{T-1}_{t=1}q_{\gamma_{1t},\gamma_{2t}}(V_{t})$, where $\boldsymbol{\gamma}_{1}=(\gamma_{11},\gamma_{12},\cdots,\gamma_{1(T-1)})$, $\boldsymbol{\gamma}_{2}=(\gamma_{21},\gamma_{22},\cdots,\\ \gamma_{2(T-1)})$, $q_{\gamma_{1t},\gamma_{2t}}(V_{t})$ is Beta Distribution with parameters $(\gamma_{1t},\gamma_{2t})$.

• •

$q_{\boldsymbol{\Phi}}(\mathbf{Z})=\prod^{n}_{i=1}q_{\boldsymbol{\phi}_{i}}(Z_{i})$; where $\mathbf{Z}=(Z_{1},Z_{2},\cdots,Z_{n})$, $\boldsymbol{\Phi}=(\boldsymbol{\phi}_{1},\boldsymbol{\phi}_{2},\cdots,\boldsymbol{\phi}_{n})$,$\boldsymbol{\phi}_{i}=(\phi_{i,1},\phi_{i,2},\cdots,\phi_{i,T})$ , $\phi_{i,t}=q(Z_{i}=t)$, $q_{\boldsymbol{\phi}_{i}}(Z_{i})$ is Multinomial distribution with parameters $\boldsymbol{\phi}_{i}$.

The algorithm is summarized in Algorithm 1. Via iterating these steps we could update the variational parameters. After convergence of $\boldsymbol{\Phi}$, $\mathbf{p}$, $\mathbf{m}$, $\boldsymbol{\tau}$, $\boldsymbol{\gamma}_{1}$ and $\boldsymbol{\gamma}_{2}$, we get an approximation of the posterior by plugging in these estimated parameters. The parameters we are interested in are $\boldsymbol{\Phi}$,$\mathbf{p}$ and $\mathbf{m}$. (In the algorithm $\|\cdot\|_{\infty,\infty}$ means the element-wise maximum absolute value; $\text{logit}(x)=(1+\exp(-x))^{-1}$.)

3.2 Posterior Computation

Given approximate posterior estimates $\hat{\mathbf{p}}$, $\hat{\mathbf{m}}$, $\hat{\boldsymbol{\Phi}}$, we construct a MAP estimator of $G$. Remind that $\hat{m}_{t}$ is the nonzero cluster center; $\hat{p}_{t}$ is the probability mass of zero of component indexed by $t$; each entry $\hat{\phi}_{it}$ of $\hat{\boldsymbol{\Phi}}$ is the posterior probability of $Z_{i}$ belonging to the cluster $t$. The approximate posterior distribution of $\eta^{\ast}_{Z_{i}}$ is

[TABLE]

The most probable posterior assignment of $\eta^{\ast}_{Z_{i}}$ based on above posterior is denoted as $\hat{\eta}^{\ast}_{Z_{i}}$. MAP estimate of cluster weights including zero clusters is $\tilde{w}_{t}=\#\{k:\hat{\eta}^{\ast}_{Z_{i}}=\hat{m}_{t}\}/n$ and $\tilde{w}_{0}=\#\{k:\hat{\eta}^{\ast}_{Z_{i}}=0\}/n$. Then the estimated prior is

[TABLE]

Plugging in $\hat{G}_{\text{MAP}}$, data generation process we are considering could be written as ($\textsf{N}^{\kappa}$ denotes normal likelihood to the power of $\kappa$):

[TABLE]

This is different from the original data generating process 4 but could be regarded as an approximation of 4 by setting $\sigma^{2}=0$. The reason why we use this approximation is that instead of using MCMC, we could compute a closed-form posterior distribution of $\boldsymbol{\theta}$ given $X^{n}$ using Bayes Formula: since $X_{i}\sim\textsf{N}(\theta_{i},1)$ and $\theta_{i}\sim\hat{G}_{\text{MAP}}$, the posterior distribution of $\theta_{i}$ given $X^{n}$ has the form

[TABLE]

where $\hat{w}_{t}$ is the posterior weight and $\sum^{T}_{t=1}\hat{w}_{t}=1$. The final estimate of $\theta_{i}$ is the posterior mean estimator

[TABLE]

3.3 Simulation Studies

We conduct the following four simulation studies; R package for our empirical Bayes estimator is available in https://github.com/yunboouyang/VBDP; all the source code is summarized in https://github.com/yunboouyang/EBestimator. For our empirical Bayes estimator, we set the maximal number of clusters to be $T=10$. We set fractional likelihood parameter $\kappa=0.99$ throughout all the simulation study. We set concentration parameter $\alpha_{0}=1$ in all of the simulation studies. For $G_{0}=w_{0}\delta_{0}+(1-w_{0})\textsf{N}(0,\sigma_{0}^{2})$, we set $w_{0}=0.01$. For simulation study 1, we set $\sigma_{0}=4$. For other simulation studies, we set $\sigma_{0}=6$. The influence of these hyper-parameters will diminish when sample size $n$ is large.

Besides our method (denoted as DP representing Dirichlet process clustering), we also consider estimators proposed by Martin and Walker (2014) (denoted as EBMW), by Koenker and Mizera (2014) (denoted as EBKM), and by Johnstone and Silverman (2004)(denoted as EBMed, we directly use functions from the EbayesThresh package by Johnstone and Silverman (2005)), as well as the Hard thresholding estimator, soft thresholding estimator, SURE estimator (using waveThresh package) and FDR estimator with parameters $q=0.01,0.1,4$.

Experiment 1

In the first simulation study, we take sample $X^{n}$ of dimension $n=200$ from the normal mean model $X_{i}\sim\textsf{N}(\theta_{i},1)$. In this case, we consider the number of nonzero elements to be $s_{n}=10,20,40,80$ and the signals are fixed at values $\mu_{0}=1,3,5,7$. We compute mean squared error (MSE) and mean absolute error (MAE) based on 200 replications as two measure of performance. We summarize the results in Table 1 and Table 2.

When $\mu_{0}$ is larger than one, our DP method is the best one in terms of MSE since it is easier to figure out the cluster centers of $\theta_{i}$’s and estimate the nonzero $\theta_{i}$’s. Besides, when $s_{n}$ is large, our DP method has the best performance among all since the more nonzero $\theta_{i}$’s we have, the easier we could estimate the cluster centers. When both $\mu_{0}$ and $\theta_{i}$ are small, which is a tough case since high dimensional noise might make it much more difficult to detect weak signals, our method still has the comparable performance among all methods. Our method has better performance when measured by MSE than that in MAE. Other methods such as EBKM and EBmed also have comparable good performance.

Experiment 2

In the second simulation study we increase the data dimension to 500. We are interested in the case when $s_{n}=25,50,100$ and all nonzero elements are fixed at $\mu_{0}=3,4,5$. This is a more challenging case since the signal is not very strong. MSE and MAE are shown in Table 3 and Table 4 respectively.

In terms of MSE, our DP method is consistently the best one among all the different configurations. The performance of EBKM is comparable. One drawback of EBKM is that the computational cost is high since solving a high dimensional optimization problem is needed. However for DP, we use a computationally efficient variational inference algorithm, which dramatically saves the computational time.

Experiment 3

In the third simulation study we maintain all the features of the second experiment except that non-zero elements of the true mean vector are now generated from standard Gaussian distribution centered at the original values. All different $\theta_{i}$’s form a data cloud which is centered at the distinct non-zero values in simulation study 2.

Even though our DP method is not as good as EBKM and EBMed in some configurations, our DP method has quite similar performance even it’s more difficult to estimate the clustering centers in Simulation Study 3 than Simulation Study 2.

Experiment 4

In the fourth example, we consider $1000$-dimensional vector estimation, with the first 10 entries of $\boldsymbol{\theta}^{\ast}$ equal 10, the next 90 entries equal $A$, and the remaining 900 entries equal 0. We consider a range of $A$, from $A=2$ to $A=7$. DP method will automatically estimate 3 cluster centers instead of 2. Average MSE and MAE are recorded in Table 7 and Table 8.

In some configurations, DP method is the best one among all the state-of-the-art methods. Although in some cases EBKM and EBMed performs better, MSE and MAE of our DP method are still comparable.

4 Conclusions and Discussions

Empirical Bayes method introduced in Martin and Walker (2014) does not have good performance in simulation study since a prior centered at $X_{i}$ may introduce noise in posterior estimation, whereas our nonparametric Bayes based clustering method will decrease the noise but capture the general pattern: the estimate of $\theta_{i}$ is adaptively shrunk toward the mean of nearby data points. Therefore our DP method outperforms the empirical Bayes estimator based on Martin and Walker (2014).

Theoretical results shows with proper choice of parameters in the prior, our posterior mean estimator achieves asymptotically minimax rate. In the implementation we propose a fast variational inference method which approximates posterior distribution, which is more efficient than Empirical Bayes method proposed by Koenker (2014) in terms of computation time. To our knowledge, this is a first work connecting high dimensional sparse vector estimation with clustering. Our nonparametric Bayesian estimator could be applied to high dimensional classification, feature selection, hypothesis testing and nonparametric function estimation. Possible extension in implementation includes using other well-studied clustering methods to estimate prior and comparing their performance.

Proofs

Proof for Lemma 1

Denote $\mathbf{w}=(w_{1},w_{2},\cdots,w_{T})$ and $\mathbf{m}=(m_{1},m_{2},\cdots,m_{T})$. Write $D_{n}$ in terms of the conditional prior $(\theta_{1},\theta_{2},\cdots,\theta_{n})|w\sim\Pi_{w,\mathbf{w},\mathbf{m}}$ and the marginal prior $w\sim\Pi$. Therefore

[TABLE]

where $p_{\theta_{i}}(X_{i})=\frac{1}{\sqrt{2\pi}}\exp(-\frac{(X_{i}-\theta_{i})^{2}}{2})$. For given $w$, the inner expectation involves an average over all configurations of the indicators $(1_{\theta_{1}=0},1_{\theta_{2}=0},\cdots,1_{\theta_{n}=0})$. This average is clearly larger than just the case where the indicators exactly match up with the support $S^{\ast}$ of $\theta^{\ast}$, times the probability of that configuration, that is,

[TABLE]

For each $i\in S^{\ast}$ and each $1\leq t\leq T$,

[TABLE]

Therefore we have

[TABLE]

Since

[TABLE]

We have $\prod_{i\in S^{\ast}}w_{t_{i}}\geq(Cs_{n}/n)^{s_{n}}=\exp(\log C\cdot s_{n}-\varepsilon_{n})$. Since $\sum_{i\in S^{\ast}}(X_{i}-m_{t_{i}})^{2}\leq 2\sum_{i\in S^{\ast}}(\theta_{i}^{\ast}-m_{t_{i}})^{2}+2\sum_{i\in S^{\ast}}(X_{i}-\theta_{i}^{\ast})^{2}$ and $\sum_{i\in S^{\ast}}(X_{i}-\theta_{i}^{\ast})^{2}\sim\chi^{2}_{s_{n}}$, by concentration inequality, we have

[TABLE]

Therefore $\sum_{i\in S^{\ast}}(X_{i}-\theta_{i}^{\ast})^{2}=o_{p}(\varepsilon_{n})$. Hence

[TABLE]

According to Law of Large Numbers,

[TABLE]

Therefore

[TABLE]

Using the same trick in the proof of Lemma 1 in Martin and Walker (2014), $\int^{1}_{0}w^{n-s_{n}}(1-w)^{s_{n}}\Pi(dw)$ is lower bounded by $\frac{\alpha}{1+\alpha}\exp[-2\varepsilon_{n}-\alpha s_{n}+o(s_{n})]$ as $n$ is sufficiently large. Putting these pieces together, we obtain

[TABLE]

4.1 Proof for Theorem 1

The main aim of the proof is to show the numerator, for sets $A_{n}$ away from $\boldsymbol{\theta}^{\ast}$, is not too large. Let $N_{n}$ be the numerator for $Q_{n}(A_{M\varepsilon_{n}})$, i.e.,

[TABLE]

Taking expectation of $N_{n}$, with respect to $P_{\theta^{\ast}}$, we get

[TABLE]

Since $\int_{\mathbb{R}}\{\frac{p_{\theta_{i}}(x_{i})}{p_{\theta_{i}^{\ast}}(x_{i})}\}^{\kappa}p_{\theta^{\ast}_{i}}(x_{i})dx_{i}=\exp\{-\frac{\kappa(1-\kappa)}{2}(\theta_{i}-\theta_{i}^{\ast})^{2}\}$, we have

[TABLE]

Since $A_{M\varepsilon_{n}}=\{\boldsymbol{\theta}|~{}~{}\|\boldsymbol{\theta}-\boldsymbol{\theta}^{\ast}\|^{2}>M\varepsilon_{n}\}$, we have

[TABLE]

Let $\frac{\kappa(1-\kappa)}{2}=c$. Next, take $M$ such that $cM>3$, and then take $K\in(3,cM)$, then using Markov Inequality we get $P_{\theta^{\ast}}(N_{n}>e^{-K\varepsilon_{n}})\leq e^{-(cM-K)\varepsilon_{n}}$. This upper bound has a finite sum over $n\geq 1$, so the Borel-Catelli lemma gives that $N_{n}\leq e^{-K\varepsilon_{n}}$ with probability 1 for all large $n$. Therefore we get

[TABLE]

Since $s_{n}=o(\varepsilon_{n})$, $Q_{n}(A_{M\varepsilon_{n}})\rightarrow 0$ as $n\rightarrow\infty$.

Variational Inference Algorithm Derivation

We will derive the variational inference algorithm for Dirichlet process mixture model. $\alpha_{0}$, $T$, $w_{0}$, $\sigma_{0}^{2}$ and the data vector $X^{n}$ is given in advance. The data generating process is summarized in 7. We treat $\mathbf{Z}$ as latent variables and $\mathbf{V},\boldsymbol{\eta^{\ast}},\boldsymbol{\xi}$ as parameters. Posterior distribution of all the parameters and latent variables is proportional to

[TABLE]

Recall that under the fully factorized variational assumption, we have

[TABLE]

Define $P(Z_{i}=t)=\phi_{i,t}$. First we find the optimal form of $q(\boldsymbol{\eta^{\ast}},\boldsymbol{\xi})$, which satisfies

[TABLE]

where $\log q(\xi_{t},\eta^{\ast}_{t})=1_{\xi_{t}=1}(\log w_{0}+\log\delta_{0}(\eta^{\ast}_{t})-\frac{\sum^{n}_{i=1}\phi_{i,t}X_{i}^{2}}{2})+1_{\xi_{t}=0}(\log(1-w_{0})-\log\sqrt{2\pi\sigma_{0}^{2}}-\frac{(\eta^{\ast}_{t})^{2}}{2\sigma_{0}^{2}}-\frac{\sum^{n}_{i=1}\phi_{i,t}(X_{i}-\eta^{\ast}_{t})^{2}}{2})+\text{const}$. Therefore the optimal form of $q_{\mathbf{p},\mathbf{m},\boldsymbol{\tau}}(\boldsymbol{\eta^{\ast}},\boldsymbol{\xi})$ is fully factorized across different clusters:

[TABLE]

In order to determine the updating formula for $p_{t},m_{t},\tau_{t}$, we use Method of Undetermined Coefficients. Suppose $q(\xi_{t},\eta^{\ast}_{t})$ is a mixture of a point mass of zero and normal distribution,

[TABLE]

therefore

[TABLE]

Even though there’s a normalizing constant, but the difference between multipliers of $1_{\xi_{t}=1}$ and $1_{\xi_{t}=0}$ is invariant with respect to the constant. Therefore we have the following equation:

[TABLE]

which holds for any $\eta^{\ast}_{t}\in\mathds{R}$. The solutions are given as follows:

[TABLE]

Next we deal with the optimal form for $q(\mathbf{V})$, which satisfies

[TABLE]

where $\log q(V_{1})=\sum^{n}_{i=1}\phi_{i,1}\cdot\log V_{1}+(\alpha_{0}-1+\sum^{T}_{t=2}\sum^{n}_{i=1}\phi_{i,t})\log(1-V_{1})+\text{const}$, $\log q(V_{2})=\sum^{n}_{i=1}\phi_{i,2}\log V_{2}+(\alpha_{0}-1+\sum^{T}_{t=3}\sum^{n}_{i=1}\phi_{i,t})\cdot\log(1-V_{2})+\text{const}$, $\cdots$, $\log q(V_{T-1})=\sum^{n}_{i=1}\phi_{i,T-1}\cdot\log V_{T-1}+(\alpha_{0}-1+\sum^{n}_{i=1}\phi_{n,T})\log(1-V_{T-1})+\text{const}$. Thus we proved

[TABLE]

Besides, $V_{1}$ follows Beta Distribution with parameters $(\gamma_{11},\gamma_{21})=(\sum^{n}_{i=1}\phi_{k,1}+1,\alpha_{0}+\sum^{T}_{t=2}\sum^{n}_{i=1}\phi_{i,t})$, $V_{2}$ follows Beta Distribution with parameters $(\gamma_{12},\gamma_{22})=(\sum^{n}_{i=1}\phi_{k,2}+1,\alpha_{0}+\sum^{T}_{t=3}\sum^{n}_{i=1}\phi_{i,t})$,$\cdots$, $V_{T-1}$ follows Beta Distribution with parameters $(\gamma_{T-1,2},\gamma_{T-1,2})=(\sum^{n}_{i=1}\phi_{k,T-1}+1,\alpha_{0}+\sum^{n}_{i=1}\phi_{i,t})$.

Finally, we deal with $q(\mathbf{Z})$. We have the following optimal form

[TABLE]

therefore we proved $q_{\boldsymbol{\Phi}}(\mathbf{Z})=\prod^{n}_{i=1}q_{\boldsymbol{\phi}_{i}}(Z_{i})$. Since $\mathds{E}_{\xi_{t},\eta^{\ast}_{t}}(X_{i}-\eta^{\ast}_{t})^{2}=X_{i}^{2}-2(1-p_{t})m_{t}X_{i}+(1-p_{t})(m_{t}^{2}+\tau^{2}_{t})$, we have

[TABLE]

Therefore $q(Z_{i})$ is the probability mass function of Multinomial Distribution. Once we fix $i$, $\phi_{i,t}\propto\exp(S_{t})$, where $S_{t}=\exp[\mathds{E}_{\gamma_{1,t},\gamma_{2,t}}(\log V_{t})+\sum^{t-1}_{i=1}\mathds{E}_{\gamma_{1,i},\gamma_{2,i}}(\log(1-V_{i}))+(1-p_{t})m_{t}X_{i}-\frac{1}{2}(1-p_{t})(m_{t}^{2}+\tau_{t}^{2})]$.