I'm doing some particle physics analysis and was hoping someone out there could give me some insight on a Gaussian-Process fit I'm trying to use to extrapolate some data.
I have data with uncertainties that I'm feeding in to the scikit-learn GaussianProcess algorithm. I'm including the uncertanties via the "nugget" argument (my implementation matches a standard example here where my "corr" is squared exponential and the "nugget" values are set to (dy/y)**2). The main concern is this: I have low absolute uncertainty (but high fractional uncertainty) at the edges of the distribution and this is producing a predicted confidence interval much larger than I expect in this region (see plot below).
The reason the uncertainties behave this way is that i'm dealing with particle physics data which is a histogram of counts of particles observed with different feature (x) values. These counts follow a Poisson distribution and thus have an uncertainty(standard deviation) of sqrt(N). So the higher count regions of the distribution have higher absolute, but lower fractional uncertainty, and vice versa for the low count regions.
I understand, as I mentioned, that the "nugget" argument in this function should have values of (fractional uncertainty)**2 when working with a squared exponential kernel. So it makes sense that if the predicted uncertainty is based on a fractional uncertainty of the input that it could be large on the edges. But I don't understand completely how this plays out in the math, and the size of the predicted uncertainty is SO much larger than the data point uncertainties on the edges that it seems wrong to me.
Can anyone comment on what's going on here? Is this behaving as expected? If so, why? Any thoughts or references to further reading on the subject would be greatly appreciated!
I'll leave you with a couple important caveats:
1) there are several data points with zero counts in the edges of the distribution. This throws a kink in the fractional uncertainty for the "nugget" because (sqrt(0)/0)**2 is not a very happy value. I made an adjustment here of just setting the nugget value for these points to 1.0, which corresponds to the value you get if this is a count of 1. I believe this is a common approximation which does affect the question at hand, but I don't think it fundamentally changes the issue.
2) The data i'm working with is actually a 2d histogram (ie, one independent variable (lets say x), another (y) and the counts as the dependent variable (z)). The plot shown is a 1d slice of the 2d data and prediction (i.e. z vs x integrated over a small range of y). I don't think this really affects the question at hand but I thought i'd mention it.