add attributes for fitting uncertainties to models by cshanahan1 · Pull Request #10552 · astropy/astropy

cshanahan1 · 2020-07-10T13:37:36Z

Fitting uncertainties are computed for models fit with LevMarLSQFitter and LinearLSQFitter. Models will have model.cov_matrix and model.stds attributes. Each model will have a parameter.std attribute as well (e.g model.slope.std).

The model attributes are set with a Covariance and StandardDeviations class, respectively, which are defined within fitting. Covariance can be indexed by two parameter names (e.g model.cov_matrix['slope', 'intercept'] for that covariance), or two integer indices corresponding to the order in model.param_names. StandardDeviations can also be indexed by parameter name or index in model.param_names. The plain numpy arrays can be accessed from these as well. The standard deviations can also be accessed from the parameters directly. These classes both have pprint method to customize how many fit parameters are printed, as well as rounding.

This does not work with units yet, but I will add that in the future.

This is a different approach to 10375, so that can be closed when this is merged.

pep8speaks · 2020-07-10T13:37:45Z

Hello @cshanahan1 👋! It looks like you've made some changes in your pull request, so I've checked the code again for style.

There are no PEP8 style issues with this pull request - thanks! 🎉

Comment last updated at 2020-09-21 14:07:05 UTC

perrygreenfield · 2020-07-10T13:51:51Z

+                par = getattr(self, param_name)
+                par.std = param_stds.pop(0)
+                setattr(self, param_name, par)
+


Should the cov_matrix setter also set the stds?

The cov_matrix setter does set the individual parameter stds (while the model.stds is set with its own setter). Do you think this should be done differently?

perrygreenfield · 2020-07-10T15:41:56Z

+                param = self.param_names[i]
+                s = param + ' ' * (longest_name - len(param)) + ' |'
+                ret_str += s + ', '.join([str(np.round(x, round_val))
+                           for x in row[0:i+1]]) + '\n'


You might consider using a format string to make this simpler. Something like
fstring = '{{0:<{0}s}} | {{1}}'.format(str(longest_name))

outside the loop, and

s = fstring.format(param, repr(np.round(row[:i+1], round_val))[7: -2]))

inside the loop.

This could be split up on more lines to be more readable (make the repr string its own line. This is a bit more complicated since part of the formatting specification is itself a variable (length of the longest parameter name)

perrygreenfield · 2020-07-10T15:49:29Z

+        else:
+            raise TypeError('Covariance must be indexed by two parameter names, or integer indicies.')
+        return(self.cov_matrix[i1][i2])
+


I'm not sure but I wonder if this handles properly the case when some of the parameters are fixed. That is, are fixed parameters represented in the covariance matrix? I wouldn't think so, but I'm guessing the supplied parameter names exclude the fixed parameters. Is that so?

Fixed and tied parameters are excluded from the list of fitted parameters. For example, a gaussian with one fixed parameter will have a cov_matrix of shape (2, 2), instead of (3, 3).

@perrygreenfield that should be how it works, but since i'm currently using model.param_names, which contains all
parameters including fixed + tied parameters, it doesn't work in it's current state. I'll change that.

I think this does work as expected now.

perrygreenfield · 2020-07-10T15:50:18Z

+            param = self.param_names[i]
+            s = param + ' ' * (longest_name - len(param) + 1)
+            ret_str += s + str(np.round(std, round_val)) + '\n'
+        return(ret_str)


A similar simplification using format strings can be done here too.

perrygreenfield

Basically only one significant comment that looks like an error, one comment about why a nested definition is necessary, and one about why 2D models aren't included. The first one should be addressed. The latter two don't need to hold up merging though the 2D issue should either be explained, or added as an issue if it is a valid concern.

perrygreenfield · 2020-08-25T17:40:34Z

+        masked = False or hasattr(y, 'mask')
+        msa = model.model_set_axis
+        # check if invertable. if not, can't calc covariance. 
+        def _is_invertable(m):


Why is this defined here? I don't see any dependence on the enclosing method's context.

cshanahan1 · 2020-09-03T13:37:21Z

@perrygreenfield a few notes

I think computing fitting uncertainties should be optional, with the default being not to (to not mess up anyones existing code that might be slowed / bogged down by switching this on, though I haven't actually timed it to see what that is). I don't currently have it this way but I would imagine it would involve adding a keyword argument to each fitter, and for those where it's not implemented raise an error if someone tries to set the switch to True? What do you think?

PG - I agree it should be optional. No argument there. I suspect it could be made an part of the base Fitter class so that you don't have to add it to all. Those that can't do it won't use it (e.g., setting an attribute of the fitter and then using it if it can handle them)

I'm not sure the linear fitter will handle masked/fitting with outlier removal for 2D models correctly. Part of the computation involves the number of points used for fitting. For the masked/outlier removal case, this will differ from the length of input x. For this reason, I create a masked array always for x and if none of the values are masked the mask will reflect that. Then, using x.count instead of len(x) give the number of points actually used in the fitting. This is straightforward for the 1D case.

PG - can you point to the specific code you are concerned about? I should take a look at it.

However for 2D, instead of the number of points it seems to be the number of rows in the 2d x/y passed in? I don't understand the math behind this - I was using the LevMarLSQ fitter, where the numpy fitter actually returns the covariance matrix, as a reference for the 'correct' value and realized I had to do it this way instead of the total number of points in the input array. For example, if x & y are both 10x10 arrays, the computation requires 10 rather than 100 to return the correct result and match LevMarLSQ. I'm not sure how this should work with 2d arrays that have various values masked.

I'm not sure the repr/str is actually useful, especially for model sets. I found myself accessing the plain numpy array instead. I also don't know if I like that the Covariance class has a pprint method to customize how many parameters are displayed/to what precision. I've seen methods named this in python before, is it reserved for something that follows a specific convention? Should printing options be controlled another way? Another issue is that I need to figure out how to account for the size of the window/notebook cell in printing. In instances when you have many parameters printing out covariances it will run off the line and look strange.

PG - don't bother optimizing for the window size now. That can be left to the user. If it is too much of a hassle at this point, just dump it all and we can worry about providing options later.
Thanks for reviewing!

nden · 2020-09-13T14:54:27Z

@cshanahan1 Please add a change log entry.

nden · 2020-09-16T17:48:12Z

This is pretty much ready for merging. @eteq or @karllark Are you able to review it?
It would be good if this can be tested against known results. If anyone has test cases and time please give it a try.

karllark · 2020-09-16T19:40:00Z

I can have a look by the end of next week. Is that soon enough?

nden · 2020-09-16T20:24:50Z

@karllark When you have time, much appreciated! I may merge it which will make it easier to test, but you can still leave comments.

nden · 2020-09-22T01:41:14Z

🎉

saimn · 2020-09-23T20:35:14Z

                model_lhs = lhs[good]
                model_rhs = model_rhs[good][..., np.newaxis]
-
+                a = model_lhs


I ran into some conflicts with this PR and #10687, and I do not understand how this line can work if there are multiple models. a should be different for each model isn't it ? cc @nden @cshanahan1
Also uncertanties is used everywhere thought that should be uncertainties, right?

I also spotted invertable -> invertible.

ahhh whoops! i am way too reliant on my spellchecker. I will fix these typos and take a look at the issue involving multiple models.

Thank you! It would be good to figure that out so we can merge the other PR. I've only just seen this and am not clear enough on what it's doing to offer any immediate insight.

@cshanahan1 - I already opened a PR for the typos (#10764), and can add the invertable correction there (thanks @jehturner ).

astropy-bot Bot added the modeling label Jul 10, 2020

pllim added the Enhancement label Jul 10, 2020

pllim added this to the v4.2 milestone Jul 10, 2020

pllim requested review from nden and perrygreenfield July 10, 2020 13:40

perrygreenfield reviewed Jul 10, 2020

View reviewed changes

nden mentioned this pull request Jul 19, 2020

fix fitting test #10582

Closed

saimn mentioned this pull request Jul 28, 2020

Improving astropy.modeling.fitting by returning parameter errors and final statistics value #10607

Open

cshanahan1 force-pushed the fitting_uncertanties_class branch from f7c175e to 50a4c6e Compare July 30, 2020 13:49

perrygreenfield reviewed Aug 27, 2020

View reviewed changes

cshanahan1 force-pushed the fitting_uncertanties_class branch 2 times, most recently from 36f09f7 to 9f791f5 Compare September 3, 2020 04:55

nden mentioned this pull request Sep 7, 2020

Fitting uncertainties on model / parameters #10375

Closed

cshanahan1 force-pushed the fitting_uncertanties_class branch 3 times, most recently from 8a6f3f8 to f801e12 Compare September 11, 2020 20:54

cshanahan1 force-pushed the fitting_uncertanties_class branch 2 times, most recently from db1706c to 417f79f Compare September 15, 2020 04:47

fitting uncertanties

19ed935

cshanahan1 force-pushed the fitting_uncertanties_class branch from 417f79f to 19ed935 Compare September 21, 2020 04:03

cshanahan1 added 2 commits September 21, 2020 09:59

tests for indexing

0d8a11f

trying this again..

ab2c5db

perrygreenfield approved these changes Sep 21, 2020

View reviewed changes

perrygreenfield merged commit 4a596d2 into astropy:master Sep 21, 2020

saimn added a commit to saimn/astropy that referenced this pull request Sep 23, 2020

Add back line from astropy#10552 after rebase

2f3ffc0

saimn reviewed Sep 23, 2020

View reviewed changes

saimn mentioned this pull request Sep 23, 2020

Fix typos in modeling #10764

Merged

saimn added a commit to saimn/astropy that referenced this pull request Oct 5, 2020

Add back line from astropy#10552 after rebase

113f725

saimn added a commit to saimn/astropy that referenced this pull request Oct 26, 2020

Add back line from astropy#10552 after rebase

f9fddf8

namurphy pushed a commit to namurphy/astropy that referenced this pull request Feb 21, 2021

Add back line from astropy#10552 after rebase

e23e4eb

alexmaragko mentioned this pull request Jul 13, 2021

Estimate feature strengths and uncertainties PAHFIT/pahfit#141

Closed

karllark mentioned this pull request Sep 10, 2021

Documentation needed for fit parameter uncertainties #12160

Open

Conversation

cshanahan1 commented Jul 10, 2020

Uh oh!

pep8speaks commented Jul 10, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Comment last updated at 2020-09-21 14:07:05 UTC

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

perrygreenfield left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

cshanahan1 commented Sep 3, 2020 • edited by perrygreenfield Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

nden commented Sep 13, 2020

Uh oh!

nden commented Sep 16, 2020

Uh oh!

karllark commented Sep 16, 2020

Uh oh!

nden commented Sep 16, 2020

Uh oh!

nden commented Sep 22, 2020

Uh oh!

saimn Sep 23, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

8 participants

pep8speaks commented Jul 10, 2020 •

edited

Loading

cshanahan1 commented Sep 3, 2020 •

edited by perrygreenfield

Loading

saimn Sep 23, 2020 •

edited

Loading