Gaussian1dModel.stddev is incorrectly used by astrofrog · Pull Request #1106 · astropy/astropy

astrofrog · 2013-05-17T21:48:11Z

While trying to fit some data, I noticed that the Gaussian1DModel is wrong:

    def eval(self, x, params):
        return params[0] * np.exp((-(1/(params[2]**2)) * 
                                                (x-params[1])**2))

it should be:

    def eval(self, x, params):
        return params[0] * np.exp((-(0.5/(params[2]**2)) * 
                                                (x-params[1])**2))

embray · 2013-05-17T21:23:39Z

+1 Good catch.

astrofrog · 2013-05-17T21:27:10Z

Shall I attach code to this issue?

embray · 2013-05-17T21:35:53Z

Can't hurt.

astrofrog · 2013-05-17T21:44:19Z

Weird, I'm getting a 'forbidden' error when trying to attach my fix (in my modeling/fix-gaussian-stddev branch) to this issue. Have you had that before?

astrofrog · 2013-05-17T21:48:33Z

Ah, turns out I had to update hub to the latest version.

embray · 2013-05-17T21:52:24Z

Yeah I dunno, I still just use @eteq's issue2pr.py script.

embray · 2013-05-17T22:03:59Z

👍 This is a bit less prone to embarrassing misreading I think 😉

astrofrog · 2013-05-17T22:19:43Z

@nden - should I just update the expected values for the tests?

nden · 2013-05-17T23:03:28Z

@astrofrog Yes, sure. Thank you.

astrofrog · 2013-05-18T09:58:06Z

@nden - I fixed all the failing tests except one which I don't understand. Do you think you could have a look to see if you can figure out what is going on?

astrofrog · 2013-05-18T10:00:46Z

I think the Gaussian2DModel has a similar issue, so looking into that

astrofrog · 2013-05-18T10:10:26Z

I tidied up and fixed the equation for the 2D Gaussian - however, I would appreciate if someone could double check the equations, since there are no tests for the numerical accuracy of the 2D Gaussian.

astrofrog · 2013-05-18T10:14:53Z

By the way, I think there's some optimizations that could be made for the calculations of the Gaussians, but I'll leave that to a separate issue (for example, avoiding recomputing the square of a parameter several times)

embray · 2013-05-18T14:37:28Z

Indeed, there are ways lots of these models could be optimized by breaking them up into separate terms that only need to be recomputed when one of the parameters is changed. But that is indeed something that needs to be tackled on a case by case basis probably, along with some profiling. Also almost all of the functions should eventually be redone to be written in Cython (and at some point I still want to explore using things like numexpr and/or Numba in Astropy but again that's another issue).

astrofrog · 2013-05-18T14:41:15Z

I actually think some of the functions are so generic that they should live outside astropy.modeling, for example 1 and 2-d Gaussians. As you can see from this PR, the 1D Gaussian equation appears quite a few times, so would benefit from being factored out.

embray · 2013-05-18T14:59:52Z

Indeed. That would be a necessary step anyways to providing multiple alternative implementations.

nden · 2013-05-18T16:03:11Z

@astrofrog The two failing tests are the only ones that use/test scipy. I can't reproduce the error locally. How do I find out the version of scipy running on Travis? Is it in travis.yml - scipy 0.12?

astrofrog · 2013-05-18T16:05:54Z

@nden - the version of scipy used on travis is given in the header of the test output:

platform linux2 -- Python 2.7.3 -- pytest-2.3.3
Running tests with Astropy version 0.3.dev1536.
Running tests in astropy.
Platform: Linux-2.6.32-042stab061.2-x86_64-with-Ubuntu-12.04-precise
Executable: /home/travis/virtualenv/python2.7/bin/python
Full Python Version: 
2.7.3 (default, Aug  1 2012, 05:14:39) 
[GCC 4.6.3]
Numpy: 1.7.1
Scipy: 0.12.0
Matplotlib: not available
h5py: 2.1.3

So it's using 0.12. I'm also using 0.12 locally.

adrn · 2013-05-19T21:50:56Z

There is another bug here with the order of param_names. This function should be:

def eval(self, x, y, p): a = 0.5 * ((np.cos(p[5]) / p[3]) ** 2 + (np.sin(p[5]) / p[4]) ** 2) b = 0.5 * (np.cos(p[5]) * np.sin(p[5]) * (1. / p[3] ** 2 - 1. / p[4] ** 2)) c = 0.5 * ((np.sin(p[5]) / p[3]) ** 2 + (np.cos(p[5]) / p[4]) ** 2) return p[0] * np.exp(-(a * (x - p[1]) ** 2 + b * (x - p[1]) * (y - p[2]) + c * (y - p[2]) ** 2))

The eval was getting stddev's and mean's wrong.

Thanks - this should be fixed

nden · 2013-05-20T17:57:23Z

@astrofrog @cdeil I've ran variations of the failing test. The result of the SLSQP fitter depends on the level of noise added to the signal. For example, if sigma in the test is defined as sigma = 8 * np.ones_like(self.xdata) the test passes. I would say it is not a very reliable algorithm. But it handles constraints internally (and actually the other similar test test_LSQ_SLSQP_cons passes. What do you think? Should the noise be decreased?

cdeil · 2013-05-20T20:13:23Z

Yes, change the test in whatever way necessary to achieve absolutely reproducible results in the unit tests.

Generally we don't want to have any test with borderline behaviour (e.g. different on 32 and 64 bit machines)
or test too many digits in fit results and then have some astropy.modeling unit tests fail on some user machines.

astrofrog · 2013-05-21T06:31:39Z

I'll try and change sigma as part of this PR to get the test passing.

astrofrog · 2013-05-21T07:07:30Z

If the tests pass, this is ready for a review.

cdeil · 2013-05-21T07:13:00Z

Please put stddev before fwhm so that Gaussian1DModel(1, 2, 3) will set the parameters corresponding to param_names, i.e. stddev=3. I think this is less confusing.

…ng for input paramters

astrofrog · 2013-05-21T08:40:27Z

@cdeil - should theta be the 6th argument for GaussianModel2D?

cdeil · 2013-05-21T08:55:59Z

IMO having the model param_names as the first model constructor arguments in the same order would be a good rule, so +1 to moving theta before the fwhm parameters, even if fwhm is related to sigma.

astrofrog · 2013-05-21T09:16:50Z

@cdeil - done

adrn · 2013-05-21T13:08:26Z

Maybe we can change the eval() in Gaussian2D to be a little clearer -- a bit verbose, but maybe easier to understand?

amp,x_mean,y_mean,x_stddev,y_stddev,theta = p

a = 0.5 * ((np.cos(theta) / x_stddev) ** 2 +
           (np.sin(theta) / y_stddev) ** 2)
b = 0.5 * (np.cos(theta) * np.sin(theta) *
           (1. / x_stddev ** 2 - 1. / y_stddev ** 2))
c = 0.5 * ((np.sin(theta) / x_stddev) ** 2 +
           (np.cos(theta) / y_stddev) ** 2)

return amp * np.exp(-(a * (x - x_mean) ** 2 +
                       b * (x - x_mean) * (y - y_mean) +
                       c * (y - y_mean) ** 2))

astrofrog · 2013-05-21T22:14:05Z

@nden - what do you think about @adrn's last suggestion?

nden · 2013-05-22T00:15:16Z

@astrofrog @adrn I timed it both ways fitting an 100x100 array and @astrofrog 's implementation is about 10% faster. May be we should sacrifice readability in favor of performance here.

nden · 2013-05-22T00:20:41Z

This test is written like a class. A function at the module level is probably sufficient, something like:

def test_gaussian2D(): m = models.Gaussian2DModel(2., 3., 4., x_stddev=1., y_stddev=5., theta=30.) assert m(3., 4.) == 2.

astrofrog · 2013-05-22T06:16:09Z

@nden - ok, sounds good. If you are happy with this, and the tests pass, then this is good to go.

Gaussian1dModel.stddev is incorrectly used This also fixed points 1 and 2 in #1077

ghost assigned nden May 17, 2013

embray reviewed May 17, 2013
View reviewed changes

astrofrog added 5 commits May 18, 2013 12:10

Fixed Gaussian1dModel stddev

31209a4

Removed redundant parentheses

dc426a1

Fixed some of the Gaussian1DModel tests

a435ae8

Fixed a couple of tests using Gaussian1DModel

2deca89

Tidied up and fixed equation for Gaussian2DModel

ec97c1e

Fixed expression for derivatives

b931633

cdeil mentioned this pull request May 19, 2013

models: bounds / tied issues #1064

Closed

adrn reviewed May 19, 2013
View reviewed changes

astrofrog added 2 commits May 21, 2013 08:53

Reduce noise to make tests of non-linear fitters more robust

6aef10d

Fixed order of parameters in Gaussian2DModel.eval

7a60197

cdeil reviewed May 21, 2013
View reviewed changes

Swap order of arguments for fwhm and stddev, and do more error checki…

164e77d

…ng for input paramters

Moved theta to be 6th argument for Gaussian2DModel

5cfb39a

nden reviewed May 22, 2013
View reviewed changes

Simplified Gaussian2DModel test

5de9a89

nden added a commit that referenced this pull request May 22, 2013

Merge pull request #1106 from astrofrog/modeling/fix-gaussian-stddev

24fac9d

Gaussian1dModel.stddev is incorrectly used This also fixed points 1 and 2 in #1077

nden merged commit 24fac9d into astropy:master May 22, 2013

astrofrog deleted the modeling/fix-gaussian-stddev branch July 5, 2016 15:50

Conversation

astrofrog commented May 17, 2013

Uh oh!

embray commented May 17, 2013

Uh oh!

astrofrog commented May 17, 2013

Uh oh!

embray commented May 17, 2013

Uh oh!

astrofrog commented May 17, 2013

Uh oh!

astrofrog commented May 17, 2013

Uh oh!

embray commented May 17, 2013

Uh oh!

embray May 17, 2013

Choose a reason for hiding this comment

Uh oh!

astrofrog commented May 17, 2013

Uh oh!

nden commented May 17, 2013

Uh oh!

astrofrog commented May 18, 2013

Uh oh!

astrofrog commented May 18, 2013

Uh oh!

astrofrog commented May 18, 2013

Uh oh!

astrofrog commented May 18, 2013

Uh oh!

embray commented May 18, 2013

Uh oh!

astrofrog commented May 18, 2013

Uh oh!

embray commented May 18, 2013

Uh oh!

nden commented May 18, 2013

Uh oh!

astrofrog commented May 18, 2013

Uh oh!

adrn May 19, 2013

Choose a reason for hiding this comment

Uh oh!

astrofrog May 21, 2013

Choose a reason for hiding this comment

Uh oh!

nden commented May 20, 2013

Uh oh!

cdeil commented May 20, 2013

Uh oh!

astrofrog commented May 21, 2013

Uh oh!

astrofrog commented May 21, 2013

Uh oh!

cdeil May 21, 2013

Choose a reason for hiding this comment

Uh oh!

astrofrog May 21, 2013

Choose a reason for hiding this comment

Uh oh!

astrofrog commented May 21, 2013

Uh oh!

cdeil commented May 21, 2013

Uh oh!

astrofrog commented May 21, 2013

Uh oh!

adrn commented May 21, 2013

Uh oh!

astrofrog commented May 21, 2013

Uh oh!

nden commented May 22, 2013

Uh oh!

nden May 22, 2013

Choose a reason for hiding this comment

Uh oh!

astrofrog May 22, 2013

Choose a reason for hiding this comment

Uh oh!

astrofrog commented May 22, 2013

Uh oh!