| // Copyright ©2016 The Gonum Authors. All rights reserved. |
| // Use of this source code is governed by a BSD-style |
| // license that can be found in the LICENSE file. |
| |
| package distmv |
| |
| import ( |
| "math" |
| "sort" |
| |
| "golang.org/x/exp/rand" |
| "golang.org/x/tools/container/intsets" |
| |
| "gonum.org/v1/gonum/floats" |
| "gonum.org/v1/gonum/mat" |
| "gonum.org/v1/gonum/stat" |
| "gonum.org/v1/gonum/stat/distuv" |
| ) |
| |
| // StudentsT is a multivariate Student's T distribution. It is a distribution over |
| // ℝ^n with the probability density |
| // p(y) = (Γ((ν+n)/2) / Γ(ν/2)) * (νπ)^(-n/2) * |Ʃ|^(-1/2) * |
| // (1 + 1/ν * (y-μ)ᵀ * Ʃ^-1 * (y-μ))^(-(ν+n)/2) |
| // where ν is a scalar greater than 2, μ is a vector in ℝ^n, and Ʃ is an n×n |
| // symmetric positive definite matrix. |
| // |
| // In this distribution, ν sets the spread of the distribution, similar to |
| // the degrees of freedom in a univariate Student's T distribution. As ν → ∞, |
| // the distribution approaches a multi-variate normal distribution. |
| // μ is the mean of the distribution, and the covariance is ν/(ν-2)*Ʃ. |
| // |
| // See https://en.wikipedia.org/wiki/Student%27s_t-distribution and |
| // http://users.isy.liu.se/en/rt/roth/student.pdf for more information. |
| type StudentsT struct { |
| nu float64 |
| mu []float64 |
| // If src is altered, rnd must be updated. |
| src rand.Source |
| rnd *rand.Rand |
| |
| sigma mat.SymDense // only stored if needed |
| |
| chol mat.Cholesky |
| lower mat.TriDense |
| logSqrtDet float64 |
| dim int |
| } |
| |
| // NewStudentsT creates a new StudentsT with the given nu, mu, and sigma |
| // parameters. |
| // |
| // NewStudentsT panics if len(mu) == 0, or if len(mu) != sigma.Symmetric(). If |
| // the covariance matrix is not positive-definite, nil is returned and ok is false. |
| func NewStudentsT(mu []float64, sigma mat.Symmetric, nu float64, src rand.Source) (dist *StudentsT, ok bool) { |
| if len(mu) == 0 { |
| panic(badZeroDimension) |
| } |
| dim := sigma.Symmetric() |
| if dim != len(mu) { |
| panic(badSizeMismatch) |
| } |
| |
| s := &StudentsT{ |
| nu: nu, |
| mu: make([]float64, dim), |
| dim: dim, |
| src: src, |
| } |
| if src != nil { |
| s.rnd = rand.New(src) |
| } |
| copy(s.mu, mu) |
| |
| ok = s.chol.Factorize(sigma) |
| if !ok { |
| return nil, false |
| } |
| s.sigma = *mat.NewSymDense(dim, nil) |
| s.sigma.CopySym(sigma) |
| s.chol.LTo(&s.lower) |
| s.logSqrtDet = 0.5 * s.chol.LogDet() |
| return s, true |
| } |
| |
| // ConditionStudentsT returns the Student's T distribution that is the receiver |
| // conditioned on the input evidence, and the success of the operation. |
| // The returned Student's T has dimension |
| // n - len(observed), where n is the dimension of the original receiver. |
| // The dimension order is preserved during conditioning, so if the value |
| // of dimension 1 is observed, the returned normal represents dimensions {0, 2, ...} |
| // of the original Student's T distribution. |
| // |
| // ok indicates whether there was a failure during the update. If ok is false |
| // the operation failed and dist is not usable. |
| // Mathematically this is impossible, but can occur with finite precision arithmetic. |
| func (s *StudentsT) ConditionStudentsT(observed []int, values []float64, src rand.Source) (dist *StudentsT, ok bool) { |
| if len(observed) == 0 { |
| panic("studentst: no observed value") |
| } |
| if len(observed) != len(values) { |
| panic(badInputLength) |
| } |
| |
| for _, v := range observed { |
| if v < 0 || v >= s.dim { |
| panic("studentst: observed value out of bounds") |
| } |
| } |
| |
| newNu, newMean, newSigma := studentsTConditional(observed, values, s.nu, s.mu, &s.sigma) |
| if newMean == nil { |
| return nil, false |
| } |
| |
| return NewStudentsT(newMean, newSigma, newNu, src) |
| |
| } |
| |
| // studentsTConditional updates a Student's T distribution based on the observed samples |
| // (see documentation for the public function). The Gaussian conditional update |
| // is treated as a special case when nu == math.Inf(1). |
| func studentsTConditional(observed []int, values []float64, nu float64, mu []float64, sigma mat.Symmetric) (newNu float64, newMean []float64, newSigma *mat.SymDense) { |
| dim := len(mu) |
| ob := len(observed) |
| |
| unobserved := findUnob(observed, dim) |
| |
| unob := len(unobserved) |
| if unob == 0 { |
| panic("stat: all dimensions observed") |
| } |
| |
| mu1 := make([]float64, unob) |
| for i, v := range unobserved { |
| mu1[i] = mu[v] |
| } |
| mu2 := make([]float64, ob) // really v - mu2 |
| for i, v := range observed { |
| mu2[i] = values[i] - mu[v] |
| } |
| |
| var sigma11, sigma22 mat.SymDense |
| sigma11.SubsetSym(sigma, unobserved) |
| sigma22.SubsetSym(sigma, observed) |
| |
| sigma21 := mat.NewDense(ob, unob, nil) |
| for i, r := range observed { |
| for j, c := range unobserved { |
| v := sigma.At(r, c) |
| sigma21.Set(i, j, v) |
| } |
| } |
| |
| var chol mat.Cholesky |
| ok := chol.Factorize(&sigma22) |
| if !ok { |
| return math.NaN(), nil, nil |
| } |
| |
| // Compute mu_1 + sigma_{2,1}ᵀ * sigma_{2,2}^-1 (v - mu_2). |
| v := mat.NewVecDense(ob, mu2) |
| var tmp, tmp2 mat.VecDense |
| err := chol.SolveVecTo(&tmp, v) |
| if err != nil { |
| return math.NaN(), nil, nil |
| } |
| tmp2.MulVec(sigma21.T(), &tmp) |
| |
| for i := range mu1 { |
| mu1[i] += tmp2.At(i, 0) |
| } |
| |
| // Compute tmp4 = sigma_{2,1}ᵀ * sigma_{2,2}^-1 * sigma_{2,1}. |
| // TODO(btracey): Should this be a method of SymDense? |
| var tmp3, tmp4 mat.Dense |
| err = chol.SolveTo(&tmp3, sigma21) |
| if err != nil { |
| return math.NaN(), nil, nil |
| } |
| tmp4.Mul(sigma21.T(), &tmp3) |
| |
| // Compute sigma_{1,1} - tmp4 |
| // TODO(btracey): If tmp4 can constructed with a method, then this can be |
| // replaced with SubSym. |
| for i := 0; i < len(unobserved); i++ { |
| for j := i; j < len(unobserved); j++ { |
| v := sigma11.At(i, j) |
| sigma11.SetSym(i, j, v-tmp4.At(i, j)) |
| } |
| } |
| |
| // The computed variables are accurate for a Normal. |
| if math.IsInf(nu, 1) { |
| return nu, mu1, &sigma11 |
| } |
| |
| // Compute beta = (v - mu_2)ᵀ * sigma_{2,2}^-1 * (v - mu_2)ᵀ |
| beta := mat.Dot(v, &tmp) |
| |
| // Scale the covariance matrix |
| sigma11.ScaleSym((nu+beta)/(nu+float64(ob)), &sigma11) |
| |
| return nu + float64(ob), mu1, &sigma11 |
| } |
| |
| // findUnob returns the unobserved variables (the complementary set to observed). |
| // findUnob panics if any value repeated in observed. |
| func findUnob(observed []int, dim int) (unobserved []int) { |
| var setOb intsets.Sparse |
| for _, v := range observed { |
| setOb.Insert(v) |
| } |
| var setAll intsets.Sparse |
| for i := 0; i < dim; i++ { |
| setAll.Insert(i) |
| } |
| var setUnob intsets.Sparse |
| setUnob.Difference(&setAll, &setOb) |
| unobserved = setUnob.AppendTo(nil) |
| sort.Ints(unobserved) |
| return unobserved |
| } |
| |
| // CovarianceMatrix calculates the covariance matrix of the distribution, |
| // storing the result in dst. Upon return, the value at element {i, j} of the |
| // covariance matrix is equal to the covariance of the i^th and j^th variables. |
| // covariance(i, j) = E[(x_i - E[x_i])(x_j - E[x_j])] |
| // If the dst matrix is empty it will be resized to the correct dimensions, |
| // otherwise dst must match the dimension of the receiver or CovarianceMatrix |
| // will panic. |
| func (st *StudentsT) CovarianceMatrix(dst *mat.SymDense) { |
| if dst.IsEmpty() { |
| *dst = *(dst.GrowSym(st.dim).(*mat.SymDense)) |
| } else if dst.Symmetric() != st.dim { |
| panic("studentst: input matrix size mismatch") |
| } |
| dst.CopySym(&st.sigma) |
| dst.ScaleSym(st.nu/(st.nu-2), dst) |
| } |
| |
| // Dim returns the dimension of the distribution. |
| func (s *StudentsT) Dim() int { |
| return s.dim |
| } |
| |
| // LogProb computes the log of the pdf of the point x. |
| func (s *StudentsT) LogProb(y []float64) float64 { |
| if len(y) != s.dim { |
| panic(badInputLength) |
| } |
| |
| nu := s.nu |
| n := float64(s.dim) |
| lg1, _ := math.Lgamma((nu + n) / 2) |
| lg2, _ := math.Lgamma(nu / 2) |
| |
| t1 := lg1 - lg2 - n/2*math.Log(nu*math.Pi) - s.logSqrtDet |
| |
| mahal := stat.Mahalanobis(mat.NewVecDense(len(y), y), mat.NewVecDense(len(s.mu), s.mu), &s.chol) |
| mahal *= mahal |
| return t1 - ((nu+n)/2)*math.Log(1+mahal/nu) |
| } |
| |
| // MarginalStudentsT returns the marginal distribution of the given input variables, |
| // and the success of the operation. |
| // That is, MarginalStudentsT returns |
| // p(x_i) = \int_{x_o} p(x_i | x_o) p(x_o) dx_o |
| // where x_i are the dimensions in the input, and x_o are the remaining dimensions. |
| // See https://en.wikipedia.org/wiki/Marginal_distribution for more information. |
| // |
| // The input src is passed to the created StudentsT. |
| // |
| // ok indicates whether there was a failure during the marginalization. If ok is false |
| // the operation failed and dist is not usable. |
| // Mathematically this is impossible, but can occur with finite precision arithmetic. |
| func (s *StudentsT) MarginalStudentsT(vars []int, src rand.Source) (dist *StudentsT, ok bool) { |
| newMean := make([]float64, len(vars)) |
| for i, v := range vars { |
| newMean[i] = s.mu[v] |
| } |
| var newSigma mat.SymDense |
| newSigma.SubsetSym(&s.sigma, vars) |
| return NewStudentsT(newMean, &newSigma, s.nu, src) |
| } |
| |
| // MarginalStudentsTSingle returns the marginal distribution of the given input variable. |
| // That is, MarginalStudentsTSingle returns |
| // p(x_i) = \int_{x_o} p(x_i | x_o) p(x_o) dx_o |
| // where i is the input index, and x_o are the remaining dimensions. |
| // See https://en.wikipedia.org/wiki/Marginal_distribution for more information. |
| // |
| // The input src is passed to the call to NewStudentsT. |
| func (s *StudentsT) MarginalStudentsTSingle(i int, src rand.Source) distuv.StudentsT { |
| return distuv.StudentsT{ |
| Mu: s.mu[i], |
| Sigma: math.Sqrt(s.sigma.At(i, i)), |
| Nu: s.nu, |
| Src: src, |
| } |
| } |
| |
| // TODO(btracey): Implement marginal single. Need to modify univariate StudentsT |
| // to be three-parameter. |
| |
| // Mean returns the mean of the probability distribution at x. If the |
| // input argument is nil, a new slice will be allocated, otherwise the result |
| // will be put in-place into the receiver. |
| func (s *StudentsT) Mean(x []float64) []float64 { |
| x = reuseAs(x, s.dim) |
| copy(x, s.mu) |
| return x |
| } |
| |
| // Nu returns the degrees of freedom parameter of the distribution. |
| func (s *StudentsT) Nu() float64 { |
| return s.nu |
| } |
| |
| // Prob computes the value of the probability density function at x. |
| func (s *StudentsT) Prob(y []float64) float64 { |
| return math.Exp(s.LogProb(y)) |
| } |
| |
| // 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. |
| func (s *StudentsT) Rand(x []float64) []float64 { |
| // If Y is distributed according to N(0,Sigma), and U is chi^2 with |
| // parameter ν, then |
| // X = mu + Y * sqrt(nu / U) |
| // X is distributed according to this distribution. |
| |
| // Generate Y. |
| x = reuseAs(x, s.dim) |
| tmp := make([]float64, s.dim) |
| if s.rnd == nil { |
| for i := range x { |
| tmp[i] = rand.NormFloat64() |
| } |
| } else { |
| for i := range x { |
| tmp[i] = s.rnd.NormFloat64() |
| } |
| } |
| xVec := mat.NewVecDense(s.dim, x) |
| tmpVec := mat.NewVecDense(s.dim, tmp) |
| xVec.MulVec(&s.lower, tmpVec) |
| |
| u := distuv.ChiSquared{K: s.nu, Src: s.src}.Rand() |
| floats.Scale(math.Sqrt(s.nu/u), x) |
| |
| floats.Add(x, s.mu) |
| return x |
| } |