Negative standard deviation from Gaussian fit

[ysBach]

It has already been pointed out in #1105.
I have the following 17 by 17 array called test, and I did a Gaussian 2D fitting and got negative x_stddev:

test = [[-54.33, 13.81, -34.55, 8.95, -143.71, -0.81, 59.25, -14.78, -204.9, -30.87, -124.39, 123.53, 70.81, -109.48, -106.77, 35.64, 18.29], [-126.19, -89.13, 63.13, 50.74, 61.83, 19.06, 65.7, 77.94, 117.14, 139.37, 52.57, 236.04, 100.56, 242.28, -180.62, 154.02, -8.03], [91.43, 96.45, -118.59, -174.58, -116.49, 80.11, -86.81, 14.62, 79.26, 7.56, 54.99, 260.13, -136.42, -20.77, -77.55, 174.52, 134.41], [33.88, 7.63, 43.54, 70.99, 69.87, 33.97, 273.75, 176.66, 201.94, 336.34, 340.54, 163.77, -156.22, 21.49, -148.41, 94.88, 42.55], [82.28, 177.67, 26.81, 17.66, 47.81, -31.18, 353.23, 589.11, 553.27, 242.35, 444.12, 186.02, 140.73, 75.2, -87.98, -18.23, 166.74], [113.09, -37.01, 134.23, 71.89, 107.88, 198.69, 273.88, 626.63, 551.8, 547.61, 580.35, 337.8, 139.8, 157.64, -1.67, -26.99, 37.35], [106.47, 31.97, 84.99, -125.79, 195.0, 493.65, 861.89, 908.31, 803.9, 781.01, 532.59, 404.67, 115.18, 111.11, 28.08, 122.05, -58.36], [183.62, 45.22, 40.89, 111.58, 425.81, 321.53, 545.09, 866.02, 784.78, 731.35, 609.01, 405.41, -19.65, 71.2, -140.5, 144.07, 25.24], [137.13, -86.95, 15.39, 180.14, 353.23, 699.01, 1033.8, 1014.49, 814.11, 647.68, 461.03, 249.76, 94.8, 41.17, -1.16, 183.76, 188.19], [35.39, 26.92, 198.53, -37.78, 638.93, 624.41, 816.04, 867.28, 697.0, 491.56, 378.21, -18.46, -65.76, 98.1, 12.41, -102.18, 119.05], [190.73, 125.82, 311.45, 369.34, 554.39, 454.37, 755.7, 736.61, 542.43, 188.24, 214.86, 217.91, 7.91, 27.46, -172.14, -82.36, -80.31], [-55.39, 80.18, 267.19, 274.2, 169.53, 327.04, 488.15, 437.53, 225.38, 220.94, 4.01, -92.07, 39.68, 57.22, 144.66, 100.06, 34.96], [130.47, -4.23, 46.3, 101.49, 115.01, 217.38, 249.83, 115.9, 87.36, 105.81, -47.86, -9.94, -82.28, 144.45, 83.44, 23.49, 183.9], [-110.38, -115.98, 245.46, 103.51, 255.43, 163.47, 56.52, 33.82, -33.26, -111.29, 88.08, 193.2, -100.68, 15.44, 86.32, -26.44, -194.1], [109.36, 96.01, -124.89, -16.4, 84.37, 114.87, -65.65, -58.52, -23.22, 42.61, 144.91, -209.84, 110.29, 66.37, -117.85, -147.73, -122.51], [10.94, 45.98, 118.12, -46.53, -72.14, -74.22, 21.22, 0.39, 86.03, 23.97, -45.42, 12.05, -168.61, 27.79, 61.81, 84.07, 28.79], [46.61, -104.11, 56.71, -90.85, -16.51, -66.45, -141.34, 0.96, 58.08, 285.29, -61.41, -9.01, -323.38, 58.35, 80.14, -101.22, 145.65]]
g_init = Gaussian2D(x_mean=8, y_mean=8)
fitter = LevMarLSQFitter()
y, x = np.mgrid[:17,:17]
g_fit = fitter(g_init, x, y, test)
#  <Gaussian2D(
#amplitude=984.7694929790363, 
#x_mean=7.198391516587464, 
#y_mean=7.49720660088511,
#x_stddev=-1.9840185107597297, 
#y_stddev=3.1840618351417307, 
#theta=-18.144136589252792)>

So here I have two questions:

• Can I just take absolute value of it?
• Why do theta has so large absolute value, not constrained in (-pi, +pi)?
  • I didn't show it here, but if I constrain theta to some arbitrary value, such as (-tt1, +tt2), I get only -tt1 or +tt2, i.e., the boundary values. It seems like the fitter could not find the local minimum well, unless unconstrained. Is there any reason for this?

[pllim]

Possible fix in #6019. Here is the result using that PR:

>>> g_fit
<Gaussian2D(amplitude=984.7681025920282, x_mean=7.198391178949762, y_mean=7.497210019182043,
x_stddev=1.984021909876468, y_stddev=3.184065373194362, theta=16.41338035215769)>

[ysBach]

@pllim Thanks for your effort :)
But I am curious why the actual results of all values are slightly different from mine, e.g.,

amplitude = 984.7694929790363 in my result but
amplitude = 984.7681025920282 in the comment's result.

[pllim]

What machine do you have? I used RHEL6 64-bit with Python 3.5.3, Numpy 1.12, and SciPy 0.18.

[ysBach]

I am using Ubuntu 16.04.2 LTS (64-bit), Python 3.6.1, Numpy 1.12.1, SciPy 0.19.0.
Maybe SciPy gives the difference..?

[pllim]

I updated SciPy to 0.19.0 but still gets 984.7681025920282. However, does the difference matter? It is only 0.0001%. Although I do need to tweak "rtol" a bit to make the test pass:

result = 984.7681025920282
ans = 984.7694929790363
assert_allclose(result, ans, rtol=1.5e-6)

[pllim]

To answer your original questions:

1. Yes, it appears abs(x_stddev) is sufficient (without new constraints from my PR).
2. I don't know. It does seem weird, as Gaussian2D doc clearly mentioned theta should be (-90, 90) degrees, i.e., (-pi/2, pi/2) radians. But when I impose this as default constraint, the fitted values don't make sense (see below).

<Gaussian2D(amplitude=154.0972375365179,
x_mean=-11.952893612287438, y_mean=-7.442634264703583,
x_stddev=102.29185380715673, y_stddev=32.1563499422225,
theta=1.5707963267948966)>

[ysBach]

@pllim

However, does the difference matter?

I cannot think a case this small difference is a big matter, but I was just a bit curious why the discrepancy happens, though it's small. So I asked just in case you have a clear reason for it :)

I don't know. It does seem weird, as Gaussian2D doc clearly mentioned theta should be (-90, 90) degrees. But when I impose this is default constraint, the fitted values don't make sense (see below).

I encountered similarly peculiar results when I was using bounding option for fitting different data (It was Moffat2D). I still have no clue for such strange fitted values..

[pllim]

FWIW even if theta is out of bounds, the fitted result still looks decent. If only there is some documentation on how to convert crazy looking theta back to some reasonable value; e.g.:

-18.1441 rad (-1039.5804 deg) -> magical formula -> 0.7054 rad (40.4196 deg) ?

Here are the visualized results from my PR (your test data, fitted Gaussian, residuals):

[ysBach]

@pllim

If only there is some documentation on how to convert crazy looking theta back to some reasonable value

I have encountered that theta issue before, and I used the following for correcting it to [-pi,+pi]:

sig_x = np.abs(g_fit.x_stddev.value)
sig_y = np.abs(g_fit.y_stddev.value)
theta = g_fit.theta.value % (2*np.pi)
if theta > np.pi:
    theta = theta - np.pi
if sig_x < sig_y:
    theta = theta - np.pi/2

The second if was introduced because the x axis does not necessarily the semi-major axis in astropy fitting. Then regardless of sig_x > sig_y or sig_x < sig_y, the resulting theta from this is the angle from x-axis to the semi-major axis.

[pllim]

I could possibly add a new property to Gaussian2D that would return a theta that makes sense using your formula but I don't know of a good name to call it. theta_within_bounds? nice_theta? theta_that_is_not_crazy? @nden?