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

     1  // Copyright ©2018 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  	"gonum.org/v1/gonum/mathext"
    13  )
    14  
    15  // InverseGamma implements the inverse gamma distribution, a two-parameter
    16  // continuous distribution with support over the positive real numbers. The
    17  // inverse gamma distribution is the same as the distribution of the reciprocal
    18  // of a gamma distributed random variable.
    19  //
    20  // The inverse gamma distribution has density function
    21  //
    22  //	β^α / Γ(α) x^(-α-1)e^(-β/x)
    23  //
    24  // For more information, see https://en.wikipedia.org/wiki/Inverse-gamma_distribution
    25  type InverseGamma struct {
    26  	// Alpha is the shape parameter of the distribution. Alpha must be greater than 0.
    27  	Alpha float64
    28  	// Beta is the scale parameter of the distribution. Beta must be greater than 0.
    29  	Beta float64
    30  
    31  	Src rand.Source
    32  }
    33  
    34  // CDF computes the value of the cumulative distribution function at x.
    35  func (g InverseGamma) CDF(x float64) float64 {
    36  	if x < 0 {
    37  		return 0
    38  	}
    39  	// TODO(btracey): Replace this with a direct call to the upper regularized
    40  	// gamma function if mathext gets it.
    41  	//return 1 - mathext.GammaInc(g.Alpha, g.Beta/x)
    42  	return mathext.GammaIncRegComp(g.Alpha, g.Beta/x)
    43  }
    44  
    45  // ExKurtosis returns the excess kurtosis of the distribution.
    46  func (g InverseGamma) ExKurtosis() float64 {
    47  	if g.Alpha <= 4 {
    48  		return math.Inf(1)
    49  	}
    50  	return (30*g.Alpha - 66) / (g.Alpha - 3) / (g.Alpha - 4)
    51  }
    52  
    53  // LogProb computes the natural logarithm of the value of the probability
    54  // density function at x.
    55  func (g InverseGamma) LogProb(x float64) float64 {
    56  	if x <= 0 {
    57  		return math.Inf(-1)
    58  	}
    59  	a := g.Alpha
    60  	b := g.Beta
    61  	lg, _ := math.Lgamma(a)
    62  	return a*math.Log(b) - lg + (-a-1)*math.Log(x) - b/x
    63  }
    64  
    65  // Mean returns the mean of the probability distribution.
    66  func (g InverseGamma) Mean() float64 {
    67  	if g.Alpha <= 1 {
    68  		return math.Inf(1)
    69  	}
    70  	return g.Beta / (g.Alpha - 1)
    71  }
    72  
    73  // Mode returns the mode of the distribution.
    74  func (g InverseGamma) Mode() float64 {
    75  	return g.Beta / (g.Alpha + 1)
    76  }
    77  
    78  // NumParameters returns the number of parameters in the distribution.
    79  func (InverseGamma) NumParameters() int {
    80  	return 2
    81  }
    82  
    83  // Prob computes the value of the probability density function at x.
    84  func (g InverseGamma) Prob(x float64) float64 {
    85  	return math.Exp(g.LogProb(x))
    86  }
    87  
    88  // Quantile returns the inverse of the cumulative distribution function.
    89  func (g InverseGamma) Quantile(p float64) float64 {
    90  	if p < 0 || 1 < p {
    91  		panic(badPercentile)
    92  	}
    93  	return (1 / (mathext.GammaIncRegCompInv(g.Alpha, p))) * g.Beta
    94  }
    95  
    96  // Rand returns a random sample drawn from the distribution.
    97  //
    98  // Rand panics if either alpha or beta is <= 0.
    99  func (g InverseGamma) Rand() float64 {
   100  	// TODO(btracey): See if there is a more direct way to sample.
   101  	return 1 / Gamma(g).Rand()
   102  }
   103  
   104  // Survival returns the survival function (complementary CDF) at x.
   105  func (g InverseGamma) Survival(x float64) float64 {
   106  	if x < 0 {
   107  		return 1
   108  	}
   109  	return mathext.GammaIncReg(g.Alpha, g.Beta/x)
   110  }
   111  
   112  // StdDev returns the standard deviation of the probability distribution.
   113  func (g InverseGamma) StdDev() float64 {
   114  	return math.Sqrt(g.Variance())
   115  }
   116  
   117  // Variance returns the variance of the probability distribution.
   118  func (g InverseGamma) Variance() float64 {
   119  	if g.Alpha <= 2 {
   120  		return math.Inf(1)
   121  	}
   122  	v := g.Beta / (g.Alpha - 1)
   123  	return v * v / (g.Alpha - 2)
   124  }