Fix documentation and permute LHC

sample/example_burnin_test.go

@@ -18,7 +18,7 @@ func ExampleMetropolisHastings_burnin() {
 	n := 1000    // The number of samples to generate.
 	burnin := 50 // Number of samples to ignore at the start.
 	var initial float64
-	// target is the distribution from which we would like to sample
+	// target is the distribution from which we would like to sample.
 	target := dist.Weibull{K: 5, Lambda: 0.5}
 	// proposal is the proposal distribution. Here, we are choosing
 	// a tight Gaussian distribution around the current location. In
@@ -30,6 +30,6 @@ func ExampleMetropolisHastings_burnin() {
 	samples := make([]float64, n+burnin)
 	MetropolisHastings(samples, initial, target, proposal, nil)
 
-	// Remove the initial samples through slicing
+	// Remove the initial samples through slicing.
 	samples = samples[burnin:]
 }

sample/example_rate_test.go

@@ -10,7 +10,7 @@ func max(a, b int) int {
 }
 
 func ExampleMetropolisHastings_samplingRate() {
-	// See Burnin example for a description of these quantities
+	// See Burnin example for a description of these quantities.
 	n := 1000
 	burnin := 300
 	var initial float64

sample/sample.go

@@ -27,13 +27,18 @@ var (
 // for easy generation from the unit interval.
 func LatinHypercube(samples []float64, q dist.Quantiler, src *rand.Rand) {
 	n := len(samples)
-	f64 := rand.Float64
+	var perm []int
+	var f64 func() float64
 	if src != nil {
 		f64 = src.Float64
+		perm = src.Perm(n)
+	} else {
+		f64 = rand.Float64
+		perm = rand.Perm(n)
 	}
 	for i := range samples {
 		v := f64()/float64(n) + float64(i)/float64(n)
-		samples[i] = q.Quantile(v)
+		samples[perm[i]] = q.Quantile(v)
 	}
 }
 

sample/sample_test.go

@@ -5,6 +5,7 @@ package sample
 
 import (
 	"math"
+	"sort"
 	"testing"
 
 	"github.com/gonum/stat"
@@ -25,6 +26,7 @@ func TestLatinHypercube(t *testing.T) {
 			dist.Normal{Mu: 5, Sigma: 3},
 		} {
 			LatinHypercube(samples, dist, nil)
+			sort.Float64s(samples)
 			for i, v := range samples {
 				p := dist.CDF(v)
 				if p < float64(i)/float64(nSamples) || p > float64(i+1)/float64(nSamples) {