Moment2 (FWHM) Images Do Not Agree With Spectral Line FWHM at Positions Near Sigma Clip Limits

[jmangum]

While verifying moment2 (FWHM) image values by comparing spectra toward specific positions in an input spectral line cube with their corresponding FWHM values in the moment2 image, noted significant disagreement near the signal-to-noise clip limit (i.e. 3-sigma). Attached is an example where the FWHM of the line (HCN 1-0) appears to be about 75 km/s, but the FWHM value from the moment2 image is about 25 km/s. This difference is rather systematic along the boundaries of the FWHM image.

[keflavich]

the test we need to do is run the FWHM estimation (i.e., moment 2 calculation) on Gaussians cutting at higher and higher thresholds. In principle it should stay constant and close to the correct value until there are only a few pixels left.

[jmangum]

ok. Why is this necessary? Note that we are using the "old" sigma-clipping approach, so the masking that you developed previously, which looked to be producing reasonable results for moment0. Is there something special about the moment2 calculation, and I guess msubcube.linewidth_fwhm(), that works poorly on this calculation for individual pixels?

[keflavich]

Is there something special about the moment2 calculation, and I guess msubcube.linewidth_fwhm(), that works poorly on this calculation for individual pixels?

That's the question I want to answer

[jmangum]

ok. So is this specifically related to msubcube.linewidth_fwhm()?

[jmangum]

Is this issue associated with #9 ?

[keflavich]

OK, so, in short: moment 2 is a systematically biased indicator of linewidth.

Demonstration:

import numpy as np
xaxis = np.linspace(-50, 150, 100)

sigma = 10.
center = 50.
synth_data = np.exp(-(xaxis-center)**2/(sigma**2 * 2.))

# Add noise
stddev = 0.1
noise = np.random.randn(xaxis.size)*stddev
error = stddev*np.ones_like(synth_data)
data = noise+synth_data

# Examine some test thresholds
threshold = np.linspace(-0.1, 0.9)

# Define our FWHM function
def fwhm(x,y,sel=None):
    if sel is not None:
        x=x[sel]
        y=y[sel]
    m1 = (x*y).sum() / y.sum()
    m2 = ((x-m1)**2*y).sum() / y.sum()
    return np.sqrt(m2 * (8*np.log(2)))

# calculate FWHM as a function of threshold
fwhms = [fwhm(xaxis, data, data > th) for th in threshold]

# plot it
pl.plot(threshold, u.Quantity(fwhms))
pl.axhline(sigma * np.sqrt(8*np.log(2)), linestyle='--', color='k')
pl.ylabel("FWHM"); pl.xlabel("Data Threshold")

result, where black dashed line is the right answer:

This is significantly more biased than I thought it was. I swear I've used unbiased estimators of the FWHM before, though, so maybe I'm forgetting something here.

Sample data is this:

[jmangum]

I am afraid that I have never used a moment2 calculation to derive a FWHM for a biased or clipped measurement.

It looks to me like no threshold (threshold = 0) badly overestimates the actual FWHM. Not sure I understand this behaviour. Looking at your fwhm function, I don't follow why you threshold both the x and y axes? For the moment calculations aren't we just thresholding on intensity (y-axis)?

[keflavich]

In the absence of noise, the calculation is correct. But it seems to be extremely sensitive to noise and threshold.

legend is stddev:

re: fwhm function - we only threshold on the data, on the y-axis. x[sel] is strictly necessary, you have to have arrays of the same length.

[jmangum]

Would like to understand why noise and threshold affect FWHM via moment2, but do need to find a path forward for getting FWHM from spectral cubes. How can we possibly use a moment2-based FWHM when noise, let alone threshold, corrupts the result? Or can you just not threshold and get a reliable FWHM from a moment2 measurement of noisy data?

[keflavich]

there are some s/n tests to do, but my experience in general is that not masking is much worse in a moderate s/n regime. In high s/n, not masking should be fine.

There are some alternatives to explore. I'm going to do some reading. I've used moment2 estimates of width for a long time, but I think all of my use cases have been as inputs to fitting a Gaussian. That's my only explanation for how I failed to realize this bias. There must be literature on this. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2464285/ is an example, but my first pass over that suggests that fitting might be involved - which we want to avoid, because any fitting approach is numerically unstable.

[jmangum]

Would it make any sense to instead calculate a FWZI by completely bypassing a moment2 calculation and doing this calculation on a pixel-by-pixel basis on the input image cube? One could then estimate an equivalent Gaussian FWHM from that FWZI if desired.

[keflavich]

I would expect FWZI to be much less well-behaved - "ZI" is poorly defined. But calculating the FWHM by stepping down from a peak, e.g. with https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.peak_widths.html, may work.