-
Notifications
You must be signed in to change notification settings - Fork 1.1k
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
(complementary) cumulative distribution for use in JointDistributionSequential #501
Comments
I am not sure there is a similar
|
Thanks a lot @junpenglao, I hadn't thought of that workaround! |
Hi @junpenglao sorry for the delay, I couldn't work on this until something else was solved :-) As a check if what I'm doing is correct, could you please take a look at the following code (please excuse the clumsy translation to Python ;-))? Also, in case you'd say it's correct but could easily be made more elegant (in a language-independent way I mean, not the Python syntax), I'd be happy for any suggestions. The model is meant to represent durations (censored and uncensored), stored in the variable Just wondering, would you have a better idea for dealing with rates for an exponential? Many thanks!
|
I would use matrix operation and juggling the shape a bit to cut down some lines: # place holder data for simulation.
data_array = np.hstack([np.ones((1000, 1)),
np.random.randn(1000, 6)])
def model(data_array):
return tfd.JointDistributionSequential(
[
tfd.Sample(
tfd.Normal(loc=0., scale=1.),
sample_shape=7),
lambda betas:
tfd.Independent(
tfd.Exponential(
rate = 1/tf.transpose(
tf.matmul(tf.cast(data_array, dtype=betas.dtype),
tf.transpose(betas)))
), reinterpreted_batch_ndims = 1)
]
)
unconstraining_bijectors = [
tfb.Exp(),
tfb.Identity()
]
def get_exponential_lccdf(betas, data_array, check_time):
e = tfd.Independent(
tfd.Exponential(
rate = 1/tf.transpose(
tf.matmul(tf.cast(data_array, dtype=betas.dtype),
tf.transpose(betas)))
), reinterpreted_batch_ndims = 1)
cum_prob = e.cdf(check_time)
return (1- cum_prob)
# construct model for uncensored data only
m = model(data_array)
# this function is written to allow for easy comparison of resulting parameters
# dependent on whether we use the potential or not
def get_log_prob (model_data, data_array, censored_data = None):
def log_prob(betas):
lp = m.log_prob([betas, model_data.check_time])
potential = get_exponential_lccdf(betas, censored_data) if censored_data is not None else 0
return (lp + potential)
return log_prob |
That is a lot nicer, thank you!! |
Hi Junpeng, thanks again :-) I'd have a final question (see below for complete code with shape annotations)... With the above refactoring, I had to transpose the target array(s), e.g.
otherwise I'd get an error during logprob calculation while sampling. To be precise, this would still work:
but then during sampling I'd get
Honestly I don't really understand, as the output shape is not affected by the refactoring...
|
hi @junpenglao sorry for bothering again, I'd have an additional, more urgent question...
instead. But when I try that (removing the bijector, although that doesn't make a difference) the model does not work at all (all parameter samples keep their initial value, no steps accepted, in spite of trying different step sizes and initial values). |
Hi there sorry about the slow response - but regarding your more pressing issue, the |
Hi @junpenglao thank you!! Whatever I did, I haven't been able to make the Even with minimal step sizes of I've extracted the Stan code from
Would you have an opinion which of the TFP models
that conforms to, in spirit? (The Thanks again! |
Yep that's precisely the root cause - no matter how small you set the step size if the initial return value is nan HMC is not going to work. The one with bijector actually transform the initial value to reasonable range first - in another word, to make it work for the new model, you need to set the initial value to something reasonable, and you should keep the bijector as your model still contains bound variable and HMC only really works on variable defined on the real line. |
Hi @junpenglao thanks again :-) I've hardcoded some values in the snippet below (this is actually the complete censored portion, 103 rows, and it's really not that different from the uncensored one...) In the below code, I've varied initial values and step sizes but I never get a single HMC step accepted... Many thanks in any case!
|
Hi, I decided to publish the post with the previous model that worked fine: https://blogs.rstudio.com/tensorflow/posts/2019-07-31-censored-data/. I think that should be ok too. In case you still had ideas reg. the above, or the shape thing #501 (comment), I would of course be happy to hear them, but I know you have things to work on so feel free to close :-) Thanks again for your help! |
(Close this for now but please feel free to follow up if you have another question!) |
Thanks again! I had tried all of the above, also in combination of 2, but never all of them together. Updated the post, - should be acceptable to statisticians (as opposed to DL people who are used to alchemy ;-)) now as well :-) Thanks again! |
Glad it helped! There are quite a bit of preprocessing/preconditioning in brms (and many other similar libraries) so that the inference could be ran automatically, and when you are using low level libraries like TFP or pymc3 you would need to spell out these boilerplate codes explicitly - it's not always trivial :-) |
Yeah! Here my main error in thinking was that, as I had first used normalization in the non-exp model but noticed that TFP seemed to be handling the substantial differences in scale between the predictors without problems, the exp model should work fine without as well, given that I experimented with standard deviations of 0.001, and was using 0 for initial values... But using all three together made the difference :-) |
Hi,
I would like to construct a simple model, using
JointDistributionSequential
, for durations ,where part of the data is censored.The non-censored part can be modeled using
tfd.Exponential
, while the censored data should be modeled using the (exponential) complementary cumulative distribution (exponential_lccdf
in Stan, see https://mc-stan.org/docs/2_19/functions-reference/exponential-distribution.html).Would someone have a suggestion how to do this? I suppose that in the exponential case, as the ccdf is simple to construct, one would have to wrap the formula
exp(- lambda * x)
in aTransformedDistribution
somehow, or in a subclass ofDistribution
?Many thanks for any help!
(BTW there's also a PyMC3 example https://docs.pymc.io/notebooks/bayes_param_survival_pymc3.html where they wrap the ccdf in a Potential class, but I don't think there is a comparable object in TFP...?)
The text was updated successfully, but these errors were encountered: