gonum.org/v1/gonum@v0.14.0/stat/distuv/pareto.go (about)

     1  // Copyright ©2017 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package distuv
     6  
     7  import (
     8  	"math"
     9  
    10  	"golang.org/x/exp/rand"
    11  )
    12  
    13  // Pareto implements the Pareto (Type I) distribution, a one parameter distribution
    14  // with support above the scale parameter.
    15  //
    16  // The density function is given by
    17  //
    18  //	(α x_m^{α})/(x^{α+1}) for x >= x_m.
    19  //
    20  // For more information, see https://en.wikipedia.org/wiki/Pareto_distribution.
    21  type Pareto struct {
    22  	// Xm is the scale parameter.
    23  	// Xm must be greater than 0.
    24  	Xm float64
    25  
    26  	// Alpha is the shape parameter.
    27  	// Alpha must be greater than 0.
    28  	Alpha float64
    29  
    30  	Src rand.Source
    31  }
    32  
    33  // CDF computes the value of the cumulative density function at x.
    34  func (p Pareto) CDF(x float64) float64 {
    35  	if x < p.Xm {
    36  		return 0
    37  	}
    38  	return -math.Expm1(p.Alpha * math.Log(p.Xm/x))
    39  }
    40  
    41  // Entropy returns the differential entropy of the distribution.
    42  func (p Pareto) Entropy() float64 {
    43  	return math.Log(p.Xm) - math.Log(p.Alpha) + (1 + 1/p.Alpha)
    44  }
    45  
    46  // ExKurtosis returns the excess kurtosis of the distribution.
    47  func (p Pareto) ExKurtosis() float64 {
    48  	if p.Alpha <= 4 {
    49  		return math.NaN()
    50  	}
    51  	return 6 * (p.Alpha*p.Alpha*p.Alpha + p.Alpha*p.Alpha - 6*p.Alpha - 2) / (p.Alpha * (p.Alpha - 3) * (p.Alpha - 4))
    52  
    53  }
    54  
    55  // LogProb computes the natural logarithm of the value of the probability
    56  // density function at x.
    57  func (p Pareto) LogProb(x float64) float64 {
    58  	if x < p.Xm {
    59  		return math.Inf(-1)
    60  	}
    61  	return math.Log(p.Alpha) + p.Alpha*math.Log(p.Xm) - (p.Alpha+1)*math.Log(x)
    62  }
    63  
    64  // Mean returns the mean of the probability distribution.
    65  func (p Pareto) Mean() float64 {
    66  	if p.Alpha <= 1 {
    67  		return math.Inf(1)
    68  	}
    69  	return p.Alpha * p.Xm / (p.Alpha - 1)
    70  }
    71  
    72  // Median returns the median of the pareto distribution.
    73  func (p Pareto) Median() float64 {
    74  	return p.Quantile(0.5)
    75  }
    76  
    77  // Mode returns the mode of the distribution.
    78  func (p Pareto) Mode() float64 {
    79  	return p.Xm
    80  }
    81  
    82  // NumParameters returns the number of parameters in the distribution.
    83  func (p Pareto) NumParameters() int {
    84  	return 2
    85  }
    86  
    87  // Prob computes the value of the probability density function at x.
    88  func (p Pareto) Prob(x float64) float64 {
    89  	return math.Exp(p.LogProb(x))
    90  }
    91  
    92  // Quantile returns the inverse of the cumulative probability distribution.
    93  func (p Pareto) Quantile(prob float64) float64 {
    94  	if prob < 0 || 1 < prob {
    95  		panic(badPercentile)
    96  	}
    97  	return p.Xm / math.Pow(1-prob, 1/p.Alpha)
    98  }
    99  
   100  // Rand returns a random sample drawn from the distribution.
   101  func (p Pareto) Rand() float64 {
   102  	var rnd float64
   103  	if p.Src == nil {
   104  		rnd = rand.ExpFloat64()
   105  	} else {
   106  		rnd = rand.New(p.Src).ExpFloat64()
   107  	}
   108  	return p.Xm * math.Exp(rnd/p.Alpha)
   109  }
   110  
   111  // StdDev returns the standard deviation of the probability distribution.
   112  func (p Pareto) StdDev() float64 {
   113  	return math.Sqrt(p.Variance())
   114  }
   115  
   116  // Survival returns the survival function (complementary CDF) at x.
   117  func (p Pareto) Survival(x float64) float64 {
   118  	if x < p.Xm {
   119  		return 1
   120  	}
   121  	return math.Pow(p.Xm/x, p.Alpha)
   122  }
   123  
   124  // Variance returns the variance of the probability distribution.
   125  func (p Pareto) Variance() float64 {
   126  	if p.Alpha <= 2 {
   127  		return math.Inf(1)
   128  	}
   129  	am1 := p.Alpha - 1
   130  	return p.Xm * p.Xm * p.Alpha / (am1 * am1 * (p.Alpha - 2))
   131  }