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.