check if RelaxedBernoulli (supporting rsample) would do what you want. For more advanced use, modeling uncertainty, sample Bernoulli parameter from a trainable Beta distribution:
pseudo-code is something like:
P = Beta(softplus(raw_a), softplus(raw_b)) # raw_a = nn.Parameter(…)
p = P.rsample()
B = RelaxedBernoulli(probs=p)
b = B.rsample()
loglik = P.log_prob(p) + B.log_prob(b) #for proper probabilistic models