# Trouble sampling a model with the exponent of a random variable

Hello,

I’m unable to recover the parameters of a model whenever I have a non-unity exponent of a random variable. Following is the code

``````gamma = 5.
beta=np.sqrt(1-1/gamma**2)
angles = np.random.normal(np.pi/4,np.pi/6,100)
alpha=1.
delta = (1/gamma * 1/(1-beta*np.cos(angles)))**alpha + np.random.normal(0,0.1,100)

with pm.Model() as model:

BoundedNormal = pm.Bound(pm.Normal, lower=1.25)
gamma_prior=BoundedNormal('gama',3,5)
alpha_prior=pm.Uniform('alpha',0,2)
s=pm.Normal('noise',0,1)

beta_val=np.sqrt(1-1/gamma_prior**2)
obs_val=np.power(1/gamma_prior * 1/(1-beta_val*np.cos(angles)),alpha)

obs = pm.Normal('observation',mu=obs_val+s, observed=delta,testval=0.1)
trace = pm.sample(5000,return_inferencedata=True,tune=5000)
``````

When I set alpha=1, I’m able to recover the parameters with high accuracy. However, when I set alpha to, for example, 0.4, all I get is a wide and flat posterior on gamma. Should I be using a different sampler? Please help me with this issue.

Thanks,
Karthik

I don’t think it has anything to do with the unity/non-unity. Your data is just more informative about `gamma` for higher values of `alpha`

``````alphas = (0.4, 0.8, 1, 1.2, 1.6)
gammas = (1, 3, 5, 7)
_, ax = plt.subplots(len(deltas), len(gammas), figsize=(20, 16))
for row, alpha in enumerate(deltas):
for col, gamma in enumerate(gammas):
beta=np.sqrt(1-1/gamma**2)
angles = np.random.normal(np.pi/4,np.pi/6,1000)
delta = (1/gamma * 1/(1-beta*np.cos(angles)))**alpha + np.random.normal(0,0.1,1000)
ax[row, col].hist(delta, ec='k', bins=20)
ax[row, col].set_title(f'alpha={alpha}, gamma={gamma}')

plt.tight_layout()
``````

Accordingly, if I generate data with `alpha=1.5` I get an even more narrow posterior for `gamma` than with `alpha=1.0`.

2 Likes

The `alpha` parameter on the other hand seems completely unlearnable

Ah, that makes sense now. Thanks very much!