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.

[jmangum]

I would think that ZI down to a threshold would be well-behaved, but a FWHM determined by stepping down from the peak would also give a usable linewidth.

[jmangum]

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.

Not a problem I think, but scipy.signal.peak_widths does not play nicely with NaN values.

[keflavich]

we can come up with an alternative implementation, it's a pretty simple algorithm.

[keflavich]

https://ui.adsabs.harvard.edu/#abs/2005ApJ...623..826R/abstract appendix B gives the systematic correction for moment 2:

[jmangum]

Nice find! This looks like a correction that could be applied on a per-pixel basis. I guess the question now is whether to implement this correction or to implement the scipy.signal.peak_widths algorithm? I have to admit that I kind of like the scipy.signal.peak_widths algorithm option as it gets away from assuming Gaussian profiles (which is a generally poor assumption for galaxies).

[keflavich]

I would recommend going with the per-pixel correction approach, actually. Even though the line profiles may not be intrinsically Gaussian, it's probably a better extrapolation than any other profile. We can try the peak_widths approach, but I suspect it will have much worse behavior in the presence of noise. Of course, the Right Thing is to do a systematic comparison and calibrate on simulated data.

[jmangum]

ok. Your call. I should post a feature request to develop a simulated data set for testing like what you suggest.

[keflavich]

Here's the correction implemented:

from astropy import units as u
import numpy as np
import pylab as pl

def correction_factor(P, xref):
    """
    What exactly is xref?  xmax is roughly where the truncation happens
    """
    from scipy.integrate import quad
    xmax = np.sqrt(2*np.log(P))
    # split up the four terms
    t1 = np.array([quad(lambda x: x**4 * np.exp(x**2 / 2), 0, xm)[0] for xm in xref])
    t2 = np.array([quad(lambda x: x**2 * np.exp(x**2 / 2), 0, xm)[0] for xm in xref])
    t3 = np.array([quad(lambda x: x**4 * np.exp(x**2 / 2), 0, xm)[0] for xm in xmax])
    t4 = np.array([quad(lambda x: x**2 * np.exp(x**2 / 2), 0, xm)[0] for xm in xmax])
    return t1/t2*(t3/t4)**-1

# 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)))

pl.clf()
for stddev in (0, 0.001, 0.01):
    import numpy as np
    xaxis = np.linspace(-50, 150, 200)

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

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

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

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

    # plot it
    L, = pl.plot(threshold, u.Quantity(fwhms), label=stddev)

    # plot with the correction factor
    # what's xref?
    pl.plot(threshold, u.Quantity(fwhms)*correction_factor(Tmax/Tclip, [2]*len(threshold))**0.5, linestyle=':', color=L.get_color())

    pl.axhline(sigma * np.sqrt(8*np.log(2)), linestyle='--', color='k')
    pl.ylabel("FWHM"); pl.xlabel("Data Threshold")
pl.legend(loc='best')

I'm not entirely sure why it doesn't seem to work that well close to zero threshold

[jmangum]

The Rosolowsky and Blitz reference was for 2-sigma clip, which I interpreted as the reason for the "1/2" factors in each of the exp terms of f(P), and also for xref. Or did I misinterpret that part of the description above? Otherwise I don't see how the the specific Tclip chosen is incorporated into the correction.

It appears to me that the correction is not very good for most thresholds. Rosolowsky and Blitz did not say how good the correction makes the resultant FWHM, but implied that it should be perfect.

[keflavich]

the 1/2's in the exponents come from the definition of the Gaussian. But I also interpreted the xref=2 as being 2-sigma. However, when I tried setting xref=n sigma, the results are completely wrong - off the scale bad, I won't bother with a plot.

The tclip comes in via the P value which goes into xmax.

I agree, it doesn't look like this correction works, which I think hints that I've implemented it incorrectly. I therefore summon @low-sky to see if he can correct my error.

(one could also use the empirical correction that is plotted above)

[jmangum]

As we do not seem to have a solution to this issue, it might be prudent to warn users that the moment2 derived using this library is not usable, or even to disable the moment2 calculation?