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
5package distuv
6
7import (
8	"math"
9
10	"golang.org/x/exp/rand"
11
12	"gonum.org/v1/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
23type 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.
30func (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.
37func (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.
46func (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.
53func (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.
63func (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.
71func (f F) NumParameters() int {
72	return 2
73}
74
75// Prob computes the value of the probability density function at x.
76func (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.
81func (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.
90func (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.
99func (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.
111func (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.
119func (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.
126func (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}
134