8

Let's say I want to compare samples taken from Poisson distributions that have different values of lambda (rate). If the samples were from a normal distribution, I could convert each observation into a z-score based on the mean and standard deviation of the distribution that it was taken from. Then I could compare the z-scores to a common distribution.

How would I do this for samples taken from a Poisson distribution?

The specific use case is that I have a whole bunch of samples generated by many processes that have different poisson rates. The question I am asking is "Is the class of processes Poisson, or are they generated in a significantly more structured way?" I do not have enough data from any single process to answer this question, but I have more than enough data when I pool across all the processes.

Here is a more in depth explanation: Place cells are neurons that fire action potentials preferentially when an animal is in a certain location.
I can calculate the rate that a given cell fires in a given location by dividing the number of spikes fired by the amount of time spent in that location. FR(x)
I have observations of how many spikes were actually fired by a particular cell during a particular time interval at a particular location.
If I had enough data, I could compare my distribution of spike counts at a particular location in a particular time interval to the poisson distribution.
However, this is impossible, because the time interval is different for every observation.

Thus, I would like a standardized statistic that I could translate each spike count observation into given the spike rate, time interval, and location corresponding to that observation. I could then pool these standardized statistics and compare them to a single standard distribution.

honi
  • 151
  • 1
  • 9
  • 2
    A z-score is, *by definition*, a data value that has been recentered and scaled by the sample mean and sample standard deviation. Doing anything else would not be a z-score. Focusing on what a z-score is and how to compute it might be beside the point and perhaps even detrimental to your analysis. What comparison do you want to do and why? Are you attempting to test the hypothesis that the underlying rates are equal, for instance? – whuber Mar 08 '15 at 18:12
  • Yes I am looking for a parallel to a z-score, not the z-score itself. I am trying to ask whether a group of samples are or are not poisson distributed. The difficulty is that the samples are taken from different distributions. – honi Mar 08 '15 at 18:13
  • 2
    A z-score, or some analog of it, won't help you decide whether the data have a Poisson distribution. For that you would conduct a goodness of fit test--and it won't involve anything like a z-score in this situation. One issue to clear up is the nature of your data: because you refer to *processes,* it is possible you have more than a bunch of independent counts. Perhaps you also have the individual times between events? – whuber Mar 08 '15 at 18:25
  • a z-score analog would help because i need to "normalize" all my data in order to conduct a goodness of fit test. – honi Mar 08 '15 at 18:26
  • the individual times between events is not relevant in this case. do you want me to explain the entire experiment for you? – honi Mar 08 '15 at 18:27
  • 1
    I am not arguing: I am trying to steer you away from doing something wrong or inferior to other available solutions. You do *not* need to normalize your data; in fact, since Poisson data are discrete counts, normalizing them will make it almost impossible to conduct a suitable test. The individual times can be relevant because they provide more information than the counts alone. If you have the times, they can be exploited to improve your procedures, but if you don't have them, that's still ok too. – whuber Mar 08 '15 at 18:29
  • yes i meant "normalize" them so that they can all be compared together to a single distribution, just as you would "normalize" your samples taken from a normal distribution into z-scores so that you can conduct a goodness of fit test against a single normal distribution with mean 0, std dev 1. – honi Mar 08 '15 at 18:31
  • this question is similar to mine: http://stats.stackexchange.com/questions/70767/non-parametric-z-score-for-comparison-of-outliers-extreme-values-across-set-of-n?rq=1 – honi Mar 08 '15 at 18:33
  • 1
    What you might like to know is that different distributional assumptions require qualitatively different goodness-of-fit tests. Using a z-score can sometimes help when the underlying distribution is not discrete (specifically, when it is a location-scale family, which the Poisson is not), but GoF tests of discrete data usually proceed in a completely different way. Do you want to test whether all your distributions are Poisson--but potentially with different rates--or do you want to test whether they are all Poisson *with the same rate*? – whuber Mar 08 '15 at 18:37
  • whether they are all poisson with different rates – honi Mar 08 '15 at 18:39
  • 1
    OK. I would like to suggest, then, that replacing your first paragraph (which currently discusses z-scores) with this more specific request for such a GoF test would increase your chances of getting good answers. – whuber Mar 08 '15 at 18:40
  • is this better? – honi Mar 08 '15 at 18:42

3 Answers3

6

Assuming that the typical number of spikes is large, I would suggest to use a Variance stabilizing transformation. For Poisson distribution, it goes as follows: let $X \sim \mathcal P(\lambda)$, and let $Y \sim 2 \sqrt X$. A first order approximation gives $$Y \simeq 2 \sqrt\lambda + {1\over \sqrt\lambda}(X-\lambda),$$ from which we get $\text{var}(Y) \simeq {1\over \lambda} \times \lambda = 1$. Moreover for $\lambda$ large $Y$ is approximately normal $Y \sim \mathcal N(2\sqrt\lambda,1)$.

You can use $Y - 2\sqrt\lambda = 2(\sqrt X - \sqrt\lambda)$ as $z$-score.

Elvis
  • 11,870
  • 36
  • 56
  • 1
    Since the OP wishes to test whether processes might be Poisson, it seems unlikely that they know the parameter $\lambda$! – whuber Mar 15 '15 at 20:57
  • It seems that I misunderstood the question... – Elvis Mar 15 '15 at 23:04
  • I do know lambda! I know the number of spikes that occur at each location and how much time was spent at that location. The question is are those spikes evenly distributed over trials or not? This is seems to be exactly what I was looking for. – honi Mar 17 '15 at 21:34
  • What do you mean by "typical number of spikes is large"? typical mu? – honi Mar 17 '15 at 21:35
  • Large... $\lambda$, yes. I just corrected the mixture between $\mu$ and $\lambda$ in my answer... – Elvis Mar 18 '15 at 01:00
  • I am very unsettled by your question. I don’t know much about neurology, but I assumed that the (excited) neurons would fire very regularly — and that they cannot fire twice too fast, because of some "refueling time". If this is true, the process is not Poisson. However for non excited neurons the process might be (roughly) Poisson... (afaik) – Elvis Mar 18 '15 at 01:05
  • 1
    many neurons do not fire regularly at all. these neurons that i am studying have even been shown to be "overdispersed", that they fire more irregularly than poisson. – honi Mar 18 '15 at 01:57
  • 1
    the "refueling" time as you call it, or refractory period as the field calls it, does not affect the poisson-ness of neural firing that much as most neurons fire at a rate that is much lower than the length of time in the refractory period. eg. normal rate of 20 Hz and refractory period of 2 ms. – honi Mar 18 '15 at 01:59
  • "refractory period" — I knew it! but I was unable to find this back... What is overdispersed? Is that the number of spikes fired in a certain amount of time, measured several times? I mean $n$ measures $X_1, \dots, X_n$ of number of spikes during time $t$; if the process is Poisson with rate $\lambda$, each $X_i \sim \mathcal P(\lambda t)$; and the empirical variance is significantly higher than the empirical mean, whereas if the model where correct they should be near to each other? Is that the kind of thing you want to test? – Elvis Mar 18 '15 at 09:24
  • If you want to test this, can you show us a few "typical" values? – Elvis Mar 18 '15 at 09:33
  • 2
    I can't comment on the OP's application, but I've called this metric [a Poisson z-score](http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2551472) as well. That was for an application more in line with control charting though. It is buried in a footnote, but I did some simulations and this works ok with fairly overdispersed data for $\lambda$ as low as 1 and 5, the distribution just gets spread out (e.g. the 2.5th and the 97.5th percentiles are at 3 instead of 2). – Andy W Mar 19 '15 at 16:10
3

Is the reason that you want to test for Poisson distribution, that you think the spikes are generated by a Poisson process (that is, the instantaneous probability of seeing a spike is constant over time)? If so, and you know the length between each spike, consider testing for a Poisson process instead. This amounts to testing that the gaps between spikes follow an exponential distribution, with a fixed rate parameter at each location. This allows you to aggregate all the different events at each location, which hopefully gives you a large enough sample to have a reasonably powerful test. (Of course, few processes are actually Poisson, so it might be more informative to look at how much the distribution deviates from exponential and in what way, but that's a separate question.)

If you don't have between-spike timings, you may have to do something more complicated. You could, for instance, find the maximum-likelihood Poisson rate for a location, and compare your real data to a simulation from the resulting model in various ways to see how different it looks. As far as I know there are no off-the-shelf tests for this case though.

Ben Kuhn
  • 5,373
  • 1
  • 16
  • 27
  • the time intervals over which the firing rate is constant are too small for me to have inter-spike intervals. – honi Mar 17 '15 at 21:38
1

This is similar to Ben's question, but are you just trying to test whether a multi-variate Poisson distribution fits this data? Your major problem here is that you're having to estimate the poisson mean for each cell from the data that you're trying to get a p-value/likelihood score/whatever for. Now strictly speaking this may not be a problem, but if you don't have enough data for any of these individually, you're not going to have enough data for them collectively without some additional restriction on your hypothesis (which could come in a variety of forms).

Here's an example of the model you're vaguely describing. Let $i$ be each cell and $Y_i$ be the count within each cell.

$$ Y_i\sim Poisson(\lambda_i)\\ \lambda_i=\mu+\alpha_i $$

This is a general Poisson regression where the number of parameters is equal to the number of observations. Unless the ratio of those two approaches 0, your estimates for the number of parameters are not statistically consistent and your likelihood ratio won't tell you very much.

Now the model Ben Describes is a little better if you have the time interval data: here let $Y_ij$ by the time between $(j-1)th$ and $j$th observations:

$$ Y_{ij}\sim Exp(\lambda_i)\\ \lambda_i=\mu+\alpha_i $$

This is a little bit better, and if the ratio is small enough you may be able to test this via something like a Q-Q Plot on the probabilities of each time (which again, credit to Ben, is roughly what he's describing). If the number of locations is roughly the log of the number of observations, you may be in business there, but that QQ-Plot is still going to have some real bias, especially toward making the big outliers look more normal.

Some real improvements would be if there is any kind of additional knowledge you could impose on the problem. Like for example if you could say that sensors in nearby locations should have similar rates of firing, then you could leverage that knowledge by imposing a linear model on lat/long or maybe a thin-plate spline if you're not sure how the spatial relationship should look. Or maybe if the sensors of a certain "type" and you could say that sensors of the same "type" should have similar rates, and you could again make it linear and do an Anova test, or make a random effects model or something like that.

Anyway, there are a lot of ways to go with this if you have some kind of additional knowledge, or if your # of locations is $O(\log n)$. HTH

Mike Nute
  • 1,309
  • 7
  • 10
  • the time intervals over which the firing rate is constant are too small for me to have inter-spike intervals. – honi Mar 17 '15 at 21:39