depth_points = pm.Uniform(“depth_points”, lower=0, upper=max_depth, shape=num_layers)
I want to generate an increasing sequence of numbers. Is there any way to do this? I’d appreciate it.
depth_point_offsets = pm.Uniform(“depth_point_offsets ”, lower=0, upper=max_depth, shape=num_layers)
depth_points = pm.Determinstic('depth_points', depth_point_offsets.cumsum())
You can do arbitrary transformations of RVs after you’ve made them, so in this case you could just take their cumulative sum. If you want to enforce a maximum, you could pump the depth_points though some bounded function (like inverse_logit).
1 Like
THanks so much!
You should think about what prior you want over the sequence. And whether it’s upper bounded by max_depth here. If it’s upper bounded by max_depth, then you can’t just take uniform points in (0, max_depth) and cumulative sum without risk of going over.
An alternative that respects n upper bound is to take a simplex, apply cumulative sum, then scale by max_depth. That will imply a prior based on the prior for the simplex. You can apply a Jacobian correction for the cumulative sum and scaling to adjust it back to uniform over the increasing sequences bounded by 0 and max_depth and then apply whatever prior you want. For example, if it’s cutpoints in an ordinal logistic regression, you might want to apply a zero-avoiding prior to the differences in the sequence.
3 Likes
Thanks so much. I will try this way.
To be clear, @bob-carpenter 's suggestion is related to the logp computation for a hypothetical “ordered bounded uniform” prior. Your could implement that in PyMC as a custom prior, but it’s not typically the way we suggest you work. In PyMC you are declaring a generative forward model (as opposed to the logp/density model that is declared in STAN). The implied logp is then automatically worked out for you by inspecting the forward graph (that generates the random draws). So the following is perfectly valid and works out to similar (but not identical! Mine starts with an unconstrained uniform latent quantity and transforms into a non uniform constrained prior whereas Bob’s suggestion starts with a constrained non uniform quantity and transform it into a constrained uniform prior):
depth_points_unsorted = pm.Uniform("depth_points_unsorted", lower=0, upper=max_depth, shape=num_layers)
depth_points = pm.Deterministic('depth_points', pt.sort(depth_points_unsorted), dims=['layer'])
You should not, in general, ever be handling jacobian adjustments yourself or working out log probabilities in PyMC. That said, if you’re really interested in that approach, you can use built-in transformations to help you. For example, we have the ordered transformation that automatically handles the transformation (and jacobian adjustment) that Bob is talking about. These transformations can be chained together, so you could also apply an interval transformation to force the values to be in a certain range (this is actually always quietly applied for you in the background when you use pm.Uniform). The disadvantage of doing this is that you lose access to prior and posterior predictive sampling, because the forward model is not aware of the transformations (which only act on the logp graph). You also have to be more careful about initial values for your priors, as we will see below.
For comparison, here are the resulting priors under both approaches. First mine:
import pymc as pm
import arviz as az
import matplotlib.pyplot as plt
import pytensor.tensor as pt
from pymc.distributions.transforms import Ordered
import numpy as np
def plot_layers(idata, group='prior'):
fig, ax = plt.subplots(figsize=(14, 4), dpi=144)
for layer in coords['layer']:
az.plot_posterior(idata[group],
var_names=['depth_points'],
ax=ax,
coords={'layer':layer},
hdi_prob='hide',
point_estimate=None,
c=plt.color_sequences['tab10'][layer],
label=f'Layer {layer}')
ax.legend()
plt.show()
max_depth = 10
num_layers = 5
coords = {'layer':range(num_layers)}
with pm.Model(coords=coords) as m:
depth_points_unsorted = pm.Uniform("depth_points_unsorted", lower=0, upper=max_depth, dims=['layer'])
depth_points = pm.Deterministic('depth_points', pt.sort(depth_points_unsorted), dims=['layer'])
idata = pm.sample_prior_predictive()
plot_layers(idata)
And using a transform (with automatic Jacobian adjustment):
with pm.Model(coords=coords) as m:
depth_points = pm.Uniform("depth_points", lower=0, upper=max_depth, dims=['layer'],
transform=Ordered(),
# You have to manually specify that the initial values satisfy the ordering constraint, otherwise you will get a "-inf logp at the initial value" error
initval=np.linspace(1e-4, max_depth - 1e-4, num_layers))
# Cannot use sample_prior_predictive here, you will just get back untransformed uniform draws if you do
idata = pm.sample()
plot_layers(idata, group='posterior')
As you can see, the results are very similar.
2 Likes
To clarify, I’m not sure the sort is a good strategy because NUTS may not love the discontinuities in the gradient (if it ever finds them after tuning). A scaled cumsum of Uniform/Dirichlet may be just fine.
In general it’s hard to find a forward constraining transform that’s equivalent to parameter transform+jacobian. For instance the sort won’t give the same results if the uniform prior is not IID.
Either way it can be hard to think about these sorts of priors, specially when they mix dimensions. We are assigning them some well understood prior densities but the implied parameters can be hard to grok.
Like if you assign a normal density to a log transformed variable it will certainly never behave very “normally” :). And similarly, ordered variables will often not look anything like the original independent densities (and they don’t here).
I usually prefer to do the forward transform like @jessegrabowski illustrated (also because it’s much more natural in PyMC) and rely on prior predictive to understand its implications after constraining. I’m too lazy to start from the constrained prior and work backwards. It also may matter little if the likelihood swamps the prior anyway.
The forward approach is also how I would go about generating fake data, so the model doesn’t end up looking too absurd to me (there’s a circularity here, if I used mcmc to generate fake data I might be of the opposite opinion)
2 Likes
@ricardoV94 @jessegrabowski Thanks sooooo much for help
3 posts were split to a new topic: Separation between data and parameters
It’s a privilege to have you here – thanks for sharing your knowledge. The work from the Stan team has been instrumental in the development of PyMC and I hope there’s ways we can contribute back.
Feel free to open a thread to ask any dev questions