github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/stat/distuv/alphastable.go (about)

     1  // Copyright ©2020 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  // AlphaStable represents an α-stable distribution with four parameters.
    14  // See https://en.wikipedia.org/wiki/Stable_distribution for more information.
    15  type AlphaStable struct {
    16  	// Alpha is the stability parameter.
    17  	// It is valid within the range 0 < α ≤ 2.
    18  	Alpha float64
    19  	// Beta is the skewness parameter.
    20  	// It is valid within the range -1 ≤ β ≤ 1.
    21  	Beta float64
    22  	// C is the scale parameter.
    23  	// It is valid when positive.
    24  	C float64
    25  	// Mu is the location parameter.
    26  	Mu  float64
    27  	Src rand.Source
    28  }
    29  
    30  // ExKurtosis returns the excess kurtosis of the distribution.
    31  // ExKurtosis returns NaN when Alpha != 2.
    32  func (a AlphaStable) ExKurtosis() float64 {
    33  	if a.Alpha == 2 {
    34  		return 0
    35  	}
    36  	return math.NaN()
    37  }
    38  
    39  // Mean returns the mean of the probability distribution.
    40  // Mean returns NaN when Alpha <= 1.
    41  func (a AlphaStable) Mean() float64 {
    42  	if a.Alpha > 1 {
    43  		return a.Mu
    44  	}
    45  	return math.NaN()
    46  }
    47  
    48  // Median returns the median of the distribution.
    49  // Median panics when Beta != 0, because then the mode is not analytically
    50  // expressible.
    51  func (a AlphaStable) Median() float64 {
    52  	if a.Beta == 0 {
    53  		return a.Mu
    54  	}
    55  	panic("distuv: cannot compute Median for Beta != 0")
    56  }
    57  
    58  // Mode returns the mode of the distribution.
    59  // Mode panics when Beta != 0, because then the mode is not analytically
    60  // expressible.
    61  func (a AlphaStable) Mode() float64 {
    62  	if a.Beta == 0 {
    63  		return a.Mu
    64  	}
    65  	panic("distuv: cannot compute Mode for Beta != 0")
    66  }
    67  
    68  // NumParameters returns the number of parameters in the distribution.
    69  func (a AlphaStable) NumParameters() int {
    70  	return 4
    71  }
    72  
    73  // Rand returns a random sample drawn from the distribution.
    74  func (a AlphaStable) Rand() float64 {
    75  	// From https://en.wikipedia.org/wiki/Stable_distribution#Simulation_of_stable_variables
    76  	const halfPi = math.Pi / 2
    77  	u := Uniform{-halfPi, halfPi, a.Src}.Rand()
    78  	w := Exponential{1, a.Src}.Rand()
    79  	if a.Alpha == 1 {
    80  		f := halfPi + a.Beta*u
    81  		x := (f*math.Tan(u) - a.Beta*math.Log(halfPi*w*math.Cos(u)/f)) / halfPi
    82  		return a.C*(x+a.Beta*math.Log(a.C)/halfPi) + a.Mu
    83  	}
    84  	zeta := -a.Beta * math.Tan(halfPi*a.Alpha)
    85  	xi := math.Atan(-zeta) / a.Alpha
    86  	f := a.Alpha * (u + xi)
    87  	g := math.Sqrt(1+zeta*zeta) * math.Pow(math.Cos(u-f)/w, 1-a.Alpha) / math.Cos(u)
    88  	x := math.Pow(g, 1/a.Alpha) * math.Sin(f)
    89  	return a.C*x + a.Mu
    90  }
    91  
    92  // Skewness returns the skewness of the distribution.
    93  // Skewness returns NaN when Alpha != 2.
    94  func (a AlphaStable) Skewness() float64 {
    95  	if a.Alpha == 2 {
    96  		return 0
    97  	}
    98  	return math.NaN()
    99  }
   100  
   101  // StdDev returns the standard deviation of the probability distribution.
   102  func (a AlphaStable) StdDev() float64 {
   103  	return math.Sqrt(a.Variance())
   104  }
   105  
   106  // Variance returns the variance of the probability distribution.
   107  // Variance returns +Inf when Alpha != 2.
   108  func (a AlphaStable) Variance() float64 {
   109  	if a.Alpha == 2 {
   110  		return 2 * a.C * a.C
   111  	}
   112  	return math.Inf(1)
   113  }