How do I compute a likelyhood from a model?

[zenon]

Hi,

say I have a model

m0 = @model X begin
    N = length(X)
    a ~ Beta(1, 1)
    X ~ Bernoulli(a) |> iid(N)
end

now, for m0(X = [false, false, false, true, true]), I want to compute the likelyhood with a = 0.5, is there a Soss way to do that?

(I see there is a method called likelihood, however, it returns a model, not a number.)

Thank you, z.

[millerjoey]

Hi! Your model is specified in a different style. You don't need to pass the data as an argument to the model; think of the argument as simply allowing your Model to define a family of distributions. Note that this is different from how models are written in Turing.

So to write your Beta-Bernoulli model, you can let your argument be the number of observations:

julia> m0 = @model N begin
       a ~ Beta(1,1)
       X ~ Bernoulli(a) |> iid(N)
       end;

julia> rand(m0(N=3))
(a = 0.626776184502901, X = Bool[1, 1, 1])

Then, the (log) likelihood of X=[1,1,0] and a=0.5 is

julia> logpdf(m0(N=3), (a=0.5, X = [1, 1, 0],))
-2.0794415416798357

[zenon]

Hello Joey,

thank you for the response. I still don't understand what's wrong with my first model. And the fact that

 logpdf(m0(X = [1, 1, 0]), (a=0.8,))

gives zero for ANY value of of a between 0 and 1 for that model doesn't help :-)

Can you give me a bit more about how to think about this? Pointer to what I should read is of course welcome!

Thank you and kind greetings, z.

PS.: Result of logpdf(m1(N=3), (a=0.5, X = [1, 1, 0],)) seems not to vary when I change only the value for N. Here I stop understanding again.

[cscherrer]

A model is a random generative process. An easy way to check a model is to see what code is generated for rand and logpdf. In this case,

julia> sourceRand(m0)
:(_rng->begin
          a = rand(_rng, Beta(1, 1))
          N = length(X)
          (N = N, a = a, X = X)
      end)

julia> sourceLogpdf(m0)
quote
    _ℓ = 0.0
    _ℓ += logpdf(Beta(1, 1), a)
    N = length(X)
    return _ℓ
end

Randomness of X is not taken into account. Because it's passed as a parameter, Soss considers it to be a fixed value. Soss is very different from Turing (and many PPLs) in this way. But by making models "function-like", we get some nice composability benefits.

Compare with @millerjoey 's suggestion (I'll call it m1):

julia> sourceRand(m1)
:(_rng->begin
          a = rand(_rng, Beta(1, 1))
          X = rand(_rng, iid(N, Bernoulli(a)))
          (a = a, X = X)
      end)

julia> sourceLogpdf(m1)
quote
    _ℓ = 0.0
    _ℓ += logpdf(Beta(1, 1), a)
    _ℓ += logpdf(Bernoulli(a) |> iid(N), X)
    return _ℓ
end

[zenon]

Thank you Chad, I will have a deeper look here.

I still don't understand why

 m2 = @model N begin
     a ~ Beta(1,1)
     X ~ Bernoulli(a) |> iid(N)
 end;

where sourceLogpdf(m2) gives

 quote
   _ℓ = 0.0
   _ℓ += logpdf(Beta(1, 1), a)
   _ℓ += logpdf(Bernoulli(a) |> iid(N), X)
   return _ℓ
 end

doesnt change the result of logpdf when I change N.

 logpdf(m2(N=3), (a=0.5, X = [1, 1, 0],))

and

 logpdf(m2(N=5), (a=0.5, X = [1, 1, 0],))

both give -2.0794415416798357

Kind greetings, z.

[cscherrer]

Ahh, good catch. That's because the dimension is currently only used for generation, when there's no value known. The code for this is in iid.jl:

function Distributions.logpdf(d::iid,x)
    s = zero(Float64)
    Δs(xj) = logpdf(d.dist, xj)

    @inbounds @simd for j in eachindex(x)
        s += Δs(x[j])
    end
    s
end