sampleuv

type IIDer ¶

type IIDer struct {
	Dist distuv.Rander
}

IIDer generates a set of independently and identically distributed samples from the input distribution.

func (iid IIDer) Sample(batch []float64)

type Importance ¶

type Importance struct {
	Target   distuv.LogProber
	Proposal distuv.RandLogProber
}

Importance is a type for performing importance sampling using the given Target and Proposal distributions.

Importance sampling is a variance reduction technique where samples are generated from a proposal distribution, q(x), instead of the target distribution p(x). This allows relatively unlikely samples in p(x) to be generated more frequently.

The importance sampling weight at x is given by p(x)/q(x). To reduce variance, a good proposal distribution will bound this sampling weight. This implies the support of q(x) should be at least as broad as p(x), and q(x) should be "fatter tailed" than p(x).

func (l Importance) SampleWeighted(batch, weights []float64)

type LatinHypercube ¶

type LatinHypercube struct {
	Q   distuv.Quantiler
	Src rand.Source
}

LatinHypercube is a type for sampling using Latin hypercube sampling from the given distribution. If src is not nil, it will be used to generate random numbers, otherwise rand.Float64 will be used.

Latin hypercube sampling divides the cumulative distribution function into equally spaced bins and guarantees that one sample is generated per bin. Within each bin, the location is randomly sampled. The distuv.UnitUniform variable can be used for easy sampling from the unit hypercube.

func (l LatinHypercube) Sample(batch []float64)

type MHProposal ¶

type MHProposal interface {
	// ConditionalDist returns the probability of the first argument conditioned on
	// being at the second argument
	//  p(x|y)
	ConditionalLogProb(x, y float64) (prob float64)

	// ConditionalRand generates a new random location conditioned being at the
	// location y.
	ConditionalRand(y float64) (x float64)
}

MHProposal defines a proposal distribution for Metropolis Hastings.

type MetropolisHastings ¶

type MetropolisHastings struct {
	Initial  float64
	Target   distuv.LogProber
	Proposal MHProposal
	Src      rand.Source

	BurnIn int
	Rate   int
}

MetropolisHastings is a type for generating samples using the Metropolis Hastings algorithm (http://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm), with the given target and proposal distributions, starting at the location specified by Initial. If src != nil, it will be used to generate random numbers, otherwise rand.Float64 will be used.

Metropolis-Hastings is a Markov-chain Monte Carlo algorithm that generates samples according to the distribution specified by target using the Markov chain implicitly defined by the proposal distribution. At each iteration, a proposal point is generated randomly from the current location. This proposal point is accepted with probability

p = min(1, (target(new) * proposal(current|new)) / (target(current) * proposal(new|current)))

If the new location is accepted, it becomes the new current location. If it is rejected, the current location remains. This is the sample stored in batch, ignoring BurnIn and Rate (discussed below).

The samples in Metropolis Hastings are correlated with one another through the Markov chain. As a result, the initial value can have a significant influence on the early samples, and so, typically, the first samples generated by the chain are ignored. This is known as "burn-in", and the number of samples ignored at the beginning is specified by BurnIn. The proper BurnIn value will depend on the mixing time of the Markov chain defined by the target and proposal distributions.

Many choose to have a sampling "rate" where a number of samples are ignored in between each kept sample. This helps decorrelate the samples from one another, but also reduces the number of available samples. This value is specified by Rate. If Rate is 0 it is defaulted to 1 (keep every sample).

The initial value is NOT changed during calls to Sample.

package main

import (
	"gonum.org/v1/gonum/stat/distuv"
	"gonum.org/v1/gonum/stat/sampleuv"
)

type ProposalDist struct {
	Sigma float64
}

func (p ProposalDist) ConditionalRand(y float64) float64 {
	return distuv.Normal{Mu: y, Sigma: p.Sigma}.Rand()
}

func (p ProposalDist) ConditionalLogProb(x, y float64) float64 {
	return distuv.Normal{Mu: y, Sigma: p.Sigma}.LogProb(x)
}

func main() {
	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 := distuv.Weibull{K: 5, Lambda: 0.5}
	// proposal is the proposal distribution. Here, we are choosing
	// a tight Gaussian distribution around the current location. In
	// typical problems, if Sigma is too small, it takes a lot of samples
	// to move around the distribution. If Sigma is too large, it can be hard
	// to find acceptable samples.
	proposal := ProposalDist{Sigma: 0.2}

	samples := make([]float64, n)
	mh := sampleuv.MetropolisHastings{Initial: initial, Target: target, Proposal: proposal, BurnIn: burnin}
	mh.Sample(samples)
}

Output:

// See Burnin example for a description of these quantities.
n := 1000
burnin := 300
var initial float64
target := distuv.Weibull{K: 5, Lambda: 0.5}
proposal := ProposalDist{Sigma: 0.2}

// Successive samples are correlated with one another through the
// Markov Chain defined by the proposal distribution. One may use
// a sampling rate to decrease the correlation in the samples for
// an increase in computation cost. The rate parameter specifies
// that for every accepted sample stored in `samples`, rate - 1 accepted
// samples are not stored in `samples`.
rate := 50

mh := sampleuv.MetropolisHastings{
	Initial:  initial,
	Target:   target,
	Proposal: proposal,
	BurnIn:   burnin,
	Rate:     rate,
}

samples := make([]float64, n)
mh.Sample(samples)

Output:

func (m MetropolisHastings) Sample(batch []float64)

type Rejection ¶

type Rejection struct {
	C        float64
	Target   distuv.LogProber
	Proposal distuv.RandLogProber
	Src      rand.Source
	// contains filtered or unexported fields
}

Rejection is a type for sampling using the rejection sampling algorithm.

Rejection sampling generates points from the target distribution by using the proposal distribution. At each step of the algorithm, the proposed point is accepted with probability

p = target(x) / (proposal(x) * c)

where target(x) is the probability of the point according to the target distribution and proposal(x) is the probability according to the proposal distribution. The constant c must be chosen such that target(x) < proposal(x) * c for all x. The expected number of proposed samples is len(samples) * c.

The number of proposed locations during sampling can be found with a call to Proposed. If there was an error during sampling, all elements of samples are set to NaN and the error can be accessed with the Err method. If src != nil, it will be used to generate random numbers, otherwise rand.Float64 will be used.

Target may return the true (log of) the probability of the location, or it may return a value that is proportional to the probability (logprob + constant). This is useful for cases where the probability distribution is only known up to a normalization constant.

func (r *Rejection) Err() error

func (r *Rejection) Proposed() int

func (r *Rejection) Sample(batch []float64)

type SampleUniformWeighted ¶

type SampleUniformWeighted struct {
	Sampler
}

SampleUniformWeighted wraps a Sampler type to create a WeightedSampler where all weights are equal.

func (w SampleUniformWeighted) SampleWeighted(batch, weights []float64)

type Sampler ¶

type Sampler interface {
	Sample(batch []float64)
}

Sampler generates a batch of samples according to the rule specified by the implementing type. The number of samples generated is equal to len(batch), and the samples are stored in-place into the input.

type Weighted ¶

type Weighted struct {
	// contains filtered or unexported fields
}

Weighted provides sampling without replacement from a collection of items with non-uniform probability.

func NewWeighted(w []float64, src rand.Source) Weighted

func (s Weighted) Len() int

func (s Weighted) Reweight(idx int, w float64)

func (s Weighted) ReweightAll(w []float64)

func (s Weighted) Take() (idx int, ok bool)

type WeightedSampler ¶

type WeightedSampler interface {
	SampleWeighted(batch, weights []float64)
}

WeightedSampler generates a batch of samples and their relative weights according to the rule specified by the implementing type. The number of samples generated is equal to len(batch), and the samples and weights are stored in-place into the inputs. The length of weights must equal len(batch), otherwise SampleWeighted will panic.