1.1 Background and Summary

Consider the problem of modeling a dataset of numbers accidents in a population of insured persons, or the incidence of an illness in a population. One can consider a spectrum of levels of granularity in describing the data and the corresponding population, from histograms, to a single distribution from a parametric family, to a bimodal distribution, to a mixture of two or more distributions, to modeling the population at the individual level via compound (predictive) distributions. Examples included are a dataset of employee days ill (a discrete variable) and a dataset of family expenditures on groceries (a continuous variable). Fits to the data obtained by various levels of granularity are compared using the Bayesian Information Criterion.

1.2 Levels of “Granularity”

In discussing the level at which to analyze a dataset, it will be shown how to go from individuals to histograms and modes to clusters or mixture components, back to individuals. The various levels proceed from histograms to clusters to predictive distributions

That is, one can go from modes to sub-populations back to individuals. These choices might be called choosing the level of “granularity” at which to analyze the dataset. Methods for various levels of granularity include those indicated in the table.

There will be two extended analyses of datasets. These datasets are neither new nor large nor high-dimensional but hopefully will still be found to be interesting, particularly from the point of view of this paper. These two datasets are:

• •

  Expenditure in a week on fruits and vegetables for 60 English families (Connor & Morrell, Statistics in Theory & Practice)

• •

  Days ill in a year for 50 miners (hypothetical data from Kenkel, Statistics for Management & Economics)

2.1 Extreme values

But, if one does fit a Poisson to the dataset, omitting, say, the upper two values 17 and 18 days, what conclusions might be drawn? Would the 17 and 18 be considered to be particularly unusual? That is to say, later, having fit a distribution, we can assess the probability of such extreme observations.

For now, we note that the mean of the other 50-2 = 48 observations is 5.1 days; the variance, 17.47, still very different values, so a single Poisson would not provide a good fit even after omitting the two largest, possibly outlying, observations. Note that also, besides the gap at 15 and 16 days, there is a gap at 10 days, nobody having been ill that number of days. Perhaps a mixture of distributions would make more sense.

Figure 1:

Histogram of days ill

2.2 Poisson mixture model

Consequently, a mixture of two Poissons was fit. The mixture model has p.m.f. (probability mass function) p(x) = π₁ p₁(x) + π₂ p₂(x), where p₁(·) is the p.m.f. of a Poisson distribution with parameter λ₁ and p₂(·) is the p.m.f. of a Poisson distribution with parameter λ₂.

Parameter estimates for the Poisson mixture model. Starting values for iterative estimation were obtained by a Student’s t method of clustering into two groups (Sclove 2016). In this method, each cut-point is tried.

In more detail: Obtain the order statistic, resulting from ordering the the observations {x₁, x₂, … , x_n} as

Candidate cut-off points c₁, c₂, … , c_n−1 are chosen in between the order statistics:

For example, one can take c_j to be the midpoint c_j = (x_(j) + x₍j + 1))/2. One then computes two-sample t for each clustering, that is, for each cut-point. The best clustering into two clusters is the one which gives the largest value of |t|, or, equivalently, of t². This is because t² compares the within-groups sum of squares and the between-groups sum of squares.

Pooled t assumes equal variances in the two clusters. Unpooled t avoids this assumption. Not assuming equal variances, one can consider the standardized difference between means as an objective function, that is, one can use an unpooled two-sample t as the criterion, maximizing it over cut points. The unpooled t statistic is given by t2=(x¯1−x¯2)2/(s12/n1+s22/n2) , that is, |t|=|x¯1−x¯2|/s12/n1+s22/n2 .

The value of the two-sample, unequal variance Student’s t is computed for each cut-point, and the cut-point giving the largest t² is chosen.

Iteration. The means of the resulting clusters are taken as tentative estimates of the distributional means; the cluster relative frequencies, as tentative estimates of the mixing probabilities. These estimates were 2.1 days, 8.9 days, with cluster frequencies of 17 out of 50 and 33 out of 50, that is, mixing probabilities of .34, .66. By doing a grid search in the vicinity of these initial estimates to maximize the mixture-model likelihood, the following estimates were obtained: λˆ1=2.77 days , πˆ1=.40 , λˆ2=9.12 days , πˆ2=.60 . Also, taking these as starting points for the EM (Expectation-Maximization) algorithm, only slightly different estimates were obtained; these were λˆ1=2.84 days , λˆ2=9.20 days , πˆ1=.41 , πˆ2=.59 .

Extreme value assessment with the fitted Poisson mixture. With the fitted Poisson mixture, the estimate of Pr{X > 16} is .007, fairly small. So perhaps results of 17 or more days could and should be considered outliers.

2.3 Comparison of models by model-selection criteria

The two fits, by histogram and by Poisson mixture, were compared by means of model-selection criteria. Given K alternative models, indexed by k = 1, 2, … , K, penalized-likelihood model-selection criteria are smaller-is-better criteria that can be written in the form

where m_k is the number of free parameters used in fitting Model k, LL_k is the log maximum likelihood of Model k, and a(n) = ln n for BIC (Bayesian Information Criterion; Schwarz 1978) and a(n) = 2 for all n for AIC (Akaike’s Information Criterion; Akaike 1973, 1974; Kashyap 1982; Sakamoto 1992). That is, for k = 1, 2, … , K alternative models, and

The number of parameters for the Poisson mixture is 2 means plus 2 mixing probabilities, less 1 because the probabilities must add to 1. That is 3 free parameters for the Poisson mixture. The number of parameters for the histogram, scored by the multinomial distribution with 17 categories (0 through 18, but 15 and 16 are missing), less 1 because the multinomial probabilities must add to 1, leaving 16 free parameters.

The results are in the next table. The histogram wins by a bit according to AIC, but the Poisson mixture wins by a wide margin according to BIC. To see this, note that BIC is derived (Schwarz 1978) as the first terms in the Taylor series expansion of (-2 times) the posterior probability of Model k, Pr(Model k | data) = pp_k, say. That is,

To calculate the posterior probabilities, one subtracts a large constant from each, divides by 2, exponeniates the negative of this, and sums these, dividing by the sum to normalize.

Different bin widths. If one makes several histograms, with different bin widths, how should the likelihood for histograms be computed? Given data points x₁, x₂, … , x_n, the likelihood for a given p.m.f. p(·) is

Here p(x_i) is the p.m.f. at the data point x_i. (For continuous data, we would write the p.d.f., f(x_i). ) But in the context of histograms what we can take p(x_i) to be?

Denote the number of bins by J. Let the bin width be denoted by h. This is an increment along the x-axis.

Let the bins be indexed by j, j = 1, 2, … , J. The class limits are x₀, x₀ + h, x₀ + 2h, … , x₀ + Jh. The class intervals (bins) are [x₀, x₀ + h), [x₀ + h, x₀ + 2h), … , [x₀ + (J − 1)h, x₀ + Jh). In the present application, x₀ = 0. Now, let j(x_i) denote the bin containing x_i and n_j(xi) be the frequency in that bin. To approximate f(x_i), motivated by f(x₁) ≈ [F(x₂) − F(x₁)]/(x₂ − x₁) = [F(x₁ + h) − F(x₁)]/[(x_i + h) − x₁] = [F(x₁ + h) − F(x₁)]/h, write

and

That is, the concept is that probability density is probability per unit length along the x axis. Thus the likelihood is

Note that p(x1,x2, … , xn)=∏i=1npj(xi)=∏j=1Jpjnj is a multinomial p.m.f. with probabilities p_j and frequencies n_j for the J categories. The maximized likelihood L is this multinomial (with p_j estimated as n_j/n), divided by hⁿ, which may be viewed as an adjustment to the likelihood due to the bin width h. In computing the likelihood, the probability density is to be used, where “density” is probability / bin width. Note that with a continuous variable we would compute probability density as f(x_i)/h, that is, f(x_i)/(Lebesgue measure of the bin interval), whereas with a discrete variable we are really computing probability density as p(x_i)/h, where now h is the counting measure of the bin interval.

Varying bin widths. In the case of non-constant bin widths, with a bin-width of h_j for the j-th interval, take the probability density at x_i to be n_j(i)/h_j(i), where h_j(i) is the width of the interval in which x_i falls and n_j(i) (short for n_j(xi) ) is the frequency (count) in that interval. The likelihood is

According to AIC, the histogram with varying bin widths wins, the Poisson mixture coming in second. According to BIC (and, equivalentlly, posterior probability), the Poisson mixture scores the best, by far. But the point is not just which model wins, but that such a comparison, comparing histograms on the one hand with fitted distributions on the other, can be made.

Levels of granularity, cont’d. Perhaps another level of granularity is approached by predictive distributions, which may be viewed as getting to the individual level of granularity. Predictive distributions may be viewed in the light of compound distributions resulting from a prior distribution on the parameter at the individual level. From the viewpoint of modern statistics, a predictive distribution is merely the marginal distribution of the observable random variable, having integrated out the prior on the parameter. (Details to follow.)

The Yule-Greenwood model approaches modeling at the individual level, stating that each individual may have his or her own accident rate λ and so is an example of a compound model. In terms of granularity, the Yule-Greenwood model is a classical example at the level of the individual in that it employs a Poisson model for each individual’s accident rate λ and then puts a (Gamma) distribution over the population of values of λ. The model is the Gamma-Poisson model (sometimes called the Poisson-Gamma model) and is a prime example of a compound model. The Gamma is a conjugate prior distribution for the Poisson, meaning that the posterior distribution of λ is also a member of the Gamma family. We discuss this further below; first, however, we fit histograms and mixtures to a dataset with continuous data.

3.1 Fitting a Gamma distribution

The two-parameter Gamma p.d.f. is

The mean is mβ. The variance is mβ². Method-of-moments estimates are, for the scale parameter β = σ²/μ, βˆ=s2/x¯=0.2081/1.022=0.2035 . and for the shape parameter m = μ/β, so mˆ=x¯/βˆ=1.022/0.2035=5.0246 .

3.2 Fitting mixture models

3.3 Comparison of models by model-selection criteria

The table summarizes the results. According to AIC, the ranking is: Gamma; Gaussian mixture; Gamma mixture; histogram with bin width .2, histogram with bin width .1. According to BIC, the ranking is: Gaussian mixture, Gamma, Gamma mixture, histogram with bin width .2, histogram with bin width .1 The Gamma mixture has only a small posterior probability.

5.1 Probability functions for the Gamma-Poisson model

In the Gamma-Poisson model, the distribution of X is Poisson with parameter usually called λ. The p.m.f. is

The mean and variance are both equal to λ.

Such a distribution can be considered, say, for the number of accidents per individual per year. For the days ill dataset (days ill in a year for a sample of n = 50 miners), we have fit a single Poisson (with mean 6.58 days per year). We looked at histograms and observed bimodality. Further, the fact that the sample variance of 19.06 was considerably larger than the sample mean was a hint of inadequacy of a single Poisson. As discussed above, a mixture of Poissons was fitted, with mixing probabilities about .6 and .4 and means about 3 days and 9 days. A finer level of granularity would be obtained by saying that each person has his own value of λ and putting a distributon on these over the population.

5.2 Gamma family of distributions

A Gamma distribution could be a good choice. The choice of the Gamma family is non-restrictive in that the family can achieve a wide variety of shapes. The single-parameter gamma has a shape parameter m; the two-parameter gamma family has, in addition, a scale parameter, β. (The reciprocal of β is the rate parameter, so-called because it is the rate, or intensity, of the associated Poisson process.) A Gamma distribution with parameter m, has p.d.f. f(λ) = Const. λ^m−1 e^−λ, λ > 0. The constant is 1/Γ(m). More generally, the two-parameter Gamma can be used: the p.d.f. is

The mean is mβ. The variance is mβ².

The Negative Exponential family of distributions. The special case of m = 1 in the Gamma family gives the negative exponential family of distributions. So the

The mean is β. The variance is β².

5.3 Development of the Gamma-Poisson model

Putting a population distribution over a parameter can be a very helpful way of modeling. The resulting model is called a compound model. In a compound model, the random variable X is considered as the result of sampling that yields an individual and that individual’s value of a parameter, and then the individual’s value of X is observed, from a distribution with that value of the parameter. Note that a compound model can be viewed as an infinite mixture model.

In this discussion, focus is on a couple of particular compound models, the Gamma-Poisson, and later, the Beta-Binomial.

The Yule-Greenwood model, from a modern viewpoint, is an application of the Gamma-Poisson model to a financial, in fact, actuarial, situation. It is in terms of a model for accident rates in a population. Suppose that the yearly number of accidents of any given individual i in a population is distributed according to a Poisson distribution with parameter λ_i accidents per year. (This is count data, similar to the days ill data.) Then the probability that individual i, with parameter value λ_i, has exactly k accidents in a year, k = 0, 1, 2, … , is

Some individuals are more accident prone (have a higher accident rate) than others, so different individuals have different values of λ. A distribution can be put on λ to deal with this. This is the Yule-Greenwood model, dating from 1920; a precursor of the Predictive Distributions of the new Predictive Analytics, predating even Abraham Wald (1950) as a founder of modern mathematical statistics and decision theory and Jimmie Savage (1954) as a founder of modern Bayesian Statistics.

The standard choice of a prior distribution on λ is a Gamma distribution. The Gamma family is a conjugate family to the Poisson, meaning that the prior and posterior distributions of λ are both in the Gamma family.

The joint distribution of X and Λ. The joint probability function of X and Λ is

The expressions for the Gamma and Poisson are put into this. That is, the weight assigned to p_X||Λ(x|λ) is f_Λ(λ).

The joint probability function is used to obtain

• •

  the marginal distribution of X, by integrating out λ, and

• •

  then the posterior distribution of Λ given x, by dividing the joint probability function by the marginal probability mass function of X.

Putting in the expressions for the Gamma and Poisson, it is seen that the marginal distribution of X, the number of accidents that a randomly selected individual has in a year, is of the form

When the prior is Gamma and the conditional is Poisson, this marginal distribution can be shown to be negative binomial. Its parameters are m and p = 1/(1 + β).

In the Bayesian model, the parameter of the conditional distribution of X, say θ, is treated as a random variable Θ.

In the Gamma-Poisson model, θ is the Poisson parameter λ.

The conditional distribution of X given that Θ = θ is Poisson(λ). The probability mass function is

The joint p.d.f. of X and Λ can be written as p_X||Λ(x|λ) which is f_Λ(λ)(x, λ) = f_Λ p_X|Λ(x|λ.

As mentioned above, from this, the posterior distribution of Λ, that is, the distribution of Λ given x, can be computed, and the marginal distribution of X can be computed. Marginal maximum likelihood estimation can be applied to obtain estimates of the remaining parameters.

Posterior distribution of Λ. Analogous to Pr(B|A) = Pr(A ∩ B) / Pr(A), the p.d.f. of the posterior distribution is the joint p.d.f., divided by the marginal p.d.f. of X :

This will turn out to be a Gamma distribution, that is, it is in the same family as the prior. The Gamma is a conjugate prior for the Poisson.

Marginal distribution of X. In Predictive Analytics, the marginal distribution of X is computed as a model of a future observation or observations of X.

The marginal distribution of X is obtained by integrating the joint distribution with respect to the parameter. Note that this computation combines information, by weighting the conditional distribution of X given λ with the prior on λ. This computation of the p.d.f. is, as stated above, fX(x)=∫0∞t f(x|λ) fΛ(λ) dλ .

Moments. The mean of the marginal distribution of X is mq/p = m β. The variance of the marginal distribution of X is mq/p² = m β(1 + β).

5.4 Empirical Bayes estimation

Empirical Bayes estimation, at least in the present context, means estimating the parameters of the prior using observations from the marginal distribution.

The hyperparameters in terms of the moments of the marginal. The parameters of the prior are called hyperparameters. In this case, they are λ and β. Suppose we solve for them in terms of the first two moments of the marginal.

Estimating the prior parameters from the marginal. Estimates of the prior parameters m and β can be obtained by, for example, taking the expressions for the hyperparameters m and β in terms of the first two raw moments and plugging in estimates m1′ and m2′ . Given a sample X₁, X₂, … , X_n, we have m1′=X¯=∑i=1NXi/N and m2′=∑i=1nXi2/n .

5.5 Application to the days ill dataset

Returning to the days ill dataset for n = 50 miners, the days ill in a year ranged from 0 to 18; the distribution seems to be bimodal.

The p.m.f. of the Negative Binomial distribution with parameters m and p is where k is the number of trials in excess of m required to get m Heads. In the Gamma-Poisson model, the marginal distribution of X is Negative Binomial with parameters with parameters m and p = 1/(1 + β).

Given that the true mean of the marginal (“predictive distribution”) Negative Binomial is μ = mq/p = mβ and the true variance is σ² = mq/p² = mβ(1 + β), and the sample mean x¯=6.58 and the sample variance s² = 19.07, one can set up two equations and solve for method of moments estimates of the hyperparameters m and β in the Gamma prior for λ.

The equations are [1] : m β = 6.58; [2] : m β(1 + β) = 19.07.

Putting [1] in [2] gives 6.58(1 + β) = 19.07, 1 + β = 19.07/6.58 ≈ 2.898, betaˆ≈1.898 . Then m ≈ 6.58/β = 6.58/1.898 ≈ 3.467. Now, μ = m(1 − p)/p = m/p − m, μ + m = m/p, p = m/(μ + m) or, estimating p = 3.467/(6.58 + 3.467) = 3.467/10.05 = 0.345. So now we have estimates of the hyperparameters.

To estimate the mean and variance of the Gamma prior, one can proceed as follows. The mean of the prior is mβ, estimated as 6.58 days ill per year. The variance of the prior is mβ², estimated as 6.58(1.898) ≈ 12.49. The standard deviation is thus estimated as 12.49≈3.03 days ill per year.

Maximum likelihood estimates are not in closed form but numerical values for them could be obtained by numerical maximization of the likelihood function. It is helpful to use the method of moments as a quick and simple method to get an idea of the values of the parameters.

The table includes the marginal negative binomial with m = 3 and p = .344 in the comparison.

According to AIC, the histogram with varying bin widths still wins, the Negative Binomial coming in second. According to BIC (and, equivalentlly, posterior probability), the Negative Binomial scores the best, by far. This Negative Binomial is unimodal with a mode of .115 at 3 days. Because it is unimodal, it perhaps does not capture the flavor of the original data, which is reflected better by the Poisson mixture.

6.1 Beta-Binomial model

Another compound model is the Beta-Binomial model.In this model, the conditional distribution of X given p is Binomial(n,p). The prior on p is Beta(α, β). The posterior distribution of p given x is Beta(α + x, β + n − x). It is as if there had been a first round of α + β trials, with α successes, followed by a second round of n trials, with x successes. Method of Moments estimates of the parameters of the prior can be relatively easilty obtained. So can the Bayes estimates.

6.2 Normal-Normal model

We have considered the Gamma-Poisson model and, briefly, another prominent compound model, the Beta-Binomial model, with a Beta prior on the Binomial success probability parameters. Still another compound model is the Normal-Normal model.

In the Normal-Normal model, X is distributed according to 𝒩(μ, σ²), the Gaussian distribution with mean μ and variance σ². The prior on μ can be taken to be 𝒩(μ0,σ02) , or perhaps a Gaussian with a different mean if there is some particular reason to do this.

The posterior distribution is again Normal. That is, the Normal family is the conjugate family for the Normal distribution. The marginal distribution is also Normal, with mean ℰ[X] = ℰ[ℰ[X|μ] = ℰ[μ] = μ₀. The variance of the marginal distribution is the mean of the conditional variance plus the variance of the conditional mean, 𝒱[X]=ℰ[σ2]+𝒱[ℰ[μ]=σ2+𝒱[μ]=σ2+σ02 . These two terms are the “components of variance”. The decomposition of the variance can be obtained also by doing the requisite algebra on the product of the prior and conditional. This model is similar to a Random Effects (Model II) model in ANOVA. The group effects (group mean minus overall mean) are considered as a sample from 𝒩(0,σ02) .