How do I make discrete predictions for discrete data that appears to have a normal distribution?

Question

I have discrete data that comes from the distribution of a discrete random variable Y. The data appears to follow the normal distribution (ie if I make a kde plot it looks like a normal distribution). I've conducted Bayesian inference and assumed that the mean mu follows a normal prior. Thus, the posterior of mu also follows a normal distribution.

Now I'm interested in conducting predictions for future samples from the pdf of Y given the samples I've already seen, but I'm unsure how to proceed. It doesn't make sense to use a normal distribution for prediction because all samples can only take non-negative Interger values. To me it doesn't seem to make sense to use a Poisson because the variance doesn't equal the mean nor does it make sense to use a Binomial because the value that the sample can take is theoretically unbounded (but in practice it is often bounded by ~25).

Answer 1

Not a lot of info to go off from, but a few things to note:

• OP mentions the marginal distribution is roughly symmetric and discrete. I've surmised this from their mention of a KDE plot.
• OP mentions "To me it doesn't seem to make sense to use a Poisson because the variance doesn't equal the mean". Now, since OP mentions examining the marginal distribution I will assume OP has computed the sample mean and variance of the outcome and noticed they are not equal. The size of the discrepancy is not mentioned, but let's assume for now they are so different that sampling variability could not reasonably explain the difference.
• OP doesn't mention any covariates so I'm inclined to think they are only looking at the marginal distribution.

There is no mention of what kind of predictions OP would like to make. Let's assume for now that they are interested mostly in the posterior and they could make predictions of many types from there.

Perhaps OP's data looks like

set.seed(0)
n <- 250
x <- rnorm(n)
y <- rpois(n, exp(5 + 0.03*x))
hist(y)

The data are roughly symmetric, discrete, and the mean is not equal to the variance. Can we model this with Bayesian statistics?

Sure, why not. Let's assume a Gamma prior on the rate parameter and use a Poisson likelihood. I'll assume OP has enough data so that the prior is negligible. Our model will be

$$ y_i \sim \mbox{Poisson}(\lambda_i) $$ $$ \lambda_i \sim \mbox{Gamma}(1, 1) $$

Because the prior is conjugate, we have a posterior we can easily sample from (its a negative binomial distribution, but let's pretend we don't know that for a moment). I will simulate 1000 datasets from the posterior predictive distribution and compare them to my y (again, pretending we don't know the x I've used).

Here are some histograms from the posterior predictive, real data in dark blue. The replications look fine, nothing too different

Additionally, we can look at the means and variances of the generated data. While the data are more variable than the replicates, it could be argued that the posterior predictive generates data that looks like the real data.

Anyway, I'm not saying OP should use this model. I'm trying to show OP how you can use Bayesian statistics with a discrete likelihood even though the mean and variance are not equal. I'm under the impression OP is being too restrictive with their modelling assumptions.

Answer 2

This sounds to me like an over-dispersion problem.
1. For the Poisson, you can "prove" than the distribution approximates the Normal given the right values of $\lambda$ and $n$. 2. It's quite common (I think McCullough and Nelder said safer to assume as default) that count data IRL have variance bigger than mean, and there is a whole world of over-dispersion models to use. You can use a form of the Negative Binomial (for example) to model your data which respects the discrete non-negative property of your data, but doesn't restrict you to having variance equal to the mean. One way of thinking about the Negative Binomial is a Poisson with a the location parameter $\lambda$ assumed to follow a Gamma distribution. I think (more by accident than design) I've tended to assume a log-Normal distribution for $\lambda$ whenever I've done this kind of modelling, but the nice thing about the Poisson-Gamma in a Bayesian framework is you get a nice simple (well, analytically tractable) posterior distribution. That sounds perfect for your problem.

Answer 3

If they are integer-valued samples, then the probability of getting each sample value $V$ would logically be the total density of the Gaussian prior you've fitted over $(V-0.5; V+0.5)$.

Anything greater - you're falling into the next bucket, less - previous one, as you draw the samples.

Since this is Gauss, you can just calculate a difference between two cdf at the two endpoints of the interval even in, like, Excel or whatever.