diff --git a/stat/distmv/normal.go b/stat/distmv/normal.go
index b4c10a165..8270e8037 100644
--- a/stat/distmv/normal.go
+++ b/stat/distmv/normal.go
@@ -277,14 +277,13 @@ func (n *Normal) Rand(x []float64) []float64 {
 	return x
 }
 
-// ScoreInput returns the score function with respect to the input of the
-// distribution at the location specified by x. The score function is the
-// derivative of the log-probability at the input.
-//  ∂/∂x log(p(x))
-// If deriv is nil, a new slice will be allocated and returned. If deriv is of
+// ScoreInput returns the gradient of the log-probability with respect to the
+// input x. That is, ScoreInput computes
+//  ∇_x log(p(x))
+// If score is nil, a new slice will be allocated and returned. If score is of
 // length the dimension of Normal, then the result will be put in-place into deriv.
 // If neither of these is true, ScoreInput will panic.
-func (n *Normal) ScoreInput(deriv, x []float64) []float64 {
+func (n *Normal) ScoreInput(score, x []float64) []float64 {
 	// Normal log probability is
 	//  c - 0.5*(x-μ)' Σ^-1 (x-μ).
 	// So the derivative is just
@@ -292,21 +291,21 @@ func (n *Normal) ScoreInput(deriv, x []float64) []float64 {
 	if len(x) != n.Dim() {
 		panic(badInputLength)
 	}
-	if deriv == nil {
-		deriv = make([]float64, len(x))
+	if score == nil {
+		score = make([]float64, len(x))
 	}
-	if len(deriv) != len(x) {
+	if len(score) != len(x) {
 		panic(badSizeMismatch)
 	}
 	tmp := make([]float64, len(x))
 	copy(tmp, x)
 	floats.Sub(tmp, n.mu)
 
-	dv := mat.NewVector(len(deriv), deriv)
+	dv := mat.NewVector(len(score), score)
 	dt := mat.NewVector(len(tmp), tmp)
 	dv.SolveCholeskyVec(&n.chol, dt)
-	floats.Scale(-1, deriv)
-	return deriv
+	floats.Scale(-1, score)
+	return score
 }
 
 // SetMean changes the mean of the normal distribution. SetMean panics if len(mu)
diff --git a/stat/distmv/normal_test.go b/stat/distmv/normal_test.go
index 36203ce88..64d33bc94 100644
--- a/stat/distmv/normal_test.go
+++ b/stat/distmv/normal_test.go
@@ -561,13 +561,13 @@ func TestNormalScoreInput(t *testing.T) {
 		}
 		x := make([]float64, len(test.x))
 		copy(x, test.x)
-		deriv := normal.ScoreInput(nil, x)
+		score := normal.ScoreInput(nil, x)
 		if !floats.Equal(x, test.x) {
 			t.Errorf("x modified during call to ScoreInput")
 		}
-		derivFD := fd.Gradient(nil, normal.LogProb, x, nil)
-		if !floats.EqualApprox(deriv, derivFD, 1e-4) {
-			t.Errorf("Case %d: derivative mismatch. Got %v, want %v", cas, deriv, derivFD)
+		scoreFD := fd.Gradient(nil, normal.LogProb, x, nil)
+		if !floats.EqualApprox(score, scoreFD, 1e-4) {
+			t.Errorf("Case %d: derivative mismatch. Got %v, want %v", cas, score, scoreFD)
 		}
 	}
 }