add normal quantile
MODIFY
distmv/normal.go
distmv/normal.go
@@ -247,6 +247,31 @@ func (n *Normal) Prob(x []float64) float64 {
return math.Exp(n.LogProb(x))
}
+// Quantile returns the multi-dimensional inverse cumulative distribution function.
+// If x is nil, a new slice will be allocated and returned. If x is non-nil,
+// len(x) must equal len(p) and the quantile will be stored in-place into x.
+// All of the values of p must be between 0 and 1, inclusive, or Quantile will panic.
+func (n *Normal) Quantile(x, p []float64) []float64 {
+ dim := n.Dim()
+ if len(p) != dim {
+ panic(badInputLength)
+ }
+ if x == nil {
+ x = make([]float64, dim)
+ }
+ if len(x) != len(p) {
+ panic(badInputLength)
+ }
+
+ // Transform to a standard normal and then transform to a multivariate Gaussian.
+ tmp := make([]float64, len(x))
+ for i, v := range p {
+ tmp[i] = distuv.UnitNormal.Quantile(v)
+ }
+ n.TransformNormal(x, tmp)
+ return x
+}
+
// Rand generates a random number according to the distributon.
// If the input slice is nil, new memory is allocated, otherwise the result is stored
// in place.

MODIFY
distmv/normal_test.go
distmv/normal_test.go
@@ -126,6 +126,50 @@ func TestNormRand(t *testing.T) {
}
}
+func TestNormalQuantile(t *testing.T) {
+ for _, test := range []struct {
+ mean []float64
+ cov []float64
+ }{
+ {
+ mean: []float64{6, 7},
+ cov: []float64{
+ 5, 0.9,
+ 0.9, 2,
+ },
+ },
+ } {
+ dim := len(test.mean)
+ cov := mat64.NewSymDense(dim, test.cov)
+ n, ok := NewNormal(test.mean, cov, nil)
+ if !ok {
+ t.Errorf("bad covariance matrix")
+ }
+
+ nSamples := 1000000
+ rnd := rand.New(rand.NewSource(1))
+ samps := mat64.NewDense(nSamples, dim, nil)
+ tmp := make([]float64, dim)
+ for i := 0; i < nSamples; i++ {
+ for j := range tmp {
+ tmp[j] = rnd.Float64()
+ }
+ n.Quantile(samps.RawRowView(i), tmp)
+ }
+ estMean := make([]float64, dim)
+ for i := range estMean {
+ estMean[i] = stat.Mean(mat64.Col(nil, i, samps), nil)
+ }
+ if !floats.EqualApprox(estMean, test.mean, 1e-2) {
+ t.Errorf("Mean mismatch: want: %v, got %v", test.mean, estMean)
+ }
+ estCov := stat.CovarianceMatrix(nil, samps, nil)
+ if !mat64.EqualApprox(estCov, cov, 1e-2) {
+ t.Errorf("Cov mismatch: want: %v, got %v", cov, estCov)
+ }
+ }
+}
+
func TestConditionNormal(t *testing.T) {
// Uncorrelated values shouldn't influence the updated values.
for _, test := range []struct {

add normal quantile