Sensible Defaults

[willtebbutt]

While the example in the README is important in that it exposes all of the useful bits of AdvancedHMC, it would be nice if it were possible to define a couple of sensible defaults. e.g.

vanilla_hmc_sampler = HMC(step_size, n_steps)
reasonably_sensible_nuts_sampler = NUTS(step_size)

I think that this kind of thing would be a helpful addition to the package, as it makes it more straightforward for a new user to get up and running with their first sampler / for an expert to prototype stuff before optimising.

In particular, it would make writing examples for Stheno.jls documentation quite a lot less verbose.

edit:

I think the above would likely tie a sampler and an adapter together. The sampling API would be something like

samples, stats = sample(h, vanilla_hmc_sampler, n_samples_to_draw; progress=true)

and by default it would maybe spend n_samples_to_draw steps burning in / adapting, then draw the samples.

[scheidan]

I'd like to bump this issue. I'm planning to use AdvancedHMC for a teaching exercise (outside of Turing) so really need a a much simpler, user-friendly interface.

For my own use I wrote this wrapper function below (more or less directly from the Readme). Is something like this of scope for this package? It would make it much easier for users not using Turing to profit from your work.

# NUTS sampler as used by STAN
Based on documentation of AdvancedHMC.

```Julia
stanHMC(lp::Function, ∇lp::Function,  θ_init,;
        n_samples::Int=1000, n_adapts::Int=n_samples÷2)
```

## Arguments
- `lp`:  log density to sample from (up to a constant)
- `∇lp`: function that computes gradient
- `n_samples::Int=1000`: number of samples
- `n_adapts::Int=n_samples÷2`: length of adaptation

### Return Value
A named tuple with fields:
- `samples`: array containing the samples
- `stats`: contains various statistics form the sampler. See AdvancedHMC documentation.
"""
function advancedHMC_NUTS(lp::Function, ∇lp::Function,  θ_init,;
                          n_samples::Int=1000, n_adapts::Int=n_samples÷2)

    # Define a Hamiltonian system
    D = length(θ_init)   # number of parameters
    metric = DiagEuclideanMetric(D)

    # choose AD framework or provide a function manually
    hamiltonian = Hamiltonian(metric, lp, θr -> (lp(θr), ∇lp(θr)))

    # Define a leapfrog solver, with initial step size chosen heuristically
    initial_ϵ = find_good_stepsize(hamiltonian, θ_init)
    integrator = Leapfrog(initial_ϵ)

    # Define an HMC sampler, with the following components
    #   - multinomial sampling scheme,
    #   - generalised No-U-Turn criteria, and
    #   - windowed adaption for step-size and diagonal mass matrix
    proposal = NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)
    adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator))

    # -- run sampler
    samples, stats = sample(hamiltonian, proposal, θ_init, n_samples,
                            adaptor, n_adapts; progress=true)

    return (samples=samples, stats=stats)

end

[yebai]

@JaimeRZP's recent work likely fixes this. See e.g. #323