Possible problem with pscale adjustment in mcmc()

[vandalt]

Hi,

I'm not sure if this is a bug or a problem with how I set my model up, but I faced a problem with the MCMC scale of some paramters. It seems to be related to the lines that set pscales in the MCMC code.

When using a big value (order of 10^6 or 10^4) for parameters like tp or gamma, the scale value is 0.1*val which, I think, can make the parameters go out of the prior bounds (or far from their expected value).

I tested if changing mcmcscale to 10^-5 manually for these parameters worked and it does. If this is the expected behaviour and that I did not make a mistake when setting up my model, I feel like it would be useful to document it somewhere as initializing the parameters with a 10% gaussian "ball" can span a lot of parameter values and surprise users. Personally, I typically use a scale around 1e-5 except in special cases. I saw that the default for values set to 0 is 1e-5. I feel like this would be a good default value in general: it might cause the chains to converge a bit slower, but it would not cause the initial values to be outside the prior (except for very narrow priors, but I guess a warning could be raised in this case).

Let me know what you think, I'm happy to contribute if you think that this requires a modification to the code :)

[bjfultn]

The default scale for tc is 0.1 days, but it looks like I forgot to include tp in that if statement. I don't think we've run into this problem before because its generally not good to fit in the tp (synth) bases since tp and w go undefined at e=0. However, there are some use cases where it is desired and/or necessary. Why don't you include tp in the if statement for tc with your PR for the other issue.

            if post.vector.vector[par][2] == 0:
                if names[i].startswith('per'):
                    pscale = np.abs(val * 1e-5*np.log10(val))
                elif names[i].startswith('logper'):
                    pscale = np.abs(1e-5 * val)
                elif names[i].startswith('tc'):
                    pscale = 0.1
                elif val == 0:
                    pscale = .00001
                else:
                    pscale = np.abs(0.10 * val)

[bjfultn]

Also, its always best to subtract the mean or median from the RVs so that they are roughly centered around zero. dvdt and curve fit better this way and some little things in the plotting code might not work well with large gamma terms.

[vandalt]

I don't think we've run into this problem before because its generally not good to fit in the tp (synth) bases since tp and w go undefined at e=0.

Agreed! I usually don't use it either (which is why I think I never ran into the issue before), but we're using RadVel in an assignment for a grad course and we had to use the synth basis. Thank you for the confirmation, I'll include tp in the if statement in my other PR.

Also, its always best to subtract the mean or median from the RVs so that they are roughly centered around zero.

Thank you, I'll make sure to subtract the mean or median before fitting from now on. I was giving gamma as an example here, but I suspect this could happen with other parameters as well, especially with the custom model functionality where a parameter with an arbitrary name could have the same problem as tp or tc (large value vs. narrow prior). Are there many parameters in the RV or GP models that require 0.1 * np.abs(value) or do you think it would be a viable option to use 0.1 or 1e-5 as the default scale for all parameters except per and logper (and maybe other RV parameters if required) ? Otherwise maybe documenting the default pscale behaviour somewhere could help make users aware of possible problems (and then they can simply set mcmcscale manually for their parameters ? Let me know if you think I should add one of these suggestions to the PR mentioned above.

Thanks!

EDIT: Added suggestion about documentation