12,826,440 members (24,967 online)

Stats

40K views
15 bookmarked
Posted 3 Feb 2010

Finding Probability Distribution Parameters from Percentiles

, 26 Jun 2014 BSD
How to determine the parameters of a probability distribution given two percentile constraints
 ```from scipy import stats, optimize from math import log, pow from functools import partial def normal_parameters(x1, p1, x2, p2): "Find parameters for a normal random variable X so that P(X < x1) = p1 and P(X < x2) = p2." denom = stats.norm.ppf(p2) - stats.norm.ppf(p1) sigma = (x2 - x1) / denom mu = (x1*stats.norm.ppf(p2) - x2*stats.norm.ppf(p1)) / denom return (mu, sigma) def cauchy_parameters(x1, p1, x2, p2): "Find parameters for a Cauchy random variable X so that P(X < x1) = p1 and P(X < x2) = p2." denom = stats.cauchy.ppf(p2) - stats.cauchy.ppf(p1) sigma = (x2 - x1) / denom mu = (x1*stats.cauchy.ppf(p2) - x2*stats.cauchy.ppf(p1)) / denom return (mu, sigma) def weibull_parameters(x1, p1, x2, p2): "Find parameters for a Weibull random variable X so that P(X < x1) = p1 and P(X < x2) = p2." gamma = (log(-log(1-p2)) - log(-log(1-p1))) / (log(x2) - log(x1)) beta = x1 / pow(-log(1-p1), 1.0/gamma) return (gamma, beta) def gamma_parameters(x1, p1, x2, p2): "Find parameters for a gamma random variable X so that P(X < x1) = p1 and P(X < x2) = p2." # Standardize so that x1 < x2 and p1 < p2 if p1 > p2: (p1, p2) = (p2, p1) (x1, x2) = (x2, x1) # function to find roots of for gamma distribution parameters def objective(alpha): return stats.gamma.ppf(p2, alpha) / stats.gamma.ppf(p1, alpha) - x2/x1 # The objective function we're wanting to find a root of is decreasing. # We need to find an inverval over which is goes from positive to negative. left = right = 1.0 while objective(left) < 0.0: left /= 2 while objective(right) > 0.0: right *= 2 alpha = optimize.bisect(objective, left, right) beta = x1 / stats.gamma.ppf(p1, alpha) return (alpha, beta) def inverse_gamma_parameters(x1, p1, x2, p2): (shape, scale) = gamma_parameters(1/x1, 1-p1, 1/x2, 1-p2) return (shape, 1/scale) def beta_parameters(x1, p1, x2, p2): "Find parameters for a beta random variable X so that P(X < x1) = p1 and P(X < x2) = p2." def square(x): return x*x def objective(v): (a, b) = v temp = square( stats.beta.cdf(x1, a, b) - p1 ) temp += square( stats.beta.cdf(x2, a, b) - p2 ) return temp # arbitrary initial guess of (3, 3) for parameters xopt = optimize.fmin(objective, (3, 3)) return (xopt[0], xopt[1]) def test(): #-------------------------------------------------------------------------- # Normal distribution # Make up a problem. p1 = 0.2 p2 = 0.8 mu = -3.4 sigma = 5.6 x1 = stats.norm.ppf(p1, mu, sigma) x2 = stats.norm.ppf(p2, mu, sigma) print "Testing normal distribution" print "True parameters: ", mu, sigma # Now solve the problem to make sure you get the original parameters back. (mu, sigma) = normal_parameters(x1, p1, x2, p2) print "Calculated parameters: ", mu, sigma print #-------------------------------------------------------------------------- # Cauchy distribution # Make up a problem. p1 = 0.1 p2 = 0.75 mu = -3.4 sigma = 5.6 x1 = stats.cauchy.ppf(p1, mu, sigma) x2 = stats.cauchy.ppf(p2, mu, sigma) print "Testing Cauchy distribution" print "True parameters: ", mu, sigma # Now solve the problem to make sure you get the original parameters back. (mu, sigma) = cauchy_parameters(x1, p1, x2, p2) print "Calculated parameters: ", mu, sigma print #-------------------------------------------------------------------------- # Weibull distribution # Make up a problem. p1 = 0.2 p2 = 0.95 gamma = 2.1 beta = 3.7 x1 = stats.exponweib.ppf(p1, 1, gamma, scale=beta) x2 = stats.exponweib.ppf(p2, 1, gamma, scale=beta) print "Testing Weibull distribution" print "True parameters: ", gamma, beta # Now solve the problem to make sure you get the original parameters back. (gamma, beta) = weibull_parameters(x1, p1, x2, p2) print "Calculated parameters: ", gamma, beta print #-------------------------------------------------------------------------- # Gamma distribution # Make up a problem. p1 = 0.2 p2 = 0.87 alpha = 0.4 beta = 0.86 x1 = stats.gamma.ppf(p1, alpha, scale=beta) x2 = stats.gamma.ppf(p2, alpha, scale=beta) print "Testing gamma distribution" print "True parameters: ", alpha, beta # Now solve the problem to make sure you get the original parameters back. (alpha, beta) = gamma_parameters(x1, p1, x2, p2) print "Calculated parameters: ", alpha, beta print #-------------------------------------------------------------------------- # Inverse gamma distribution # Make up a problem. p1 = 0.01 p2 = 0.38 alpha = 3.4 beta = 5.6 x1 = stats.invgamma.ppf(p1, alpha, scale=beta) x2 = stats.invgamma.ppf(p2, alpha, scale=beta) print "Testing inverse gamma distribution" print "True parameters: ", alpha, beta # Now solve the problem to make sure you get the original parameters back. (alpha, beta) = inverse_gamma_parameters(x1, p1, x2, p2) print "Calculated parameters: ", alpha, beta print #-------------------------------------------------------------------------- # Beta distribution # Make up a problem. p1 = 0.1 p2 = 0.8 a = 4.7 b = 2.8 x1 = stats.beta.ppf(p1, a, b) x2 = stats.beta.ppf(p2, a, b) print "Testing beta distribution" print "True parameters: ", a, b # Now solve the problem to make sure you get the original parameters back. (a, b) = beta_parameters(x1, p1, x2, p2) print "Calculated parameters: ", a, b print if __name__ == "__main__": test() ```

If a file you wish to view isn't highlighted, and is a text file (not binary), please let us know and we'll add colourisation support for it.

Share

 Singular Value Consulting United States
I am an independent consultant in software development and applied mathematics. I help companies learn from their data to make better decisions.

Check out my blog or send me a note.

You may also be interested in...

 Pro Pro