github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/stat/distuv/f.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  	"github.com/jingcheng-WU/gonum/mathext"
    13  )
    14  
    15  // F implements the F-distribution, a two-parameter continuous distribution
    16  // with support over the positive real numbers.
    17  //
    18  // The F-distribution has density function
    19  //  sqrt(((d1*x)^d1) * d2^d2 / ((d1*x+d2)^(d1+d2))) / (x * B(d1/2,d2/2))
    20  // where B is the beta function.
    21  //
    22  // For more information, see https://en.wikipedia.org/wiki/F-distribution
    23  type F struct {
    24  	D1  float64 // Degrees of freedom for the numerator
    25  	D2  float64 // Degrees of freedom for the denominator
    26  	Src rand.Source
    27  }
    28  
    29  // CDF computes the value of the cumulative density function at x.
    30  func (f F) CDF(x float64) float64 {
    31  	return mathext.RegIncBeta(f.D1/2, f.D2/2, f.D1*x/(f.D1*x+f.D2))
    32  }
    33  
    34  // ExKurtosis returns the excess kurtosis of the distribution.
    35  //
    36  // ExKurtosis returns NaN if the D2 parameter is less or equal to 8.
    37  func (f F) ExKurtosis() float64 {
    38  	if f.D2 <= 8 {
    39  		return math.NaN()
    40  	}
    41  	return (12 / (f.D2 - 6)) * ((5*f.D2-22)/(f.D2-8) + ((f.D2-4)/f.D1)*((f.D2-2)/(f.D2-8))*((f.D2-2)/(f.D1+f.D2-2)))
    42  }
    43  
    44  // LogProb computes the natural logarithm of the value of the probability
    45  // density function at x.
    46  func (f F) LogProb(x float64) float64 {
    47  	return 0.5*(f.D1*math.Log(f.D1*x)+f.D2*math.Log(f.D2)-(f.D1+f.D2)*math.Log(f.D1*x+f.D2)) - math.Log(x) - mathext.Lbeta(f.D1/2, f.D2/2)
    48  }
    49  
    50  // Mean returns the mean of the probability distribution.
    51  //
    52  // Mean returns NaN if the D2 parameter is less than or equal to 2.
    53  func (f F) Mean() float64 {
    54  	if f.D2 <= 2 {
    55  		return math.NaN()
    56  	}
    57  	return f.D2 / (f.D2 - 2)
    58  }
    59  
    60  // Mode returns the mode of the distribution.
    61  //
    62  // Mode returns NaN if the D1 parameter is less than or equal to 2.
    63  func (f F) Mode() float64 {
    64  	if f.D1 <= 2 {
    65  		return math.NaN()
    66  	}
    67  	return ((f.D1 - 2) / f.D1) * (f.D2 / (f.D2 + 2))
    68  }
    69  
    70  // NumParameters returns the number of parameters in the distribution.
    71  func (f F) NumParameters() int {
    72  	return 2
    73  }
    74  
    75  // Prob computes the value of the probability density function at x.
    76  func (f F) Prob(x float64) float64 {
    77  	return math.Exp(f.LogProb(x))
    78  }
    79  
    80  // Quantile returns the inverse of the cumulative distribution function.
    81  func (f F) Quantile(p float64) float64 {
    82  	if p < 0 || p > 1 {
    83  		panic(badPercentile)
    84  	}
    85  	y := mathext.InvRegIncBeta(0.5*f.D1, 0.5*f.D2, p)
    86  	return f.D2 * y / (f.D1 * (1 - y))
    87  }
    88  
    89  // Rand returns a random sample drawn from the distribution.
    90  func (f F) Rand() float64 {
    91  	u1 := ChiSquared{f.D1, f.Src}.Rand()
    92  	u2 := ChiSquared{f.D2, f.Src}.Rand()
    93  	return (u1 / f.D1) / (u2 / f.D2)
    94  }
    95  
    96  // Skewness returns the skewness of the distribution.
    97  //
    98  // Skewness returns NaN if the D2 parameter is less than or equal to 6.
    99  func (f F) Skewness() float64 {
   100  	if f.D2 <= 6 {
   101  		return math.NaN()
   102  	}
   103  	num := (2*f.D1 + f.D2 - 2) * math.Sqrt(8*(f.D2-4))
   104  	den := (f.D2 - 6) * math.Sqrt(f.D1*(f.D1+f.D2-2))
   105  	return num / den
   106  }
   107  
   108  // StdDev returns the standard deviation of the probability distribution.
   109  //
   110  // StdDev returns NaN if the D2 parameter is less than or equal to 4.
   111  func (f F) StdDev() float64 {
   112  	if f.D2 <= 4 {
   113  		return math.NaN()
   114  	}
   115  	return math.Sqrt(f.Variance())
   116  }
   117  
   118  // Survival returns the survival function (complementary CDF) at x.
   119  func (f F) Survival(x float64) float64 {
   120  	return 1 - f.CDF(x)
   121  }
   122  
   123  // Variance returns the variance of the probability distribution.
   124  //
   125  // Variance returns NaN if the D2 parameter is less than or equal to 4.
   126  func (f F) Variance() float64 {
   127  	if f.D2 <= 4 {
   128  		return math.NaN()
   129  	}
   130  	num := 2 * f.D2 * f.D2 * (f.D1 + f.D2 - 2)
   131  	den := f.D1 * (f.D2 - 2) * (f.D2 - 2) * (f.D2 - 4)
   132  	return num / den
   133  }