Specifying model with threshold / ifelse likelihood in Turing / Soss / Tilde

Specific Domains Probabilistic programming
, ,
1

Hi everyone,

I have a model which is a simplified version of the one presented in Mulder & Hoff (2021) and looks something like this:

\begin{align} \mathbf{z}_{i} &\sim \textrm{MvNormal}(\boldsymbol{\theta}_{i},\mathbf{I}) \\[0.8ex] \boldsymbol{\theta}_{ir} &= \mathbf{x}_{ir}^\top \boldsymbol{\beta} \\[0.8ex] c_{i} \mid \mathbf{z}_{i} &\sim \textrm{truncate} \Big( \textrm{Normal}(0, \sigma_c), -\infty,\max(\textbf{z}_{i}) \Big) \\[0.8ex] y_{ir} &= \begin{cases} 1 &\textrm{if } z_{ir} > c_{i} \\ 0 & \textrm{otherwise} \end{cases} \end{align}

In the paper, they discuss a Gibbs sampler implementation but I would like to specify this in a Julia PPL. Hereā€™s how far I got in Turing (probably very non-optimized):

@model function MRREM(X, y)
	N, A, K = size(X)
	
	Ī² ~ filldist(Normal(0, 1), K)
	Ļƒc ~ Exponential()
	
	Īø = [@views X[i,:,:] * Ī² for i in 1:N]

	z ~ arraydist([MvNormal(t, I) for t in Īø])
	zmax = maximum(z, dims=1) |> vec
	c ~ arraydist([Truncated(Normal(0, Ļƒc), -Inf, zm) for zm in zmax])
	# y ~ ???
end

let 
	N, A, k = 100, 20, 3
	X = rand(N, A, k)
	MRREM(X, missing) |> rand
end

But I donā€™t really know how to specify the ifelse for the y_{ir} in a PPL or how to transform this into something that you can sample from / evaluate the logdensity of. Does someone know what to do here?

1 Like
2

Hi @Jakob ,

Itā€™s pretty typical for PPLs to have difficulty ā€œobservingā€ deterministic functions of a distribution. I think the standard way to do this would be to integrate out c, so each y[i] follows a Bernoulli distribution as a function of Ļƒc and zmax[i].

You could write this in the model itself, or (in Soss or Tilde) build a new measure to express this, then use that measure in the model. The advantage of this is that you can test the measure independently and do any performance doing in isolation.

2 Likes
3

Thanks, @cscherrer. That sounded like more math than Iā€™m comfortable with but luckily I found that further down in the paper they give a specification of the probability of the y_{ir}'s as:

Pr(y_{ir}=1 |\mathbf{z_i},\sigma_c)=\frac{\Phi\big(\frac{z_{ir}}{\sigma_c}\big)}{\Phi \Big(\frac{\max(\textbf{z}_{i})}{\sigma_c}\Big)}

So I guess I can just work in terms of that.

Btw., I remember seeing something about truncated measures in MeasureTheory but couldnā€™t find it anymore. Is there something I overlooked or is there anything planned?

1 Like
4

Great! Yes, without working through it explicitly thatā€™s how Iā€™d expect the result to look.

Thereā€™s a start on it here:

Itā€™s not hard in principle, just a matter of working through some corner cases.

1 Like
5

Ah great, Iā€™m looking forward to giving the model a go in Soss / Tilde when that lands.

6

Please post the code, if you get something working.