Solved – Chi squared test for goodness of fit

chi-squared-testcurve fittinggoodness of fithypothesis testingscipy

I have a signal coming from a measurement device (e.g., Volt vs., time) and perform a fit to the signal. I would like to do a goodness of fit test like the chi-squared test. However, I am not sure whether it is correct to apply it in my case.

In most online-tutorials (e.g. here http://www.stat.yale.edu/Courses/1997-98/101/chigf.htm), the first step of the chi-squared test is to construct the relative residuals as follows: $\sum (O_i – E_i)^2/E_i$ where $O_i$ would be the observed value and $E_i$ the expected value. This is most often done with data which are merely counts, e.g., counts in bins.

However, on the other hand I know that often in science the relative residuals are done in the following way: $\sum (O_i – E_i)^2/\sigma_i^2$ where $\sigma_i$ denotes the error of data point $i$.

I have two questions:

Is the first version only true when the noise on the data is Poisson, and thus $\sigma_i^2 = E_i$.
Since the examples I found are always for counts in some categories, is the chi-squared test even the right one to measure the goodness of a fit for, e.g., a voltage vs. time signal (or more mathematically: I don't draw multiple random variables according to some distribution and categorize them into bins)?

Edit:
Some more information about the fit and my data. For every timestep dt I have exactly one voltage value and I fit (using scipy curve_fit) a hyperbolic curve of the form
$f(x) = \frac{mx}{k+x}$ to it. I know the error of the data points. Here is a plot of the points used for the fit (red) and the fitted curve (blue). The way I would think to apply the test is calculating $\frac{1}{n_{\textrm{points}}-2}\sum (y_i – f(y_i))^2/\sigma_i^2$ (in the first denominator -2 as we have two parameters in the fitting function). And to my knowledge this should be around 1. As for the precise test, I have to say that I don't know which chi square distribution to use (the number of degrees of freedom) to make the hypothesis test.

Best Answer

Is the first version only true when the noise on the data is Poisson, and thus $σ^2_i=E_i$.

Not quite; for example it works for the multinomial (see Pearson 1900 ); $E_i$ is no longer the standard deviation but the dependence between cells in the multinomial exactly compensates for it; see also the test of independence.

Since the examples I found are always for counts in some categories, is the chi-squared test even the right one to measure the goodness of a fit for, e.g., a voltage vs. time signal

Under some very particular assumptions perhaps, including conditionally independent Gaussian response and known $\sigma_i$. I often see it applied where it clearly doesn't apply (e.g. where there's substantial observation error in the $x$'s and its a situation where an errors-in-variables model would apply, or where the supposedly-known $\sigma$ values are clearly inconsistent with the spread of the data around the fit).

To my recollection it doesn't actually apply to nonlinear models when estimating parameters (except approximately/asymptotically).

Related Solutions

Solved – Goodness of fit for 2D histograms

OK, I've extensively revised this answer. I think rather than binning your data and comparing counts in each bin, the suggestion I'd buried in my original answer of fitting a 2d kernel density estimate and comparing them is a much better idea. Even better, there is a function kde.test() in Tarn Duong's ks package for R that does this easy as pie.

Check the documentation for kde.test for more details and the arguments you can tweak. But basically it does pretty much exactly what you want. The p value it returns is the probability of generating the two sets of data you are comparing under the null hypothesis that they were being generated from the same distribution . So the higher the p-value, the better the fit between A and B. See my example below where this easily picks up that B1 and A are different, but that B2 and A are plausibly the same (which is how they were generated).

# generate some data that at least looks a bit similar
generate <- function(n, displ=1, perturb=1){
    BV <- rnorm(n, 1*displ, 0.4*perturb)
    UB <- -2*displ + BV + exp(rnorm(n,0,.3*perturb))
    data.frame(BV, UB)
}
set.seed(100)
A <- generate(300)
B1 <- generate(500, 0.9, 1.2)
B2 <- generate(100, 1, 1)
AandB <- rbind(A,B1, B2)
AandB$type <- rep(c("A", "B1", "B2"), c(300,500,100))

# plot
p <- ggplot(AandB, aes(x=BV, y=UB)) + facet_grid(~type) + 
    geom_smooth() +     scale_y_reverse() + theme_grey(9)
win.graph(7,3)
p +geom_point(size=.7)

enter image description here

> library(ks)
> kde.test(x1=as.matrix(A), x2=as.matrix(B1))$pvalue
[1] 2.213532e-05
> kde.test(x1=as.matrix(A), x2=as.matrix(B2))$pvalue
[1] 0.5769637

MY ORIGINAL ANSWER BELOW, IS KEPT ONLY BECAUSE THERE ARE NOW LINKS TO IT FROM ELSEWHERE WHICH WON'T MAKE SENSE

First, there may be other ways of going about this.

Justel et al have put forward a multivariate extension of the Kolmogorov-Smirnov test of goodness of fit which I think could be used in your case, to test how well each set of modelled data fit to the original. I couldn't find an implementation of this (eg in R) but maybe I didn't look hard enough.

Alternatively, there may be a way to do this by fitting a copula to both the original data and to each set of modelled data, and then comparing those models. There are implementations of this approach in R and other places but I'm not especially familiar with them so have not tried.

But to address your question directly, the approach you have taken is a reasonable one. Several points suggest themselves:

Unless your data set is bigger than it looks I think a 100 x 100 grid is too many bins. Intuitively, I can imagine you concluding the various sets of data are more dissimilar than they are just because of the precision of your bins means you have lots of bins with low numbers of points in them, even when the data density is high. However this in the end is a matter of judgement. I would certainly check your results with different approaches to binning.
Once you have done your binning and you have converted your data into (in effect) a contingency table of counts with two columns and number of rows equal to number of bins (10,000 in your case), you have a standard problem of comparing two columns of counts. Either a Chi square test or fitting some kind of Poisson model would work but as you say there is awkwardness because of the large number of zero counts. Either of those models are normally fit by minimising the sum of squares of the difference, weighted by the inverse of the expected number of counts; when this approaches zero it can cause problems.

Edit - the rest of this answer I now no longer believe to be an appropriate approach.

I think Fisher's exact test may not be useful or appropriate in this situation, where the marginal totals of rows in the cross-tab are not fixed. It will give a plausible answer but I find it difficult to reconcile its use with its original derivation from experimental design. I'm leaving the original answer here so the comments and follow up question make sense. Plus, there might still be a way of answering the OP's desired approach of binning the data and comparing the bins by some test based on the mean absolute or squared differences. Such an approach would still use the $n_g\times 2$ cross-tab referred to below and test for independence ie looking for a result where column A had the same proportions as column B.

I suspect that a solution to the above problem would be to use Fisher's exact test, applying it to your $n_g \times 2$ cross tab, where $n_g$ is the total number of bins. Although a complete calculation is not likely to be practical because of the number of rows in your table, you can get good estimates of the p-value using Monte Carlo simulation (the R implementation of Fisher's test gives this as an option for tables that are bigger than 2 x 2 and I suspect so do other packages). These p-values are the probability that the second set of data (from one of your models) have the same distribution through your bins as the original. Hence the higher the p-value, the better the fit.

I simulated some data to look a bit like yours and found that this approach was quite effective at identifying which of my "B" sets of data were generated from the same process as "A" and which were slightly different. Certainly more effective than the naked eye.

With this approach testing the independence of variables in a $n_g \times 2$ contingency table, it doesn't matter that the number of points in A is different to those in B (although note that it is a problem if you use just the sum of absolute differences or squared differences, as you originally propose). However, it does matter that each of your versions of B has a different number of points. Basically, larger B data sets will have a tendency to return lower p-values. I can think of several possible solutions to this problem. 1. You could reduce all of your B sets of data to the same size (the size of the smallest of your B sets), by taking a random sample of that size from all the B sets that are bigger than that size. 2. You could first fit a two dimensional kernel density estimate to each of your B sets, and then simulate data from that estimate that is equal sizes. 3. you could use some kind of simulation to work out the relationship of p-values to size and use that to "correct" the p-values you get from the above procedure so they are comparable. There are probably other alternatives too. Which one you do will depend on how the B data were generated, how different the sizes are, etc.

Hope that helps.

Solved – Appropriate goodness of fit measure for signal with unknown errors

If you estimate the curve parameters (i.e. $m$ and $k$) using least squares, then you are implicitly using the root-mean-squared error as the misfit metric (i.e. objective to be minimized).

In general for a regression problem you hypothesize a model of the form $$y=f_\theta(x)+\epsilon$$ where $(x,y)$ is the observed data, $f_\theta$ is a function depending on unknown parameters $\theta$, and $\epsilon$ is an unknown pointwise error (with expected value zero). Generally the parameters $\theta$ are estimated by minimizing some function of the residuals, $r=y-f_\theta(x)$.

In the case of least squares, this is $E(r)=\overline{(r^2)}$, so the "RMSE" corresponds to an estimated standard deviation of the error term (computed over the sampled residuals). If the errors $\epsilon$ are normally distributed, then least squares is a maximum-likelihood estimator of the parameters.

In your case, the errors appear to have some outliers (around $t=0.5$), which will inflate the error-scale estimate (and possibly bias the parameter estimate). You could mitigate this by using a robust estimate for the residual-dispersion (either as a post-processing, or as part of a robust regression scheme). In any case, while a single-number summary is convenient, it is always good to also inspect the residual distribution (both vs. time, and the bulk PDF).

For parameter uncertainty, a simple approach would be bootstrapping. This will also have the benefit of showing the impact of any outliers (as these will be less likely to be included in bootstrap samples).

Best Answer

Related Solutions

Solved – Goodness of fit for 2D histograms

Solved – Appropriate goodness of fit measure for signal with unknown errors

Related Question