Golang math.Erf函数代码示例

// https://stackoverflow.com/questions/5971830/need-code-for-inverse-error-function
func erfinv(y float64) float64 {
	if y < -1.0 || y > 1.0 {
		panic("invalid input")
	}

	var (
		a = [4]float64{0.886226899, -1.645349621, 0.914624893, -0.140543331}
		b = [4]float64{-2.118377725, 1.442710462, -0.329097515, 0.012229801}
		c = [4]float64{-1.970840454, -1.624906493, 3.429567803, 1.641345311}
		d = [2]float64{3.543889200, 1.637067800}
	)

	const y0 = 0.7
	var x, z float64

	if math.Abs(y) == 1.0 {
		x = -y * math.Log(0.0)
	} else if y < -y0 {
		z = math.Sqrt(-math.Log((1.0 + y) / 2.0))
		x = -(((c[3]*z+c[2])*z+c[1])*z + c[0]) / ((d[1]*z+d[0])*z + 1.0)
	} else {
		if y < y0 {
			z = y * y
			x = y * (((a[3]*z+a[2])*z+a[1])*z + a[0]) / ((((b[3]*z+b[3])*z+b[1])*z+b[0])*z + 1.0)
		} else {
			z = math.Sqrt(-math.Log((1.0 - y) / 2.0))
			x = (((c[3]*z+c[2])*z+c[1])*z + c[0]) / ((d[1]*z+d[0])*z + 1.0)
		}
		x = x - (math.Erf(x)-y)/(2.0/math.SqrtPi*math.Exp(-x*x))
		x = x - (math.Erf(x)-y)/(2.0/math.SqrtPi*math.Exp(-x*x))
	}

	return x
}

// Black-Scholes (European put and call options)
// C Theoretical call premium (non-dividend paying stock)
// c = sn(d1) - ke^(-rt)N(d2)
// d1 = ln(s/k) + (r + s^2/2)t
// d2 = d1- st^1/2
// v = stock value
// k = Stock strike price
// s = Spot price
// t = time to expire in years
// r = risk free rate
// v = volitilaty (sigma)
// e math.E 2.7183
// putcall = "c" for a call or "p" for a put
func BS(s, k, t, r, v float64, putcall string) float64 {
	d1 := (math.Log(s/k) + ((r + (math.Pow(v, 2) / 2)) * t)) / (v * math.Sqrt(t))
	d2 := d1 - (v * math.Sqrt(t))
	if putcall == call {
		return math.Erf(d1)*s - (math.Erf(d2) * k * math.Pow(math.E, (-1*r*t)))
	}
	if putcall == put {
		return k*math.Pow(math.E, (-1*r*t)) - s + ((math.Erf(-1*d2) * k * math.Pow(math.E, (-1*r*t))) - (math.Erf(-1*d1) * s))
	}
	panic(NOOR)
}

func Test_2dinteg02(tst *testing.T) {

	//verbose()
	chk.PrintTitle("2dinteg02. bidimensional integral")

	// Γ(1/4, 1)
	gamma_1div4_1 := 0.2462555291934987088744974330686081384629028737277219

	x := utl.LinSpace(0, 1, 11)
	y := utl.LinSpace(0, 1, 11)
	m, n := len(x), len(y)
	f := la.MatAlloc(m, n)
	for i := 0; i < m; i++ {
		for j := 0; j < n; j++ {
			f[i][j] = 8.0 * math.Exp(-math.Pow(x[i], 2)-math.Pow(y[j], 4))
		}
	}
	dx, dy := x[1]-x[0], y[1]-y[0]
	Vt := Trapz2D(dx, dy, f)
	Vs := Simps2D(dx, dy, f)
	Vc := math.Sqrt(math.Pi) * math.Erf(1) * (math.Gamma(1.0/4.0) - gamma_1div4_1)
	io.Pforan("Vt = %v\n", Vt)
	io.Pforan("Vs = %v\n", Vs)
	io.Pfgreen("Vc = %v\n", Vc)
	chk.Scalar(tst, "Vt", 0.0114830435645548, Vt, Vc)
	chk.Scalar(tst, "Vs", 1e-4, Vs, Vc)

}

func main() {
	for true {
		r := bufio.NewReader(os.Stdin)
		s, err := r.ReadString('\n')
		if err == os.EOF {
			break
		}
		s = strings.TrimRight(s, "\n")
		a := strings.Split(s, " ")
		f := a[0]
		x, err := strconv.Atof64(a[1])
		switch f {
		case "erf":
			fmt.Println(math.Erf(x))
		case "expm1":
			fmt.Println(math.Expm1(x))
		case "phi":
			fmt.Println(phi.Phi(x))
		case "NormalCDFInverse":
			fmt.Println(normal_cdf_inverse.NormalCDFInverse(x))
		case "Gamma":
			fmt.Println(math.Gamma(x))
		case "LogGamma":
			r, _ := math.Lgamma(x)
			fmt.Println(r)
		case "LogFactorial":
			fmt.Println(log_factorial.LogFactorial(int(x)))
		default:
			fmt.Println("Unknown function: " + f)
			return
		}
	}
}

func avalancheTest1(t *testing.T, k Key) {
	const REP = 100000
	r := rand.New(rand.NewSource(1234))
	n := k.bits()

	// grid[i][j] is a count of whether flipping
	// input bit i affects output bit j.
	grid := make([][hashSize]int, n)

	for z := 0; z < REP; z++ {
		// pick a random key, hash it
		k.random(r)
		h := k.hash()

		// flip each bit, hash & compare the results
		for i := 0; i < n; i++ {
			k.flipBit(i)
			d := h ^ k.hash()
			k.flipBit(i)

			// record the effects of that bit flip
			g := &grid[i]
			for j := 0; j < hashSize; j++ {
				g[j] += int(d & 1)
				d >>= 1
			}
		}
	}

	// Each entry in the grid should be about REP/2.
	// More precisely, we did N = k.bits() * hashSize experiments where
	// each is the sum of REP coin flips.  We want to find bounds on the
	// sum of coin flips such that a truly random experiment would have
	// all sums inside those bounds with 99% probability.
	N := n * hashSize
	var c float64
	// find c such that Prob(mean-c*stddev < x < mean+c*stddev)^N > .99
	for c = 0.0; math.Pow(math.Erf(c/math.Sqrt(2)), float64(N)) < .99; c += .1 {
	}
	c *= 2.0 // allowed slack - we don't need to be perfectly random
	mean := .5 * REP
	stddev := .5 * math.Sqrt(REP)
	low := int(mean - c*stddev)
	high := int(mean + c*stddev)
	for i := 0; i < n; i++ {
		for j := 0; j < hashSize; j++ {
			x := grid[i][j]
			if x < low || x > high {
				t.Errorf("bad bias for %s bit %d -> bit %d: %d/%d\n", k.name(), i, j, x, REP)
			}
		}
	}
}

func pValueFromNormalDifferences(correct []float64, estimate []float64) ([]float64, error) {
	if len(correct) != len(estimate) {
		return nil, fmt.Errorf("p-value calculation requires two lists of equal size")
	}
	deviations := []float64{}
	for i := range correct {
		if math.IsInf(correct[i], 0) || math.IsNaN(correct[i]) {
			continue
		}
		if math.IsInf(estimate[i], 0) || math.IsNaN(estimate[i]) {
			continue
		}
		deviations = append(deviations, estimate[i]-correct[i])
	}
	meandev := 0.0
	for _, dev := range deviations {
		meandev += dev
	}
	meandev /= float64(len(deviations))
	stddev := 0.0
	for _, dev := range deviations {
		stddev += (meandev - dev) * (meandev - dev)
	}
	stddev /= float64(len(deviations)) - 1
	stddev = math.Sqrt(stddev) // Now we have a good estimate for the population standard deviation.
	pvalues := make([]float64, len(estimate))
	for i := range pvalues {
		if math.IsInf(correct[i], 0) || math.IsNaN(correct[i]) {
			pvalues[i] = math.NaN()
			continue
		}
		if math.IsInf(estimate[i], 0) || math.IsNaN(estimate[i]) {
			pvalues[i] = math.NaN()
			continue
		}
		difference := estimate[i] - correct[i]
		tvalue := (difference - meandev) / stddev
		pvalue := 1 - math.Erf(math.Abs(tvalue)/math.Sqrt2)
		pvalues[i] = pvalue
	}
	return pvalues, nil
}

//calculates the attempt frequency of a particle
func attemptf2(p *particle) float64 {
	volume := cube(p.r) * 4 / 3. * math.Pi
	hk := 2. * Ku1 / (p.msat * mu0)

	gprime := Alpha * gamma0 * mu0 / (1. + (Alpha * Alpha))

	delta := Ku1 * volume / (kb * Temp)

	bx, by, bz := B_ext(T)
	bextvector := vector{bx / mu0, by / mu0, bz / mu0}
	hoverhk := math.Abs(bextvector.dot(p.u_anis)) / hk
	if math.Signbit(bextvector.dot(p.m)) {
		hoverhk *= -1.
	}

	postpostfactor := 1. / (math.Erf(math.Sqrt(delta) * (1 - hoverhk)))

	return gprime * hk / math.Sqrt(delta*math.Pi) * postpostfactor * math.Exp(-delta)

}

// CDF calculates the cumulative distribution function for any standard
// distribution.
func CDF(x, mu, sigma float64) float64 {
	return 0.5 * (1 + math.Erf((x-mu)/(sigma*math.Sqrt2)))
}

//Unit Normal CDF
func (u *UnitNormal) CDF(x float64) float64 {
	return 0.5 * (1.0 + math.Erf(x/(math.Sqrt2)))
}

func (n *Normal) CalcCDF(x float64) (float64, error) {
	return 0.5 * (1.0 + math.Erf((x-n.Mean)/(n.StdDev*math.Sqrt2))), nil
}

// Cumulative Distribution Function for the Normal distribution
func Normal_CDF(μ, σ float64) func(x float64) float64 {
	return func(x float64) float64 { return ((1.0 / 2.0) * (1 + math.Erf((x-μ)/(σ*math.Sqrt2)))) }
}

// CDF computes the value of the cumulative density function at x.
func (n Normal) CDF(x float64) float64 {
	return 0.5 * (1 + math.Erf((x-n.Mu)/(n.Sigma*math.Sqrt2)))
}

func GaussCumulativeTo(x float64) float64 {
	return math.Erf(x/math.Sqrt2)/2 + 0.5
}

/* Cumulative Distribution Function for the Normal distribution */
func Normal_CDF(mu, sigma float64) func(x float64) float64 {
	return func(x float64) float64 {
		return (0.5 * (1 - math.Erf((mu-x)/(sigma*math.Sqrt2))))
	}
}

// Cdf computes the cumulative probability function @ x
func (o DistLogNormal) Cdf(x float64) float64 {
	if x < TOLMINLOG {
		return 0
	}
	return (1.0 + math.Erf((math.Log(x)-o.M)/(o.S*math.Sqrt2))) / 2.0
}