Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Repo info
Activity
    Alex Lubbock
    @alubbock
    I would recommend switching to Python 3.x wherever possible, since 2.7 isn't developed anymore and therefore PySB is dropping support for it too. If you have Python 3.6 working with Cython, then I wouldn't think it's worth your time trying to fix weave on Python 2.7.
    saubhagya-gatech
    @saubhagya-gatech
    Dear Forum members, I have a quick question about PyDream. I will attach the trace plots for the context.
    trace_plots.png
    PyDream did something weird at 10,000 steps as I see a huge spike in all the parameters. I am trying to understand what happened there.
    saubhagya-gatech
    @saubhagya-gatech
    my apologies spike at 12000 is a concatenating error. However, at every 2000 steps, I see bigs jumps. is that because of adapted gammas?
    Oscar Ortega
    @ortega2247
    Hi Saubhagya, are you plotting iterations vs. logp?
    saubhagya-gatech
    @saubhagya-gatech
    Hi Oscar, no, its iterations vs parameter value. I will guess at every 2000 steps, it is flagging chain outliers and putting them back close to other chains. But it seems like every now and then different chains are venturing into a region with high value of first parameter. I wonder if there is a parameter to tune this as that region might be acceptable region. I am not sure though.
    converged chains look something like this for 13 parameters. You can see clearly in first 3 parameters that the outlier chain is pulled back at every 2000 steps.
    Trace_Plots.png
    Oscar Ortega
    @ortega2247
    Hi Saubhagya, have you checked if your chains have converged? PyDREAM outputs a file with GelmanRubin_iteration in its name, this is the result of applying the Gelman-Rubin statistics to your chains. You can assume that your chains have converged if all the values are less than 1.2. Regarding the venturing of chains, there is parameter named snooker that you can pass to the run_dream function. This snooker parameter takes a value between 0-1 and correspond to the probability of doing a snooker update, the default value is set to 0.1. If you want to forego the snooker update you can set it to 0.
    snooker_update.png
    ^ This is the schematic of how a snooker update works
    For more information about the theory behind pydream you can also check: https://escholarship.org/uc/item/6j55p5kh
    saubhagya-gatech
    @saubhagya-gatech
    Thank you Oscar. Yes, my chains are converged with GR below 1.2 for all parameters. I was slightly skeptical because that huge shift is happening right at every 2000 steps. I wonder if it is happening due automatically adapting Gamma and Crossover probabilities? My number of iteration is 2000, so after every 2000 steps, Gammas and Crossover probabilities are adapted. I followed your Robertson problem example. After every 2000 steps, chains are set to last step of previous set of iterations using
    starts = [sampled_params[chain][-1, :] for chain in range(nchains)]
    saubhagya-gatech
    @saubhagya-gatech
    So, I did some checks and found that that chains which are "outlier" exploring in the region away from rest of the chains are pulled back into their initial position after each 2000 iterations. Is there any such check or feature which will not allow these outlier chains to continue from their last step but put them back to initial location of the "iteration set"?
    Oscar Ortega
    @ortega2247
    how are you restarting the pydream runs? PyDREAM should be loading the history of sampled parameters and shouldn't start from the initial positions, they should start from the latest parameter values sampled. Also, if you could send us a reproducible example of how you are getting these results, I could take a look at this
    saubhagya-gatech
    @saubhagya-gatech
    Hi Oscar, I guess this is the main confusion for me. How do we restart pydream? There are two arguments that deals with it -- "starts" and "restart". As I understand, either we can set "restart" true, the history file will be automatically loaded using model name from kwargs, OR we can give start position as the last step of previous run. So, if we have "restart" true and also have model name given in kwargs, we do not need to provide "start" argument. However, pydream still gives error: "Restart run specified but no start positions given."
    I looked at the sources code:
    if restart:
    if start == None:
    raise Exception('Restart run specified but no start positions given.')
    if 'model_name' not in kwargs:
    raise Exception('Restart run specified but no model name to load history and crossover value files from given.')
    we need either of the two conditions satisfied: either argument is pass to "start" or we have model name ink wargs. It seems the code above needs both conditions to be satisfied.
    I am using this example from the example directory to test restart. I have made iterations equal to 100 so that I can check automatic restart.
    saubhagya-gatech
    @saubhagya-gatech
    My apologies for bad formatting in chat. I will be careful next time. I am new to Gitter. I appreciate your help.

    In the example, work around is used by providing start position before each set of iterations as:

    if np.any(GR > 1.2):
    starts = [sampled_params[chain][-1, :] for chain in range(nchains)]
    while not converged:

    However, since it is provided before the while statement, the start position is never updated after first set of interations, 2000 iterations in my case, hence I see a huge shift after 4000, 6000, 8000 ....., iterations as they are set back to location at the end of 2000. For now I can move the start assignment after the while statement, but it would be nice to automate this through history file. I will be happy to provide more inputs if needed.

    Oscar Ortega
    @ortega2247

    Hi Saubhagya, I am sorry for the confusion. So, if you want to restart a pydream run, you must pass values for the start and model_name parameters. The start values should be the last parameter values sampled in your previous pydream run. The model_name should be the same model name you used in the previous pydream run. The model_name parameter is used to load the history of all parameters sampled, the crossover probabilities and gamma probabilities. Pydream can't use the history of parameters sampled to set the start argument because the history file doesn't have information about the chain that originated the sampled parameters, that's why it is necessary to pass the start argument manually.

    Regarding the Robertson example, you are right. There is an error in that example, the starts variable should be passed to the run_pydream function, it should read something like this:

            sampled_params, log_ps = run_dream(sampled_parameter_names, likelihood, niterations=niterations,
                                               nchains=nchains, multitry=False, gamma_levels=4, adapt_gamma=True, start=starts
                                               history_thin=1, model_name='robertson_nopysb_dreamzs_5chain', verbose=True, restart=True)
    saubhagya-gatech
    @saubhagya-gatech
    Hi Oscar, that helps a lot. Thank you very much for the clarification. I appreciate your assistance.
    Just another correction in the Robertson example, < starts = [sampled_params[chain][-1, :] for chain in range(nchains)] > should be after while statement. Currently it is before the while statement.
    if np.any(GR > 1.2):
           while not converged:
           starts = [sampled_params[chain][-1, :] for chain in range(nchains)]
    Oscar Ortega
    @ortega2247
    yes, thanks for pointing that out
    saubhagya-gatech
    @saubhagya-gatech
    Hi Oscar, another quick question. Does PyDream has the outlier detection functionality that is available in the original DREAM?
    Oscar Ortega
    @ortega2247
    PyDREAM implements the DREAM(zs) algorithm, which samples the parameter space based on past states. Thus, quoting Vrugt: "outlier chains do not need forceful treatment. Such chains can always sample their own past and with a periodic value of gamma=1 jump directly to the mode of the target."
    saubhagya-gatech
    @saubhagya-gatech
    Ah ok, that makes sense. Thanks a lot!
    Oscar Ortega
    @ortega2247
    Happy to help!
    saubhagya-gatech
    @saubhagya-gatech
    Good Morning,
    A questions regarding history file saving in PyDream. When we restart PyDream, prehistory files appended or over-written? I wonder if there is any difference between running a straight 50,000 steps run vs running 25,000 steps and then restarting the run from 25,000th step and same model name?
    Oscar Ortega
    @ortega2247
    yep, pydream should append the values to the history chains, crossover and gamma probabilities. If you want to save the chains sampled parameters and their correspondent logps you have to do it manually:
            for chain in range(len(sampled_params)):
                np.save('robertson_nopysb_dreamzs_5chain_sampled_params_chain_' + str(chain) + '_' + str(total_iterations),
                            sampled_params[chain])
                np.save('robertson_nopysb_dreamzs_5chain_logps_chain_' + str(chain) + '_' + str(total_iterations),
                            log_ps[chain])
    it should be the same to run 50,000 steps in one run or two sequential runs of 25,000 steps
    saubhagya-gatech
    @saubhagya-gatech
    Ah thanks. I have a question about convergence check. It is very likely I have a parameter space which is Bimodal, and the inferior mode is sampled less frequently by chains. Hence, GR factor goes below 1.2 and comes back up when a new chain ventures into inferior mode. I have picture of trace plot and GR factors attached. Does it mean that I just need to run chains for a much longer time, when all chains have sampled the inferior modes ad hence all have similar variances? Or I can assume that target distribution has been sampled the first time I get GR <1.2? There is a mention in the PyDream documentation.
    "also note that the presence of multiple modes in the posterior distributions will cause Gelman Rubin to fail, possibly falsely."
    Trace_plots.png
    GR Plot.png
    Oscar Ortega
    @ortega2247
    Yes, you are right. If the posterior distribution is bimodal the Gelman-Rubin diagnostics would fail. I don't know what to do in those cases, but I found this paper where the authors developed a modified Gelman-Rubin test that works for multimodal distributions: https://arxiv.org/pdf/1812.09384.pdf You could try to apply it to your results
    saubhagya-gatech
    @saubhagya-gatech
    Thanks. I will try that. I am having little bit f hard time understanding how gamma is adapted. Since, its default value is False, it appears that it is not a key functionality of DREAM(ZS). If I keep gamma adaptation disabled, I do not see any argument to pass gamma_0 value. I was exploring this option to improve the acceptance rate by reducing reducing the gamma_0 value. Or, is increasing Gamma level a better way to do this?
    Oscar Ortega
    @ortega2247
    You can set adapt_gamma=True if you want the gamma probabilities to be adapted. Otherwise, the gamma probabilities are defined: gamma_probabilities = [1/float(gamma_levels) for i in range(gamma_levels)], so if you want to include lower gamma probabilities you can increase the integer passed to the gamma_levels argument
    vw-liane
    @vw-liane

    Hi there! I hope you and yours are well.

    New person to Gitter & PySB here. I have some questions after going through the tutorial at https://pysb.readthedocs.io/en/latest/tutorial.html.

    1) When I run the cmd visualization command, I don't get errors, but neither do I get any visual output. Would you know what could be causing this or things I should try to debug? I can post the screenshots (how do I attach something?).

    2) Is there a spot for me to learn more about the "Model calibration" and "Modules" sections? They are empty at the moment.

    3) I saw multiple typos in the tutorial. Can I send you somehow the locations of the typos? Sorry about that; spelling is one of those small things that bugs me.

    Thank you for your time and I look forward to hearing from you.

    Alex Lubbock
    @alubbock
    Hi @vw-liane,
    1) What operating system are you using? Unfortunately, I think Gitter doesn't a way to upload images directly. You'd have to upload it somewhere like imgur and paste the link here.
    2) For model calibration, PyDREAM has PySB integration (here is a publication on that). There's more detail on the PySB modules here. We'll update the tutorial with a few links in the next release. Thanks for bringing it to my attention.
    3) We'd welcome your suggestions for improvements and fixes to the documentation. You could either open an issue on Github describing the issues, or even better, if you're comfortable with Github you could make the edits to the tutorial directly and open a pull request.
    vw-liane
    @vw-liane

    Hi @alubbock ,
    Sorry for my delay; life is wild.
    1) I am on Windows 10 64 bit. For now I will work around the visualization aspect until I explore that deeper in my project. Thank you for the tip on using imgur!
    2) I'm checking out PyDREAM as well! Thank you for that insight.
    3) Yes, I think I could do the edit and pull request! That sounds like the best option.

    When I have more thoughts/ questions, I will let you know! Thank you for your help so far.

    Alex Lubbock
    @alubbock
    @vw-liane No problem; thanks for taking the time to put together a pull request. I've just merged it.
    saubhagya-gatech
    @saubhagya-gatech

    Dear forum members,

    I am a PyDream user. I have a quick question. Is there a way to define different form of priors for different parameters?
    Currently I am using uniform priors for all my parameters. I am using following line of code:

    parameters_to_sample = SampledParam(uniform, loc=lower_limits, scale=scale)
    sampled_parameter_names = [parameters_to_sample]

    Now, I want to give uniform priors to a few parameters and gaussians to the rest. What is the best way to go about it?

    Oscar Ortega
    @ortega2247

    Hi,

    You can use any of the scipy distributions (https://docs.scipy.org/doc/scipy/reference/stats.html) as a prior in pydream.

    To use a uniform and a gaussian distribution you can do something like this:

    from scipy.stats import norm, uniform
    
    par1 = SampledParam(uniform, loc=lower_limits, scale=scale)
    par2=SampledParam(norm,  loc=mean, scale=std)
    
    sampled_parameters = [par1, par2

    Hope this is helpful

    saubhagya-gatech
    @saubhagya-gatech
    oh ok. This makes so much sense. Essentially, we have to give a list of scion distributions. Can we give distributions other than uniform or normal? Like some fitted kernel?
    image.png
    something like this:
    Oscar Ortega
    @ortega2247
    Interesting. I have never done anything like that but in theory it should be doable. Check this example from StackOverFlow to create a new scipy distribution. You just have to make sure that your distribution has a rvs and interval methods