Trouble sampling a model with the exponent of a random variable

Questions
1

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

2

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()

image
image1034×727 37.2 KB

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
3

The alpha parameter on the other hand seems completely unlearnable

4

Ah, that makes sense now. Thanks very much!