func InvCDF()

in internal/stats/dist.go [109:179]
56 lines of code
14 McCabe index (conditional complexity)

func InvCDF(dist DistCommon) func(y float64) (x float64) {
	type invCDF interface {
		InvCDF(float64) float64
	}
	if dist, ok := dist.(invCDF); ok {
		return dist.InvCDF
	}

	// Otherwise, use a numerical algorithm.
	//
	// TODO: For discrete distributions, use the step size to
	// inform this computation.
	return func(y float64) (x float64) {
		const almostInf = 1e100
		const xtol = 1e-16

		if y < 0 || y > 1 {
			return nan
		} else if y == 0 {
			l, _ := dist.Bounds()
			if dist.CDF(l) == 0 {
				// Finite support
				return l
			} else {
				// Infinite support
				return -inf
			}
		} else if y == 1 {
			_, h := dist.Bounds()
			if dist.CDF(h) == 1 {
				// Finite support
				return h
			} else {
				// Infinite support
				return inf
			}
		}

		// Find loX, hiX for which cdf(loX) < y <= cdf(hiX).
		var loX, loY, hiX, hiY float64
		x1, y1 := 0.0, dist.CDF(0)
		xdelta := 1.0
		if y1 < y {
			hiX, hiY = x1, y1
			for hiY < y && hiX != inf {
				loX, loY, hiX = hiX, hiY, hiX+xdelta
				hiY = dist.CDF(hiX)
				xdelta *= 2
			}
		} else {
			loX, loY = x1, y1
			for y <= loY && loX != -inf {
				hiX, hiY, loX = loX, loY, loX-xdelta
				loY = dist.CDF(loX)
				xdelta *= 2
			}
		}
		if loX == -inf {
			return loX
		} else if hiX == inf {
			return hiX
		}

		// Use bisection on the interval to find the smallest
		// x at which cdf(x) <= y.
		_, x = bisectBool(func(x float64) bool {
			return dist.CDF(x) < y
		}, loX, hiX, xtol)
		return
	}
}