I have a light curve, that is photon count rate versus time, of an astronomical object. These data are periodic since the source signal is periodic. I can fold the data with the period and obtain an average profile.

At this point, I want to check whether the single profiles are statistically different from the averaged one. Then, I was thinking to do a test between the average profile and the single profiles.

This test could be a 2sets KS test, or a Chi-squared test, or others. I don't know which one is better suited for my purpose.

Could you suggest/address me in this, please?

Here is a figure where the folded profile is shown in blue points, and the red line a tentative of smoothing the profile (not too relevant for this discussion)

enter image description here

Also, here a small collection of "single-profiles" compared with the model from the red line in the previous figure http: fig1 , fig2 , fgi3.

EDIT 2: here is a cumulative plot (very preliminary).

  • $\begingroup$ Dependence may be an issue. At the least, you'd want to separate out the single profile from the average of the rest (which could be done by a linear combination of the overall average and the single profile), which would remove one form of dependence (though there might still be some time-dependence issues). $\endgroup$
    – Glen_b
    Commented Apr 12, 2014 at 10:26
  • $\begingroup$ You mean, every time I want to test a single profile, I should compare it with the average one obtained without including the single I want to test, right? $\endgroup$
    – Py-ser
    Commented Apr 12, 2014 at 10:31
  • $\begingroup$ Yes, that's what I meant. Is your interest in identifying each one that's significantly different, or just that some are in fact different from others? $\endgroup$
    – Glen_b
    Commented Apr 12, 2014 at 10:58
  • $\begingroup$ Sorry, I can't understand your last question. You mean: "identify only ONE profile that is different..." ? Could you explain it better? $\endgroup$
    – Py-ser
    Commented Apr 12, 2014 at 11:18
  • $\begingroup$ I'd have thought that the use of the word 'each' conveyed that I didn't mean 'only one', but 'every individual one'. I was asking about the difference between concluding "at least one is different from at least one other" and "profiles 1, 13 and 28 each differ from the average of all others beside themselves". (Other possibilities exist, such as wanting to identify clusters, so that you might say 'profiles 1,13 and 28 fall into one group and all other profiles form the main group'.) $\endgroup$
    – Glen_b
    Commented Apr 12, 2014 at 13:46

3 Answers 3


At this point, I want to check whether the single profiles are statistically different from the averaged one. Then, I was thinking to do a test between the average profile and the single profiles.

I think which test you should use will be strongly dependent on what the profiles look like, and what you think an "average profile" should even be.

If the periodicity is exact, and the profiles appear to exist in different "classes", then maybe you can apply some sort of clustering algorithm (e.g. k-means) based on a distance metric (Euclidean, treating each profile as a vector -- but if high-dimensional, than maybe Manhattan norm is better?).

Another possibility is to compare some summary statistic applied to the profile, instead of the data directly. For example if the profiles are "similar" visually, but have different phases or magnitudes, then you might want to compute those values as an ordered pair representation that will give more robust results when you start comparing profiles.

Can you post some sample figures or data?

  • $\begingroup$ Thanks, I will edit my original post to include a figure. $\endgroup$
    – Py-ser
    Commented Apr 22, 2014 at 1:20
  • $\begingroup$ Thanks. Q: Is the number of samples (potentially) different for each profile? It looks like fig3 has fewer points than the red "average". $\endgroup$
    – Hao Ye
    Commented Apr 22, 2014 at 18:07
  • $\begingroup$ Yes, they very often vary in number of points and/or error bars. I have just put 3 profiles to give you an idea of how they are different from the folded one. $\endgroup$
    – Py-ser
    Commented Apr 23, 2014 at 0:32
  • $\begingroup$ In that case, I would just caution that any distance metric you use for similarity/dissimilarity controls for the number of sample points. (For instance, use RMSE instead of Euclidean distance.) $\endgroup$
    – Hao Ye
    Commented Apr 25, 2014 at 18:54
  • $\begingroup$ sorry, what do you mean? Are you referring to any particular test? $\endgroup$
    – Py-ser
    Commented Apr 28, 2014 at 8:02

If we ignore dependencies between time intervals, then the number of photons received during a specific time frame $[t_i, t_{i+1})$ can be modeled as a poisson distribution $\text{Pois}(\lambda_i)$. Hence, the maximum likelihood estimate of $\lambda_i$ is simply the average count during the time interval $[t_i, t_{i+1})$.

Let us denote the number of time intervals by $n$. A single profile with counts $x_1, \ldots, x_n$ is unlikely if the probability of generating it using the Poisson model with our estimated parameters is low. Since we are modeling the counts using independent Poisson variables, this probability is just the product of the point-mass-function of the different Poisson distributions. $\Pr[x_1, \ldots, x_n | \lambda_1, \ldots, \lambda_n] = \prod_{i=1}^n \frac{\lambda_i^{x_i}}{x_i!}e^{-x_i}$.

To avoid numerical issues, it is probably better to compute the log of this expression. $\log \Pr[x_1, \ldots, x_n | \lambda_1, \ldots, \lambda_n] = \sum_{i=1}^n x_i \log \lambda_i - \log (x_i !) - x_i $.

If this probability is very low then the single profile may be anomalous.

So, how can one choose a bound? The easiest way is by simulation:

  • Simulate many draws from the poisson model. i.e. take $X_i \sim \text{Pois}(\lambda_i)$ for $i=1, \ldots, n$.
  • For each draw $(X_1, \ldots, X_n)$ compute the log-likelihood as described above.

The above procedure yields scores $s_1, \ldots, s_r$, where $r$ is the number of simulations. These scores are samples from the distribution of log-likelihood scores under the null hypothesis. Thus an $\alpha=5\%$ level bound can be estimated by taking the 5-th percentile (from the bottom) of the set of scores $\{s_1, \ldots, s_r\}$.

As Glen_b notes, estimating the $\lambda$ values and then scoring using all the single profiles creates some nasty dependencies. It is better to estimate $\lambda_1, \ldots, \lambda_n$ using all the other profiles except the one currently being tested. However, if the total number of profiles is large, this doesn't really matter.

  • $\begingroup$ Thank you! Please, could you spend some few more words about "So, how can one choose a bound? The easiest way is to simulate many poisson draws, compute their probability and then take some low percentile of the resulting set of probabilities."? I am not sure about the meaning. $\endgroup$
    – Py-ser
    Commented Apr 30, 2014 at 3:11
  • $\begingroup$ Don't you think that neglecting the dependency of $\lambda_{i+1}$ from $\lambda_i$ could be a problem? If I follow correctly, you assume here that time intervals are all independent, but we observe here a rise and then a fall in the count of photons. Py-ser, I understand there is no model at all for the photon emissions, am I right? $\endgroup$
    – Zag
    Commented Apr 30, 2014 at 8:55
  • $\begingroup$ @Py-ser - I've edited for clarity. Can you understand it now? $\endgroup$ Commented Apr 30, 2014 at 14:48
  • $\begingroup$ @Zad - yes, this assumes independence and I've added a caveat in my answer. This may or may not make sense, but much more details and graphs about the problem domain are needed. $\endgroup$ Commented Apr 30, 2014 at 14:52
  • 1
    $\begingroup$ @Zag, no, there is no such a model in astronomy, the receiving of photons is considered poissonian. $\endgroup$
    – Py-ser
    Commented May 2, 2014 at 1:14

If I understand correctly (and please correct me if I am wrong), we have a stochastic process, that exhibits a well-defined "seasonality" (the orbital phase cycle -OPC). We count photons as the orbital phase evolves (hopefully in the same time intervals and of equal length). Then we chop the process into pieces, each piece of exactly the same length (one full OPC), we "stack" the pieces together, and calculate their average... and we want to test whether each piece (individual profile) is "statistically different" to this average you calculated.

The problem is, we don't know why the OP wants to do that. I am writing this to qualify what will come below -because it may prove useless to the OP.

When we treat a stochastic process like this, then, in order for such averaging calculations to be meaningful, it must be the case that the process exhibits stationarity and/or ergodicity (here in a seasonal/periodic sense, namely, that the process repeats itself stochastically in each OPC). Otherwise, this average profile is meaningless.

So if what I am describing is a correct translation of the case in question, then, testing single sets of process realizations (individual profiles) against the average, is not the way to go. One should test first for stationarity of the process, taking into account its periodic nature.

A simple way to do this is the following: Denote the instances of counts inside OPC $i$ by $C_{1i},C_{2i},...,C_{Ti}$ where $i=1,...,N$. Then take the same count-instance from each OPC and create a time series

$$\{C_{11}, C_{12},...,C_{1N}\},\;\; \text{...,} \;\;\{C_{T1}, C_{T2},...,C_{TN}\}$$

These time series have "by-passed" the periodic nature of the process. If these time series are stationary, then the whole process is seasonally stationary, repeating itself in each OPC.

This approach will also permit to detect perhaps other trends, or temporary phenomena, in your data. Even by visual inspection one could detect shifts or increased variability, and then, if these phenomena happen for the same bunch of $i$'s inside each time series (i.e. for the same OPC's), we will have found a period in the evolution of the process where something has happened, out of the "ordinary".

Assume now we found to our satisfaction that the process is OPC-stationary. Then creating the average profile is meaningful:

$$\{\bar C_{1\cdot}, \bar C_{2\cdot},..., \bar C_{T\cdot}\}$$

Since we have accepted stationarity, this is a series that includes "good approximations" to the expected value of the process in each count instance inside an OPC. Each instance inside an OPC may be governed by a different distribution. The average profile contains estimates of the mean values of these distributions.

Can we now go back to the original question? Can we compare the individual profiles against the average? What we would be comparing here? As said, the average profile is a series of mean values of underlying distributions. Each individual profile is a series containing realizations of these same distributions (per count instance). And we do not test whether a realization of a random variable is "statistically different" from its own expected value, do we?

Let alone the fact that by the look of it, each count instance inside an OPC is not characterized by the same distribution (mean values certainly look different). So a whole OPC cannot be viewed as data coming from the same distribution (it may come from the same family of course). Then, tests that were created to compare two sets of values each coming from a distribution, in order to test whether they come from the same distribution after all, are not relevant here, since each set of values separately does not come from one distribution to begin with. This makes questionable any such comparison between two individual profiles also.

Another answer provided a likelihood approach that is valid in that direction, through the use of a joint distribution of Poissons that permits differences in moments (of course it is valid only if the process is OPC stationary, and so averaging to obtain estimates for the $\lambda_i$s is valid).

But measuring somehow directly the "overall distance" between the average profile and each individual profile, if it tells us something, is how strong was the variability contained in this specific individual profile compared to the average profile. It would not be a "reject/not reject" statistical test of distributional similarity.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.