California-Planet-Search / radvel

General Toolkit for Modeling Radial Velocity Data
http://radvel.readthedocs.io
MIT License
58 stars 52 forks source link

jitter has same value as gamma in example and in execution #99

Closed bpholden closed 6 years ago

bpholden commented 6 years ago

I maybe doing something wrong, but the readthedocs example shows the same behavior.

If you follow the Advanced Usage example, you have the value for jitter and for gamma converging to the same value. In the example, you find:

mod = initialize_model() like = radvel.likelihood.RVLikelihood(mod, t, vel, errvel) like.params['gamma'].value = 0.1 like.params['jit'].value = 1.0 but when

print(like) jit 1 True gamma 1 True

I find the same behavior in code that is modeled on the Advanced Usage example. In the case below you see the best fitting values below:

jit 39.5228 True gamma 39.5228 True

bjfultn commented 6 years ago

I can't reproduce this issue in the K2-24 jupyter notebook example.

Can you confirm that you are running the latest version of RadVel (1.0.0)?

Do you see this issue when you run through the K2-24 notebook? Or only when you are trying to reproduce it for your own application?

bpholden commented 6 years ago

1) This issue appeared when I upgraded to 1.0.0

2) Does K2-24 have the exact same jitter and gamma value in your notebook? As that is what the example in the following webpage shows: http://radvel.readthedocs.io/en/latest/tutorials.html#k2-24-fitting-mcmc

res = optimize.minimize( post.neglogprob_array, # objective function is negative log likelihood post.get_vary_params(), # initial variable parameters method='Powell', # Nelder-Mead also works )

plot_results(like) # plot best fit model print(post) parameter value vary per1 20.8853 False tc1 2072.79 False secosw1 0.01 False sesinw1 0.01 False logk1 1.56037 True per2 42.363 False tc2 2082.63 False secosw2 0.01 False sesinw2 0.01 False logk2 1.80937 True dvdt -0.0364432 True curv -0.00182455 True jit 2.62376 True gamma 2.62376 True

bpholden commented 6 years ago

The previous example was the purely circular orbit in the documentation. The eccentric orbit has the exact same behavior, the jitter and gamma are the same value:

res = optimize.minimize( post.neglogprob_array, post.get_vary_params(), method='Nelder-Mead',)

plot_results(like) print post

parameter value vary per1 20.8853 False tc1 2072.79 False secosw1 0.0821103 True sesinw1 0.402699 True logk1 1.15535 True per2 42.363 False tc2 2082.63 False secosw2 -0.187327 True sesinw2 -0.71348 True logk2 2.25931 True dvdt -0.0189353 True curv -0.00144404 True jit 3.46874 True gamma 3.46874 True

bjfultn commented 6 years ago

OK, the example in the notebook is slightly different than the advanced usage tutorial. They should be the same so I'll fix this.

Regardless, this is a very sneaky bug...

TLDR: These lines from the example:

like.params['gamma'].value = 1.0
like.params['jit'].value = 1.0

should actually be:

like.params['gamma'] = radvel.Parameter(value=1.0)
like.params['jit'] = radvel.Parameter(value=1.0)

The explanation: The likelihood object gives default values of radvel.Parameter(value=np.nan) to gamma and jitter as it is initialized. The problem is that it assigns a single instance of radvel.Parameter(value=np.nan) to both jitter and gamma. Since our radvel.Parameter object is mutable then updating its value attribute actually updates the value attributes for everything that was set to be equal to that particular radvel.Parameter instance. So, since initially like.params['gamma'] == radvel.Parameter(value=np.nan) == like.params['jit'] then if like.params['gamma'].value is set to 1.0 then like.params['jit'].value == 1.0 must also be true.

I've fixed this in version 1.0.1 (along with the docs) which will be available shortly, but you should update your code to explicitly set gamma and jit to their own radvel.Parameter instances as well since this is the more "Pythonic" way to do it.

bjfultn commented 6 years ago

@sblunt @petigura this "bug" actually uncovered an unintended feature. I think that if, say, we wanted to force "w1 == w2" I believe you can do:

params['w1'] = radvel.Parameters(value=XX)
params['w2'] = params['w1']

and I believe that will force this constraint in the fitting.

bpholden commented 6 years ago

Thanks BJ, and I should have noticed that!

I would like to say that the 1.0.0 API is very nice, much clearer than the previous versions.

On Wed, Nov 8, 2017 at 9:48 AM BJ Fulton notifications@github.com wrote:

@sblunt https://github.com/sblunt @petigura https://github.com/petigura this "bug" actually uncovered an unintended feature. I think that if, say, we wanted to force "w1 == w2" I believe you can do:

params['w1'] = radvel.Parameters(value=XX) params['w2'] = params['w1']

and I believe that will force this constraint in the fitting.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/California-Planet-Search/radvel/issues/99#issuecomment-342898199, or mute the thread https://github.com/notifications/unsubscribe-auth/ABKTmtFUvfgUikbB_2fPq6Mx-5LutJykks5s0elRgaJpZM4QVqup .