There’s a subtle issue I’ve mentioned before involving `rand`

and `logpdf`

that has gotten in the way of composability. I’m thinking in terms of Soss and MeasureTheory.jl, but I thought I’d post here since it may be helpful to others.

First, the problem. There’s an implicit law in Distributions.jl that for any distribution `d`

,

```
logpdf(d, rand(d))
```

should be valid.

Now, suppose we have a Soss model

```
m = @model N begin
μ ~ Normal()
σ ~ HalfNormal()
x ~ Normal(μ,σ) |> iid(N)
return x
end
```

So `x = rand(m(N=10))`

would just give us a vector of 10 values, and evaluating `logpdf(m(N=10), x)`

would require integration (by MCMC or other methods).

This is sometimes what we want, but often it’s not. The problem is that the internal state is discarded.

In other contexts (StatsModels, Turing, etc) we have a `sample`

that’s similar in some ways to `rand`

, but might have some additional context included. Also, in MeasureTheory.jl I’ve been working in terms of `logdensity`

instead of `logpdf`

, because many densities won’t be *probability* densities.

So currently I’m thinking we could keep `logpdf(d, rand(d))`

as is and add a similar law for `logdensity(d, sample(d))`

, where the return value of `sample`

includes that of `rand`

in addition to a representation of the internal state. For example, here it could be a named tuple over `(μ,σ,x)`

. Then for PPL, we can work in terms of `sample`

except for the special case where we want to intentionally discard the internal state.

Any thoughts/concerns?