add implementation for gamma distribution

This commit is contained in:
NunoSempere 2024-06-09 23:01:49 +02:00
parent 5f1188c516
commit d57300e8ce

View File

@ -47,6 +47,58 @@ func Sample_to(low float64, high float64, r Src) float64 {
return math.Exp(Sample_normal_from_90_ci(loglow, loghigh, r)) return math.Exp(Sample_normal_from_90_ci(loglow, loghigh, r))
} }
func Sample_gamma(alpha float64, r Src) float64 {
// a simple method for generating gamma variables, marsaglia and wan tsang, 2001
// https://dl.acm.org/doi/pdf/10.1145/358407.358414
// see also the references/ folder
// note that the wikipedia page for the gamma distribution includes a scaling parameter
// k or beta
// https://en.wikipedia.org/wiki/gamma_distribution
// such that gamma_k(alpha, k) = k * gamma(alpha)
// or gamma_beta(alpha, beta) = gamma(alpha) / beta
// so far i have not needed to use this, and thus the second parameter is by default 1.
if alpha >= 1 {
var d, c, x, v, u float64
d = alpha - 1.0/3.0
c = 1.0 / math.Sqrt(9.0*d)
OuterLoop:
for {
InnerLoop:
for {
x = Sample_unit_normal(r)
v = 1.0 + c*x
if v > 0.0 {
break InnerLoop
}
}
v = v * v * v
u = Sample_unit_uniform(r)
if u < 1.0-0.0331*(x*x*x*x) { // Condition 1
// the 0.0331 doesn't inspire much confidence
// however, this isn't the whole story
// by knowing that Condition 1 implies condition 2
// we realize that this is just a way of making the algorithm faster
// i.e., of not using the logarithms
return d * v
}
if math.Log(u) < 0.5*(x*x)+d*(1.0-v+math.Log(v)) { // Condition 2
return d * v
}
}
} else {
return Sample_gamma(1.0+alpha, r) * math.Pow(Sample_unit_uniform(r), 1.0/alpha)
}
}
func Sample_mixture(fs []func64, weights []float64, r Src) float64 { func Sample_mixture(fs []func64, weights []float64, r Src) float64 {
// fmt.Println("weights initially: ", weights) // fmt.Println("weights initially: ", weights)