Bernoulli stuck at 1e-6

I try to use Hamiltonian Monte Carlo to fit a model consists of a single Bernoulli distribution.
The model can get stuck at p=1e-6. Replacing Bernoulli by RelaxedBernoulli doesn't solve the problem.

If I use tensorflow-probability to fit a more complex model that contains Bernoulli distributions, would the parameters all get stuck at 0.0 and 1.0?

Result:

[[1.e-06]
 [1.e-06]
 [1.e-06]
 ...
 [1.e-06]
 [1.e-06]
 [1.e-06]]
mean:0.0000  stddev:0.0000  acceptance:0.0000

Expected result:
I expect the parameter p should be 0.3 instead of 1e-6

Test:

import numpy as np
import tensorflow as tf
import tensorflow_probability as tfp

tfd = tfp.distributions
from functools import partial


def logp_func(observations, p):
    b = tfd.RelaxedBernoulli(np.float32(0.5), probs=p)
    return tf.reduce_sum(b.log_prob(observations))


def main():
    obs = np.random.choice(2, p=[0.7, 0.3], size=1000).astype(np.float32)
    observations = tf.constant(np.array(obs))
    logp = partial(logp_func, observations)

    step_size = tf.get_variable(
        name='step_size',
        initializer=1.,
        use_resource=True,  # For TFE compatibility.
        trainable=False)

    num_results = int(1e4)
    num_burnin_steps = int(1e3)
    hmc = tfp.mcmc.HamiltonianMonteCarlo(
        target_log_prob_fn=logp,
        num_leapfrog_steps=3,
        step_size=step_size,
        step_size_update_fn=tfp.mcmc.make_simple_step_size_update_policy(
            num_adaptation_steps=int(num_burnin_steps * 0.8)))

    samples, kernel_results = tfp.mcmc.sample_chain(
        num_results=num_results,
        num_burnin_steps=num_burnin_steps,
        current_state=np.array([1e-6], dtype=np.float32),
        kernel=hmc)

    # Initialize all constructed variables.
    init_op = tf.global_variables_initializer()

    with tf.Session() as sess:
        init_op.run()
        samples_, kernel_results_ = sess.run([samples, kernel_results])
        print(samples_)
        print('mean:{:.4f}  stddev:{:.4f}  acceptance:{:.4f}'.format(
            samples_.mean(), samples_.std(),
            kernel_results_.is_accepted.mean()))


main()

Start by running the provided Python reproduction with the single Bernoulli and RelaxedBernoulli cases. Inspect the HMC setup, especially sample_chain, current_state, step-size adaptation, and the distribution log probability, then compare the observed acceptance and samples with the expected p≈0.3. Done means identifying whether the behavior is expected or correcting it and adding a regression test.