github.com/gopherd/gonum@v0.0.4/stat/sampleuv/sample_test.go (about)

     1  // Copyright ©2015 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package sampleuv
     6  
     7  import (
     8  	"sort"
     9  	"testing"
    10  
    11  	"github.com/gopherd/gonum/floats/scalar"
    12  	"github.com/gopherd/gonum/stat"
    13  	"github.com/gopherd/gonum/stat/distuv"
    14  )
    15  
    16  const tol = 1e-2
    17  
    18  type lhDist interface {
    19  	Quantile(float64) float64
    20  	CDF(float64) float64
    21  }
    22  
    23  func TestLatinHypercube(t *testing.T) {
    24  	for _, nSamples := range []int{1, 2, 5, 10, 20} {
    25  		samples := make([]float64, nSamples)
    26  		for _, dist := range []lhDist{
    27  			distuv.Uniform{Min: 0, Max: 1},
    28  			distuv.Uniform{Min: 0, Max: 10},
    29  			distuv.Normal{Mu: 5, Sigma: 3},
    30  		} {
    31  			LatinHypercube{Q: dist}.Sample(samples)
    32  			sort.Float64s(samples)
    33  			for i, v := range samples {
    34  				p := dist.CDF(v)
    35  				if p < float64(i)/float64(nSamples) || p > float64(i+1)/float64(nSamples) {
    36  					t.Errorf("probability out of bounds")
    37  				}
    38  			}
    39  		}
    40  	}
    41  }
    42  
    43  func TestImportance(t *testing.T) {
    44  	// Test by finding the expected value of a Normal.
    45  	trueMean := 3.0
    46  	target := distuv.Normal{Mu: trueMean, Sigma: 2}
    47  	proposal := distuv.Normal{Mu: 0, Sigma: 5}
    48  	nSamples := 100000
    49  	x := make([]float64, nSamples)
    50  	weights := make([]float64, nSamples)
    51  	Importance{Target: target, Proposal: proposal}.SampleWeighted(x, weights)
    52  	ev := stat.Mean(x, weights)
    53  	if !scalar.EqualWithinAbsOrRel(ev, trueMean, tol, tol) {
    54  		t.Errorf("Mean mismatch: Want %v, got %v", trueMean, ev)
    55  	}
    56  }
    57  
    58  func TestRejection(t *testing.T) {
    59  	// Test by finding the expected value of a Normal.
    60  	trueMean := 3.0
    61  	target := distuv.Normal{Mu: trueMean, Sigma: 2}
    62  	proposal := distuv.Normal{Mu: 0, Sigma: 5}
    63  
    64  	nSamples := 20000
    65  	x := make([]float64, nSamples)
    66  	r := &Rejection{Target: target, Proposal: proposal, C: 100}
    67  	r.Sample(x)
    68  	ev := stat.Mean(x, nil)
    69  	if !scalar.EqualWithinAbsOrRel(ev, trueMean, tol, tol) {
    70  		t.Errorf("Mean mismatch: Want %v, got %v", trueMean, ev)
    71  	}
    72  }
    73  
    74  type condNorm struct {
    75  	Sigma float64
    76  }
    77  
    78  func (c condNorm) ConditionalRand(y float64) float64 {
    79  	return distuv.Normal{Mu: y, Sigma: c.Sigma}.Rand()
    80  }
    81  
    82  func (c condNorm) ConditionalLogProb(x, y float64) float64 {
    83  	return distuv.Normal{Mu: y, Sigma: c.Sigma}.LogProb(x)
    84  }
    85  
    86  func TestMetropolisHastings(t *testing.T) {
    87  	// Test by finding the expected value of a Normal.
    88  	trueMean := 3.0
    89  	target := distuv.Normal{Mu: trueMean, Sigma: 2}
    90  	proposal := condNorm{Sigma: 5}
    91  
    92  	burnin := 500
    93  	nSamples := 100000 + burnin
    94  	x := make([]float64, nSamples)
    95  	mh := MetropolisHastings{
    96  		Initial:  100,
    97  		Target:   target,
    98  		Proposal: proposal,
    99  
   100  		BurnIn: burnin,
   101  	}
   102  	mh.Sample(x)
   103  
   104  	ev := stat.Mean(x, nil)
   105  	if !scalar.EqualWithinAbsOrRel(ev, trueMean, tol, tol) {
   106  		t.Errorf("Mean mismatch: Want %v, got %v", trueMean, ev)
   107  	}
   108  }