stat/distuv: correct Gamma Mode() and LogProb(0)/Prob(0) for alpha <= 1

stat/distuv/gamma.go

@@ -47,12 +47,15 @@ func (g Gamma) ExKurtosis() float64 {
 // LogProb computes the natural logarithm of the value of the probability
 // density function at x.
 func (g Gamma) LogProb(x float64) float64 {
-	if x <= 0 {
+	if x < 0 {
 		return math.Inf(-1)
 	}
 	a := g.Alpha
 	b := g.Beta
 	lg, _ := math.Lgamma(a)
+	if a == 1 {
+		return math.Log(b) - lg - b*x
+	}
 	return a*math.Log(b) - lg + (a-1)*math.Log(x) - b*x
 }
 
@@ -63,11 +66,11 @@ func (g Gamma) Mean() float64 {
 
 // Mode returns the mode of the normal distribution.
 //
-// The mode is NaN in the special case where the Alpha (shape) parameter
+// The mode is 0 in the special case where the Alpha (shape) parameter
 // is less than 1.
 func (g Gamma) Mode() float64 {
 	if g.Alpha < 1 {
-		return math.NaN()
+		return 0
 	}
 	return (g.Alpha - 1) / g.Beta
 }

stat/distuv/gamma_test.go

@@ -83,12 +83,20 @@ func testGamma(t *testing.T, f Gamma, i int) {
 	checkProbContinuous(t, i, x, 0, math.Inf(1), f, quadTol)
 	checkQuantileCDFSurvival(t, i, x, f, 5e-2)
 	if f.Alpha < 1 {
-		if !math.IsNaN(f.Mode()) {
-			t.Errorf("Expected NaN mode for alpha < 1, got %v", f.Mode())
+		if f.Mode() != 0 {
+			t.Errorf("Expected 0 mode for alpha < 1, got %v", f.Mode())
+		}
+		if !math.IsInf(f.Prob(0), 1) {
+			t.Errorf("Expected pdf(0) to be +Infinity for alpha < 1, got %v", f.Prob(0))
 		}
 	} else {
 		checkMode(t, i, x, f, 2e-1, 1)
 	}
+	if f.Alpha == 1 {
+		if math.IsInf(f.LogProb(0), 0) {
+			t.Errorf("Expected log(pdf(0)) to be finite for alpha = 1, got %v", f.LogProb(0))
+		}
+	}
 	cdfNegX := f.CDF(-0.0001)
 	if cdfNegX != 0 {
 		t.Errorf("Expected zero CDF for negative argument, got %v", cdfNegX)