def wakefield_bf(beta,se,p,w = 0.1): # calculate Wakefield bayes factor z = stats.norm.isf(p/2) r = w/(se**2 + w) bf = math.sqrt(1-r) * math.exp(z**2/2*r) return bf