Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
package main
import (
"flag"
"log"
"math"
)
// constants and inverse error function
// found from http://libit.sourceforge.net/math_8c.html
const (
erfinv_a3 = -0.140543331
erfinv_a2 = 0.914624893
erfinv_a1 = -1.645349621
erfinv_a0 = 0.886226899
erfinv_b4 = 0.012229801
erfinv_b3 = -0.329097515
erfinv_b2 = 1.442710462
erfinv_b1 = -2.118377725
erfinv_b0 = 1
erfinv_c3 = 1.641345311
erfinv_c2 = 3.429567803
erfinv_c1 = -1.62490649
erfinv_c0 = -1.970840454
erfinv_d2 = 1.637067800
erfinv_d1 = 3.543889200
erfinv_d0 = 1
PositiveSign = 1.0
NegativeSign = -1.0
)
func InvErf(x float64) float64 {
var sign, r, y float64
if x < -1 || x > 1 {
return math.NaN()
}
if x == 0 {
return 0;
}
if x > 0 {
sign = PositiveSign
} else {
sign = NegativeSign
}
if x < 0.7 {
r = x * (((erfinv_a3 * x * x + erfinv_a2) * x * x + erfinv_a1) * x * x + erfinv_a0)
r /= (((erfinv_b4 * x * x + erfinv_b3) * x * x + erfinv_b2) * x * x + erfinv_b1) * x * x + erfinv_b0
} else {
y = math.Sqrt(-math.Log((1 - x) / 2))
r = (((erfinv_c3 * y + erfinv_c2) * y + erfinv_c1) * y + erfinv_c0)
r /= ((erfinv_d2 * y + erfinv_d1) * y + erfinv_d0)
}
r = r * sign
x = x * sign
r -= (math.Erf(r) - x) / (2 / math.Sqrt(math.Pi) * math.Exp(-r * r))
r -= (math.Erf(r) - x) / (2 / math.Sqrt(math.Pi) * math.Exp(-r * r))
return r
}
type Simulation struct {
Mean float64
Stdev float64
Probability float64
}
func (sim Simulation) InvNorm() float64 {
var sign float64
if sim.Probability >= 0.5 {
sign = PositiveSign
} else {
sign = NegativeSign
}
return sim.Mean + (sim.Stdev * sign * math.Sqrt(2) * InvErf((2*sim.Probability) - 1))
}
func (sim Simulation) InvLogNorm() float64 {
return math.Exp(sim.InvNorm())
}
func main() {
prob := flag.Float64("probability", 0.5, "probability from 0.0 to 1.0")
mean := flag.Float64("mean", 0, "mean (average) of the distribution")
stdev := flag.Float64("stdev", 1, "standard deviation of the distribution")
flag.Parse()
sim := Simulation{*mean, *stdev, *prob}
log.Println("compare your answers against Google Sheets'")
log.Printf("NORMINV(%v, %v, %v) = %v", sim.Probability, sim.Mean, sim.Stdev, sim.InvNorm())
log.Printf("LOGINV(%v, %v, %v) = %v", sim.Probability, sim.Mean, sim.Stdev, sim.InvLogNorm())
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment