Something wrong with Gaussian fitting?

[hopehhchen]

See the data and the codes @catherinezucker shared. HC thinks it's due to mishandling of units.

[hopehhchen]

Or the background subtraction?

[catherinezucker]

This also applies to the Plummer fitting. It keeps returning Rflat=0.000 and P=1.000. @hopehhchen do you have an ETA for when this bug will be fixed? I still need to determine the new widths for my paper, and then I need the bug-free version for the updated radfil tutorial I want to get to Bernie ASAP. Thanks!

[hopehhchen]

@catherinezucker, I will commit a bug-free version before the end of Thursday (8/10). Thanks for letting me know the problem.

[catherinezucker]

Sounds great, @hopehhchen . Thanks!

[hopehhchen]

@catherinezucker, your comment must be a lucky charm of some kind. I just found out that not assigning the bounds parameter in the Plummer fit would fix the problem (d41aef7f3a1c5f2b369af7a59ba51df557d5a8a4). Please try to break the current code with your unlikeliest filaments. I'll check on my part why the bounds parameter introduces unexpected behaviors. Thanks!

[catherinezucker]

Thanks, Hope! That definitely helped. Except now R_flat is going negative. And at least for the test filament, the answers are still offset by 25% from the original values. I think it might have something to do with the background subtraction. For reference, here are the original bounds and guesses I used with lmfit:

    #Plummer model
    if model=="Plummer":
        mod=Model(plummer)
        params = Parameters()
        params.add('N_0', value=np.max(self.mastery[include]),min=0.0)
        params.add('R_flat', value=np.median(self.distpc)/2.0,min=0.0)
        params.add('p', value=2.0,min=0.0)

    #Gaussian model
    if model=="Gaussian":
        mod=Model(gaussian)
        params=Parameters()
        params.add('amp',value=np.max(self.mastery[include]),min=0)
        params.add('wid',value=np.median(self.distpc)/3.0,min=0)

[catherinezucker]

Actually, no I don't think it's background. I think you can solve this by just replicating the guesses and at least the lower bounds I have.

[hopehhchen]

Thanks for the suggestion, @catherinezucker! I just set all bounds parameters to (0., np.inf), and it did solve the issue. Please test again and see if you still get numbers that don't make sense.

[catherinezucker]

Unfortunately it is failing for Musca now (R_flat goes back to 0.000).

[catherinezucker]

Can you see if you can recreate the tutorial fits for Musca with the same binning, samp_int, etc? And then I'll go through and redo all my filaments, quantify how different the best fits are between the two radfil versions, and then see if any additional changes need to be made. Does this sound OK @hopehhchen ?

[hopehhchen]

@catherinezucker, I think the problem is the weird behavior of boundary conditions in fitting a Plummer. I have updated the code by setting no bounds and made sure that R_flat is positively defined by taking the absolute value of R_flat after fitting. This is mathematically sound, because the target function to optimize should be symmetric for R_flat around 0.

Tests on B5, Filament 1, and Musca produce reasonable values. For Musca, I got Guassian and Plummer fits with each parameter within 1% of the value in the original RadFil_Tutorial.ipynb, using exactly the same setup parameters.

Please test the codes again to see if this actually solves the problem. Thanks!

[hopehhchen]

@catherinezucker, could you confirm that the fitting is good for your many filaments?

[catherinezucker]

Yeah, it looks really good so far, I think! I am trying to wrap this up today or tomorrow before I leave, so we will be able to close this issue soon. Will keep you updated.

On Fri, Aug 18, 2017 at 11:48 AM, Hope Chen notifications@github.com wrote:

@catherinezucker https://github.com/catherinezucker, could you confirm that the fitting is good for your many filaments?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/catherinezucker/radfil/issues/9#issuecomment-323390058, or mute the thread https://github.com/notifications/unsubscribe-auth/ARfXSUCjeueO2Gtmsyq_nyMyH-efo-Goks5sZbJRgaJpZM4Ov4tz .

[hopehhchen]

Cool! We can meet after our meeting with Bernie if that helps. Thanks, @catherinezucker!