# Applying a variance-stabilizing transform to a fitted function (rather than data)

## Outline

I'm working with data corrupted by a mixed Poisson-Gaussian noise model (for example with images gathered in astronomy or electron microscopy), and have been using the generalized Anscombe data transformation described in the paper:

"Optimal inversion of the generalized Anscombe transformation for Poisson-Gaussian noise" http://www.cs.tut.fi/~foi/papers/OptGenAnscombeInverse-doublecolumn-preprint.pdf doi:10.1109/TIP.2012.2202675

The forward transform is:

$$f_{\sigma}\left ( z \right ) = \left\{\begin{matrix} 2\sqrt{z+\frac{3}{8}+\sigma^{2}}, & z > -\frac{3}{8}-\sigma^{2}\\ 0, & z \leq -\frac{3}{8}-\sigma^{2} \end{matrix}\right.$$

where $z$ is the data and $\sigma$ represents the level of Gaussian noise added to the data.

The inverse transform is more complex (as explained in the paper), but essentially I'm using tabulated values and interpolation to apply the inverse transform to my data.

## The Problem

My problem lies in applying the transform to the parameters of a function that describes the data.

What I'd like to do is fit a curve to a noisy dataset, and I would like to see if prior application of this transform is better than a direct fit to the data.

In my scenario the model I'd like to fit is a Gaussian curve of the form

$$f(x)=A \exp \left ( -\frac{(B-x)^{2}}{2C^{2}} \right ) + D$$

The forward transform of the noisy data is straightforward enough, and fitting the model to this data is also fine.

The problem arises in transforming the model back to the original data range. The parameters $A$ and $D$ are simple, because the inverse transform (as described in the paper) can be used to return them to the original range. Similarly, parameter $B$ is simple to deal with, as it won't change.

But parameter $C$ is problematic. As you can see from the images, scaling the curve with the transformation affects both its height and its width, and I'm unsure as to how to proceed in transforming the model back to get the parameters.

For example, in the 1st image below, the original value is $C = 10$, but in the model that I've fitted to the transformed data (2nd image), $C=13.2$.

My current method is to fit the curve and use the model to generate a "clean" dataset. I then transform this data back with the inverse transform, and then fit again.

This is what I've done to get the 3rd image, where the "Original curve" represents the data before it was corrupted by the noise. You can see that it's a pretty good match.

My question is: Is it possible to avoid this second fitting step and transform the model directly, and so obtain the correct value for the width, $C$?

I have to fit more complex models than this (think a sum of Gaussian peaks) and as the number of parameters increases, the time penalty of having to fit twice become problematic.

## Ideas

1. My initial thought was to try and use the full-width-half-maximum (FWHM) of the peak to see if that was a possible route, but I hit a dead-end.

2. An alternative might be to try and come up with a table for the inverse transform of $C$, and then interpolate the values accordingly.

• Could you please explain what the "Original curve" in the bottom plot represents?
– whuber
Aug 22, 2014 at 16:08
• Of course - the "original curve" is what the data looked like before it was corrupted by noise. I realise I could have made this clearer. Have updated the image Aug 22, 2014 at 16:09
• The phrasing of this question suggests that the model as fit to the transformed data should be a Gaussian, but that is not the case. This inconsistency makes rather a puzzle of figuring out what you're trying to do. If I cover up 95% of the question and jump from the beginning to the end, ignoring everything in between, it seems to be asking "how do I estimate the width of a Gaussian signal that has Poisson noise?" Is this a correct interpretation of what you are looking for?
– whuber
Sep 17, 2014 at 19:59
• @whuber Both curves (i.e. untransformed and transformed) have the same form $f(x)$, but with different values of $A, B, C, D$, was my understanding. Because the Anscomeb tranform is a data transform, rescaling the background level $D$ and the amplitude $A$ are straight forward, but for $C$ that is not the case. It's converting the width back after the reverse transformation that's the problem. Sep 18, 2014 at 7:13
• My understanding is that the OP wants to fit a nonlinear parametrized mean value function $f$. The last comment suggests that the OP has already tried a direct maximum-likelihood approach with the given error model, and wants to compare this with a transformation based approach. Here the same parametrized family of functions $f$ is fitted to the transformed data, and an inverse transformation is sought to get the fit back on the original scale. I would expect the direct approach to be statistically more efficient under the error model.
– NRH
Sep 18, 2014 at 8:14