# SampleChains

**URL:** <https://discourse.julialang.org/t/samplechains/55937>\
**Category:** Probabilistic Programming\
**Created:** [February 24, 2021, 3:24pm UTC](https://discourse.julialang.org/t/samplechains/55937 "2021-02-24T15:24:32Z")\
**Posts on this page:** 8\
**Page:** 2

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [February 24, 2021, 11:42pm UTC](https://discourse.julialang.org/t/samplechains/55937/21 "2021-02-24T23:42:08Z")

</div>

Also I think it’s very useful to be able to solve problems like “I want a estimate of \mathbb{E}[f(X)] with a maximum standard error of \varepsilon.” If samples are expensive enough, we’d want to track the standard error and stop sampling when we’re close enough.

If a chain knows how to get more samples and we have, say, an `OnlineStat` for the standard error, this seems straightforward.

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [February 25, 2021, 2:50pm UTC](https://discourse.julialang.org/t/samplechains/55937/22 "2021-02-25T14:50:02Z")

</div>

I’m also just realizing AbstractMCMC probably won’t handle weighted observations, which we’ll need for VI. I think it’s best for AbstractMCMC to only manage the sampling logic, produce an iterator, allow this to be called as needed, and pass one sample at a time to whatever package is managing the chains.

---

<div class="post-metadata">

**Author:** ![cpfiffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cpfiffer/32/208747_2.png) [@cpfiffer](https://discourse.julialang.org/u/cpfiffer)\
**Post date:** [February 25, 2021, 4:49pm UTC](https://discourse.julialang.org/t/samplechains/55937/23 "2021-02-25T16:49:23Z")

</div>

> [@cscherrer](#):
>
> But how does it know where the sampler is? In Soss the `model` will get you the `logdensity` very quickly, but other things would need to be recomputed. If it’s all stored in the `chains` , that seems equivalent to what I’m doing now.

That state tells you where the sampler is?

> [@cscherrer](#):
>
> Also I think it’s very useful to be able to solve problems like “I want a estimate of E[f(X)]\mathbb{E}[f(X)] with a maximum standard error of ε\varepsilon .” If samples are expensive enough, we’d want to track the standard error and stop sampling when we’re close enough.
> 
> If a chain knows how to get more samples and we have, say, an `OnlineStat` for the standard error, this seems straightforward.

Again this is just just a difference in opinion on which thing does what. You want your chain to be more than a data structure and to have knowledge of the outside world – how samples work, how to get them, etc. Our philiosophy has been that it is just a chain. It holds the samples and maybe does some analysis. It should do that well and have _no idea how to get more samples_. It seems seems really bizarre to me to add that kind of logic to it.

> [@cscherrer](#):
>
> I’m also just realizing AbstractMCMC probably won’t handle weighted observations, which we’ll need for VI. I think it’s best for AbstractMCMC to only manage the sampling logic, produce an iterator, allow this to be called as needed, and pass one sample at a time to whatever package is managing the chains.

That is just not true – it is exceptionally general. You can do absolutely whatever you want within the framework of AbstractMCMC. If you want to weight, build that into your inference algorithm.

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [February 25, 2021, 6:30pm UTC](https://discourse.julialang.org/t/samplechains/55937/24 "2021-02-25T18:30:51Z")

</div>

> [@cpfiffer](#):
>
> That state tells you where the sampler is?

Ok, I thought by “state” you meant the state for an iterator. It sounds like you mean something a lot more general, like almost everything you’d have in an iterator, but you should have to jump through some more hoops to get to an actual iterator. I’m just having a lot of trouble seeing any advantage to this.

> [@cpfiffer](#):
>
> Our philiosophy has been that it is just a chain. It holds the samples and maybe does some analysis.

If it holds a state that can tell you where to find a sampler, that’s a lot more than you’re describing here.

> [@cpfiffer](#):
>
> It should do that well and have _no idea how to get more samples_ . It seems seems really bizarre to me to add that kind of logic to it.

Again, it doesn’t have that logic. It gets it from other places, and just holds it for the convenience of the end user.

I understand this is a difference of opinion. But there are strong advantages to the approach I’m suggesting. I hope I’ve made them clear; please let me know if I haven’t. It would be really helpful to hear any disadvantages, so we can get beyond our current opinions. I’m willing to rethink this if need be, I just haven’t yet heard a good argument against it.

> [@cpfiffer](#):
>
> You can do absolutely whatever you want within the framework of AbstractMCMC. If you want to weight, build that into your inference algorithm.

Thanks, that’s really interesting. I’ll start taking a closer look and try to get a better feel for how things can connect.

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [February 25, 2021, 6:41pm UTC](https://discourse.julialang.org/t/samplechains/55937/25 "2021-02-25T18:41:40Z")

</div>

I do think it’s critical to keep things very simple for end-users. So getting samples and adding some more should be no harder than this:

```julia
julia> using SampleChainsDynamicHMC, TransformVariables

julia> function ℓ(nt)
               z = nt.x/nt.σ
               return -z^2 - nt.σ - log(nt.σ)
       end
ℓ (generic function with 1 method)

julia> t = as((x=asℝ, σ=asℝ₊))
TransformVariables.TransformTuple{NamedTuple{(:x, :σ), Tuple{TransformVariables.Identity, TransformVariables.ShiftedExp{true, Float64}}}}((x = asℝ, σ = asℝ₊), 2)

julia> chain = initialize!(DynamicHMCChain, ℓ, t)
1-element Chain with schema (x = Float64, σ = Float64)
(x = -0.17±0.0, σ = 0.43±0.0)

julia> drawsamples!(chain, 9)
10-element Chain with schema (x = Float64, σ = Float64)
(x = 0.05±0.24, σ = 0.512±0.3)

julia> drawsamples!(chain, 90)
100-element Chain with schema (x = Float64, σ = Float64)
(x = -0.01±0.32, σ = 0.588±0.37)

```

Or with multiple chains:

```julia
julia> chains = initialize!(4, DynamicHMCChain, ℓ, t)
4-element MultiChain with 4 chains and schema (x = Float64, σ = Float64)
(x = 0.02±0.4, σ = 0.83±0.6)

julia> drawsamples!(chains, 9)
40-element MultiChain with 4 chains and schema (x = Float64, σ = Float64)
(x = -0.04±0.5, σ = 0.82±0.54)

julia> drawsamples!(chains, 90)
400-element MultiChain with 4 chains and schema (x = Float64, σ = Float64)
(x = -0.14±0.62, σ = 0.92±0.75)

```

---

<div class="post-metadata">

**Author:** ![cpfiffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cpfiffer/32/208747_2.png) [@cpfiffer](https://discourse.julialang.org/u/cpfiffer)\
**Post date:** [February 25, 2021, 6:50pm UTC](https://discourse.julialang.org/t/samplechains/55937/26 "2021-02-25T18:50:18Z")

</div>

Perhaps I understand what you’re doing finally. We have exactly the same logic as you:

```julia
chain = sample(model, sampler, 1000)

chain = resume(chain, 1000)

```

But in this case `chain` just has an empty bucket to store the model and the state. The sampler and the model decide what to put in the empty bucket. It’s equivalent to

```julia
chain = resume(model, sampler, chain, 1000)

```

None of this seems any more complicated for the end-user than what you are doing. More importantly, none of this is precluded by AbstractMCMC or anything – the only thing that happens here is that the model and the sampler take what they need from the chain in order to resume, rather than storing an iterator directly.

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [February 25, 2021, 7:15pm UTC](https://discourse.julialang.org/t/samplechains/55937/27 "2021-02-25T19:15:14Z")

</div>

> [@cpfiffer](#):
>
> We have exactly the same logic as you:
> 
> ```julia
> chain = sample(model, sampler, 1000)
> 
> chain = resume(chain, 1000)
> 
> ```
> 
> But in this case `chain` just has an empty bucket to store the model and the state. The sampler and the model decide what to put in the empty bucket.

Ok, yes I think that’s the same. It’s very confusing to call `chain` an empty bucket, or to think of it as only holding the samples. The “state” in your case is equivalent to an “iterator” as I’ve been talking about it. They’re exactly the same, otherwise `chain = resume(chain, 1000)` would be impossible.

Does `resume(chain, 1000)` allocate new memory and copy everything over?

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [March 4, 2021, 4:28pm UTC](https://discourse.julialang.org/t/samplechains/55937/28 "2021-03-04T16:28:51Z")

</div>

This is now registered:

> **[GitHub - cscherrer/SampleChainsDynamicHMC.jl](https://github.com/cscherrer/SampleChainsDynamicHMC.jl)**
>
> Contribute to cscherrer/SampleChainsDynamicHMC.jl development by creating an account on GitHub.

My goals here are

- Make this very easy to use from DynamicHMC directly
- Have similarly easy-to-use implementations for other back-ends
- Connect to PPLs like Soss with very little extra code
- Work toward a connection with AbstractMCMC.jl

I’m still thinking through that last point. The biggest benefit seems to be some nice-looking code for running chains in parallel. Ideally, a one-time effort to interface to AbstractMCMC would mean getting samplers “for free” down the road. I’m skeptical it would really work out this way, but I hope I might be wrong about that.

Anyway, I don’t think this will happen without some help. AbstractMCMC seems designed with `MCMCChains.jl` in mind. Even if it’s possible in principle to go outside this, the direct approach would seem to involve a lot of extra allocation and copying.

One very minor point: AbstractMCMC has `chain_type` very far down the list of arguments. In other Julia functions that take a `::Type` argument that will become the return type, this is one of the first arguments. The more common approach makes it more prominent, and also makes it easy to dispatch as

```julia
mynewmethod(::Type{WhateverType}, args...) = someotherfunction(args...) 

```

[Previous page](https://discourse.julialang.org/t/samplechains/55937.md?page=1)