What would be a good way to get sample mean of joint distribution?

It appears the .mean() method is not implemented for customized joint distribution? We may use sampling method to estimate but it does not appear stable. Can anyone suggest an effective method?

For example, we take the guide book distribution:

joint = tfd.JointDistributionNamedAutoBatched(dict(
e= tfd.Exponential(rate=[100, 120]),
g=lambda e: tfd.Gamma(concentration=e[0], rate=e[1]),
n= tfd.Normal(loc=0, scale=2.),
m=lambda n, g: tfd.Normal(loc=n, scale=g),
x=lambda m: tfd.Sample(tfd.Bernoulli(logits=m), 12),
))

We would expect joint.mean() to give marginal mean of each component, but it is not implemented per error message. We use sampling method as follow:

joint.mean() not implemented

small_sample=joint.sample(2000)
large_sample=joint.sample(3000)
for k in small_sample.keys():
print('Smaller sample mean',k,tf.reduce_mean(small_sample[k],axis=0).numpy())

print()

for k in large_sample.keys():
print('Larger sample mean',k,tf.reduce_mean(large_sample[k],axis=0).numpy())

However, estimates from different sample sizes are quite different for sequential variables (e.g., n and m):
Smaller sample mean e [0.00983032 0.00844856]
Smaller sample mean g 9.612298
Smaller sample mean n -0.028275907
Smaller sample mean m 4.339757
Smaller sample mean x [0 0 0 0 0 0 0 0 0 0 0 0]

Larger sample mean e [0.00985759 0.00873891]
Larger sample mean g 12.751243
Larger sample mean n -0.0056907604
Larger sample mean m -14.086281
Larger sample mean x [0 0 0 0 0 0 0 0 0 0 0 0]

Are there any better method to get marginal means?

Thanks,
Xiaolong

The issue provides a JointDistributionNamedAutoBatched example and sample-mean code but names no repository files or tests. Start by reproducing the .mean() failure with the supplied distribution and determine the intended marginal-mean behavior for sequential variables; done requires an agreed implementation or documented supported method, with coverage for the example.