Exact Minimax Predictive Density for Sparse Count Data

This paper discusses predictive densities under the Kullback--Leibler loss in high-dimensional sparse count data models. In particular, Poisson sequence models under sparsity constraints are discussed. Sparsity in count data implies zero-inflation or quasi zero-inflation, that is, situations where there exists an excess of zeros or near-zero counts. We investigate the exact asymptotically minimax Kullback--Leibler risks in sparse and quasi-sparse Poisson sequence models. We also provide a class of Bayes predictive densities that attain exact asymptotic minimaxity. To this end, we employ spike-and-slab priors with improper slab priors and introduce the notion of predictive specification of the scale in the mixture of improper priors. For application, we provide exact asymptotically minimax predictive densities that are adaptive to an unknown sparsity, and discuss the performance of the proposed Bayes predictive densities in settings where current observations are missing completely at random. The simulation studies as well as applications to real data demonstrate the efficiency of the proposed Bayes predictive densities.

Introduction
Predictive density is a probability density of future observations on the basis of current observations. It is used not only to estimate future observations but also to quantify their uncertainty. It has a wide range of application in statistics, information theory, and machine learning. The simplest class of predictive densities is the class of plug-in predictive densities. A plug-in predictive density is constructed by substituting an estimator into an unknown parameter of a statistical model. Another class of predictive densities is the class of Bayes predictive densities. A Bayes predictive density is the posterior mixture of densities of future observations. There is a vast literature on predictive density within statistical models in finite dimensions; see Subsection 1.2 for the literature review. One of the important findings in the literature is that one can construct the Bayes predictive density better than a plug-in predictive density under the Kullback-Leibler loss; refer to [28]. Conversely, little is known about predictive density within statistical models in high dimensions. Important works in this direction are [42,41]. In [42,41], several predictive densities (including a Bayes predictive density) superior to all plug-in predictive densities even in sparse high-dimensional Gaussian models are constructed.
The aim of this paper is to construct an efficient predictive density for highdimensional sparse count data. The efficiency of a predictive density is measured by the supremum of the Kullback-Leibler risk under sparsity constraints. In this paper, two types of sparsity are discussed: exact and quasi sparsity. Exact sparsity in count data means that there exhibits an excess of zeros. Quasi sparsity means that there exhibits an excess of near-zero counts. Exact sparsity is captured by the 0 -norm, and quasi sparsity is captured by the number of near-zero counts. See Subsection 1.1 for the formulation.
The motivation for analyzing sparse count data is well-known. In analyzing high-dimensional count data, there often exhibits exact sparsity and quasi sparsity corresponding to inflation of zeros and near-zero counts, respectively. Data with an overabundance of zeros include agriculture [21], environmental sciences [1], and manufacturing [34]. Data with an overabundance of near-zero counts include DNA sequencing and terrorist attacks [12]. Another example (Japanese crime statistics) is presented in Section 4.

Problem setting and contributions
The main results are summarized with the problem formulation ahead. Let X i (i = 1, 2, . . . , n) be a current observation independently distributed according to Po(rθ i ), and let Y i (i = 1, 2, . . . , n) be a future observation independently distributed according to Po(θ i ), where θ = (θ 1 , . . . , θ n ) is an unknown parameter and r is a known constant. Constant r represents the ratio of the mean of the i-th (i = 1, . . . , n) current observation to that of the i-th future observation. By a sufficiency and transformation reduction, this constant represents the ratio of sample sizes of current observations to those of future observations. Suppose that X = (X 1 , . . . , X n ) and Y = (Y 1 , . . . , Y n ) are independent. The densities of X and Y with parameter θ are denoted by p(x | θ) and q(y | θ), respectively: Two types of parameter spaces are of interest: the exact and the quasi sparse parameter spaces. Given s ∈ (0, n), the exact sparse parameter space is described as Θ[s] := {θ ∈ R n + : θ 0 s}, where · 0 is the 0 -norm given by θ 0 := #{i : θ i > 0}; Given s ∈ (0, n) and a threshold ε > 0, the quasi sparse parameter space is described as Θ[s, ε] := {θ ∈ R n + : N (θ, ε) s}, where N (θ, ε) := #{i : θ i > ε}, ε > 0. A threshold ε determines whether the parameter value of each coordinate is near-zero or not.
The corresponding risk (expected loss) is denoted by R(θ,q) = x∈N n y∈N n p(x | θ)q(y | θ) log q(y | θ) q(y; x) . To express high-dimensional settings under sparsity constraints, we employ the high dimensional asymptotics in which n → ∞ and η n := s/n = s n /n → 0. The values of s and ε possibly depend on n and thus in what follows the dependence on n is often expressed as s = s n , ε = ε n , and η = η n .

The minimax Kullback-Leibler risk over Θ[s] is defined as
The main theoretical contributions of this paper are summarized as follows: (i) The exact asymptotic minimax risks R(Θ[s n ]) and R(Θ[s n , ε n ]) over Θ[s n ] and Θ[s n , ε n ] are identified in Theorem 2.1. Constant r is shown to be the only key parameter in describing the ratio of the predictive minimax risks to the estimative minimax risks as shown in Proposition 2.1. This constant highlights the interesting parallel between Gaussian and Poisson decision theories as discussed in Subsection 1.2; (ii) A class of Bayes predictive densities attaining the exact asymptotic minimaxity is presented in Theorem 2.2; (iii) A class of exact asymptotically minimax predictive densities that are adaptive to an unknown sparsity is also presented in Theorem 2.3; (iv) The proposed predictive densities are shown to attain the exact asymptotic minimaxity in the settings where current observations are missing completely at random (MCAR) in Section 3. Consideration for adaptation and for missing observations are important in application. In the two real-data examples in Section 4, there is no prior knowledge of s n , and there exhibit missing values of observations; hence we take care of adaptation and missing observations. Practical effectiveness of the proposed Bayes predictive densities is examined by both simulation studies and applications to real data in Section 4. These studies show that the proposed Bayes predictive densities are effective in senses of not only point prediction but also predictive uncertainty quantification.
The proposed class of predictive densities builds upon spike-and-slab prior distributions with improper slab priors. The proposed class involves two ideas. The first idea is to control the tail behavior of a slab prior. Interestingly, spikeand-slab prior distributions with slab priors having exponential tails do not yield asymptotically exact minimax predictive densities as Proposition 2.2 indicates. To obtain predictive densities that are not only asymptotically minimax but are also easily implemented by the exact sampling, spike-and-slab priors with improper slab priors are used. The second idea is to specify the scale of an improper slab prior so as to minimize an upper bound of the minimax risks. In the objective Bayesian literature, it is well-known that the scale of improper priors having the different dimensions within their mixture affects the posterior distribution unlike a single improper prior; see [23,35,13]. The scale is uniquely chosen to minimize an upper bound of supremum risks and is dependent on the constant r. See Subsection 2.2 for more details.

Literature review
There is a rich literature on constructing predictive densities in fixed finite dimensions. In the literature, Bayes predictive densities are showed to dominate plug-in predictive densities. Studies of Bayes predictive densities date back to [2,43,3,44]. The first quantitative comparison of Bayes and plug-in predictive densities in a wide class of parametric models is [28]. [28] showed that there exists a Bayes predictive density that dominates a plug-in predictive density under the Kullback-Leibler loss, employing asymptotic expansions of Bayes predictive densities; see also [22] for asymptotic expansions of Bayes predictive densities. Minimax Bayes predictive densities for unconstrained parameter spaces are studied in [36,4]. These predictive densities under parametric constraints are studied in [17,33,37]. Shrinkage priors for Bayes predictive densities under Gaussian models are investigated in [29,19,27,40]; see also [26,5] for the cases where the variances are unknown. Shrinkage priors for Bayes predictive densities under Poisson models are developed in [30,32]. The cases under α-divergence losses are covered by [11,46,39,50].
Relatively little is known about constructing predictive densities in high dimensions. [42,41] elegantly construct an asymptotically minimax predictive density for sparse Gaussian models. [47] obtained an asymptotically minimax predictive density for nonparametric Gaussian regression models under Sobolev constraints; thereafter, [49] obtained an adaptive minimax predictive density for these models. See also [48]. All above results employ Gaussian likelihood and the corresponding results for count data have been not known.
Poisson models deserve study in their own right as prototypical count data modeling. Poisson models exhibit several correspondences to Gaussian models. [9,24,25,38] find the correspondence in estimation of means using the rescaled squared loss defined as [20,30,31,32] find the correspondence in prediction using the Kullback-Leibler loss. In particular, [25,38] find it in the exact asymptotic minimaxity under ellipsoidal and rectangle constraints in high-dimensional Poisson models using the re-scaled squared loss. In spite of the elegant correspondence in [25,38], the re-scaled squared loss is not compatible with sparsity: the loss diverges if θ i = 0 andθ i (X) = 0 for at least one index i.
Employing the Kullback-Leibler divergence, this paper presents the results of exact asymptotic minimaxity in both estimation and prediction within sparse Poisson models, which are cleanly parallel to the result for sparse Gaussian models by [42]; see Subsection 2.2 for detailed discussions. Here we present a new strategy in constructing a predictive density, which covers several new topics in predictive density under sparsity constraints: Quasi-sparsity; The adaptation to sparsity; Missing Completely At Random (MCAR). See real-data examples in Subsection 4.2 for the importance of these topics.
Relatively scarse are theoretical studies of zero-inflated or quasi zero-inflated Poisson models in high dimensions in spite of their importance. [12] constructs local-global shrinkage priors for high-dimensional quasi zero-inflated Poisson models, providing theoretical properties of the shrinkage factors and of the multiple testing statistics. In our own opinion, a prior should be constructed according to the context in which the prior is used. Our priors adopt the information of a predictive setting where they are used, i.e., the constant r. This indicates our priors are suitable in prediction. We confirm in Section 4 that our priors broadly outperform their priors in predictive density. In contrast, we consider that their priors would be more suitable than our priors in multiple testing or in interpreting shrinkage factors.

Notations
Additional notations are summarized here. The notation a n ∼ b n signifies that a n /b n converges to 1 as n goes to infinity. The notation a n b n signifies that a n /b n converges to a constant as n goes to infinity. The notation O(a n ) indicates a term of which the absolute value divided by a n is bounded for a large n. The notation o(a n ) indicates a term of which the absolute value divided by a n goes to zero in n.
For a function f : indicates the expectation of f (X, Y ) with respect to p(x | θ)q(y | θ). Likewise, for a function g : N → R, the expectation E λ [g(X 1 )] indicates the expectation of g(X 1 ) with respect to Po(λ).

Organization
The rest of this paper is organized as follows. Section 2 presents theoretical results in predictive density for sparse Poisson models without missing observations. Theorems 2.1-2.3 are the main theorems. Theorem 2.1 provides exact asymptotic minimax risks; Theorem 2.2 presents a class of Bayes predictive densities attaining the exact asymptotic minimaxity; Theorem 2.3 presents a class of adaptive predictive densities. Section 3 provides the description of problem settings and theoretical results in predictive density for sparse Poisson models with MCAR settings. Section 4 reports simulation studies and applications to real data. Section 5 presents the outlines of Theorems 2.1 and 3.1. All proofs of theorems and propositions for Section 2 are provided in Appendix A. All proofs of theorems and propositions for Section 3 are provided in Appendix B.

Main results
Results for sparse Poisson models without missing observations are provided in order; the precise description of the exact asymptotic minimax risk, and the construction of Bayes predictive densities that attain exact asymptotic minimaxity. Discussions on comparison between prediction and estimation and on adaptation are provided in the subsequent subsection.
(a) For the exact sparse parameter space Θ[s n ], the following holds: (b) For the quasi sparse parameter space Θ[s n , ε n ] with any shrinking sequence ε n > 0 such that ε n = o(η n ), the following holds: The implication of this theorem in comparison with the estimative minimax risk is discussed in Subsection 2.2. The proof outline of this theorem is provided in Section 5. The full proof is provided in Appendices A.1-A.2.
The second theorem provides a class of Bayes predictive densities that attain the exact asymptotic minimaxity with the knowledge of sparsity s n . For h > 0 and κ > 0, let Π[h, κ] be an improper prior of the form where δ 0 is the Dirac measure centered at 0. The corresponding Bayes predictive density is well-defined and can be written as Here, the coordinate-wise marginal distribution of q Π[h,κ] is just a zero-inflated negative binomial distribution, and sampling from q Π[h,κ] is easy. Let Theorem 2.2. Fix r ∈ (0, ∞) and κ > 0. Fix also a sequence s n ∈ (0, n) such that η n = s n /n = o(1). The predictive density q Π[Lηn,κ] with L > 0 and κ > 0 satisfies and thus is asymptotically minimax: two asymptotic equalities hold as n → ∞, where ε n > 0 is a shrinking sequence such that ε n = o(η n ), and Υ 1 and Υ 2 represent terms that are independent of L or that are O(s n η n ); Furthermore, L * minimizes the right hand sides in (1) with respect to L.
The idea behind the proposed class with the scale L * is discussed in Subsection 2.2. The proof is provided in Appendix A.3.
(a) The predictive density q Π[L * η n ,κ] is adaptive in the exact minimax sense on the class of exact sparse parameter spaces: for any sequence s n ∈ (0, n) such that inf n s n > 0 and s n = o(n 1/2 ), we have The predictive density q Π[L * η n ,κ] is adaptive in the exact minimax sense on the class of quasi sparse parameter spaces: for any shrinking sequence ε n > 0 such that ε n = o(η n ) and for any sequence s n ∈ (0, n) such that inf n s n > 0 and s n = o(n 1/2 ), we have The main ingredient of the proof is discussed in Subsection 2.2. The proof is provided in Appendix A.4.

Discussions
Several discussions are provided in order. All proofs of propositions are provided in Appendix A.5.

Prediction and estimation
On the basis of Theorem 2.1, minimax risks for estimation and prediction are compared. The minimax risk E(Θ[s n ]) restricted to plug-in predictive densities is provided in the following proposition.
hold as n → 0, where s n ∈ (0, n) is any sequence such that η n = s n /n = o(1), and ε n > 0 is any shrinking sequence such that ε n = o(η n ). According to Theorem 2.1 and Proposition 2.1, the exact asymptotic constants of predictive and estimative minimax risks (R(Θ[s n ]) and E(Θ[s n ])) are different with r, and the rates of convergence of predictive and estimative minimax risks are identical with n. The constant r is the key parameter in describing the exact constants of the risks: the exact constant of predictive minimax risk increases as r decreases but is bounded by 1, while that of estimative minimax risk grows to infinity as r decreases. A similar phenomenon occurs in sparse Gaussian models as found in [42]. However, the exact constants of predictive minimax risk are quite different in Poisson and Gaussian models, and they are not intuitively derivable. Figures 1 and 2 show the comparisons of the exact constants of minimax risks for sparse Poisson and Gaussian models. The vertical line indicates values of the risks and the horizontal line indicates values of r. They show the similarity of the behavior with respect to r of minimax risks in Poisson and Gaussian cases. An interesting observation in comparison of Poisson and Gaussian cases is that the exact constants of predictive minimax risks in both cases get closer to 1 as r approaches to 0.

Spike-and-slab priors with improper slab priors
For Theorem 2.2, an idea behind the construction of Π[L * η n , κ] with L * defined before Theorem 2.2 is discussed. We employ "spike-and-slab" priors with improper slab priors. Spike-and-slab priors are widely used in the literature of Bayesian sparse modelings; see [18,45,7,8] for the recent developments. As discussed in the literature, the tail behavior of the slab prior strongly affects the performance, and thus the choice of a slab prior is substantially important. We provide two propositions indicating the effects of the tail behavior on the performance.
The first proposition shows an example in which spike-and-slab priors with slab priors having exponential tails suffer from the sub-optimality. Let s n ∈ (0, n) be a sequence such that η n = s n /n = o(1). For µ 0, ν = (ν 1 , . . . , ν n ) ∈ R n + , and k = (k 1 , . . . , k n ) ∈ R n , let where Ga(·; ν 1 , k 1 ) is Gamma distribution with inverse scale parameter ν 1 and shape parameter k 1 . The sub-optimality is discussed in the exact sparse parameter space.
Proposition 2.2. Fix ν ∈ R n + and a shrinking sequence µ = O(η n ). Then, the following holds: Note that the sup-optimality of spike-and-slab priors with bounded uniform priors in sparse Gaussian models is reported in [41]. Bounded uniform priors have exponential tails.
It is of theoretical interest to reveal a condition on the tail of a slab prior under which the Bayes predictive density based on a spike-and-slab prior is asymptotically minimax. The second proposition shows that a polynomial decay of the tail of a slab prior is allowed; thus there exists a proper prior distribution of which the posterior mixture attains exact asymptotic minimaxity. Let s n ∈ (0, n) be a sequence such that η n = s n /n = o(1). Let Proposition 2.3. Fix r ∈ (0, ∞). Then, the following holds: Here we shall address the predictive specification of the scale L in Π[Lη n , κ], i.e., L * . As the former part of Theorem 2.2 indicates, there is a practical ambiguity in the theoretical result that Π[h, κ] yields the exact aymptotically minimax predictive density as long as h = Lη n with L > 0. Our choice of L * avoids this ambiguity. The upper bound (1) is obtained by applying just one inequality sup θ {f (θ) + g(θ)} sup θ f (θ) + sup θ g(θ) with f, g to the supremum risk; hence the scale L * is close to the second-order optimal among the class Π[Lη n , κ]. This also shows that the scale of improper priors within their mixture can be specified by the predictive setting (characterized by r) in which priors are used. Our idea here is relevant to [23]: In [23], the scale of improper priors within their mixture is determined to yield log-posterior probabilities that coincide with log maximum likelihood plus an Akaike factor; see also the appendix of [35]. In this light, [23] and this paper indicate that the specifications of the scale can be done from a predictive viewpoint.

Adaptation
Discussions for Theorem 2.3 are provided.
First, the key idea of Theorem 2.3 is to evaluate a point-wise performance of an empirical Bayes predictive density with an estimator of the sparsity pluggedin. When a convenient estimatorŝ n for s n is available, pluggingη n :=ŝ n /n into η n of Π[L * η n , κ] yields the empirical Bayes predictive density q Π[L * η n ,κ] . The next proposition (Proposition 2.4) provides a point-wise performance guarantee for q Π[L * η n,κ] with an estimatorŝ n that satisfies the following condition. (i)ŝ n γ a.s. for some γ > 0; (ii) There exists δ > 0 for which we have Proposition 2.4. Fix r ∈ (0, ∞) and κ > 0. Suppose that Condition 2.1 holds. There exists a positive constant c depending only on r, κ, and γ such that (a) for θ ∈ Θ[s n ] with any sequence s n ∈ (0, n), we have with any sequence s n ∈ (0, n) and any sequence ε n > 0, we have The estimator max{1, #{i : Second, Condition 2.1 is discussed. Consider the O(s n log s n )-term in Proposition 2.4. This term can be improved provided that δ < 1. However, the assumption δ < 1 is considered to be too stringent for θ near 0, since for i = 1, . . . , n, it is intrinsically difficult to test whether θ i > 0 or θ i = 0 for small θ i ; Hence, paying the cost of adding the O(s n log s n )-term, we obtain Proposition 2.4 for any δ > 0. Conversely, Theorem 2.3 requires an assumption on s n . The condition on s n in Theorem 2.3 is discussed in Section 6.

Predictive Density for Sparse Poisson models with
Missing-Completely-At-Random settings In this section, a series of results discussed in Section 2 is generalized into settings with MCAR settings. Prediction in sparse Poisson models with MCAR is formulated as follows. Let r i 's (i = 1, 2, . . .) be positive random variables. Given r i (i = 1, . . . , n), let X i (i = 1, 2, . . . , n) be a current observation independently distributed according to Po(r i θ i ), and let Y i (i = 1, 2, . . . , n) be a future observation independently distributed according to Po(θ i ), where θ i (i = 1, . . . , n) is an unknown parameter. Suppose that X = (X 1 , . . . , X n ) and Y = (Y 1 , . . . , Y n ) are independent. The true densities of X and Y with parameter θ given r i s are denoted by p(x | θ) and q(y | θ), that is For the sequel theoretical analyses, we assume that r i s are independent and identically distributed according to a sampling distribution G, and make the following condition on G. Let E G be the expectation with respect to G.
Condition 3.1 (i) is usual. Condition 3.1 (ii) excludes any distribution G that is highly concentrated around 0 and it is not stringent. Consider a longitudinal situation in which X i (i = 1, . . . , n) is obtained as the sum of {X i,j : j = 1, . . . , r i }, where r i (i = 1, . . . , n) represents the sample size in the i-th coordinate, and for each i, X i,j (j = 1, . . . , r i ) follows Po(θ i ). Condition 3.1 implies that for each coordinate there exists at least one observation: r i 1. Also, for verifying Condition 3.1 (ii) in application, refer to Section 4.

Main results
In this section, results for sparse Poisson models with MCAR are presented. Notations for the theorems are introduced ahead. Fix an infinite sequence The following theorems provide theoretical results for settings with MCAR. The first theorem describes the exact minimax risk. The second theorem states that when G satisfies Condition 3.1, q Π[Lηn,κ] with 0 < κ 1 attains exact asymptotic minimaxity regardless of a sampling distribution G of r i s. Here let The third theorem provides adaptive predictive densities. Letŝ n := max{1, #{i : X i 1, i = 1, . . . , n}} and letη n :=ŝ n /n.
is a sequence such that η n = s n /n = o(1), and ε n > 0 is any shrinking sequence such that ε n = o(η n ).
hold, where s n ∈ (0, n) is a sequence such that η n = s n /n = o(1), and ε n > 0 is any shrinking sequence such that ε n = o(η n ). (a) For exact sparse parameter spaces, the following asymptotic equality holds: for any sequence s n ∈ (0, n) such that inf n s n > 0 and s n = o(n 1/2 ), we have (b) For quasi sparse parameter spaces, the following asymptotic equality holds: for any shrinking sequence ε n > 0 such that ε n = o(η n ) and for any sequence s n ∈ (0, n) such that inf n s n > 0 and s n = o(n 1/2 ), we have The proof outline of Theorem 3.1 is provided in Section 5. The full proofs of all above theorems are provided in Appendices B.1-B.3.

Discussions
In this subsection, by considering the following two settings, we compare Theorem 3.1 with Theorem 2.1: (A) G is the Gamma distribution with shape parameter r/l and scale parameter l for 0 < l 1 and r 2; (B) G is the distribution of 1 + S, where S follows the binomial distribution with trial number N and success probability p.
In Setting (A), the mean remains r for any l and the variance is rl, which means that as l → 0, G is weakly convergent to the Dirac measure δ r centered at r corresponding to a non-MCAR case. In Setting (B), G is weakly convergent to the Dirac measure δ 1 centered at 1 as p → 0, and G is weakly convergent to the Dirac measure δ 1+N as p → 1.  In Figure 3, the vertical line indicates exact constants and the horizontal line indicates values of l. The blue line denotes the case with r = 2, the green line denotes the case with r = 4, and the red line denotes the case with r = 6. In Figure 4, the vertical line indicates exact constants and the horizontal line indicates values of p. The blue line denotes the case with N = 10, the green line denotes the case with N = 20, and the red line denotes the case with N = 30. Figures 3 and 4 show the exact constants of predictive minimax risks in Settings (A) and (B), respectively. According to Figure 3, the constant gets larger as the variance of G increases. The constant in a MCAR case approaches to that in a non-MCAR case in the limit l → 0. According to Figure 4, the constant gets smaller as the missing probability 1 − p gets smaller. Further, the numerical result in Setting (B) is consistent to the results in [15,16] for the literature of nonparametric regression in the presence of missing observations. Theorems 1 and 2 in [15] provide tight lower and upper bounds on mean integrated squared errors (MISE) in nonparametric regression with predictors missing at random. Those theorems also provide an exact asymptotically minimax estimator for MCAR cases and show that if predictors are MCAR, the minimax MISE gets smaller as the missing probability approaches zero.
Next, we provide a proposition useful for comparing the estimation and prediction in cases with MCAR. This proposition provides the same comparison as in Subsection 2.2.  (1), and ε n > 0 is any shrinking sequence such that ε n = o(η n ).

Simulation studies
In this section, we present simulation studies to compare the performance of various predictive densities.

Sparse Poisson model
Consider a sparse Poisson model without MCAR settings. Parameter θ and observations X and Y are drawn from respectively. Here, • ν 1 , . . . , ν n are independent samples from the Gamma distribution with a shape parameter 10 and a scale parameter 1; • S is drawn from the uniform distribution on all subsets having exactly s elements; • ν 1 , . . . , ν n and S are independent.
Here for a subset J ⊂ {1, . . . , n}, e J indicates the vector whose i-th component is 1 if i ∈ J and 0 otherwise. We examine two cases (n, s, r) = (200, 5, 1) and (200,5,20), and generate 500 current observations X's and 500 future observations Y 's. See Appendix C for the results with different choices of (n, s, r).
We compare the following four predictive densities: • The proposed predictive density based on Π[L * η n , κ]; • The Bayes predictive density based on the shrinkage prior proposed in [30]; • The Bayes predictive density based on the Gauss hypergeometric prior proposed in [12]; • The plug-in predictive density based on an 1 -penalized estimator.
For the proposed predictive density with Π[L * η n , κ], an estimatorŝ n for sparsity is the simple estimator described before Theorem 2.3. The second predictive density is shown in [30] to dominate the Bayes predictive density based on the Jeffreys prior. This predictive density has a hyper-parameter β and in simulation studies it is fixed to be 1. The third predictive density employs the global-local prior proposed in [12] and the specification of the hyper-parameters follows the online support pages the authors provide.
The performance of predictive densities is evaluated using the following three measures: • the mean of the 1 distance ( n i=1 |u i − v i | for u, v ∈ R n ) between the mean of a predictive density and a future observation, • the predictive log likelihood, that is, the log of the value of a predictive density at sampled Y and X, and • the (empirical) coverage probability of Y on the basis of the joint 90%prediction set constructed by a predictive density. Tables 1 and 2 show the results of the comparison. The following abbreviations are used. The Bayes predictive density proposed in [12] is abbreviated to GH. The Bayes predictive density proposed in [30] is abbreviated to K04. The plug-in density based on an 1 -penalized estimator with regularization parameter rλ is abbreviated to 1 (λ). The abbreviation 1 distance represents a mean Table 1 Comparison of predictive densities without MCAR and with (n, s, r) = (200, 5, 1): the 1 distance, PLL, and 90%CP represent the mean 1 distance, the predictive log likelihood, and the empirical coverage probability based on a 90%-prediction set, respectively. For each result, the averaged value is followed by the corresponding standard deviation. Underlines indicate the best performance. gard to the empirical coverage probabilities, the predictive densities based on Π[L * η n , 0.1] and the Gauss hypergeometric prior give the empirical coverage probabilities of Y that are relatively close to the nominal level. The prediction set of the plug-in predictive density based on the 1 -penalized estimator is too narrow to cover future observations. This is mainly because for this plug-in predictive density, an 1 -penalized estimator returns zero for a coordinate at which the current observation is zero and most of the marginal predictiev intervals degenerate into zero. This degeneracy also induces the divergence of a predictive log likelihood value of the plug-in predictive density based on an 1 -penalized estimator.

Sparse Poisson model with Missing-Completely-At-Random settings
Consider a sparse Poisson model with MCAR. Parameter θ and observations X and Y are drawn in the following way: where ν 1 , . . . , ν n and S follow the same distributions as those in the previous subsection. In addition, r i − 1 (i = 1, . . . , n) are independently drawn from the binomial distribution Bi(m, p) with the parameters (m, p) either (1, 0.9) or (10, 0.9). We set (n, s) = (200, 5), and generate 500 current observations X's and 500 future observations Y 's. We compare the following four predictive densities: • The proposed predictive density based on Π[Lη n , κ]; • The Bayes predictive density based on the shrinkage prior proposed in [32]; • The Bayes predictive density based on the Gauss hypergeometric prior proposed in [12]; • The plug-in predictive density based on an 1 -penalized estimator.
An estimatorŝ n is determined in the same manner as in the previous subsection. In [32], the second predictive density is shown to dominate the Bayes predictive density based on the Jeffreys prior in the case where the numbers of observations are coordinate-wise different. In simulation studies, each hyper parameter β i of the second predictive density is fixed to be 1.
In comparing the performance, we use the weighted 1 distance between the mean of the predictive density and a future observation. Here the the weighted 1 distance between u, v ∈ R n is given by For the construction of prediction sets, we also use this weighted 1 distance. Table 3 Comparison of predictive densities with MCAR and with (n, s, m, p) = (200, 5, 1, 0.9): the W-1 distance, PLL, and 90%CP represent the weighted mean 1 distance, the predictive log likelihood, and the empirical coverage probability based on a 90%-prediction set, respectively.
For each result, the averaged value is followed by the corresponding standard deviation. Underlines indicate the best performance.  Table 4 Comparison of predictive densities with MCAR and with (n, s, m, p) = (200, 5, 10, 0.9): the W-1 distance, PLL, and 90%CP represent the weighted mean 1 distance, the predictive log likelihood, and the empirical coverage probability based on a 90%-prediction set, respectively.
For each result, the averaged value is followed by the corresponding standard deviation. Underlines indicate the best performance.  Tables 3 and 4 show the results. In addition to the abbreviations used in the previous subsection, the abbreviation W-1 distance represents a mean weighted 1 distance. We see that the weighted 1 distances by the proposed predictive density based on Π[Lη n , 0.1] and the plug-in predictive density based on the 1 -penalized estimator are in the smallest level of all predictive densities compared here. As is the case without MCAR, the proposed predictive density with Π[Lη n , 0.1] broadly has the coverage probability that is relatively close to the nominal level, whereas the other predictive densities (with one expecption) do not.
In conclusion, regardless of whether there exists a missing of observations or not, these simulation studies show that the proposed predictive density is highly effective not only in point prediction but also in uncertainty quantification for future observations. Also, they suggest that the proposed predictive density with κ = 0.1 is stable than that with κ = 1.0.

Applications to Real Data
In this section, we apply our methods to two real datasets; Japanese crime data and exome sequencing data.

Pickpockets in Tokyo Prefecture
We discuss the validity of the proposed predictive densities empirically using Japanese crime data. A motivation for this analysis comes from the importance of taking measures against future crimes by utilizing past crime data.
We apply our methods to crime data from an official database called the number of crimes in Tokyo by type and town [14]. This database reports the total numbers of crimes in Tokyo Prefecture. They are classified by a town and also by the type of crimes.
We use pickpocket data from 2012 to the first half of 2018 at 978 towns in eight wards (Bunkyo Ward, Chiyoda Ward, Chuo Ward, Edogawa Ward, Koto Ward, Minato Ward, Sumida Ward, and Taito Ward). Figure 5 shows total counts of pickpockets from 2012 to 2017 for all towns in the wards. The scale of the pickpocket occurrences in each town is expressed by a gradation of colors: there have occurred more pickpockets in a deeper-colored town over 6 years. There have occurred no pickpockets in white-colored towns. As seen from are used as future observations. Since the counts in the first half of 2018 would be considered as the half of the total counts in 2018, in general, the ratio r of sample sizes is set as r = 12. However, some observations are missing because several towns, though in rare cases, did not report the counts. As in Subsection 4.1, we compare the proposed predictive density q Π[Lηn,κ] with the three existing predictive densities, that is, the Bayes predictive density GH based on a Gauss hypergeometric prior, the Bayes predictive density K04 based on the shrinkage prior, and the plug-in predictive density based on an   Table 5 shows a summary of comparisons. In all measures, the proposed predictive density q Π[Lηn,0.1] has the best scores.

Rare mutation rates in an oncogene
We consider an application of the proposed methods to exome sequencing data from a huge database called the Exome Aggregation Consortium (ExAC). ExAC reports the total numbers of mutant alleles in each genetic position along the whole exome, gathered from 60706 unrelated individuals. The analysis in this subsection is motivated by [12]: see Section 5 in [12]. We focus on rare allele mutations in a gene PIK3CA; For the importance of analysing rare allele mutations and the choice of the gene, see [12]. We also follow the pre-process of the data described in [12].
We apply the sparse Poisson model to the numbers of rare mutant alleles as follows. We denote by X i (i = 1 . . . , 551) the number of rare mutant alleles in the i-th position. We assume that X i is distributed according to Po(r i θ i ). Here r i (i = 1, . . . , 551) is the double number of individuals whose i-th location are sequenced, and θ i is the frequency rate common to individuals in the i-th position. The doubling is necessary because each individual has two copies of each allele. Since numerous fragments of DNA sequences are sampled and read at random, r i 's are different in general (see Appendix E); hence the data is regarded as having an MCAR structure.
Our goal is to predict the behavior of rare allele mutations under the assumption that all individuals could be sequenced at all positions and the sequencing depth is uniform across the gene; the target mutation counts Y i 's are assumed to be distributed according to ⊗ n i=1 Po(rθ i ) with r = 121412 = 2 × 60706. We compare the propose prior Π[Lη n , κ] with the two existing priors, that is, the Gauss hypergeometric prior in [12] and the shrinkage prior in [30]. For the proposed prior, κ is set to 0.1 and the estimateŝ n to be plugged in is determined by the k-means clustering with k = 2 (the resulting value ofŝ n is 17). In this study, we give the qualitative comparison using samples from predictive densities. Figure 6 shows marginal prediction intervals at a nominal level of 50%. The prediction intervals of the proposed predictive density q Π[Lηn,0.1] and of the Bayes predictive density based on GH show apparently similar behaviors: they shrink at positions whose counts are low, and they remain the scales at positions whose counts are large. The prediction intervals by K04 degenerate at most of the locations. In lower regimes of counts, there exists a dissimilarity between q Π[Lηn,0.1] and GH. Most of the medians of the intervals constructed by the predictive density based on GH are away from 0; while those constructed by q Π[Lηn,0.1] sometimes reach zero. Owing to the coordinate-wise independence, this dissimilarity is considered to reflect on that the proposed predictive density has more flexibility than GH.

Proof outlines
This section provides proof outlines of Theorems 2.1 and 3.1. Appendix A provides details of the proofs for Section 2. Appendix B provides details of the proofs for Section 3.

Outline of the proof of Theorem 2.1
The outline of the proof of Theorem 2.1 (a) is provided. The proof is divided into two parts: providing upper and lower bounds on R(Θ[s n ]).
A lower bound on R(Θ[s n ]) builds upon Bayes risk maximization based on a block-independent prior. Let Π B,ν (dθ) with ν > 0 be a block-independent prior built as follows: divide {1, 2, . . . , n} into contiguous blocks B j (j = 1, 2, . . . , s n ) with each length m n := η −1 n . In each block B j , draw (θ 1+mn(j−1) , . . . , θ mnj ) independently according to a single spike prior with spike strength ν > 0, where a single spike prior with spike strength ν > 0 is the distribution of νe I with a uniformly random index I ∈ {1, . . . , m n } and a unit length vector e i in the i-th coordinate direction. Finally, set θ i = 0 for the remaining n−m n s n components.
It will be shown that sup ν>0 B(ν) Cs n log η −1 n , from which we will have R(Θ[s n ]) Cs n log η −1 n .

An upper bound on R(Θ[s n ]) is derived by bounding the coordinate-wise
Kullback-Leibler risk of the Bayes predictive density q Π * based on Π * = Π[η n , κ] with κ > 0: where q Π * (y i | x i ) is the marginal distribution of q Π * . It will be shown that • ρ(0) = O(η n ); • sup λ>0 ρ(λ) (C + o(1))s n log η −1 n , from which we will have the following: Here, the coordinate-wise independence of the predictive density q Π * is used.
Lemma 5.1. For a prior Π of θ, ifθ Π (x; t) based on Π is strictly larger than 0 for any x ∈ N n and any t ∈ (r, 1 + r), then, the equality This formula is known in the literature, but the proof is given in Appendix D for the sake of completeness.
Lemma 5.1 reduces bounding the predictive risk of the Bayes predictive density to bounding the estimative risk of the corresponding Bayes estimator. Splitting the sample space for the control of the Bayes estimator together with Lemma 5.1, we will complete the proof; see Appendices A.1-A.2.

The main ingredient of proof of Theorem 3.1
The main ingredient of the proof of Theorem 3.1 is provided in this subsection. The next lemma is the main ingredient and describes the exact asymptotic minimax rirks for a setting with a fixed sequence {r i : i = 1, . . . , n} satisfying the following condition: holds as n → ∞ and η n = s n /n → 0; (b) for the quasi sparse parameter space Θ[s n , ε n ] with any shrinking sequence ε n such that ε n = o(η n ), the asymptotic equality R(Θ[s n , ε n ]) ∼ Cs n log(η −1 n ) holds as n → ∞ and η n → 0.
For a Bayes estimatorθ Π based on Π, let where E θ;τ is the expectation with respect to p(· | θ; τ ).

Lemma 5.3. [32]
For a prior Π of θ, ifθ Π (z; τ ) based on Π is strictly larger than 0 for any z ∈ N n and any τ ∈ [0, 1], then, the equality For completing the proof, we will split the sample space to control the behavior of the Bayes estimator as with the proof of Theorem 2.1; see Appendix B.4.1.

Discussion and conclusion
In this paper, we have investigated the exact asymptotically minimax Kullback-Leibler risks in sparse and quasi-sparse Poisson sequence models. We have constructed Bayes predictive densities attaining exact asymptotic minimaxity using spike-and-slab priors with improper priors. The scale of improper priors is specified from the the predictive point of view. For application, we have provided exact asymptotically minimax predictive densities that are adaptive to an unknown sparsity, and have discussed the performance in MCAR settings. The simulation studies as well as applications to real data have demonstrated the efficiency of the proposed densities.
For the future work, we shall discuss the condition on s n in Theorem 2.3. Consider the condition on s n in Theorem 2.3. This condition might be unavoidable if an estimator for s n is used, which is relevant to [6,10]. For sparse Gaussian models, [6] constructed an adaptive estimator for s n in the regime s n = o(n 1/2 ); while [6] removed a polynomially small proportion of non-zero components to construct an efficient estimator for s n in the regime s n ∼ n γ with γ > 1/2. [10] also found this phase transition in the estimation of non-smooth functionals. In addition, this condition on s n is conjectured to be unnecessary if another route is taken for constructing adaptive predictive densities. Either conjecture requires further theoretical developments in predictive density under sparsity constraints.

Step 2: Upper bound on R(Θ[s n ])
The explicit form ofθ Π * (·; t) is derived ahead. It depends only on whether x is 0 or not. The following two identities forθ Π * (·; t) show the form ofθ Π * : for t ∈ (r, r + 1), we havê The derivations are as follows. From the Bayes formula, we havê where we use the convention that 0 0 = 1. Substituting x 1 = 0 into equality (2) yields the first identity. Substituting x 1 1 into equality (2) yields the second identity.

A.2. Proof of Theorem 2.1 (b)
The proof follows almost the same line as in the proof of Theorem 2.1 (a).
Since N (θ, ε n ) s n for θ ∈ Θ[s n , ε n ], we have We put a prior Π * = Π[η n , κ] into Π and then we have R i (λ, q Π * ,i ) = ρ(λ). From (3), we have sup 0<λ R i (λ, q Π * ,i ) C log(n/s n ). Observe that for sufficiently small λ (depending only on r and κ), the first, the third, and the sixth terms in the right hand side of (3) are monotonically increasing with respect to λ, while the second and the forth terms of (3) are monotonically decreasing with respect to λ; hence substituting λ = ε n gives each supremum of the first, the third, and the sixth terms, while letting λ → 0 gives each supremum of the second and the forth terms. Thus using the inequality e −tλ 1 − tλ and using the rightcontinuity of λ log λ at λ = 0, we obtain sup λ εn R i (λ, q Π Π * ) = O(ε n log{η −1 n }), which completes the proof.

A.3. Proof of Theorem 2.2
We provide the proof only for Θ[s n ] and the proof for Θ[s n , ε n ] is straightforward. We begin with multiplying η n by L in Π * = Π[η n , κ] of the proof of Theorem 2.1. Then, inequalities (3) and (4) where Υ 1 is defined as the supremum of the terms except the first term in the right hand side of (3) with respect to λ. Observe that Υ 1 consists of terms independent of L or terms that are o(s n ). This completes the proof of the former part of Theorem 2.2. A simple calculus gives the latter part of Theorem 2.2 and thus we complete the proof.

A.4. Proof of Theorem 2.3
We begin with showing that δ forŝ n is bounded from above by a constant independent of n. To this end, we use the bias-variance decomposition ofŝ n −s n .
Consider the variance ofŝ n − s n . Using the inequalities we have Combining (8) with (9) shows that max{E|ŝ n /s n − 1|, E|ŝ n /s n − 1| 2 } is bounded from above by a positive constant depending only on inf n s n . Together with Proposition 2.4 (a), this completes the proof of Theorem 2.3 (a).

Proof for Case (b):
The proof is almost the same as of Case (a). Consider the bias ofŝ n − s n . Since the decomposition #{i : X i 1} = n j=1 Z j holds with Z j (j = 1, . . . , n) independent and distributed according to the Bernoulli distribution having the success probability 1 − exp(−rθ [j] ), and since the inequality θ [j] ε n holds for j = s n + 1, . . . , n, we have Consider the variance ofŝ n − s n . Using the inequality we have Combining (10) with (11) shows that max{E|ŝ n /s n − 1|, E|ŝ n /s n − 1| 2 } is bounded from above by a positive constant depending only on inf n s n . Together with Proposition 2.4 (b), this completes the proof of Theorem 2.3 (b).

A.5.1. Proof of Proposition 2.1
The proof follows almost the same line as in the proof of Theorem 2.1. We provide the proof only for Θ[s n ]. In bounding E(Θ[s n ]) below, we employ the Bayes risk maximization based on a block-independent prior Π B,ν with ν > 0.
Combining this expression ofθ ΠSS,1 (X 1 ; t) with the simple algebra that (a + b)/(c + d) b/d for any (a, b, c, d) ∈ R 4 + such that ad bc giveŝ for sufficiently large n depending only on ν 1 and for X 1 1.

A.5.3. Proof of Proposition 2.3
The proof steps are the same as in the proof of Theorem 2.2. For notational simplicity, we consider Π P as We begin with getting the following inequalities: there exist positive constants C 1 , C 2 , C 3 depending only on r such that for i = 1, . . . , n and for t ∈ (r, r + 1), we have where see the description before Lemma 5.1 for the definition ofθ ΠP (·; ·). We derive these bounds from the explicit form ofθ P (·; ·) and from the identity For λ > 0, x 1 ∈ N, and t ∈ (r, r + 1), let From the above inequalities forθ ΠP (·; ·), we have, for λ > 0 and t ∈ (r, r + 1), For λ > 0, let We will show that ρ P (0) = O(η n ) and sup λ>0 ρ P (λ) (C + o(1))s n log η −1 n . Lemma 5.1 gives t dt and thus it suffices to bound E tλ [ρ P (λ, X 1 ; t)] uniformly in t ∈ (r, r+1). Consider the decomposition Bounding the first and the second terms in the right hand side of the above decomposition is straightforward; we have Consider bounding the third term. Since log{1/(X 1 − 1)} 0 for X 1 2 and since Since we have for λ > 2 from the Bennett inequality for the Poisson distribution, and since which provides the bound on the third term and thus completes the proof.

A.5.4. Proof of Proposition 2.4 (a)
The decomposition of the Kullback-Leibler divergence gives where for θ ∈ Θ[s n ], let A := A(θ) = {i : θ i = 0}. In what follows, we will bound the second and the third terms in the rightmost side of the above equality.

A.5.5. Proof of Proposition 2.4 (b)
The proof follows almost the same line as in the proof for Proposition 2.4 (a).

B.3. Proof of Theorem 3.3
For the proof of Case (a), replacing e −rθ [j] in (8) and (9) by e −r [j] θ [j] completes the proof. For the proof of Case (b), replacing e −rθ [j] in (10) and (11)  We provide only the proof for (a) and omit the proof for (b). Following the same argument as in Subsection A.2 completes the proof for case (b).
Step 1: Lower bound on R(Θ[s n ]) A lower bound on R(Θ[s n ]) builds upon Bayes risk maximization based on a varied-spike block-independent prior. Let Π VB,ν with ν = (ν (1) , . . . , ν (sn) ) ∈ R mn × R mn × · · · × R mn × R n−mnsn be a varied-spike block-independent prior built as follows: divide {1, 2, . . . , n} into contiguous blocks B j (j = 1, 2, . . . , s n ) with each length m n := η −1 n . In each block B j , draw (θ 1+mn(j−1) , . . . , θ mnj ) independently according to a single spike prior with spike strength parameter ν (j) ∈ R sn + , where a single spike prior with spike strength parameter ν (j) ∈ R sn + is the distribution of ν (j) I e I with a uniformly random index I ∈ {1, . . . , m n } and a unit length vector e i in the i-th coordinate direction. Finally, set θ i = 0 for the remaining n−m n s n components. It is worth noting that a varied-spike block-independent prior may be different from a block-independent prior since the spike strength may be varied in each coordinate.

B.4.2. Proof of Lemma B.1
The proof is divided into two parts: deriving the minimax risk with r i 's conditioned, and taking expectation of that with respect to G. Recall that C := n i=1 C i /n. The former part of the proof essentially is completed by Lemma 5.2.
Consider the latter part of the proof. Take the expectation of C with respect to G. Since C is in (0,1), Hoeffding's inequality gives This implies that for any δ ∈ (0, 1), we have with probability higher than 1 − δ. In combining (29) with the result of the former part, we make an additional analysis on the behaviors of C 1 and C 2 in (27); see the sentence just before (27) for possible choices of these constants. Condition 3.1 assures that This implies that combining (29) with the result of the former part yields the inequality E G C + {1/(2n)} log(2/δ) + b n s n log η −1 n with probability higher than 1−δ, where b n is an o p (1) term that is independent of δ. This completes the proof.

B.4.3. Proof of Lemma B.2
The proof follows essentially the same steps as those in Subsection A.3. We provide the proof only for Θ[s n ]. For the simplicity, we conduct the proof with L replaced by 1. Following the same argument as that of Theorem 5.2, we have, for i = 1, 2, . . . , n, The above equalities forθ Π [ηn,κ] show that there exists a positive constant C 1 depending only on κ and sup i r i for which we have which completes the proof.

B.5. Proof of Proposition 3.1
The proof is a combination of the following lemma and the same argument as in Subsection B.1: hold as n → ∞ and η n → 0, where ε n is any shrinking sequence such that ε n = o(η n ).
Once we obtain Lemma B.3, the remaining part of the proof of Proposition 3.1 is easy and thus we complete the proof.

Proof of Lemma B.3:
In what follows, we prove Lemma B.3. We provide the proof only for Θ[s n ]. The proof follows almost the same line of the proof of Lemma 5.2. In bounding E(Θ[s n ]) below, we employ the Bayes risk based on the varied-spike block-independent prior Π VB,ν with ν = (1/r 1 , . . . , 1/r n ) ∈ R n + . In bounding E(Θ[s n ]) above, we use a threshold estimator. The threshold estimator used in this proof is the same as in Subsection A.5.1: for each i = 1, . . . , n, θ T,i (x i ) = η n if x i = 0 andθ T,i (x i ) = (x i + 1)/t i if x i 1.
Here for a subset J ⊂ {1, . . . , n}, e J indicates the vector of which the i-th component is 1 if i ∈ J and 0 if otherwise. We examine two cases (n, s, r) = (200,5,20) and (n, s, r) = (200, 20,20), and generate 500 current observations X's and 500 future observations Y 's. Table 6 Comparison of predictive densities without MCAR, with (n, s, r) = (200, 5,20), and with quasi-sparsity: for each result, the averaged value is followed by the corresponding standard deviation. Underlines indicate the best performance. The same abbreviations as in Table 1 are used.   20,20), and with quasi-sparsity: for each result, the averaged value is followed by the corresponding standard deviation. Underlines indicate the best performance. Tables 6 and 7 show that the performance of the proposed predictive densities does not depend on whether a parameter is exact sparse or quasi-sparse.

C.2. Effect of s
This subsection provides simulation studies highlighting the effect of s. The set-up except for s is the same as that in Subsection 4.1.1. Table 8 Comparison of predictive densities without MCAR, with (n, s, r) = (200, 50,20), and with exact sparsity: for each result, the averaged value is followed by the corresponding standard deviation. Underlines indicate the best performance.  Tables 8 and 9 find that the Bayes predictive density based on a Gauss hypergeometric prior works better in the predictive likelihood sense as the sparsity s n is relatively large.
we complete the proof. Differentiating E[log{p(Z t | θ; t)/p Π (Z t ; t)}] with respect to t yields Let e i be the unit length vector in the i-th coordinate direction (i = 1, . . . , n).
Together with the simple fact that p Π (Z t + e i ; t)/p Π (Z t ; t) =θ Π,i (Z t ; t), for any function f : N n → R) yields Similarly, the identity (d/dt) log p Π ( Combining these identities (32) and (33) with (31) gives (30), which completes the proof.  Figure 7 shows that r i s are varied but concentrated around the certain value. The summary statistics are as follows: the mean is 113783.36; the median is 119236; the standard deviation is 15888.35; the skewness is -4.02; the kurtosis is 17.42; r i s are bounded above by 121412. r i s are bounded below by 12752. From Figure 7, we see that most of locations are read from 100,000 to 120,000 times while a few locations are read less than 100,000 times. These differences can be explained by a mechanism of a sequencer: a DNA sequencer randomly reads fragments of a DNA sequence numerous times, and then reconstructs the whole sequence by combining all read fragments; hence, the numbers of individuals read at genetic locations are regarded as random variables and there occurs a difference in each value of them.


Privacy-preserving Neural Representations of Text

This article deals with adversarial attacks towards deep learning systems for Natural Language Processing (NLP), in the context of privacy protection. We study a specific type of attack: an attacker eavesdrops on the hidden representations of a neural text classifier and tries to recover information about the input text. Such scenario may arise in situations when the computation of a neural network is shared across multiple devices, e.g. some hidden representation is computed by a user’s device and sent to a cloud-based model. We measure the privacy of a hidden representation by the ability of an attacker to predict accurately specific private information from it and characterize the tradeoff between the privacy and the utility of neural representations. Finally, we propose several defense methods based on modified training objectives and show that they improve the privacy of neural representations.

Introduction
This article presents an adversarial scenario meant at characterizing the privacy of neural representations for NLP tasks, as well as defense methods designed to improve the privacy of those representations. A deep neural network constructs intermediate hidden representations to extract features from its input. Such representations are trained to predict a label, and therefore should contain useful features for the final prediction. However, they might also encode information about the input that a user wants to keep private (e.g. personal data) and can be exploited for adversarial usages.
We study a specific type of attack on neural representations: an attacker eavesdrops on the hidden representations of novel input examples (that are not in the training set) and tries to recover information about the content of the input text ( Figure 1). A typical scenario where such attacks would occur is when the computation of a deep neural net is shared between several devices (Li et al., 2017). For example, a user's device computes a representation of a textual input, and sends it a to cloudbased neural network to obtain, e.g. the topic of the text or its sentiment. The scenario is illustrated in Figure 1.
Private information can take the form of key phrases explicitly contained in the text. However, it can also be implicit. For example, demographic information about the author of a text can be predicted with above chance accuracy from linguistic cues in the text itself (Rosenthal and McKeown, 2011;Preoţiuc-Pietro et al., 2015).
Independently of its explicitness, some of this private information correlates with the output labels, and therefore will be learned by the network. In such a case, there is a tradeoff between the utility of the representation (measured by the accuracy of the network) and its privacy. It might be necessary to sacrifice some accuracy in order to satisfy privacy requirements.
However, this is not the case of all private information, since some of it is not relevant for the prediction of the text label. Still, private information might be learned incidentally. This nonintentional and incidental learning also raises privacy concerns, since an attacker with an access to the hidden representations, may exploit them to recover information about the input.
In this paper we explore the following situation: (i) a main classifier uses a deep network to predict a label from textual data; (ii) an attacker eavesdrops on the hidden layers of the network and tries to recover information about the input text of unseen examples. In contrast to previous work about neural networks and privacy (Papernot et al., 2016;Carlini et al., 2018) we do not protect the privacy of examples from the training set, but the privacy of unseen examples provided, e.g., by a user.
An example of a potential application would be a spam detection service with the following constraints: the service provider does not access verbatim emails sent to users, only their vector representations. Theses vector representations should not be usable to gather information about the user's contacts or correspondents, i.e. protect the user from profiling. This paper makes the following contributions: 1 • We propose a metric to measure the privacy of the neural representation of an input for Natural Language Processing tasks. The metric is based on the ability of an attacker to recover information about the input from the latent representation only. • We present defense methods designed against this type of attack. The methods are based on modified training objectives and lead to an improved privacy-accuracy tradeoff.

Adversarial Scenario
In the scenario we propose, each example consists of a triple (x, y, z), where x is a natural language text, y is a single label (e.g. topic or sentiment), and z is a vector of private information contained in x. Our base setting has two entities: (i) a main classifier whose role is to learn to predict y from x, (ii) an attacker who learns to predict z from the latent representation of x used by the main classifier. We illustrate this setting in Figure 1. In order to evaluate the utility and privacy of a specific model, we proceed in three phases: Phase 1. Training of the main classifier on (x, y) pairs and evaluation of its accuracy; Phase 2. Generation of a dataset of pairs (r(x), z) for the attacker, r is the representation function of the main classifier (r is defined in Section 2.1); Phase 3. Training of the attacker's network and evaluation of its performance for measuring privacy.
In the remainder of this section, we describe the main classifier (Section 2.1), and the attacker's model (Section 2.2).

Text Classifier
As our base model, we chose a standard LSTM architecture (Hochreiter and Schmidhuber, 1997) for sequence classification. LSTM-based architectures have been applied to many NLP tasks, including sentiment classification (Wang et al., 2016) and text classification (Zhou et al., 2016).
First, an LSTM encoder computes a fixed-size representation r(x) from a sequence of tokens x = (x 1 , x 2 , . . . , x n ) projected to an embedding space. We use θ r to denote the parameters used to construct r. They include the parameters of the LSTM, as well as the word embeddings. Then, the encoder output r(x) is fed as input to a feedforward network with parameters θ p that predicts the label y of the text, with a softmax output activation. In the standard setting, the model is trained to minimize the negative log-likelihood of y labels: where N is the number of training examples.

Attacker's Classifier
Once the main model has been trained, we assume that its parameters θ r and θ p are fixed. We generate a new dataset made of pairs (r(x), z(x)), where r(x) is the hidden representation used by the main model and z(x) is a vector of private categorical variables. In practice, z is a vector of binary variables, (representing e.g. demographic information about the author). In our experiments, we use the same training examples x for the main classifier and the classifier of the attacker. However, since the attacker has access to the representation function r parameterized by θ r , they can generate a dataset from any corpus containing the private variables they want to recover. In other words, it is not necessary that they have access to the original training corpus to train their classifier.
The attacker trains a second feedforward network on the new dataset {(r(x (i) ), z (i) )} i≤N . This classifier uses a sigmoid output activation to compute the probabilities of each binary variable in z: P (z|r(x); θ a ) = σ(FeedForward(r(x))).
It is trained to minimize the negative loglikelihood of z: assuming that the K variables in z are independent. Since the parameters used to construct r are fixed, the attacker only acts upon its own parameters θ a to optimize this loss.
We use the performance of the attacker's classifier as a proxy for privacy. If its accuracy is high, then an eavesdropper can easily recover information about the input document. In contrast, if its accuracy is low (i.e. close to that of a mostfrequent label baseline), then we may reasonably conclude that r does not encode enough information to reconstruct x, and mainly contains information that is useful to predict y.
In general, the performance of a single attacker does not provide sufficient evidence to conclude that the input representation r is robust to an attack. It should be robust to any type of reconstruction method. In the scope of this paper though, we only experiment with a feedforward network reconstructor, i.e. a powerful learner.
In the following sections, we propose several training method modifications aimed at obfuscating private information from the hidden representation r(x). Intuitively, the aim of these modifications is to minimize some measure of information between r and z to make the prediction of z hard. An obvious choice for that measure would be the Mutual Information (MI) between r and z. However, MI is hard to compute due to the continuous distribution of r and does not lend itself well to stochastic optimization.

Defenses Against Adversarial Attacks
In this section, we present three training methods designed as defenses against the type of attack we described in Section 2.2. The first two methods are based on two neural networks with rival objective functions (Section 3.1). The last method is meant at discouraging the model to cluster together training examples with similar private variables z (Section 3.2).

Adversarial Training
First, we propose to frame the training of the main classifier as a two-agent process: the main agent and an adversarial generator, exploiting a setting similar to Generative Adversarial Networks (GAN, Goodfellow et al., 2014). The generator learns to reconstruct examples from the hidden representation, whereas the main agent learns (i) to perform its main task (ii) to make the task of the generator difficult.
We experiment with two types of generators: a classifier that predicts the binary attributes z(x) used as a proxy for the reconstruction of x (Section 3.1.1) and a character-based language model that directly optimizes the likelihood of the training examples (Section 3.1.2).

Adversarial Classification: Multidetasking
In order not to make r(x) a good representation for reconstructing z, we make two modifications to the training setup of the main model (Phase 1): • We use a duplicate adversarial classifier, with parameters θ a , that tries to predict z from r(x). It is trained simultaneously with the main classifier. Its training examples are generated on the fly, and change overtime as the main classifier updates its own parameters. This classifier simulates an attack during training. • We modify the objective function of the main classifier to incorporate a penalty when the adversarial classifier is good at reconstructing z. In other words, the main classifier tries to update its parameters so as to confuse the duplicate attacker.
Formally, for a single data point (x, y, z), the adversarial classifier optimizes: whereas the main classifier optimizes: The first term of this equation is the log-likelihood of the y labels. The second term is designed to deceive the adversary. The hyperparameters α > 0 and β > 0 control the relative importance of both terms.
As in a GAN, the losses of both classifiers are interdependent, but their parameters are distinct: the adversary can only update θ a and the main classifier can only update θ r and θ p .
The duplicate adversarial classifier is identical to the classifier used to evaluate privacy after the main model has been trained and its parameters are fixed. However, both classifiers are completely distinct: the former is used during the training of the main model (Phase 1) to take privacy into account whereas the latter is used to evaluate the privacy of the final model (Phase 3), as is described in Section 2.

Adversarial Generation
The second type of generator we use is a characterbased LSTM language model that is trained to reconstruct full training examples. For a single example (x; y), the hidden state of the LSTM is initialized with r(x), computed by the main model. The generator optimizes: where θ is the set of parameters of the LSTM generator, x i is the i th character in the document, and C is the length of the document in number of characters. The generator has no control over r(x), and optimizes the objective only by updating its own parameters θ . Conversely, the loss of the main model is modified as follows: The first term maximizes the likelihood of the y labels whereas the second term is meant at making the reconstruction difficult by maximizing the loss of the generator. As in the loss function described in the previous section, α and β control the relative importance of both terms. Once again, the main classifier can optimize the second term only by updating θ r , since it has no control over the parameters of the adversarial generator.
A key property of this defense method is that it has no awareness of what the private variables z are. Therefore, it has the potential to protect the neural representation against an attack on any private information. From a broader perspective, the goal of this defense method is to specialize the hidden representation r(x) to the task at hand (sentiment or topic prediction) and to avoid learning anything not relevant to it.

Declustering
The last strategy we employ to make the task of the attacker harder is based on the intuition that private variables z are easier to predict from r when the main model learns implicitly to cluster examples with similar z in the same regions of the representation space.
In order to avoid such implicit clustering, we add a term to the training objective of the main model that penalizes have hidden representations r(x) and r(x ) in the same region of space. We use the following modified loss for a single example: where (x , z ) is another example sampled uniformly from the training set, α is a hyperparameter controlling the importance of the second term, and (·, ·) ∈ [0, 1] is the normalized Hamming distance.

Experiments
Our experiments are meant to characterize the privacy-utility tradeoff of neural representations on text classification tasks, and evaluating if the proposed defense methods have a positive impact on it. We first describe the datasets we used (Section 4.1) and the experimental protocol (Section 4.2), then we discuss the results (Section 4.3). We found that in the normal training regime, where no defense is taken into account, the adversary can recover private information with higher accuracy than a most frequent class baseline. Furthermore, we found that the de- fenses we implemented have a positive effect on the accuracy-privacy tradeoff.

Datasets
We experiment with two text classification tasks: sentiment analysis (Section 4.1.1) and topic classification (Section 4.1.2). The sizes of each dataset are summarized in Table 1.

Sentiment Analysis
We use the Trustpilot dataset  for sentiment analysis. This corpus contains reviews associated with a sentiment score on a five point scale, and self-reported information about the users. We use the five subcorpora corresponding to five areas (Denmark, France, Germany, United Kingdom, United States).
We filter examples containing both the birth year and gender of the author of the review and use these variables as the private information. As in previous work on this dataset (Hovy, 2015;, we bin the age of the author into two categories ('under 35' and 'over 45'). Finally, we randomly split each subcorpus into a training set (80%), a development set (10%) and a test (10%).
As an additional experimental setting, we use both demographic variables (gender and age) as input to the main model. We do so by adding two additional tokens at the beginning of the input text, one for each variable. It has been shown that those variables can be used to improve text classification (Hovy, 2015). Also, we would like to evaluate whether the attacker's task is easier when the variables to predict are explicitly in the input, compared to when these information are only potentially and implicitly in the input. In other words, this setting simulates the case where private in-formation may be used by the model to improve classification, but should not be exposed too obviously. In the rest of this section, we use RAW to denote the setting where only the raw text is used as input and +DEMO, the setting where the demographic variables are also used as input.

Topic Classification
We perform topic classification on two genres of documents: news articles and blog posts.
News article For topic classification of news article, we use two datasets: the AG news corpus 2 (Del Corso et al., 2005) and the English part of the Deutsche Welle (DW) news corpus (Pappas and Popescu-Belis, 2017).
For the AG corpus, following Zhang et al. (2015), we construct the dataset by extracting documents belonging to the four most frequent topics, and use the concatenation of the 'title' and 'description' fields as the input to the classifier. We randomly split the corpus into a training set (80%), a development set (10%) and a test set (10%). For the DW dataset, we use the 'text' field as input, and the standard split. We kept only documents belonging to the 20 most frequent topics.
The attacker tries to detect which named entities appear in the input text (each coefficient in z(x) indicates whether a specific named entity occurs in the text). For both datasets, we used the named entity recognition system from the NLTK package (Bird et al., 2009) to associate each example with the list of named entities that occur in it. We select the five most frequent named entities with type 'person', and only keep examples containing at least one of these named entities. This filtering is necessary to avoid a very unbalanced dataset (since each selected named entity appears usually in very few articles).
Blog posts We used the blog authorship corpus presented by Schler et al. (2006), a collection of blog posts associated with the age and gender of the authors, as provided by the authors themselves. Since the blog posts have no topic annotation, we ran the LDA algorithm (Blei et al., 2003) on the whole collection (with 10 topics). The LDA outputs a distribution on topics for each blog post. We selected posts with a single dominating topic (> 80%) and discarded the other posts. We binned age into two category (under 20 and over 30). We  used the age and gender of the author as the private variables. These variables have a very unbalanced distribution in the dataset, we randomly select examples to obtain uniform distributions of private variables. Finally, we split the corpus into a training set (80%), a validation set and a test set (10% each).

Protocol
Evaluation For the main task, we report a single accuracy measure. For measuring the privacy of a representation, we compute the following metrics: • For demographic variables (sentiment analysis and blog post topic classification): 1 − X, where X is the average of the accuracy of the attacker on the prediction of gender and age; • For named entities (news topic classification): 1−F , where F is an F-score computed over the set of binary variables in z that indicate the presence of named entities in the input example.
Training protocol We implemented our model using Dynet (Neubig et al., 2017 For each dataset, and each LSTM state dimension ({8, 16, 32, 64, 128}), we train the main model for 8 epochs (sentiment classification) or 16 epochs (topic classification), and select the model with the best accuracy on the development set. Then, we generate the dataset for the attacker, train the adversarial model for 16 epochs and select the model with the worst privacy on the development set (i.e. the most successful attacker).
It has to be noted that we select the models that implement defenses on their accuracy, rather than their privacy or a combination thereof. In practice, we could also base the selection strategy on a privacy budget: selecting the most accurate model with privacy above a certain threshold.

Results
This section discusses results for the sentiment analysis task (Section 4.3.1) and the topic classification task (Section 4.3.2).

Sentiment Analysis
How private are neural representations? Before discussing the effect of proposed defense methods, we motivate empirically our approach by showing that adversarial models can recover private information with reasonable accuracy when the attack is targeted towards a model that implements none of the presented defense methods.
To do so, we compare the accuracy of adversarial models to two types of baselines: • As a lower bound, we use the most frequent class baseline. • As an upper bound, we trained a classifier that can optimize the hidden representations (r) for the attacker's tasks. In other words, this baseline is trained to predict de-  Table 3: Results on the test sets of the Trustpilot dataset, +DEMO setting. Main is the accuracy on sentiment analysis. Priv. is the privacy measure (i.e. the inverse accuracy of the attacker: higher is better, see Section 4.2). The baselines are most-frequent class classifiers. The values reported for the defense methods indicate absolute differences with the standard training regime (no defense implemented) for both metrics.
mographic variables from x, as if it were the main task.
In Table 2, we compare both baselines to the best adversary in the two settings (RAW and +DEMO) among the models trained with no defenses. First of all, we observe that apart from gender on the German dataset, the trained baseline outperforms the most frequent class baseline by a wide margin (8 to 25 absolute difference). Second of all, the attacker is able to outperform the most frequent class baseline overall, even in the RAW setting. In more details, for age, the adversary is well over the baseline in all cases except US. On the other hand, gender seems harder to predict: the adversary outperforms the most frequent class baseline only in the +DEMO setting.
The same pattern is visible for the blog post dataset, also presented in the last line of Table 2: the best adversaries are 14 points over the baseline for gender and 5 points for age, i.e. almost as good as a model that can fine tune the hidden representations.
These results justify our approach, since they demonstrate that hidden representations learn private information about the input, and can be exploited to recover this information with reasonable accuracy.
Effect of defenses We report results for the main task accuracy and the representation privacy in Table 3 for the +DEMO setting and in Table 4 for the RAW setting. Recall that the privacy measure   (Priv.) is computed by 1 − X where X is the average accuracy of the attacker on gender and age predictions. When this privacy metric is higher, it is more difficult to exploit the hidden representation of the network to recover information about x. The 'Standard' columns contain the accuracy and privacy of the base model described in Section 2. The next columns present the absolute variation in accuracy and privacy for the three defense methods presented in Section 3: Multidetasking, Adversarial Generation, and Declustering. We also report for each corpus the most frequent class baseline for the main task accuracy, and the privacy of the most frequent class baselines on private variables (i.e. the upper bound for privacy). The three modified training methods designed as defenses have a positive effect on privacy. Despite a model selection based on accuracy, they lead to an improvement in privacy on all datasets, except on the France subcorpus. In most cases, we observe only a small decrease in accuracy, or even an improvement at times (e.g. multidetasking on the Germany dataset, RAW setting), thus improving the tradeoff between the utility and the privacy of the text representations.

Topic Classification
We report results on topic classification in Table 5.
News articles For the news corpora, the privacy metric is based on the F-score on the binary variables z indicating the presence or absence of a named entity in the text. First of all, we observe that defense methods that explicitly use z (i.e. multidetasking and declustering), have a very positive effect on privacy, but also a detrimental effect on the main task. We hypothesize that this is due to the strong correlations between the main task labels y and the private information z. As a result, improving the privacy of the neural representations comes at a cost in accuracy.
In contrast, the adversarial generation defense method lead to an improvement in accuracy, that is quite substantial for the DW corpus. We speculate that this is due to the secondary term in the objective function of the main model (Section 3.1.2) that helps avoiding overfitting the main task or learning spurious features.
Blog posts On the blog post dataset, the effects are smaller, which we attribute to the nature of the task of the attacker. The defense methods consistently improve privacy and, in one case, accuracy. The best effects on the tradeoff are achieved with the multidetasking and adversarial generation methods.

Discussion
The main result of our experiments is that the defenses we propose improve privacy with usually a small effect, either positive or negative, on accuracy, thus improving the tradeoff between the utility and the privacy of neural representations.
An important direction for future work is the choice of a strategy for model selection. The tradeoff between utility and privacy can be controlled in many ways. For example, the importance of both terms in the loss functions in Section 3.1 can be controlled to favor either privacy or utility. In the scope of this paper, we did not perform thorough hyperparameter tuning, but believe that doing so is important for achieving better results, since the effects of defense method can be more drastic than desired in some cases, as exemplified on the news corpora (Table 5).
Overall, we found that the multidetasking approach lead to the more stable improvements and should be preferred in most cases, since it is also the less computationnally expensive defense. On the other hand, the adversarial generation method does not require the specification of private variables, and thus is a more general approach.

Related Work
The deployment of machine learning in both academic and industrial contexts raises concerns about adversarial uses of machine learning, as well as concerns about attacks specifically targeted at these algorithms that often rely on large amounts of data, including personal data.
More generally, the framework of differential privacy (Dwork, 2006) provides privacy guarantees for the problem of releasing information without compromising confidential data, and usually involves adding noise in the released information. It has been applied to the training of deep learning models Papernot et al., 2016;Papernot et al., 2018), and Bayesian topic models (Schein et al., 2018).
The notion of privacy is particularly crucial to NLP, since it deals with textual data, oftentimes user-generated data, that contain a lot of private information. For example, textual data contain a lot of signal about authors (Hovy and Spruit, 2016). and can be leveraged to predict demographic variables (Rosenthal and McKeown, 2011;Preoţiuc-Pietro et al., 2015). Oftentimes, this information is not explicit in the text but latent and related to the usage of various linguistic traits. Our work is based on a stronger hypothesis: this latent information is still present in vectorial representations of texts, even if the representations have not been supervised by these latent variables. Li et al. (2017) study the privacy of unsupervised representations of images, and measures their privacy with the peak signal to noise ratio between an original image and its reconstruction by an attacker. They find a tradeoff between the privacy of the learned representations and the accuracy of an image classification model that uses these representations as inputs. Our setting is complementary since it is applied to NLP tasks, but explores a similar problem in the case of representations learned with a task supervision.
A related problem is the unintended memorization of private data from the training set and has been addressed by Carlini et al. (2018). They tackle this problem in the context of text generation (machine translation, language modelling).
If an attacker has access to e.g. a trained language model, they are likely to be able to generate sentences from the training set, since the language model is trained to assign high probabilities to those sentences. Such memorization is problematic when the training data contains private information and personal data. The experimental setting we explore is different from these works: we assume that the attacker has access to a hidden layer of the network and tries to recover information about an input example that is not in the training set.
In a recent study, Li et al. (2018) proposed a method based on GAN designed to improve the robustness and privacy of neural representations, applied to part-of-speech tagging and sentiment analysis. They use a training scheme with two agents similar to our multidetasking strategy (Section 3.1.1), and found that it made neural representations more robust and accurate. However, they only use a single adversary to alter the training of the main model and to evaluate the privacy of the representations, with the risk of overestimating privacy. In contrast, once the parameters of our main model are fixed, we train a new classifier from scratch to evaluate privacy.

Conclusion
We have presented an adversarial scenario and used it to measure the privacy of hidden representations in the context of two NLP tasks: sentiment analysis and topic classification of news article and blog posts. We have shown that in general, it is possible for an attacker to recover private variables with higher than chance accuracy, using only hidden representations. In order to improve the privacy of hidden representations, we have proposed defense methods based on modifications of the training objective of the main model. Empirically, the proposed defenses lead to models with a better privacy.