[SciPy-User] one-sided gauss fit -- or: how to estimate backgound noise ?

[Zachary Pincus]

> I was thinking about 1D-fitting the histogram. And I need this to be
> fully automatic so that I can apply that to many subregions of many
> images.
> I have to think about your suggestions for a while.

Well, if you don't want to estimate the mean/std of the underlying  
data but just work based on the histogram, you could use a nonlinear  
optimizer to fit a mean/std-parameterized gaussian PDF to the 1D  
histogram (with some amount of the tails chopped off). Just make the  
loss function the RMS between the data points (histogram height) and  
the fit curve at the positions of each data point? That would probably  
work too? But I'm not a statistician.

Though, now that I think on it, I seem to recall that the EM algorithm  
was originally deployed to estimate parameters of gaussians with  
censored tails. (I think the problem was: how tall was the average  
Frenchman, given distribution of heights of French soldiers and the  
knowledge that there was a height minimum for army service?) I think  
you just estimate the mean/std from the censored data, then fill in  
the censored tails with samples from the fit distribution, and then re- 
estimate the mean/std from the new data, etc. I forget exactly how one  
does this (does it work on the histogram, or the underlying data,  
e.g.) but that's the general idea.

Zach


> Thanks,
> Sebastian.
>
>
> On Fri, Jul 16, 2010 at 1:04 AM, Christoph Deil
> <Deil.Christoph at googlemail.com> wrote:
>> Hi Sebastian,
>>
>> in astronomy a method called kappa-sigma-clipping is sometimes used
>> to estimate the background level by clipping away most of the signal:
>> http://idlastro.gsfc.nasa.gov/ftp/pro/math/meanclip.pro
>> I am not aware of a python implementation, but it's just a few  
>> lines of code.
>>
>> If you can identify the background level approximately by eye,
>> e.g. by plotting a histogram of your data, you should be able to
>> just fit the tail of the Gaussian that only contains background.
>>
>> Here is my attempt at doing such a fit using  
>> scipy.stats.rv_continous.fit(),
>> similar to but not exactly what you want:
>>
>> from scipy.stats import norm, halfnorm, uniform
>> signal = - uniform.rvs(0, 3, size=10000)
>> background = norm.rvs(size=10000)
>> data = hstack((signal, background))
>> hist(data, bins=30)
>> selection = data[data>0]
>> halfnorm.fit(selection)
>> x = linspace(-3, 3, 100)
>> y = selection.sum() * halfnorm.pdf(x)/3
>> plot(x,y)
>>
>> Good luck!
>> Christoph
>>
>> On Jul 15, 2010, at 10:39 PM, Sebastian Haase wrote:
>>
>>> Hi,
>>> In image analysis one is often faced with (often unknown) background
>>> levels (offset) + (Gaussian) background noise.
>>> The overall intensity histogram of the image is in fact often  
>>> Gaussian
>>> (Bell shaped), but depending on how many (foreground) objects are
>>> present the histogram shows a positive tail of some sort.
>>>
>>> So, I just got the idea if there was a function (i.e. mathematical
>>> algorithm) that would allow to fit only the left half of a Gaussian
>>> bell curve to data points !?
>>> This would have to be done in a way that the center, the variance  
>>> (or
>>> sigma) and the peak height are free fitting parameters.
>>>
>>> Any help or ideas are appreciated,
>>> thanks
>>> Sebastian Haase
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user