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