github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/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  //  (α x_m^{α})/(x^{α+1}) for x >= x_m.
    18  //
    19  // For more information, see https://en.wikipedia.org/wiki/Pareto_distribution.
    20  type Pareto struct {
    21  	// Xm is the scale parameter.
    22  	// Xm must be greater than 0.
    23  	Xm float64
    24  
    25  	// Alpha is the shape parameter.
    26  	// Alpha must be greater than 0.
    27  	Alpha float64
    28  
    29  	Src rand.Source
    30  }
    31  
    32  // CDF computes the value of the cumulative density function at x.
    33  func (p Pareto) CDF(x float64) float64 {
    34  	if x < p.Xm {
    35  		return 0
    36  	}
    37  	return -math.Expm1(p.Alpha * math.Log(p.Xm/x))
    38  }
    39  
    40  // Entropy returns the differential entropy of the distribution.
    41  func (p Pareto) Entropy() float64 {
    42  	return math.Log(p.Xm) - math.Log(p.Alpha) + (1 + 1/p.Alpha)
    43  }
    44  
    45  // ExKurtosis returns the excess kurtosis of the distribution.
    46  func (p Pareto) ExKurtosis() float64 {
    47  	if p.Alpha <= 4 {
    48  		return math.NaN()
    49  	}
    50  	return 6 * (p.Alpha*p.Alpha*p.Alpha + p.Alpha*p.Alpha - 6*p.Alpha - 2) / (p.Alpha * (p.Alpha - 3) * (p.Alpha - 4))
    51  
    52  }
    53  
    54  // LogProb computes the natural logarithm of the value of the probability
    55  // density function at x.
    56  func (p Pareto) LogProb(x float64) float64 {
    57  	if x < p.Xm {
    58  		return math.Inf(-1)
    59  	}
    60  	return math.Log(p.Alpha) + p.Alpha*math.Log(p.Xm) - (p.Alpha+1)*math.Log(x)
    61  }
    62  
    63  // Mean returns the mean of the probability distribution.
    64  func (p Pareto) Mean() float64 {
    65  	if p.Alpha <= 1 {
    66  		return math.Inf(1)
    67  	}
    68  	return p.Alpha * p.Xm / (p.Alpha - 1)
    69  }
    70  
    71  // Median returns the median of the pareto distribution.
    72  func (p Pareto) Median() float64 {
    73  	return p.Quantile(0.5)
    74  }
    75  
    76  // Mode returns the mode of the distribution.
    77  func (p Pareto) Mode() float64 {
    78  	return p.Xm
    79  }
    80  
    81  // NumParameters returns the number of parameters in the distribution.
    82  func (p Pareto) NumParameters() int {
    83  	return 2
    84  }
    85  
    86  // Prob computes the value of the probability density function at x.
    87  func (p Pareto) Prob(x float64) float64 {
    88  	return math.Exp(p.LogProb(x))
    89  }
    90  
    91  // Quantile returns the inverse of the cumulative probability distribution.
    92  func (p Pareto) Quantile(prob float64) float64 {
    93  	if prob < 0 || 1 < prob {
    94  		panic(badPercentile)
    95  	}
    96  	return p.Xm / math.Pow(1-prob, 1/p.Alpha)
    97  }
    98  
    99  // Rand returns a random sample drawn from the distribution.
   100  func (p Pareto) Rand() float64 {
   101  	var rnd float64
   102  	if p.Src == nil {
   103  		rnd = rand.ExpFloat64()
   104  	} else {
   105  		rnd = rand.New(p.Src).ExpFloat64()
   106  	}
   107  	return math.Exp(math.Log(p.Xm) + 1/p.Alpha*rnd)
   108  }
   109  
   110  // StdDev returns the standard deviation of the probability distribution.
   111  func (p Pareto) StdDev() float64 {
   112  	return math.Sqrt(p.Variance())
   113  }
   114  
   115  // Survival returns the survival function (complementary CDF) at x.
   116  func (p Pareto) Survival(x float64) float64 {
   117  	if x < p.Xm {
   118  		return 1
   119  	}
   120  	return math.Pow(p.Xm/x, p.Alpha)
   121  }
   122  
   123  // Variance returns the variance of the probability distribution.
   124  func (p Pareto) Variance() float64 {
   125  	if p.Alpha <= 2 {
   126  		return math.Inf(1)
   127  	}
   128  	am1 := p.Alpha - 1
   129  	return p.Xm * p.Xm * p.Alpha / (am1 * am1 * (p.Alpha - 2))
   130  }