Higher order cumulants?

[mattsperanza]

Hi,

I'm not exactly sure the right place for this, but I'm thinking of implementing GaMD and would like to have a good idea if (and how) I'll be able to re-weight the data I get out of the MD simulation. I also want to make sure that I'm understanding the existing issues properly and if there are existing ways of improving the estimation. I really appreciate the availability (and readability) of these analysis scripts and example data since it helps me see how the math actually works in practice.

Based on this 2014 paper, it is concluded that it is not worth using higher order cumulants because of numerical stability. Figure 8e and 8f (Trp-cage) show that the error in the 2nd and 3rd cumulant estimations corresponds to locations where the anharmonicity is large (also near a local maxima). Then in figures S3e and S3f (alanine dipeptide), it shows that the errors also occur near maxima even in an excess of sampling.

Here are the questions I have about this:

• Does this mean that the distribution of dV values is fundamentally skewed (has higher order cumulants) near maxima, or is it just skewed because sampling those maxima is difficult?
• Is this usually not an issue?
• Would higher order moments even help if the underlying problem is truly a sampling issue?

I haven't seen anything about it in any of the more recent papers/reviews, but maybe I missed it! The main reason I'm asking is because it seems like the issues are not arising due to calculation of <dV^n>, but rather in the formula used to combine the moments into cumulants. If it were issues with <dV^n>, the Maclaurin series would have stability issues as well. I recently stumbled onto this wikipedia article, where it mentions that variance and higher order central moments have stability issues when calculated using moments.

Since the third cumulant is just the third central moment (according to this wiki article), I figured that I could try out replacing the formula in your code in terms of moments for the formula in terms of the third central moment (calculated with spicy.stats.skew()). This is the result from the 1D example given:

Maybe you have already tried this and decided the 2nd order results are good and locations where there exists anharmonicity should simply be ignored since they aren't sampled well enough. Or it is possible that this data is particularly well suited for the changes I made. As a check, I wanted to try to use higher order cumulants, but 4th order fails:

I believe because I'm required to use a combination of lower order central moments in a manner similar to how 3rd order cumulants were used prior. There could exist a stable algorithm for all cumulants, but I haven't looked too much. I'm just curious to know if this could help (in other cases, since obviously these graphs show that C2 and C3 agree pretty well?

Thanks for all the great work with GaMD its a really great idea and I'm excited to try it out!

[MiaoLab20]

Hi, Thanks for the great questions: * Does this mean that the distribution of dV values is fundamentally skewed (has higher order cumulants) near maxima, or is it just skewed because sampling those maxima is difficult? In the Improved Reweighting paper, we performed simulations using the previous aMD and aMD_DB methods. As you can see in Figs. 1 and 6, distributions of the boost potential dV values are indeed skewed. Note that these are overall distributions; we didn’t plot distributions for individual bins, like those near maxima. My guess is that near the energy barriers (maxima), the conformations are less sampled and fewer data points are populated in each bin. This leads to higher anharmonicity and larger errors. You can probably do more analysis to double check. * Is this usually not an issue? That was an issue for accurate reweighting of aMD simulations as discussed in the 2014 Improved Reweighting paper. To address the issue, GaMD was later developed to use a harmonic function to generate the boost potentials that exhibit near-Gaussian distribution, with three enhanced sampling principles applied (see https://pubs.acs.org/doi/10.1021/acs.jctc.5b00436). More recently, we have also tried to generate Gaussian boost potentials using Bayesian neural network in a new Deep Boosted MD (DBMD) method (https://pubs.acs.org/doi/10.1021/acs.jpclett.3c00926), which can further help with the reweighting. * Would higher order moments even help if the underlying problem is truly a sampling issue? If your boost potential follows a Gaussian distribution, cumulant expansion to the 2nd order (also known as “Gaussian approximation”) will be mathematically exact to approximate the reweighting factor. Higher order cumulants would be too noisy and unnecessary in the reweighting calculation. Hope this helps, Yinglong Yinglong Miao, Ph.D. Associate Professor Computational Medicine Program and Department of Pharmacology University of North Carolina - Chapel Hill Tel: 1-919-962-5696 http://miaolab.org<http://miaolab.org/> Editor-in-Chief npj Drug Discovery https://www.nature.com/npjdrugdiscov On May 10, 2025, at 12:20 PM, Matthew J. Speranza ***@***.***> wrote: [https://avatars.githubusercontent.com/u/99224182?s=20&v=4]mattsperanza created an issue (MiaoLab20/pyreweighting#7)<#7> Hi, I'm not exactly sure the right place for this, but I'm thinking of implementing GaMD and would like to have a good idea if (and how) I'll be able to re-weight the data I get out of the MD simulation. I also want to make sure that I'm understanding the existing issues properly and if there are existing ways of improving the estimation. I really appreciate the availability (and readability) of these analysis scripts and example data since it helps me see how the math actually works in practice. Based on this<https://pubs.acs.org/doi/10.1021/ct500090q> 2014 paper, it is concluded that it is not worth using higher order cumulants because of numerical stability. Figure 8e and 8f (Trp-cage) show that the error in the 2nd and 3rd cumulant estimations corresponds to locations where the anharmonicity is large (also near a local maxima). Then in figures S3e and S3f (alanine dipeptide), it shows that the errors also occur near maxima even in an excess of sampling. Here are the questions I have about this: * Does this mean that the distribution of dV values is fundamentally skewed (has higher order cumulants) near maxima, or is it just skewed because sampling those maxima is difficult? * Is this usually not an issue? * Would higher order moments even help if the underlying problem is truly a sampling issue? I haven't seen anything about it in any of the more recent papers/reviews, but maybe I missed it! The main reason I'm asking is because it seems like the issues are not arising due to calculation of <dV^n>, but rather in the formula used to combine the moments into cumulants. If it were issues with <dV^n>, the Maclaurin series would have stability issues as well. I recently stumbled onto this wikipedia article<https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance>, where it mentions that variance and higher order central moments have stability issues when calculated using moments. Since the third cumulant is just the third central moment (according to this<https://en.wikipedia.org/wiki/Cumulant> wiki article), I figured that I could try out replacing the formula in your code in terms of moments for the formula in terms of the third central moment (calculated with spicy.stats.skew()). This is the result from the 1D example given: C3_Differences.png (view on web)<https://github.com/user-attachments/assets/60b0f0ce-090b-4344-bd08-00cc72ab2083> Maybe you have already tried this and decided the 2nd order results are good and locations where there exists anharmonicity should simply be ignored since they aren't sampled well enough. Or it is possible that this data is particularly well suited for the changes I made. As a check, I wanted to try to use higher order cumulants, but 4th order fails: C4_combo.png (view on web)<https://github.com/user-attachments/assets/e2c68b2d-44fe-4bc3-9d59-8b3c8b61693b> I believe because I'm required to use a combination of lower order central moments in a manner similar to how 3rd order cumulants were used prior. There could exist a stable algorithm for all cumulants, but I haven't looked too much. I'm just curious to know if this could help (in other cases, since obviously these graphs show that C2 and C3 agree pretty well? Thanks for all the great work with GaMD its a really great idea and I'm excited to try it out! — Reply to this email directly, view it on GitHub<#7>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ARJGPFV4X4VEZX6QCFQU3PT25YRNDAVCNFSM6AAAAAB427ZWRSVHI2DSMVQWIX3LMV43ASLTON2WKOZTGA2TIMZSHA4DQMQ>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>