Golang math.Lgamma函数代码示例

/*Beta function
Special cases:
z<0: NaN
w<0: NaN
*/
func Beta(z float64, w float64) float64 {
	if z < 0 || w < 0 {
		return math.NaN()
	}
	a, _ := math.Lgamma(z)
	b, _ := math.Lgamma(w)
	c, _ := math.Lgamma(z + w)
	return math.Exp(a + b - c)
}

func lfastchoose2(n, k float64) (float64, int) {
	// mathematically the same as lfastchoose()
	// less stable typically, but useful if n-k+1 < 0
	//	r := lgammafn_sign(n-k+1, s_choose)
	r, s_choose := math.Lgamma(n - k + 1)
	p, _ := math.Lgamma(n + 1)
	q, _ := math.Lgamma(k + 1)
	return p - q - r, s_choose
}

func (sampler *Sampler) updateComponentE(targetEvent int) (lgamma float64) {
	var tPositive, tNegative, tNormalize float64
	for k := 0; k < numTop; k++ {
		tPositive, _ = math.Lgamma(float64(sampler.Model.Eventtype_histogram[k]) + sampler.eventPrior)
		tNegative, _ = math.Lgamma(float64((sampler.Model.NumESDs)-sampler.Model.Eventtype_histogram[k]) + sampler.eventPrior)
		tNormalize, _ = math.Lgamma(float64(sampler.Model.NumESDs) + 2*sampler.eventPrior)
		lgamma += ((tPositive + tNegative) - tNormalize)
	}
	return lgamma
}

func (sampler *Sampler) updateComponentP(participantID, eventID int) (lgamma float64) {
	// for each alternative participanttype
	var pPositive, pNegative, pNormalize float64
	for i := 0; i < numPar; i++ {
		pPositive, _ = math.Lgamma(float64(sampler.Model.Participanttype_eventtype_histogram[i][eventID]) + sampler.participantPrior)
		pNegative, _ = math.Lgamma(float64(sampler.Model.Eventtype_histogram[eventID]-sampler.Model.Participanttype_eventtype_histogram[i][eventID]) + sampler.participantPrior)
		pNormalize, _ = math.Lgamma(float64(sampler.Model.Eventtype_histogram[eventID]) + 2*sampler.participantPrior)
		lgamma += ((pPositive + pNegative) - pNormalize)
	}
	return
}

// Calculates Ix(a, b). Borrowed from
// NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING
// ISBN (0-521-43064-Z)
func betaInc(a, b, x float64) float64 {
	lgab, _ := math.Lgamma(a + b)
	lga, _ := math.Lgamma(a)
	lgb, _ := math.Lgamma(b)
	exp := lgab - lga - lgb + a*math.Log(x) + b*math.Log(1.0-x)
	bt := math.Exp(exp)
	if x < (a+1.0)/(a+b+2.0) {
		return bt * betaCF(a, b, x) / a
	} else {
		return 1.0 - bt*betaCF(b, a, 1.0-x)/b
	}
}

// LogGeneralizedBinomial returns the log of the generalized binomial coefficient.
// See GeneralizedBinomial for more information.
func LogGeneralizedBinomial(n, k float64) float64 {
	if n < 0 || k < 0 {
		panic(badNegInput)
	}
	if n < k {
		panic(badSetSize)
	}
	a, _ := math.Lgamma(n + 1)
	b, _ := math.Lgamma(k + 1)
	c, _ := math.Lgamma(n - k + 1)
	return a - b - c
}

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

// Lower incomplete gamma.
func lgamma(x, s float64, regularized bool) float64 {
	if x == 0 {
		return 0
	}
	if x < 0 || s <= 0 {
		return math.NaN()
	}

	if x > 1.1 && x > s {
		if regularized {
			return 1.0 - ugamma(x, s, regularized)
		}
		return math.Gamma(s) - ugamma(x, s, regularized)
	}

	var ft float64
	r := s
	c := 1.0
	pws := 1.0

	if regularized {
		logg, _ := math.Lgamma(s)
		ft = s*math.Log(x) - x - logg
	} else {
		ft = s*math.Log(x) - x
	}
	ft = math.Exp(ft)
	for c/pws > eps {
		r++
		c *= x / r
		pws += c
	}
	return pws * ft / s
}

// Upper incomplete gamma.
func ugamma(x, s float64, regularized bool) float64 {
	if x <= 1.1 || x <= s {
		if regularized {
			return 1 - lgamma(x, s, regularized)
		}
		return math.Gamma(s) - lgamma(x, s, regularized)
	}

	f := 1.0 + x - s
	C := f
	D := 0.0
	var a, b, chg float64

	for i := 1; i < 10000; i++ {
		a = float64(i) * (s - float64(i))
		b = float64(i<<1) + 1.0 + x - s
		D = b + a*D
		C = b + a/C
		D = 1.0 / D
		chg = C * D
		f *= chg
		if math.Abs(chg-1) < eps {
			break
		}
	}
	if regularized {
		logg, _ := math.Lgamma(s)
		return math.Exp(s*math.Log(x) - x - logg - math.Log(f))
	}
	return math.Exp(s*math.Log(x) - x - math.Log(f))
}

// Skewness returns the skewness of the distribution.
func (w Weibull) Skewness() float64 {
	stdDev := w.StdDev()
	firstGamma, firstGammaSign := math.Lgamma(1 + 3/w.K)
	logFirst := firstGamma + 3*(math.Log(w.Lambda)-math.Log(stdDev))
	logSecond := math.Log(3) + math.Log(w.Mean()) + 2*math.Log(stdDev) - 3*math.Log(stdDev)
	logThird := 3 * (math.Log(w.Mean()) - math.Log(stdDev))
	return float64(firstGammaSign)*math.Exp(logFirst) - math.Exp(logSecond) - math.Exp(logThird)
}

/*Normalized Gamma Function
(or Complementary Incomplete Gamma Function)

Equal to Gamma(a, x)/Gamma(a)
Evaluated by Legendre's continued fraction

Special Cases:
G(0, 0) = Infinity
G(0, positive) = 1.0
G(0, negative) = 1/(-a)
*/
func G(a float64, x float64) float64 {
	if x == 0 {
		if a == 0 {
			return math.Inf(1)
		} else if a > 0 {
			return 1.0
		} else if a < 0 {
			return 1.0 / math.Abs(a)
		}
	}

	//Evaluate Legendre's continued fraction
	//using Lentz's algorithm
	//Shift from Thomson and Barnett
	//Continued fraction:
	//http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/

	b0 := x + 1.0 - a
	C := 1.0 / (10E-30)
	D := 1.0 / b0
	if b0 == 0 {
		D = 10E30
	}
	f := D
	//numerator
	an := func(n int) float64 {
		if n == 0 {
			return 1.0
		}
		return -1.0 * float64(n) * (float64(n) - a)
	}
	//denominator
	bn := func(n int) float64 {
		return x + (float64(2*n + 1)) - a
	}
	//Lentz's algorithm until machine precision or 1,000 iterations
	for j := 1; j < 1000; j++ {
		D = bn(j) + an(j)*D
		if math.Abs(D) < 10E-20 {
			D = 10E-30
		}
		C = bn(j) + an(j)/C
		if math.Abs(C) < 10E-20 {
			C = 10E-30
		}
		D = 1.0 / D
		del := D * C
		f *= del
		if math.Abs(del-1.0) < 10E-15 {
			break
		}
	}
	lnGa, _ := math.Lgamma(a)
	return f * math.Exp(-x+a*math.Log(x)-lnGa)
}

// compute document likelihood of the events in the current esd
// all participant labelings will stay constant -> no need to compute them!
func (sampler *Sampler) documentLikelihood(label Label) float64 {
	var wordTypeFactor, wordFactor, wordNorm float64
	var typeWordTotal int
	documentLikelihood := 0.0
	// iterate over eventtypes
	for k := 0; k < numTop; k++ {
		wordFactor = 0.0
		typeWordTotal = 0
		// iterate over terms in event-vocab
		for term, histogram := range sampler.Model.Word_eventtype_histogram {
			typeWordTotal += histogram[k]
			// compute LGamma(N(word,event) + prior + udpate)
			wordTypeFactor, _ = math.Lgamma(float64(histogram[k]) + sampler.EventlmPriors[k][term])
			wordFactor += wordTypeFactor
		}
		// normalize LGamma(N(words_by_event) + V*prior + total_update)
		wordNorm, _ = math.Lgamma(float64(typeWordTotal) + sum(sampler.EventlmPriors[k]))
		documentLikelihood += (wordFactor - wordNorm)
	}
	return documentLikelihood
}

// compute document likelihood of the participant realization in question, given the proposed label
// all event doc likelihoods will stay constant w.r.t. change -> no need to compute them!
func (sampler *Sampler) documentLikelihoodP(event int, participant int, label Label) float64 {
	var wordTypeFactor, wordFactor, wordNorm float64
	var typeWordTotal /*, update*/ int
	documentLikelihood := 0.0
	// iterate over participanttypes
	for i := 0; i < numPar; i++ {
		wordFactor = 0.0
		typeWordTotal = 0
		// iterate over terms in participant vocab
		for term, histogram := range sampler.Model.Word_participanttype_histogram {
			typeWordTotal += histogram[i]
			// set 'update' according to the number of times term is present in current particip descr
			// compute LGamma(N(word,part) + prior + update)
			wordTypeFactor, _ = math.Lgamma(float64(histogram[i]) + sampler.ParticipantlmPriors[i][term])
			wordFactor += wordTypeFactor
		}
		// normalize
		wordNorm, _ = math.Lgamma(float64(typeWordTotal) + sum(sampler.ParticipantlmPriors[i]))
		documentLikelihood += (wordFactor - wordNorm)
	}
	return documentLikelihood
}

//Gamma Distribution PDF
func (g *Gamma) PDF(x float64) float64 {
	if x <= 0 {
		return math.NaN()
	}
	if math.IsInf(x, 1) {
		return 1.0
	} else if math.IsInf(x, -1) {
		return 0.0
	}
	lga, _ := math.Lgamma(g.alpha)
	logp := (g.alpha * math.Log(g.beta)) + ((g.alpha - 1.0) * math.Log(x)) - (x * g.beta) - lga
	return math.Exp(logp)
}

func TestLnGamma(t *testing.T) {
	acc := 0.0000001
	check := func(x, y float64) bool {
		if false {
			return x == y
		}
		return math.Abs(x-y) < acc
	}
	for i := 0; i < 100; i++ {
		x := NextGamma(10, 10)
		g1 := LnΓ(x)
		g2, _ := math.Lgamma(x)
		if !check(g1, g2) {
			t.Error(fmt.Sprintf("For %v: %v vs %v", x, g1, g2))
		}
	}
	//var start int64
	Seed(10)
	start := time.Now()
	for i := 0; i < 1e6; i++ {
		x := NextGamma(10, 10)
		math.Lgamma(x)
	}
	now := time.Now()
	duration2 := float64(now.Sub(start)) / 1e9

	//duration2 := float64(time.Now()-start) / 1e9
	Seed(10)
	start = time.Now()
	for i := 0; i < 1e6; i++ {
		x := NextGamma(10, 10)
		LnΓ(x)
	}
	now = time.Now()
	duration1 := float64(now.Sub(start)) / 1e9
	fmt.Printf("Mine was %f\nTheirs was %f\n", duration1, duration2)
}