Ticket #9080: trac#9080.patch

File trac#9080.patch, 10.6 KB (added by klee, 3 years ago)

  • sage/gsl/gsl_random.pxi 

    # HG changeset patch
# User Kwankyu Lee <ekwankyu@gmail.com>
# Date 1275105086 -32400
# Node ID 1636334cf6db2082da0aa02a9fb8afbb8158e6d5
# Parent  1451c00a8d44f78f01cb0dc22ebb7464cc6ad0e2
add F distribution

diff -r 1451c00a8d44 -r 1636334cf6db sage/gsl/gsl_random.pxi
    a b  
    111111  double gsl_ran_tdist ( gsl_rng * r,  double nu) 
    112112  double gsl_ran_tdist_pdf ( double x,  double nu) 
    113113   
     114  double gsl_ran_fdist ( gsl_rng * r,  double a,  double b) 
     115  double gsl_ran_fdist_pdf ( double x,  double a,  double b) 
     116 
    114117  double gsl_ran_laplace ( gsl_rng * r,  double a) 
    115118  double gsl_ran_laplace_pdf ( double x,  double a) 
    116119   
    … …  
    197200  double gsl_cdf_fdist_P ( double x,  double nu1,  double nu2) 
    198201  double gsl_cdf_fdist_Q ( double x,  double nu1,  double nu2) 
    199202   
     203  double gsl_cdf_fdist_Pinv(double x, double a, double b) 
     204  double gsl_cdf_fdist_Qinv(double x, double a, double b) 
     205 
    200206  double gsl_cdf_beta_P ( double x,  double a,  double b) 
    201207  double gsl_cdf_beta_Q ( double x,  double a,  double b) 
    202208 

  • sage/gsl/probability_distribution.pyx 

    diff -r 1451c00a8d44 -r 1636334cf6db sage/gsl/probability_distribution.pyx
    a b  
    66- ``RealDistribution``: various real-valued probability distributions. 
    77 
    88- ``SphericalDistribution``: uniformly distributed points on the 
    9   surface of an `$n-1$` sphere in `$n$` dimensional euclidean space. 
     9  surface of an `n-1` sphere in `n` dimensional euclidean space. 
    1010 
    1111- ``GeneralDiscreteDistribution``: user-defined discrete distributions. 
    1212 
    … …  
    1919- Carlo Hamalainen (2008-08): full doctest coverage, more documentation, 
    2020  GeneralDiscreteDistribution, misc fixes. 
    2121 
     22- Kwankyu Lee (2010-05-29): F-distribution support. 
     23 
    2224REFERENCES: 
    2325 
    2426    GNU gsl library, General discrete distributions 
    … …  
    2830    http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Distributions.html 
    2931""" 
    3032 
    31  
    3233############################################################################## 
    3334#         Copyright (C) 2004, 2005, 2006 Joshua Kantor <kantor.jm@gmail.com> 
    3435#  Distributed under the terms of the GNU General Public License (GPL) 
    … …  
    4647import integration 
    4748from sage.modules.free_module_element import vector 
    4849 
    49  
    5050#TODO: Add more distributions available in gsl 
    5151#available but not currently wrapped are exponential, laplace, cauchy, landau, gamma, 
    5252#gamma, beta logistic. 
    … …  
    5858    lognormal 
    5959    pareto 
    6060    t 
     61    F 
    6162    chisquared 
    6263    exppow 
    6364    weibull 
    … …  
    365366        sage: T.cum_distribution_function_inv(.5) 
    366367        2.0 
    367368  
    368     The student's t distribution has one parameter ``nu``:: 
     369    The t-distribution has one parameter ``nu``:: 
    369370 
    370371        sage: nu = 1 
    371372        sage: T = RealDistribution('t', nu) 
    … …  
    377378        0.75 
    378379        sage: T.cum_distribution_function_inv(.5) 
    379380        0.0 
     381 
     382    The F-distribution has two parameters ``nu1`` and ``nu2``::      
     383 
     384        sage: nu1 = 9; nu2 = 17                              
     385        sage: F = RealDistribution('F', [nu1,nu2]) 
     386        sage: F.get_random_element() # random 
     387        1.65211335491 
     388        sage: F.distribution_function(1) 
     389        0.669502550519 
     390        sage: F.cum_distribution_function(3.68) 
     391        0.98997177723 
     392        sage: F.cum_distribution_function_inv(0.99) 
     393        3.68224152405 
    380394  
    381     The chi squared distribution has one parameter ``nu``:: 
     395    The chi-squared distribution has one parameter ``nu``:: 
    382396 
    383397        sage: nu = 1 
    384398        sage: T = RealDistribution('chisquared', nu) 
    … …  
    537551 
    538552    def get_random_element(self): 
    539553        """ 
    540         Get a random sample the probability distribution. 
     554        Get a random sample from the probability distribution. 
    541555 
    542556        EXAMPLE:: 
    543557 
    … …  
    560574            result = gsl_ran_pareto(self.r, self.parameters[0], self.parameters[1]) 
    561575        elif self.distribution_type == t: 
    562576            result = gsl_ran_tdist(self.r, self.parameters[0]) 
     577        elif self.distribution_type == F: 
     578            result = gsl_ran_fdist(self.r, self.parameters[0], self.parameters[1]) 
    563579        elif self.distribution_type == chisquared: 
    564580            result = gsl_ran_chisq(self.r, self.parameters[0]) 
    565581        elif self.distribution_type == exppow: 
    … …  
    568584            result = gsl_ran_weibull(self.r, self.parameters[0], self.parameters[1]) 
    569585        elif self.distribution_type == beta: 
    570586            result = gsl_ran_beta(self.r, self.parameters[0], self.parameters[1]) 
     587        else: 
     588            raise TypeError, "Not a valid probability distribution" 
    571589  
    572590        return sage.rings.real_double.RDF(result) 
     591 
    573592    def set_distribution(self, name = 'uniform', parameters = []): 
    574593        """ 
    575594        This method can be called to change the current probability distribution. 
    … …  
    629648                    float(x) 
    630649                except: 
    631650                    raise TypeError, "Lognormal distributions requires real parameters" 
    632   
    633651            self.parameters = <double*>malloc(sizeof(double)*2) 
    634652            self.parameters[0] = float(parameters[0]) 
    635653            self.parameters[1] = float(parameters[1]) 
    … …  
    642660            self.parameters = <double*>malloc(sizeof(double)) 
    643661            self.parameters[0] = float(parameters) 
    644662            self.distribution_type = t 
    645   
     663        elif name == 'F': 
     664            if len(parameters) != 2: 
     665                raise TypeError, "F-distribution requires two real parameters" 
     666            try: 
     667                map(float, parameters) 
     668            except: 
     669                raise TypeError, "F-distribution requires real parameters" 
     670            self.parameters = <double *>malloc(sizeof(double)*2) 
     671            self.parameters[0] = float(parameters[0]) 
     672            self.parameters[1] = float(parameters[1]) 
     673            self.distribution_type = F        
    646674        elif name == 'chisquared': 
    647675            try: 
    648676                float(parameters) 
    … …  
    659687                    float(x) 
    660688                except: 
    661689                    raise TypeError, "exponential power distributions requires real parameters" 
    662   
    663690            self.parameters = <double*>malloc(sizeof(double)*2) 
    664691            self.parameters[0] = float(parameters[0]) 
    665692            self.parameters[1] = float(parameters[1]) 
    … …  
    671698                map(float, parameters) 
    672699            except: 
    673700                raise TypeError, "weibull distribution requires real parameters" 
    674  
    675701            self.parameters = <double *>malloc(sizeof(double)*2) 
    676702            self.parameters[0] = float(parameters[0]) 
    677703            self.parameters[1] = float(parameters[1]) 
    678704            self.distribution_type = weibull 
    679   
    680705        elif name == 'beta': 
    681706            if len(parameters) != 2: 
    682707                raise TypeError, "beta distribution requires two real parameters" 
    … …  
    684709                map(float, parameters) 
    685710            except: 
    686711                raise TypeError, "beta distribution requires real parameters" 
    687  
    688712            self.parameters = <double *>malloc(sizeof(double)*2) 
    689713            self.parameters[0] = float(parameters[0]) 
    690714            self.parameters[1] = float(parameters[1]) 
    691715            self.distribution_type = beta 
    692   
    693   
    694716        else: 
    695717            raise TypeError, "Not a valid probability distribution" 
    696718  
    697719        self.name = name 
    698720  
    699   # def _get_random_element_c(): 
     721    #def _get_random_element_c(): 
    700722  
    701  
    702723    def reset_distribution(self): 
    703724        """ 
    704725        This method resets the distribution. 
    … …  
    748769            return sage.rings.real_double.RDF(gsl_ran_pareto_pdf(x, self.parameters[0], self.parameters[1])) 
    749770        elif self.distribution_type == t: 
    750771            return sage.rings.real_double.RDF(gsl_ran_tdist_pdf(x, self.parameters[0])) 
     772        elif self.distribution_type == F: 
     773            return sage.rings.real_double.RDF(gsl_ran_fdist_pdf(x, self.parameters[0], self.parameters[1])) 
    751774        elif self.distribution_type == chisquared: 
    752775            return sage.rings.real_double.RDF(gsl_ran_chisq_pdf(x, self.parameters[0])) 
    753776        elif self.distribution_type == exppow: 
    … …  
    756779            return sage.rings.real_double.RDF(gsl_ran_weibull_pdf(x, self.parameters[0], self.parameters[1])) 
    757780        elif self.distribution_type == beta: 
    758781            return sage.rings.real_double.RDF(gsl_ran_beta_pdf(x, self.parameters[0], self.parameters[1])) 
    759   
     782        else: 
     783            raise TypeError, "Not a valid probability distribution" 
    760784 
    761785    def cum_distribution_function(self, x): 
    762786        """ 
    … …  
    781805            return sage.rings.real_double.RDF(gsl_cdf_pareto_P(x, self.parameters[0], self.parameters[1])) 
    782806        elif self.distribution_type == t: 
    783807            return sage.rings.real_double.RDF(gsl_cdf_tdist_P(x, self.parameters[0])) 
     808        elif self.distribution_type == F: 
     809            return sage.rings.real_double.RDF(gsl_cdf_fdist_P(x, self.parameters[0], self.parameters[1])) 
    784810        elif self.distribution_type == chisquared: 
    785811            return sage.rings.real_double.RDF(gsl_cdf_chisq_P(x, self.parameters[0])) 
    786812        elif self.distribution_type == exppow: 
    … …  
    789815            return sage.rings.real_double.RDF(gsl_cdf_weibull_P(x, self.parameters[0], self.parameters[1])) 
    790816        elif self.distribution_type == beta: 
    791817            return sage.rings.real_double.RDF(gsl_cdf_beta_P(x, self.parameters[0], self.parameters[1])) 
     818        else: 
     819            raise TypeError, "Not a valid probability distribution" 
    792820  
    793  
    794821    def cum_distribution_function_inv(self, x): 
    795822        """ 
    796823        Evaluate the inverse of the cumulative distribution 
    … …  
    814841            return sage.rings.real_double.RDF(gsl_cdf_pareto_Pinv(x, self.parameters[0], self.parameters[1])) 
    815842        elif self.distribution_type == t: 
    816843            return sage.rings.real_double.RDF(gsl_cdf_tdist_Pinv(x, self.parameters[0])) 
     844        elif self.distribution_type == F: 
     845            return sage.rings.real_double.RDF(gsl_cdf_fdist_Pinv(x, self.parameters[0], self.parameters[1])) 
    817846        elif self.distribution_type == chisquared: 
    818847            return sage.rings.real_double.RDF(gsl_cdf_chisq_Pinv(x, self.parameters[0])) 
    819848        elif self.distribution_type == exppow: 
    … …  
    823852            return sage.rings.real_double.RDF(gsl_cdf_weibull_Pinv(x, self.parameters[0], self.parameters[1])) 
    824853        elif self.distribution_type == beta: 
    825854            return sage.rings.real_double.RDF(gsl_cdf_beta_Pinv(x, self.parameters[0], self.parameters[1])) 
     855        else: 
     856            raise TypeError, "Not a valid probability distribution" 
    826857  
    827858    def plot(self, *args, **kwds): 
    828859        """