ramess101 / rjmc_lj_ethane Goto Github PK

@ocmadin

5/17 Notes:

• Owen is going to review the RJMC_LJ_2CLJ.py code to understand the details of what is going on and to determine ways to improve the code.

• We concluded that we will use the transition matrix with a delayed acceptance RJMC approach. We will run this by Josh and John on Wednesday.

• Transition matrix will require a way to scale epsilon, sigma, for different bond lengths and also different quadrupoles

• If we do not use a transition matrix but just use a delayed acceptance, we were thinking that we would want the step size to revert to the original guess so that it can find the optimal in 100 MCMC steps, or so. We would store the tuned step size for after the delayed acceptance criteria is applied, so that we do not need a completely new burn-in period. However, we foresee some problems if we have a non-informative prior that it will take a while to find the optimal region.

• By using a transition matrix and a delayed acceptance, we should not need to use a different step size for the model swap stage

Tasks:

• Modify RJMC_tuned to allow for delayed acceptance

• Modify posterior so that accounts for different number of parameters

• Modify AUA so that bond length is actually a parameter

• Modify LJ so that it is actually 2LJ (i.e. eps = 2*eps, bond length=0). This would only really help if we do not use a transition matrix

• Use sigma+L/2 as the property that is scaled. This would also be better if no transition matrix so that sigma is more feasible.

I will post questions/comments about the code here:

        for j in range(n_params):
    
            # Get current value for parameter j
            params = current_params.copy() # This approach updates previous param values
            
            # Propose new values
            if j == 0: #If proposing a new model
                if not i%swap_freq:
                    mod_ran = np.random.random()
                    if mod_ran < 1./3: #Use new models with equal probability
                        proposed_model = 0
                    elif mod_ran < 2./3:
                        proposed_model = 1
                    elif mod_ran < 1:
                        proposed_model = 2
                    if proposed_model != current_model:
                        model_swap_attempts += 1
                        params[0] = proposed_model
                        params[1] *= Tmatrix_eps[current_model,proposed_model]
                        params[2] *= Tmatrix_sig[current_model,proposed_model]
            else:        
params[j] = rnorm(current_params[j], prop_sd[j])

I am confused as to how the j==0 corresponds to choosing a new model. Won't this just always attempt to choose a model?

@ocmadin

I just wanted to document the changes somewhere more permanent than a slack conversation.

Copied from slack:

I (ramess101) made a few changes:

1. All the data are now found in a TRC_data/ directory with a subdirectory that is the chemical formula (C2H6, C2H4, etc.)
2. The first line of each file is now a header with, "T(K) property (units) uncertainty in property (units)"
3. I have included the TDE uncertainties (even if we never use them I thought it could be helpful)
4. I have created a Tc.txt file, rather than containing the critical temperature in the .yaml file
   TDE is ThermoData Engine. It is the software that automatically processes data. It is how we are determining which data sets to include or exclude. I only report a single Tc value (with uncertainty) rather than all the literature Tc values. This Tc value is the TDE estimated Tc which considers all experimental data (not just Tc experimental data) and finds the most thermodynamically consistent Tc value. This will be essentially an average of the experiments for the compounds of interest.
5. Previously I used REFPROP for DeltaHv (note that DeltaHv is now capitalized for the file name) because there are hardly any DeltaHv data for ethane. However, I am now only extracting real experimental DeltaHv data (this will be fairly scarce). Since we aren't actually using DeltaHv I though this would not be a problem. It was mainly for the case that we need DeltaHv experimental data down the road.
6. I created a lit_forcefields/ directory that contains the .yaml files for each compound (C2H6, C2H4, etc.)
7. I have restructured the .yaml files so that they only contain literature values for forcefield_params (epsilon, sigma, Lbond, Qpole) and the corresponding percent_deviations (rhol, Pv, gamma, Tc)
8. The literature forcefield_params are now contained within a single array where each element pertains to a different force field
9. sigma and Lbond are now in Angstrom and Qpole is in DAngstrom (to make it easier to compare with literature, we need to now divide these values by 10 when feeding into our functions that except nm or Dnm)
10. I created a Mw.txt file in TRC_data/ because I removed M_w from the yaml files

• Improve model swap acceptance for the 'All' cases

In general, model swap acceptance ratios are lower for 'All' cases than they are for 'rhol+psat' cases. We should look at using better Q proposals, model mapping techniques, and swap frequencies to improve these ratios - they will probably be different for different compounds

• Run 10^7 step RJMC for all compounds in the 'rhol+Psat' case.

Currently completed:

• C2H6 (Bayes Factor = 2.8)
• N2 (Bayes Factor = 2.4)
• O2 (Bayes Factor = 8.0)

Todo:

• C2H2 (Good data coverage except for SurfTens, only coverage from 0.62Tc-0.7Tc)

• C2H4 (Good data coverage except for SurfTens, only coverage from 0.41Tc-0.65Tc)

• Br2 (Good data coverage except for SurfTens, only coverage from 0.47Tc-0.55Tc)

• Cl2 (Only Psat data available)

• C2F4 (Good data coverage except for SurfTens, No SurfTens data available)

• C2Cl4 (rhol and Psat data only from 0.4Tc-0.65Tc, no SurfTense data available)

• F2 (Good Psat data coverage, rhol only from 0.5Tc-0.6Tc, SurfTens only from 0.45Tc-0.55Tc)

• Run 10^7 step RJMC for selected compounds in the 'All' case

• Have completed simulations for C2H6, N2, but unsure if results are valid (100% model 0 sampling for N2, uncoverged distribution for C2H6)

• Will run for O2 in the near future

• Won't be possible to do this for Cl2,C2F4,C2Cl4

• Will be able to do with modified data ranges for C2H2,C2H4,Br2,F2

This is a relatively short term plan, after completing these tasks and looking at the distributions/Bayes factors, we can determine how to move forward. If there are things you think I am missing please let me know.

@ramess101

I am copying an email correspondence I had with @jchodera a while ago:

Me:

As I interpret this equation, p(theta_j | D)/p(theta_i | D) is just the posterior distribution ratio that we normally compute with Metropolis-Hastings MCMC. The P(M_j) / P(M_i) is the prior knowledge ratio for the two models. If I assume that the prior knowledge suggests that they are equally probable this ratio will just be 1. The P(M_i | M_j) / P(M_i | M_j) is the ratio of "jump from the current model" at "pre-specified probabilities". Currently, in my code every MC step has a 50:50 probability of a model swap, therefore this ratio should be 1. The final term (T_ji/T_ij) is the one I am not as sure about. In other texts I have seen this ratio expressed as a single term, namely, the determinant of the Jacobian matrix. Currently, I am using a value of 1, but maybe this is my problem. In my specific example, I assumed that the Jacobian would be 1 because the parameters remain unchanged when going from one model to the next (since in one model the second parameter is simply ignored). However, I have come across some work where they introduce an auxiliary variable into the single parameter model so that they both have 2 parameters. Unfortunately, it is not clear to me how I would do this for my specific example.

John:

Apologies if the notation here was confusing! I realize now this is much more difficult to understand than I had hoped.

Patrick and I just worked out a more detailed version of the acceptance criteria. See the attached photo.

The likelihood ratio and prior ratio give you your posterior ratio, as you describe.

If you are not changing your parameters or putting them through a functional transformation as you change models, then the parameter proposal ratio is unity, as you describe.

And if your model jump proposal probability is symmetric, the model jump proposal probability is unit, as you describe.

Me:

Essentially, with all the new terms being equal to 1, my acceptance criterion is still just the ratio of posterior distributions. However, I am not sure if the posterior for the single parameter model should have a single prior (just for epsilon) while the two parameter model should have two priors (for both epsilon and sigma). Although it would make sense to only use 1 prior for the single parameter model, this results in an acceptance ratio of 1 for every model swap move when going from 1-parameter to 2-parameter model. So currently I am including two priors even for the single parameter model. Although this does not seem appropriate, this does result in the 1-parameter model being favored over the 2-parameter model.

John:

Going between models of different dimension is trickier. Patrick can point to some good references here.

The simplest approach in going from an [epsilon] to [epsilon,sigma] space is to simply keep the sigma variable around and use a separable prior on epsilon and sigma:

p(epsilon | M1) = p(epsilon | M2)
p(sigma | M1) = p(sigma | M2)

where you would just have one model (M1) be independent of sigma and the other model (M2) depend on both epsilon and sigma. This is essentially how the implementation in pymc would work, with a model index integer selecting which likelihood model you use, but the priors being the same for all models.

Me:

Here are my questions:

1. How should I determine T_ji(theta_i | theta_j) when the only difference between the models is the inclusion of an addition parameter? Do I need to use some sort of auxiliary parameter so they both have 2 parameters?

John:

If you are keeping the same number of parameters in each model and not changing the parameters when proposing a model change, this remains unity.
If you are changing the parameters but keeping the same number, this is not necessarily unity. If it is a deterministic transformation, the Jacobian is relevant. If it is stochastic, the probability distribution is relevant.
If you are changing the number of parameters, this gets more complicated.

Me:

2. Is proposing a model swap every step with equal probability a good approach?

John:

It depends on how many models you have. Optimizing this can in principle get you higher acceptance rates, but remember that increasing the probability of proposing a model k will also mean that the model jump proposal probability will oppose transtiions to that model.

Me:

3. Does the posterior distribution for the single parameter model only include a prior for epsilon while the two parameter model has a prior for both epsilon and sigma?

John:

If you are keeping sigma around in both models as an auxiliary variable in the model where it isn't relevant, and the prior for both parameters is the same in both models, you need to include the prior contribution in both models. The auxiliary variable will integrate out to be irrelevant in the evidence for the model in which it is not used.

Me:

4. Is there a simple way to implement this in PyMC? All I have found online for RJMC with PyMC is a GitHub code that is not completely finished.

John:

You can see an example that Patrick has coded up for GB forcefield parameterization here:

https://github.com/choderalab/gbff/blob/master/model.py#L280-L321

The model choice is a categorial variable.
All parameters have the same prior and are sampled in all models.
Only some parameters actually matter in computing the likelihood for each model, so they are free to float (within the prior distribution) if they don't appear in the likelihood.

Me:

Let me know if I need to clarify my questions and issues more.

I have been looking for a good resource for RJMC. The original article written by Green et al. in 1995 (Ref [54] in the NSF proposal) is rather difficult to digest. I would be very interested if you have any other resources to walk me through the algorithm.

John:

Hopefully Patrick can point you in the right direction!