If you had written the suggested solution using PyTensor operations, would one still be able to use gradient-based (i.e., HMC/NUTS) samplers? My understanding is that writing functions myself in pymc means I cannot use gradient-based samplers like HMC/NUTS.
I wrote the suggestion with PyTensor operations. pm.logp returns an expression that uses PyTensor operations. I could have done the same manually. If you write with PyTensor operations (which should be very similar to writing numpy code) you get gradients for free and can sample with NUTS.
Is using pm.sample() the same as sample_prior_predictive as long as I do not include data in observed=data (in my case I do not have observed data; I just want to sample a distribution)?
Yes, they should be sampling from the same distribution. The distinction is prior_predictive uses forward sampling (requiring random number generators), and sample uses mcmc (requiring logp expressions and optional gradients for those)
I’m assuming the above code will use Metropolis-Hastings or HMC; is that correct?
The example code I wrote will use NUTS.
Is there a way to specify a starting configuration where I want the sampling to start from (obviously if it is a random sampling it wouldn’t make sense, but I am using sequential Monte Carlo, so choice of starting point makes sense)? For example, let’s say I want to run the suggested code for 100 time steps so I am at a configuration C, then initialize another Metropolis-Hastings or HMC with C as a starting point. Is there a way to do that?
You can define the starting point by passing initvals to sample: pymc.sampling.mcmc.sample — PyMC dev documentation
However there’s a complication to your picture. These samples need tuning before they start giving correct draws, and during tuning they will drift from your suggested initvals. If you want more fine-grained control you have to either initialize the step samplers already tuned manually, or perhaps override their tuning routine. This is non-trivial, specially for advanced samplers like NUTS.
For sequential monte carlo, you can provide the distribution of initial particle positions via the start kwarg I think, otherwise it attempts to initialize from the prior (using forward sampling): pymc.smc.sample_smc — PyMC dev documentation