Re: regression with poisson models

[Craig Markwardt <craigmnet(at)cow.physics.wisc.edu>]

Hi Nick--

matzke@geog.ucsb.edu (Nick Matzke) writes:
> 
> I have been using the function REGRESS (IDL 5.4) to test correlations
> between images derived from two different sensors (simply using each
> pixel-pixel comparison as an observation).  This works fine, but both
> datasets appear to have poisson distributions rather than gaussian
> distributions.
... long description deleted ...

First, I wanted to note that the image fitting code you will find on
my web page is not particularly suited to your task.  As you note for
yourself, what you are really doing is a linear regression between two
datasets, and the fact that the data came from images is largely
incidental.

Nick, I think you are generally on the right track.  I believe you
will want to substitute the Poisson error in place of the usual
Gaussian error.  This is done all the time for counting statistics,
which I assume is likely to be your case too.

Now, closer to your specifics.  You are describing REGRESS in IDL 5.4,
which is still "new" in my book [i.e. I don't have regular access to
it yet :-].  It appears that major changes have occurred with all the
regression codes in IDL 5.4, and to be honest I do not think all the
changes are good.  The big change is a shift to the
MEASURE_ERROR-style keywords instead of weights.  The problem with
this is exactly what you found: you can't assign a *zero* weight
without crashing the program!

My mind actually boggles at this change.  I hope they correct course
somewhat and allow both error and weight options.  The MPFIT family of
routines does.

Workarounds, in order of difficulty:

 * use the WEIGHTS parameter instead of MEASURE_ERROR.  It is
   obsolete, but it still works.

 * remove the zero-valued entries from your data set.  They get a zero
   weight anyway, so removing them cannot hurt.
   wh = where(y GT 0)
   result = regress(x(wh), y(wh), measure_error=sqrt(y(wh)), ... )

 * use MPFITFUN with a linear function, for which WEIGHTS will always
   work, I promise.  However you don't get all those clever automatic
   statistics that you get with REGRESS.  :-)


Good luck,
Craig

-- 
--------------------------------------------------------------------------
Craig B. Markwardt, Ph.D.         EMAIL:    craigmnet@cow.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response
--------------------------------------------------------------------------