Improper initial values when supports of the prior distribution are themselves random

[bdeonovic]

While trying out my first example in Turing I ran into an issue (see discussion here: https://discourse.julialang.org/t/help-with-first-non-trivial-turing-example/38964/14)

Apparently Turing does not provide proper initial values when supports of the prior distribution are themselves random. Below is an example in which the issue arises. Note that the example may have other issues as well:

# Import Turing and Distributions.
using Turing, Distributions

n = 1000
m = 10
k = 4
theta = randn(n)
b = zeros(k,m)
for i in 1:m
    b[1,i] = randn()
    for j in 2:k
        dd = truncated(Normal(), b[j-1,i], Inf)
        b[j,i] = rand(dd)
    end
end

logit = x -> log(x / (1 - x))
invlogit = x -> exp(x)/(1 + exp(x))
y = zeros(m,n)
probs = zeros(k,m,n)
for p in 1:n
    for i in 1:m
        probs[1,i,p] = 1.0
        for j in 1:(k-1)
            Q = invlogit(theta[p] - b[j,i])
            probs[j,i,p] -= Q
            probs[j+1,i,p] = Q
        end
        y[i,p] = rand(Categorical(probs[:,i,p]))
    end
end



# Graded Response Model
@model function grm(y, n, m, k, ::Type{TC}=Array{Float64,3}, ::Type{TM}=Array{Float64,2}, ::Type{TV}=Vector{Float64}) where {TC, TM, TV}
    b = TM(undef, k, m)
    for i in 1:m
        b[1,i] ~ Normal(0,1)
        for j in 2:k
            b[j,i] ~ truncated(Normal(0,1), b[j-1,i], Inf) 
        end
    end
    probs = TC(undef, k, m, n)
    theta = TV(undef, n)
    for p in 1:n
        theta[p] ~ Normal(0,1)
        for i in 1:m
            probs[1,i,p] = 1.0
            for j in 1:(k-1)
                Q = invlogit(theta[p] - b[j,i])
                probs[j,i,p] -= Q
                probs[j+1,i,p] = Q
            end
            y[i,p] ~ Categorical(probs[:,i,p])
        end
    end
    return theta, b
end;

result = sample(grm(y, n, m, k), MH(), 1)

With error

Sampling: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:01
ERROR: ArgumentError: Categorical: the condition isprobvec(p) is not satisfied.

Here probs[:i,p] should be guarenteed to be a probability vector as long as the b[:,i] are in ascending order. If we put in some debugging print statements we find that at times the b[:,i] are not in ascending order.

[mohdibntarek]

Thanks @bdeonovic for opening this. Fixing this issue requires a non-trivial redesign of parts of DynamicPPL. I am on it though.

[mohdibntarek]

Btw this particular case can be worked around using a custom distribution and a custom bijector:

using Distributions, Bijectors

struct MyDistribution <: ContinuousMultivariateDistribution
    n::Int
end

function Base.rand(d::MyDistribution)
    b = zeros(d.n)
    b[1] = rand(Normal(0,1))
        for i in 2:d.n
        b[i] = rand(truncated(Normal(0,1), b[i-1], Inf))
    end
    return b
end
function Distributions.logpdf(d::MyDistribution, b::AbstractVector)
    l = logpdf(Normal(0,1), b[1])
    for i in 2:d.n
        l += logpdf(truncated(Normal(0,1), b[i-1], Inf), b[i])
    end
    return l
end

struct MyBijector <: Bijectors.Bijector{1} end
function (b::MyBijector)(x::AbstractVector)
    y = similar(x)
    y[1] = x[1]
    for i in 2:length(x)
        y[i] = bijector(truncated(Normal(0,1), x[i-1], Inf))(x[i])
    end
    return y
end
function (b::Inverse{<:MyBijector})(y::AbstractVector)
    x = similar(y)
    x[1] = y[1]
    for i in 2:length(y)
        x[i] = inv(bijector(truncated(Normal(0,1), x[i-1], Inf)))(y[i])
    end
    return x
end
function Bijectors.logabsdetjac(b::MyBijector, x::AbstractVector)
    l = float(zero(eltype(x)))
    for i in 2:length(x)
        l += logabsdetjac(bijector(truncated(Normal(0,1), x[i-1], Inf)), x[i])
    end
    return l
end
Bijectors.bijector(d::MyDistribution) = MyBijector()

[bdeonovic]

Thanks @mohamed82008, could you point me towards some material for how I can now use the defined custom distribution and bijector in order to sample from the posterior distribution?

[mohdibntarek]

Once you have the above definitions you can use MyDistribution as a prior in the model like any other distribution. So your model becomes:

# Graded Response Model
@model function grm(y, n, m, k, ::Type{TC}=Array{Float64,3}, ::Type{TM}=Array{Float64,2}, ::Type{TV}=Vector{Float64}) where {TC, TM, TV}
    b = TM(undef, k, m)
    for i in 1:m
        b[1,:] ~ MyDistribution(k)
    end
    probs = TC(undef, k, m, n)
    theta = TV(undef, n)
    for p in 1:n
        theta[p] ~ Normal(0,1)
        for i in 1:m
            probs[1,i,p] = 1.0
            for j in 1:(k-1)
                Q = invlogit(theta[p] - b[j,i])
                probs[j,i,p] -= Q
                probs[j+1,i,p] = Q
            end
            y[i,p] ~ Categorical(probs[:,i,p])
        end
    end
    return theta, b
end;

[bdeonovic]

Currently getting this error:

julia> result = sample(grm(y, n, m, k), MH(), 10)
ERROR: MethodError: no method matching reconstruct(::Array{MyDistribution,1}, ::Array{Float64,1})

If I try another sampler (like HMC) it doesn't budge from the initial values:

julia> result = sample(grm(y, n, m, k), HMC(0.2, 3), 10)
┌ Warning: The current proposal will be rejected due to numerical error(s).
│   isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC C:\Users\Benjamin\.julia\packages\AdvancedHMC\WJCQA\src\hamiltonian.jl:47

[mohdibntarek]

This works for me:

using Turing
using Distributions, Bijectors

struct MyDistribution <: ContinuousMultivariateDistribution
    n::Int
end

function Base.rand(d::MyDistribution)
    b = zeros(d.n)
    b[1] = rand(Normal(0,1))
        for i in 2:d.n
        b[i] = rand(truncated(Normal(0,1), b[i-1], Inf))
    end
    return b
end
function Distributions.logpdf(d::MyDistribution, b::AbstractVector)
    l = logpdf(Normal(0,1), b[1])
    for i in 2:d.n
        l += logpdf(truncated(Normal(0,1), b[i-1], Inf), b[i])
    end
    return l
end

struct MyBijector <: Bijectors.Bijector{1} end
function (b::MyBijector)(x::AbstractVector)
    y = similar(x)
    y[1] = x[1]
    for i in 2:length(x)
        y[i] = bijector(truncated(Normal(0,1), x[i-1], Inf))(x[i])
    end
    return y
end
function (b::Inverse{<:MyBijector})(y::AbstractVector)
    x = similar(y)
    x[1] = y[1]
    for i in 2:length(y)
        x[i] = inv(bijector(truncated(Normal(0,1), x[i-1], Inf)))(y[i])
    end
    return x
end
function Bijectors.logabsdetjac(b::MyBijector, x::AbstractVector)
    l = float(zero(eltype(x)))
    for i in 2:length(x)
        l += logabsdetjac(bijector(truncated(Normal(0,1), x[i-1], Inf)), x[i])
    end
    return l
end
Bijectors.bijector(d::MyDistribution) = MyBijector()

n = 10
m = 10
k = 4
theta = randn(n)
b = zeros(k,m)
for i in 1:m
    b[1,i] = randn()
    for j in 2:k
        dd = truncated(Normal(), b[j-1,i], Inf)
        b[j,i] = rand(dd)
    end
end

logit = x -> log(x / (1 - x))
invlogit = x -> exp(x)/(1 + exp(x))
y = zeros(m,n)
probs = zeros(k,m,n)
for p in 1:n
    for i in 1:m
        probs[1,i,p] = 1.0
        for j in 1:(k-1)
            Q = invlogit(theta[p] - b[j,i])
            probs[j,i,p] -= Q
            probs[j+1,i,p] = Q
        end
        y[i,p] = rand(Categorical(probs[:,i,p]))
    end
end

# Graded Response Model
@model function grm(y, n, m, k, ::Type{TC}=Array{Float64,3}, ::Type{TM}=Array{Float64,2}, ::Type{TV}=Vector{Float64}) where {TC, TM, TV}
    b = TM(undef, k, m)
    for i in 1:m
        b[:,i] ~ MyDistribution(k)
    end
    probs = TC(undef, k, m, n)
    theta = TV(undef, n)
    for p in 1:n
        theta[p] ~ Normal(0,1)
        for i in 1:m
            probs[1,i,p] = 1.0
            for j in 1:(k-1)
                Q = invlogit(theta[p] - b[j,i])
                probs[j,i,p] -= Q
                probs[j+1,i,p] = Q
            end
            y[i,p] ~ Categorical(probs[:,i,p])
        end
    end
    return theta, b
end;

result = sample(grm(y, n, m, k), HMC(0.01, 4), 100)

[bdeonovic]

Thanks a lot @mohamed82008! I've simplified the code a bit (and allowed the categories to be variable for each item. The following code is a bit shorter. Only problem is that it is a bit slow to sample. Any suggestions?

# Import Turing and Distributions.
using Turing, Distributions, Bijectors

struct MyDistribution <: ContinuousMultivariateDistribution
    n::Int
end

function Distributions.length(d::MyDistribution)
    d.n
end

function Base.rand(d::MyDistribution)
	b = zeros(d.n)
    b[1] = rand(Normal(0,1))
	for i in 2:d.n
        b[i] = rand(truncated(Normal(0,1), b[i-1], Inf))
    end
    return b
end
function Distributions._logpdf(d::MyDistribution, b::AbstractVector)
    l = logpdf(Normal(0,1), b[1])
	for i in 2:d.n
        l += logpdf(truncated(Normal(0,1), b[i-1], Inf), b[i])
    end
    return l
end

struct MyBijector <: Bijectors.Bijector{1} end
function (b::MyBijector)(x::AbstractVector)
	y = similar(x)
	y[1] = x[1]
	for i in 2:length(x)
		y[i] = bijector(truncated(Normal(0,1), x[i-1], Inf))(x[i])
	end
	return y
end
function (b::Inverse{<:MyBijector})(y::AbstractVector)
	x = similar(y)
	x[1] = y[1]
	for i in 2:length(y)
		x[i] = inv(bijector(truncated(Normal(0,1), x[i-1], Inf)))(y[i])
	end
	return x
end
function Bijectors.logabsdetjac(b::MyBijector, x::AbstractVector)
	l = float(zero(eltype(x)))
	for i in 2:length(x)
		l += logabsdetjac(bijector(truncated(Normal(0,1), x[i-1], Inf)), x[i])
	end
	return l
end
Bijectors.bijector(d::MyDistribution)= MyBijector()

n = 1000
m = 3
ks = [3,4,5]
theta = randn(n)
b = Array{Float64,1}.(undef, ks .- 1)
for i in 1:m
    b[i] .= rand(MyDistribution(ks[i]-1))
end

y = zeros(m,n)
for p in 1:n
    for i in 1:m
       y[i,p] = rand(OrderedLogistic(theta[p], b[i]))
    end
end

# Graded Response Model
@model function grm(y, n, m, ks, ::Type{TV}=Vector{Float64}) where {TV}
    b = TV.(undef, ks .- 1)
    for i in 1:m
      b[i] ~ MyDistribution(ks[i]-1)
    end
    theta = TV(undef, n)
    for p in 1:n
        theta[p] ~ Normal(0,1)
        for i in 1:m
            y[i,p] ~ OrderedLogistic(theta[p], b[i])
        end
    end
    return theta, b
end;

result = sample(grm(y, n, m, ks), HMC(0.01,4), 10);

[mohdibntarek]

Try:

using ReverseDiff, Memoization
Turing.setadbackend(:reversediff)
Turing.setrdcache(true)

before the sample line.

[mohdibntarek]

Hmm actually the caching may not work well with OrderedLogistic because it has branches inside. So ignore the Turing.setrdcache(true) and if you set it true already, set it back to false. See the correction below.

[mohdibntarek]

Alternatively, we can define a ReverseDiff custom adjoint for the logpdf of OrderedLogistic and then caching should work.

[mohdibntarek]

Hmm on third thought, caching should work fine because the branch taken by the logpdf function for each loop iteration is the same each time the function is called because y is fixed. So ignore my previous comment :)

[mohdibntarek]

You may need to call Turing.emptyrdcache() when changing y though.

[bdeonovic]

I just tried it....I tabbed over to read about something else and when I looked back the sampling was already done! I'm rerunning it out of disbelief.

[mohdibntarek]

Hahaha that's ReverseDiff for you! Check the results against other backends for a small data case before committing fully.

[mohdibntarek]

Always good to run some sanity tests :)