Jul-29-2019, 04:35 PM
Congratulations scidam your post is this weeks winner.
(Jul-23-2019, 11:10 PM)scidam Wrote: I don't understand this code too, but I made some comments, hope that helps:
def calcPrior(n, N): # this line produces a list [N/1, N/2, ... , 1], i.e. a part of harmonic sequence. # in Python3 would consist of floats, in Python2 it would consist of integers l = [N/i for i in range(1, N+1)] # this is sum of the harmonic sequence tot = sum(l) # tot is a number # each value in harmonic series was divided onto # the sum of all values; l variable is overridden l = [i/tot for i in l] # this line assumes that 0<=n<=1; e.g. if N = 10, n=0.5, then # we return l[int(10*0.5)] = l[int(5.0)] = l[5], i.e. fifth # element of the array l. int(....) just drops mantissa of a float number, e.g. # int(1.919) = 1; return l[int(n*N)] def calcL(s, initialF, finalF): N = 180 eps = 0.0005 # some magic numbers here mu = initialF # this string wouldn't be needed, if finalF was declared as mu. (this is for convenience, mu is shorter than finalF) sd = (mu*(1-mu))**0.5#/N # This is likely iterative method of finding a solution of some equation (system of two equations) for i in range(50): old_mu = mu mu = (1+s)*mu/(1+s*mu) sd = (1.0/N)*(mu*(1-mu)) + (((1+s)/(1+s*old_mu))**2)*sd # these lines could be rewritten, e.g. as # for i in range(50): # 50 is another magic number # mu, sd =(1+s)*mu/(1+s*mu), (1.0/N)*(mu*(1-mu)) + (((1+s)/(1+s*old_mu))**2)*sd # e.g. without introducing old_mu # As we found a solution of equations, we integrate probability density function of normal # distribution over interval (finalF-eps, finalF+eps) #I completely don't understand why we do so. integral of pdf is cdf, we could just use # cdf instead. eps seems to be small, so this integral is almost equal to #something like `norm(loc=mu, scale=sd).cdf(finalF)*2*eps` prob = quad(norm(loc=mu, scale=sd).pdf, finalF-eps, finalF+eps) # ok, integration is done, as a result we got some probability # Finally, we returns that probability being multiplied on [initialF*N] - index of normalized harmonic sequence return prob[0]*calcPrior(initialF, N) # to get this code worked some import should be included at the beginning of the script, # import numpy as np # from scipy.stats import norm, chi2 # from scipy.integrate import quad # function findS is not defined in this code., start variable is not defined too. p value of s=0.5 p = 1-chi2.cdf(2*(np.log10(calcL(0.5, start, 1))-np.log10(findS(0, start, 1))), 1)