0035 pre model gqs

[bob-carpenter]

This is ready to review.

I tried to flesh out how we'd actually code this and document it in the Reference Manual.

[WardBrian]

Rendered: https://github.com/stan-dev/design-docs/blob/0035-pre-model-gqs/designs/0035-pre-model-gqs.md

[nhuurre]

The Plummer (2015) paper says that naive cut algorithm does not converge to multiple imputation distribution, and I don't immediately see how a latent data block could implement the suggested tempered cut algorithm.

In my opinion, Plummer's discussion is pretty pessimistic about the utility of a latent data block:

Exact MCMC sampling from a cut model seems to require evaluating the marginal
likelihood p(Y | ϕ) which is computationally expensive
(Gelman and Meng 1998). It might be possible, as suggested
by a reviewer, that auxiliary variable methods such as the one
proposed by Møller et al. (2006) may allow this evaluation to
be skipped. Otherwise, it appears unlikely that any MCMC
method can sample exactly from the target density p∗(θ, ϕ).

[bob-carpenter]

Thanks. This I have to think about more. I'll have to go back and read Plummer's paper again!

[bob-carpenter]

I'm not convinced by Plummer's proof (p. 40). He assumes you have to satisfy a factored form of detailed balance that would be sufficient, but isn't necessary. Here's the full Bayes graphical model example from the paper, with parameters $\phi, \theta$ and data $y, z$,

$z \sim p(\cdot \mid \phi)$ and $y \sim p(\cdot \mid \phi, \theta)$.

Then BUGS lets you put a "cut" between $z, \phi$ and $y, \theta$ (at least that's how it's drawn---I like the graphic of a dashed line as "cut here").

Plummer says (section 2, right after equation (4)) that we can sample from the cut distribution using multiple imputation as

$\phi^{(n)} \sim p(\phi \mid z)$ and $\theta^{(n)} \sim p(\theta \mid \phi^{(n)}, y, z)$.

I'm pretty sure this is exactly what our latent data or cut parameters block will do. Each iteration $n$, it will sample $\phi^{(n)}$ in the cut parameters block (that gives one of the multiple imputations) and then sample `\theta^{(n)}$ as an ordinary parameter with target density $\theta^{(n)}$ for each $\phi^{(n)}$ given as above. Clearly if we sampled $\phi^{(n)}$ once and then $\theta^{(n)}$ multiple times, we'd get exactly the multiple imputation result someone might compute. So I don't see why it won't work if we only take one draw of $\theta^{(n)}$ for each $\phi^{(n)}$.

In Figure 3, Plummer shows that the "adaptive Metropllis 1D" sampler gets roughly the same answer as multiple imputation. We should also get roughly the same answer as multiple imputation if HMC is able to sample from
$p(\theta \mid \phi^{(n)}, y, z)$. The main worry we've had in the past about this is that it's a moving target for adaptation. But then not necessarily more so than other complex posteriors.

I added section at the end on functional testing and added a note saying that we'll test by making sure we can recreate the multiple imputation results in an example like Plummer's.

I may be missing something here, so I'm going to leave this open pending further comments.

[nhuurre]

What is this needed for?

[bob-carpenter]

@andrewgelman's been bugging us to do this for ages so that

a. he can write something like if (mod(iteration_number, 100)) print(...), which I think is probably not going to do what he wants because it'll print every leapfrog iteration, and

b. so we can implement something like tempering by using iteration number to modify density. Again, all sorts of issues with these algorithms.

[SteveBronder]

This feels adjacent to the main point of the proposal. As long as we don't introduce state to the model class I'm fine with this

[bob-carpenter]

I should also add a big note that this proposal does not include any mutable state in the model.

Yes, iteration number is adjacent and we could put it in a separate proposal, but it's going to be really easy to do this way, I think. I can take it out if people think it's over-reaching.

[nhuurre]

The first alternative is marginalization.
The second alternative is to implement the HMC+Gibbs algorithm in the Stan algorithms library and just allow models to declare discrete parameters.
Also we should consider special syntax for discrete parameters that may allow either marginalization or Gibbs chosen at runtime.

[bob-carpenter]

I'll clarify that I meant there isn't an existing alternative for Gibbs within Stan.

You can marginalize only when tractable, so it rules out high dimensional categorical and unbounded count data and also things like variable selection (not that those work well with Gibbs).

The original design of the model class included discrete parameters because I was thinking we could just use the entire density to update with Metropolis or with Metropolis-within-Gibbs style updates.

[bob-carpenter]

Thanks---added both clarifications as suggested. I'll leave this open for comment now rather than closing.

I think a proposal to add automatic Gibbs sampling to Stan should be another design document. We could do generic Metropolis-within-Gibbs, but I asked the PyMC devs, and they say they can compute the Markov blanket to do this more efficiently.

[SteveBronder]

Since it's just holding data that is made once I think we can use an allocator scheme instead. So we define a latent data struct inside of the model class which will also use the allocator for acquiring memory.

struct latent_data {
  int iteration_number__;  // included by default in all latent_data
  // Will hold memory owned by allocator
  Eigen::Map<Eigen::Matrix<double, -1, 1>> sens; 
  std::vector<double, stan::math::local_allocator> spec;  
};

For the allocator, we will then also make a member function in the model block that will tell us the size of the latent data struct.

std::size_t get_latent_size_bytes() {
  // Assume `sens` is size 5 and spec is size 6
  return sizeof(latent_data) + 5 * sizeof(double) + 6 * sizeof(double);
}

The allocator will look very similar to Stan's arena allocator, but instead of storing a vector of unsigned char pointers this will only store one void pointer. We will use just one void pointer so that it's easy to cast the memory inside of the allocator to the latent_data class when we need it.

struct local_allocator {
  void* mem_;
  // Must be fixed size
  local_allocator(std::size_t init_bytes) : mem_(malloc(init_bytes));
  template<typename T>
  void* allocate(std::size_t amount);

  template <typename T>
  T* cast_this() {
    return reinterpret_cast<T*>(mem_);
  }
};

Since the sizes of variables have to come from data and transformed data we can know at initialization the amount of memory we need.

stan::math::local_allocator latent_alloc(model->get_latent_size_bytes());

Then the model class will have a latent_data member function that uses the allocator as an in/out parameter

inline void
set_latent_data(stan::math::allocator& alloc) {
  latent_data* = new (alloc.allocate<latent_data>()) latent_data{
    0, 
    Eigen::Map<Eigen::Matrix<double, -1, 1>>(
      new (alloc.allocate<double>(5), 5, 1)),
    std::vector<double, stan::math::allocator>(5, alloc)
  };
  // Fill in latent_data...
  return;
}

So now we can have a new log_prob that can take in the allocator. The log_prob will call the allocators cast function to get back the pointer to the latent_data. We can either use the latent_data pointer for the rest of the code or pull out references to the members of latent_data.

auto log_prob(stan::math::allocator& alloc, ...) {
  latent_data* = alloc.cast<latent_data>();
  // Just makes a reference, very low cost
  auto&& spec = latent_data->spec;
  // Rest of the program written as normal
  auto X = Y * spec;
}

And with this scheme the model class would also have new member functions for printing latent data names and latent data sizes etc.

[bob-carpenter]

Since it's just holding data that is made once I think we can use an allocator scheme instead.

What's the advantage to doing it the way you propose? It looks very similar to me, so I'm sure I'm missing something important.

[SteveBronder]

It removes inheritance, does not need a generic factory method, and goes more into how to handle dynamically allocated objects. I don't think we need to have a base class that we upcast to the actual class. Since this is just used in the model class the rest of the stan code can just worry about passing around an allocator that holds the memory for the latent data struct.