Hi… Thanks for pointing out the error. Indeed the sum should be only on the last term and the first two terms should be multiplied by the size of the sample n. I thought by doing a sum on the entire term it should be incorporated. However, I made the changes and used NUTS instead as below.
simulated_dist = pw.Power_Law(xmin=0.5,parameters=[2.5])
test_data = simulated_dist.generate_random(1000)
with pm.Model() as powerlaw_model:
def custom_func(p1,gamma,xmin):
n = pt.tensor.shape(p1).value[0]
return (n*pt.tensor.log(gamma-1)-n*pt.tensor.log(xmin)- pt.tensor.sum(gamma*pt.tensor.log(p1/xmin)))
gamma = pm.Uniform('gamma',1,4)
xmin = pm.Uniform('xmin',0.1,1)
power_law = pm.DensityDist('power_law',gamma,xmin,logp=custom_func,observed=test_data)
with powerlaw_model:
step = pm.NUTS()
samples = pm.sample(1000,step=step,chains=3)
Getting this result now… Still not able to fit
I am wondering if I am doing something conceptually wrong here…