Gaussian Process Positive Definite Exception

[gcgibson]

I have written the following code for fitting a simple GP

using Mamba
## Model Specification
function kernel_with_noise(x,y,s2)
  dist = exp((-(x-y)^2))
  if x != y
    dist
  else
    1 + s2
  end
end

function kernel(x,y)
  dist = exp((-(x-y)^2))
end

line = Dict{Symbol, Any}(
  :x =>[1,2],
  :y => [1,2],
  :xtest => [1,2]
)
line[:xmat] = line[:x]
size_ = length(line[:x])

function posterior(xmat,xtest,y,s2)
  K = ones(length(xmat),length(xmat))
  K_ss = ones(length(xtest),length(xtest))
  K_s = ones(length(xtest),length(xmat))
  for i in 1:size(K)[1]
    for j in 1:size(K)[2]
      K[i,j] = kernel_with_noise(xmat[i],xmat[j],s2)
    end
  end
  for i in 1:size(K_ss)[1]
    for j in 1:size(K_ss)[2]
      K_ss[i,j] = kernel(xtest[i],xtest[j])
    end
  end
  for i in 1:size(K_s)[1]
    for j in 1:size(K_s)[2]
      K_s[i,j] = kernel(xtest[i],xmat[j])
    end
  end
  mu_star = K_s*inv(K)*y
  sigma_star = K_ss-transpose(K_s)*inv(K)*K_s
  mu_star,sigma_star
end

model = Model(

  y = Stochastic(1,
    (cov) ->  MvNormal(zeros(size_), cov),
    false
  ),
  ytilde = Logical(1,
    (s2,xtest,xmat) ->
    begin
      mu_star,sigma_star = posterior(xmat,xtest,line[:y],s2)
      d = MvNormal(mu_star,sigma_star)
      rand(d)
    end
  ),
  cov = Logical(2,
    (xmat,s2) ->
    begin
      x = eye(size_)
      for i in 1:size_
        for j in 1:size_
          tmp = kernel_with_noise(xmat[i],xmat[j],s2)
          x[i,j] = tmp
        end
      end
      x
    end
  ),
  s2 = Stochastic(
    () -> InverseGamma(.1, .1)
  )

)
scheme = [Slice([:s2],1.0)]

## Initial Values
inits = [
  Dict{Symbol, Any}(
    :y => line[:y],
    :s2 => rand(Gamma(.1, .1))
  )
for i in 1:1
]
setsamplers!(model, scheme)
sim1 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=1)
describe(sim1)

which throws a

Base.LinAlg.PosDefException(1)

while exploring the posterior for certain samples. Is there anyway to reject a sample based on this exception and keep going? I am not sure that is the correct way of handling bad s^2 values so if anyone has other suggestions for what to do please let me know!

[bdeonovic]

It may be easier to just write your own sampler rather than putting a lot of sampler logic in the Stochastic/Logical node definitions.

[bdeonovic]

For example see ## User-Defined Samplers in http://mambajl.readthedocs.io/en/latest/tutorial.html or the example samplers in http://mambajl.readthedocs.io/en/latest/examples/pollution.html

[gcgibson]

@bdeonovic That makes sense, but the custom samplers still require some return

Gibbs_sigma2 = Sampler([:sigma2],
  (mu, sigma2, y) ->
    begin
      a = length(y) / 2.0 + shape(sigma2.distr)
      b = sumabs2(y - mu) / 2.0 + scale(sigma2.distr)
      rand(InverseGamma(a, b))
    end
)

In this case a random draw from InvGamma. In the case of a non positive definite sample I would like to simple reject it from the chain, can I hack my sampler to do that?

EDIT: To clarify, if I wanted to put an if statement like so

Gibbs_sigma2 = Sampler([:sigma2],
  (mu, sigma2, y) ->
    begin
      a = length(y) / 2.0 + shape(sigma2.distr)
      b = sumabs2(y - mu) / 2.0 + scale(sigma2.distr)
      if matrix formed from some combination of a and b is positive semi definite:
           rand(InverseGamma(a, b))
      else
           not sure what to put here
    end
)

[bdeonovic]

You might want to check why non-positive definite matrices are being formed. Sometimes a matrix is not positive definite in julia due to round off errors. In those cases you can explicitly specify to julia that the matrix is Positive definite.

Otherwise, you can always return the current iterate if the proposed matrix is not positive definite.

[gcgibson]

How do you get access to the underlying chain to return the last accepted step, I guess you don't need it you can just return the input to your sampler?

[bdeonovic]

in the example you have

Gibbs_sigma2 = Sampler([:sigma2],
  (mu, sigma2, y) ->
    begin
      ##sigma2.value is the current value of sigma2 which you are trying to sample a new value for
      a = length(y) / 2.0 + shape(sigma2.distr)
      b = sumabs2(y - mu) / 2.0 + scale(sigma2.distr)
      if matrix formed from some combination of a and b is positive semi definite:
           rand(InverseGamma(a, b))
      else
           not sure what to put here
    end
)
``

[gcgibson]

Awesome! Thanks for your help! The matrices are PSD until s^2 gets sufficiently small

[gcgibson]

Unfortunately, I can't use Gibbs sampling because the conditional density is not available (unlike Bayesian linear regression , gaussian processes don't have an analytical distribution for s2 that I know of ). Is there anyway to get a code snippet like that into MH or NUTS?

[bdeonovic]

Your custom sampler doesn't have to be a Gibbs sampler. It can be whatever you want. For MH you would just need to calculate the the acceptance probability for example.

[gcgibson]

But then I would need a proposal dist etc right? I would basically have to re-implement MH in that code block? Gibbs seemed nice because it was just actually sampling from known dists, but I would prefer to pass a flag to the MH sampler saying treat this exception as probability 0, if possible.