JumpCount

calculations of errorbars of ΔG predictions from molecular simulations

R

nA <- 3 # number of A to B transitions nB <- 3 # number of B to A transitions alpha <- 0.05 # value of significance level (0.05 for 95 % CI) temp <- 300 # temperature in Kelvins err_top <- 8.314*temp*log(qf(p=(1-alpha/2), df1=2*nB, df2=2*nA)) err_bottom <- -8.314*temp*log(qf(p=alpha/2, df1=2*nB, df2=2*nA)) se <- 8.314*temp*sqrt(trigamma(nA) + trigamma(nB)) print(err_top/1000) # height of top errorbar in kJ/mol print(err_bottom/1000) # height of bottom errorbar in kJ/mol print(se/1000) # standard error in kJ/mol

Python
import scipy as sp from scipy.stats import f nA = 3 # number of A to B transitions nB = 3 # number of B to A transitions alpha = 0.05 # value of significance level (0.05 for 95 % CI) temp = 300.0 # temperature in Kelvins err_top = 8.314*temp*sp.log(f.ppf(q=(1-alpha/2), dfn=2*nB, dfd=2*nA)) err_bottom = -8.314*temp*sp.log(f.ppf(q=alpha/2, dfn=2*nB, dfd=2*nA)) se = 8.314*temp*sp.sqrt(sp.special.polygamma(1, nA) + sp.special.polygamma(1, nB)) print(err_top/1000.0) # height of top errorbar in kJ/mol print(err_bottom/1000.0) # height of bottom errorbar in kJ/mol print(se/1000.0) # standard error in kJ/mol

Wolfram Mathematica
nA = 3; (* number of A to B transitions *) nB = 3; (* number of B to A transitions *) alpha = 0.05; (* value of significance level (0.05 for 95 % CI) *) temp = 300.0; (* temperature in Kelvins*) errTop = 8.314*temp*Log[Quantile[FRatioDistribution[2*nB, 2*nA], 1 - alpha/2]]; errBottom = -8.314*temp*Log[Quantile[FRatioDistribution[2*nB, 2*nA], alpha/2]]; se = 8.314*temp*Sqrt[PolyGamma[1, nA] + PolyGamma[1, nB]]; Print[errTop/1000.] (* height of top errorbar in kJ/mol *) Print[errBottom/1000.] (* height of bottom errorbar in kJ/mol *) Print[se/1000.] (* standard error in kJ/mol *)

Octave
nA = 3 # number of A to B transitions nB = 3 # number of B to A transitions alpha = 0.05 # value of significance level (0.05 for 95 % CI) temp = 300.0 # temperature in Kelvins errTop = 8.314*temp*log(finv(1 - alpha/2, 2*nB, 2*nA)) errBottom = -8.314*temp*log(finv(alpha/2, 2*nB, 2*nA)) se = 8.314*temp*sqrt(psi(1, nA) + psi(1, nB)) errTop/1000. # height of top errorbar in kJ/mol errBottom/1000. # height of bottom errorbar in kJ/mol se/1000. # standard error in kJ/mol

Julia
using Distributions using SpecialFunctions nA = 3 # number of A to B transitions nB = 3 # number of B to A transitions alpha = 0.05 # value of significance level (0.05 for 95 % CI) temp = 300.0 # temperature in Kelvins d = FDist(2*nB,2*nA) err_top = 8.314*temp*log(quantile(d, 1-alpha/2.0)) err_bottom = -8.314*temp*log(quantile(d, alpha/2.0)) se = 8.314*temp*sqrt(trigamma(nA) + trigamma(nB)) print(err_top/1000.0) # height of top errorbar in kJ/mol print(err_bottom/1000.0) # height of bottom errorbar in kJ/mol print(se/1000.0) # standard error in kJ/mol

<< BACK
Support: ELIXIR CZ