Should the initial parameter samples / values be recorded for each chain

[HarrisonWilde]

Following from a discussion between myself and @mohamed82008 , he managed to set up a quick fix for passing an initial theta into a Turing model, the code looks as follows using the fix he has made a PR for here #1281:

@model gdemo(x, y) = begin
    s ~ InverseGamma(2,3)
    m ~ truncated(Normal(0.0, sqrt(s)), 0.0, 2.0)
    x ~ Normal(m, sqrt(s))
    y ~ Normal(m, sqrt(s))
end;
model = gdemo(1.0, 1.0);
varinfo = Turing.VarInfo(model);
model(varinfo, Turing.SampleFromPrior(), Turing.PriorContext((m = 0.5, s = 0.5)));
init_theta = varinfo[Turing.SampleFromPrior()]
c = sample(model, HMC(0.01, 1), 100, init_theta = init_theta)

This results in the following:

julia > c = sample(model, HMC(0.01, 1), 1000, init_theta = init_theta)
┌ Info: Using passed-in initial variable values
│   init_theta =
│    2-element Array{Float64,1}:
│     1.0
└     1.0
Sampling 100%|██████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:00
Object of type Chains, with data of type 1000×11×1 Array{Float64,3}

Iterations        = 1:1000
Thinning interval = 1
Chains            = 1
Samples per chain = 1000
internals         = acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, is_accept, log_density, lp, n_steps, nom_step_size, step_size
parameters        = m, s

2-element Array{ChainDataFrame,1}

Summary Statistics
  parameters    mean     std  naive_se    mcse      ess   r_hat
  ──────────  ──────  ──────  ────────  ──────  ───────  ──────
           m  0.9333  0.0522    0.0017  0.0164   4.0161  1.7977
           s  1.1174  0.0723    0.0023  0.0206  12.6280  1.0097

Quantiles
  parameters    2.5%   25.0%   50.0%   75.0%   97.5%
  ──────────  ──────  ──────  ──────  ──────  ──────
           m  0.8084  0.9056  0.9439  0.9709  1.0116
           s  0.9544  1.0804  1.1237  1.1693  1.2343


julia > get(c, :m)
(m = [0.997230952122332; 0.99781519383941; … ; 0.8313393573236126; 0.8277092475063401],)

So the question now is: should the initial values appear in the chain? This would also verify that passing it in is working correctly though I am fairly confident of that.

[xukai92]

Both options make sense to me.

As for MCMC, we do obtain a chain (X_0, ..., X_T) that includes the state sampled from \pi_0. So it's a reasonable option to include it in the chain.

On the other side, the indexing of the chain object actually starts from 1. So ignoring X_0 in the object aligns the indexing, i.e. c[t] refers to X_t.

I'm not sure if we can support c[0] easily but I think it's reasonable to at least provide an option to include or exlcude X_0 in c.