loginv.go

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())
}