# How to measure accuracy of prediction distribution?

I have been trying to come op with a way to test if my dataset follows a Poisson distribution. Normally one would use the DistributionFitTest function, but here is my problem. The way my data is collected, the mean of the distribution is not a constant, such that I have a dataset with entries "{$\lambda$,outcome}".

Ideally my data would look like,

dataperfect = {5# , RandomVariate[PoissonDistribution[5#]]} & /@ RandomReal[{0, 1}, 10^3];


Unfortunately I suspect there is a small deviation of my measured $\lambda$ and the true mean. So my data looks more like this,

datatest = {5#[[1]]+.5#[[2]], RandomVariate[PoissonDistribution[5#[[1]]+.5#[[2]]]]} & /@ RandomReal[{0, 1}, {10^3,2}];


My question, is there a way of determining how close my datatest to a 'pure' Poisson distribution? Many thanks in advance.

EDIT

Since there is some discussion on whether the global or local comparison to a Poisson distribution is necessary some elaboration seems to be in order.

I have a certain model, a black box, which provides simultaneously two values $\lambda'_i$ and $\lambda''_i$ for the measured outcome's $\gamma'_i$ and $\gamma''_i$ of one event $i$. The assumption is made that the $\lambda_i$'s are Poisson distributed.

Further, a reasonable assumption is that the $\lambda_i$'s are a small deviation off the 'true' distribution parameter and the goal is to asses and minimise this deviation, by tweaking the model.

Using NonlinearModelFit instead of GeneralizedLinearModelFit in the way described here seems to do the trick, giving just a bit more freedom in defining a model for the expected deviation.

• Your edit makes this a completely different question. And given that there is already an accepted answer, I suspect this current question won't get the attention it deserves. In short, I suggest you open a new question. – JimB May 20 '16 at 13:04
• If you decide to open a new question, I suggest describing a model about how the two types of observations are generated and then what parameters should be compared. It shouldn't start with a test statistic. And it might be better to be asked at Cross Validated or the mathematics site. – JimB May 20 '16 at 13:07
• @JimBaldwin Knowing what questions to ask is always half of the answer. I'll open a new question with your suggestions. And in any case thank you for helping. – user19218 May 20 '16 at 13:30
• If you can tolerate some additional comments...The outcome and predictor variables need to be clearly defined: outcomes are usually the counts and $\lambda$'s usually represent the Poisson means. Also, the uniform distribution of p-values when it does apply only applies to continuous random variables and not ever to discrete random variables such as a Poisson. That, in part, is why I suggest starting with a model that describes the deviation from a Poisson that you suspect might be in the data generation process. – JimB May 20 '16 at 13:45
• In your edit that mentions NonlinearModelFit there is a strong caveat that must be mentioned. That function is much more flexible in the fixed effects part of the model but assumes that a least squares solution is appropriate and not what is optimal for a Poisson distribution. While that is partly mitigated by using an appropriate Weight, that turns you farther away from being able to test for your original quest of "pure Poissonness". – JimB May 21 '16 at 18:12

One possibility is to use a generalized linear model for which Mathematica has the function GeneralizedLinearModelFit. For your potentially Poisson data the observed count $y_i$ given a predictor variable $x_i$ will have a Poisson distribution with mean $\lambda_i$ with $\lambda_i$ being a function of $x_i$.

$$y_i|x_i \sim Poisson(\lambda_i)$$

with maybe $\lambda_i = a + b x_i$ (an identity link) or $log{(\lambda_i)}=a+b x_i$ (a log link). One can have any number of predictor variables and there are other link functions that can be considered.

If we assume that the "x" values in datatest are fixed and known (i.e., ignoring that they are randomly selected), then from your construction of the Poisson counts we expect that the mean count is proportional to $x_i$. We would be fitting

$$\lambda_i = b x_i$$

and estimating $b$.

Here is some Mathematica code to do so:

glm = GeneralizedLinearModelFit[datatest, x, x, ExponentialFamily -> "Poisson",
DispersionEstimatorFunction -> "PearsonChiSquare"]


The IncludeConstantBasis -> False tells the function not to fit an intercept. And as pointed out in the comments by @user19218, the DispersionEstimatorFunction needs to be made explicitly to "PearsonChiSquare" as the Automatic default is to always give a value of 1.0.

Now to finally get closer to answering your original question. You can't test if your data follows a Poisson distribution (after accounting for the different values of $\lambda$) but you can certainly look for serious departures from a Poisson distribution. One aspect that can be examined is characterized by the dispersion parameter:

glm["EstimatedDispersion"]
(* 1.01864 *)


If the distributions are Poisson, then we expect the dispersion parameter to be around 1.0. For this dataset, the estimate is close to 1 so we have no evidence of a departure from a Poisson. If it is much larger than 1 (say 2 or more), then one needs to consider how the over-and-above Poisson variability arises. Possibly one can add a random effect associated with some other factor (time of day, blocks of measurements under similar conditions, etc.). But Mathematica does not yet offer (as far as I know) a generalized linear mixed model - although one could construct the necessary code relatively easily for simple models.

My bias is certainly against doing any formal test of Poissonnesss, Gaussianness, etc., as it is extremely unlikely that any dataset comes from the exact kind of distribution being tested and a large enough sample size will end in "rejection". But what you want to know (again, in my opinion) is if there are any departures from the assumptions that would affect your conclusions. That requires both statistical and subject matter knowledge.

So how close to 1 do you need to be? My statistician answer is "It depends" but if the overdisperion estimate is greater than 2, then you'll need to consider modifying the model maybe using a negative binomial or adding a random effect or accounting for excess zeros depending on the situation. To get a more definite answer I'd suggest asking that at Cross Validated.

• This was exactly the answer I was looking for, thank you! – user19218 May 18 '16 at 8:10
• What a nice answer! +1 – ciao May 20 '16 at 8:48
• @JimBaldwin could not agree more with ciao. Fantastic answer. :-) – ubpdqn May 20 '16 at 9:02
• @JimBaldwin I think I can finally finish my quest for 'pure Poissonness', it turns out you were right all along! At least with the small correction of using the option DispersionEstimatorFunction -> "PearsonChiSquare" with GeneralizedLinearModelFit such that glm["EstimatedDispersion"] is not by default 1.0 ... – user19218 May 22 '16 at 13:44
• Good. That just goes to show that I should learn to read the documentation better. I'll update the answer so that your fix isn't just in the comments. Thanks! – JimB May 22 '16 at 15:31

My question, is there a way of determining how close my datatest to a 'pure' Poisson distribution?

@JimBaldwin provided an answer for considering the whole dataset (as requested in the question.) Here is proposed a method to evaluate the closeness to Poisson distribution locally using goodness of fit tests.

Let us generate the data as shown in the question but with the perfect data being 10 times larger in order to illustrate the method. (Not needed for the method itself.)

rpars = RandomReal[{0, 1}, {10^3, 2}];

datatest = {5 #[[1]] + .5 #[[2]],
RandomVariate[PoissonDistribution[5 #[[1]] + .5 #[[2]]]]} & /@
rpars ;

dataperfect = {5 #[[1]],
RandomVariate[PoissonDistribution[5 #[[1]]]]} & /@
RandomReal[{0, 1}, {10 Length[rpars], 2}];


Here is how datatest looks like:

And here is how dataperfect looks like:

Suppose we have chosen an $x$ point, $x_0$ (see the theoretical set-up in Jim's answer). We can gather points from datatest that have x-coordinates close to $x_0$ and do a goodness of fit test using Poisson distribution with $\lambda = x_0 \pm i*h$ for some small $h > 0$ and $-k \leq i \leq k, i \in Z, k \in N$.

Here is an example for $x_0=4.5$:

{dp, r} = {4.5, 0.1};
data = Select[datatest, Abs[#[[1]] - dp] <= r &][[All, 2]];
pres = Table[{dp + h,
PearsonChiSquareTest[data, PoissonDistribution[dp + h],
"PValue"]}, {h, -0.6, 0.6, 0.1}];
ColumnForm[{Length[data],
TableForm[pres,
TableHeadings -> {Range[Length[pres]], {"x", "p-value"}}] /.
Max[pres[[All, 2]]] -> Style[Max[pres[[All, 2]]], Red]
}]


We see that the maximum p-value is not at $x_0=4.5$. (The length of the dataset given to the goodness of fit test is printed out since we have to consider the statistical power of the test.)

We can compare with the result of the same procedure applied at the same point with dataperfect:

Block[{data = Select[dataperfect, Abs[#[[1]] - dp] <= r &][[All, 2]],
pres}, pres =
Table[{dp + h,
PearsonChiSquareTest[data, PoissonDistribution[dp + h],
"PValue"]}, {h, -0.6, 0.6, 0.1}];
ColumnForm[{
Length[data],
TableForm[pres,
TableHeadings -> {Range[Length[pres]], {"x", "p-value"}}] /.
Max[pres[[All, 2]]] -> Style[Max[pres[[All, 2]]], Red]
}]
]


A couple of obvious points.

• With this method there are some assumptions about the deviations of the parameters of the underlying processes. It is assumed that small differences in $x$ would bring small differences in the corresponding $\lambda$.

• Other distributions can be used with the goodness of fit test.

• Do you have a reference for such a procedure? I'm afraid I can't see anything related to a statistical test going on. The dataperfect data is created from a known process but all that will be known is the data and not the exact process that generates the data. Losing the +.5 #[[2]] in dataperfect does not make it some sort of standard of comparison. The "x" values are assumed to be fixed and known for a generalized linear model. A random selection of the "x" values (with/without the + .5 #[[2]]) is only to obtain a dataset to plug into GeneralizedLinearModelFit. I'm not getting it. – JimB May 20 '16 at 3:52
• @JimBaldwin I am not saying we have to use dataperfect to perform the test/procedure. The procedure is to perform a goodness of fit test to points with close x-coordinates. The comparison with dataperfect was given just for illustration how that procedure performs on a known dataset that adheres to the expectations of OP. The way I see, it we use goodness of fit instead visually judging the reconstructed PDF at the $x_0$ with the expected PDF as shown here. – Anton Antonov May 20 '16 at 8:26