Click here to Skip to main content
12,301,965 members (52,804 online)
Click here to Skip to main content

Tagged as

Stats

36.8K views
399 downloads
15 bookmarked
Posted

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()

	

By viewing downloads associated with this article you agree to the Terms of Service and the article's licence.

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.

License

This article, along with any associated source code and files, is licensed under The BSD License

Share

About the Author

John D. Cook
Singular Value Consulting
United States 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...

| Advertise | Privacy | Terms of Use | Mobile
Web02 | 2.8.160530.1 | Last Updated 26 Jun 2014
Article Copyright 2010 by John D. Cook
Everything else Copyright © CodeProject, 1999-2016
Layout: fixed | fluid