Evaluating the posterior

[maximerischard]

What is the best way to evaluate the posterior for a set of parameters? For example, with

@model gdemo(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt(s))
  x[1] ~ Normal(m, sqrt(s))
  x[2] ~ Normal(m, sqrt(s))
  return s, m
end
mf = gdemo([1.5, 2.0])

how do I evaluate the unnormalised posterior, say, at s=1.2 m=0.5?

[emilemathieu]

Hi, the joint is proportional to the posterior. You can compute it the following way

@model gdemo(x,s,m) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt(s))
  x[1] ~ Normal(m, sqrt(s))
  x[2] ~ Normal(m, sqrt(s))
  return s, m
end

model = gdemo(x,1.2,0.5)
res = model()
println(res.logp)

[maximerischard]

Wonderful, thank you. Just to make sure I understand, does this create a model where x, s, and m are all assumed to be observed, right? So in order to sample from this model and evaluate the posterior, I need to create two separate models?

[xukai92]

Wonderful, thank you. Just to make sure I understand, does this create a model where x, s, and m are all assumed to be observed, right?

Yes

So in order to sample from this model and evaluate the posterior, I need to create two separate models?

If you just want to sample from the model (by some samplers) and check the posterior, in stead of evaluating the posterior of some given s and m, you only need the first model.

@model gdemo(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt(s))
  x[1] ~ Normal(m, sqrt(s))
  x[2] ~ Normal(m, sqrt(s))
  return s, m
end
mf = gdemo([1.5, 2.0])
chain = sample(mf, HMC(100, 0.2, 3)) # sample this model using HMC for 100 sampler with 0.2 leapfrog stepsize for 3 steps in each iteration
print(chain[:s])
print(chain[:m])
print(chain[:lp]) # log posterior

[maximerischard]

Thank you again, I'm starting to understand how this all works. My only problem now is that evaluating the joint using @emilemathieu's code is 77,000 times slower than hand-coding it:

using Turing
using BenchmarkTools
Turing.set_verbosity(0)

@model gdemo(x,s,m) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt(s))
  x[1] ~ Normal(m, sqrt(s))
  x[2] ~ Normal(m, sqrt(s))
  return s, m
end


function gdemo_logp(s, m)
    x = [1.5, 2.0]
    mf = gdemo(x, s, m)
    res = Base.invokelatest(mf)
    return res.logp[1]
end

function byhand_logp(s, m)
    x = [1.5, 2.0]
    sprior = InverseGamma(2,3)
    mprior = Normal(0, √s)
    likelihood = Normal(m, √s)
    lp = logpdf(sprior, s) + logpdf(mprior, m) + logpdf(likelihood, x[1]) + logpdf(likelihood, x[2])
    return lp
end
;

gdemo_logp(1.2, 0.5)

-5.338371361183928

byhand_logp(1.2, 0.5)

-5.338371361183928

@benchmark byhand_logp(1.2, 0.5)

BenchmarkTools.Trial: 
  memory estimate:  96 bytes
  allocs estimate:  1
  --------------
  minimum time:     128.658 ns (0.00% GC)
  median time:      135.709 ns (0.00% GC)
  mean time:        159.941 ns (7.77% GC)
  maximum time:     4.977 μs (94.38% GC)
  --------------
  samples:          10000
  evals/sample:     875

@benchmark gdemo_logp(1.2, 0.5)

BenchmarkTools.Trial: 
  memory estimate:  197.50 KiB
  allocs estimate:  3362
  --------------
  minimum time:     10.555 ms (0.00% GC)
  median time:      12.032 ms (0.00% GC)
  mean time:        12.424 ms (0.52% GC)
  maximum time:     25.212 ms (41.02% GC)
  --------------
  samples:          402
  evals/sample:     1

Any ideas on how to evaluate the joint in a reasonable amount of time?

[maximerischard]

I was able to get it down to “only” 100x slower than hard-coded, by manipulating the VarInfo object directly. It's not pretty, and still quite slow, so I'd welcome any suggestions.

import Turing: VarName, sym, setval!, vns, setlogp!

@model gdemox(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt(s))
  x[1] ~ Normal(m, sqrt(s))
  x[2] ~ Normal(m, sqrt(s))
  return s, m
end
mf = gdemox([1.5, 2.0])
vi = mf()

function gdemox_logp(s, m)
    names = Turing.vns(vi)
    namedict = Dict{Symbol,VarName}()
    setlogp!(vi, 0.0)
    for n in names
        namedict[sym(n)] = n
    end
    logp_before = vi.logp[1]
    setval!(vi, s, namedict[:s])
    setval!(vi, m, namedict[:m])
    mf(vi)
    return vi.logp[1]
end

gdemox_logp(1.2, 0.5)

-5.338371361183928

@benchmark gdemox_logp(1.2, 0.5)

BenchmarkTools.Trial: 
  memory estimate:  3.09 KiB
  allocs estimate:  76
  --------------
  minimum time:     12.489 μs (0.00% GC)
  median time:      14.279 μs (0.00% GC)
  mean time:        16.170 μs (2.52% GC)
  maximum time:     4.193 ms (97.26% GC)
  --------------
  samples:          10000
  evals/sample:     1

[xukai92]

This is a good question. They main reason the way you evaluate posterior is slow is because the function gdemo_logp also includes the compilation time. mf = gdemo(x, s, m) means compile the model with given observations into mf.

There is another way to do the evaluation, which takes the use of Base.invokelatest(model, vi, spl), where model is the compiled model with only observations, vi is an internal data structure we use to store model variables and spl is the sampler which can be nothing.

This can be done by:

@model gdemo2(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt(s))
  x[1] ~ Normal(m, sqrt(s))
  x[2] ~ Normal(m, sqrt(s))
  return s, m
end

mf = gdemo2(x)

# Generate a VarInfo by sampling from pior
vi = Base.invokelatest(mf)
# NOTE: you can try print(vi) to see what's inside...

function gdemo_logp2(s, m)
    vi.logp[1] = 0         # set log posterior to 0
    vi.vals[1][1] = s      # set s as you want
    vi.trans[1][1] = false # set the value of s as non-transformed
    vi.vals[1][2] = m      # set m as you want
    
    # Evalute using vi
    res = Base.invokelatest(mf, vi, nothing)
    return res.logp[1]
end

NOTE: on my machine the mean time for this is 12.998 μs which is still higher than they by_hand one on my machine which is 233.460 ns.

However as you can see we don't provide an easy-to-use interface to build vi yet as it is designed to be use internally. Also in the use case we considered, as the posterior is only evaluated internally, the compilation is only done once so it wasn't a big problem.

[xukai92]

Well you just figured it our when I was typing...

[xukai92]

100x seems a big overhead here. @yebai Do you have any idea on why and how to make it faster? Seems quite significant here.

[xukai92]

I found this is interesting. Do you know why?

function gdemo_logp(s, m)
    vi.logp = 0
    vi.logp += logpdf(InverseGamma(2,3), s) 
    vi.logp += logpdf(Normal(0, √s), m) 
    vi.logp += logpdf(Normal(m, sqrt(s)), x[1])
    vi.logp += logpdf(Normal(m, sqrt(s)), x[2])
    vi.logp
end

BenchmarkTools.Trial: 
  memory estimate:  288 bytes
  allocs estimate:  16
  --------------
  minimum time:     11.752 μs (0.00% GC)
  median time:      13.005 μs (0.00% GC)
  mean time:        13.046 μs (0.00% GC)
  maximum time:     59.533 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

But

function gdemo_logp(s, m)
    logp = 0
    logp += logpdf(InverseGamma(2,3), s) 
    logp += logpdf(Normal(0, √s), m) 
    logp += logpdf(Normal(m, sqrt(s)), x[1])
    logp += logpdf(Normal(m, sqrt(s)), x[2])
    logp
end

BenchmarkTools.Trial: 
  memory estimate:  192 bytes
  allocs estimate:  10
  --------------
  minimum time:     251.597 ns (0.00% GC)
  median time:      280.504 ns (0.00% GC)
  mean time:        310.580 ns (9.32% GC)
  maximum time:     10.900 μs (96.22% GC)
  --------------
  samples:          10000
  evals/sample:     365

It seems that the slowness is due to the use of logp filed of vi??

[xukai92]

@ChrisRackauckas Do you know why the performances in my comment above differs a lot? Is using a filed of a composite type very expensive?

[maximerischard]

That seems likely to be due to type instability. They’re both type unstable since you’re initialising to an integer (0) rather than a float. Try using the '@code_warntype' to see what the problem might be.

[xukai92]

Well I can change the 2nd one to a generic type but the 2nd actually much faster.

[ChrisRackauckas]

How is the first one parameterized? If the type is parameterized or set concrete on that field, then vi.logp = 0 does the conversion when setting so it's actually type-stable. My guess is that field was left to ride as Any. And yes the second code isn't type-stable but the + is such a small amount of the time I don't know if it makes a difference, but do make it logp = 0.0.

[xukai92]

Thanks Maxime and Chris, the results above are come from type instability issue.