Wow… indeed it gives a good approximation. Will see if I can incorporate x_{min} as well and share. Thanks 
. Also, true about the format of the likelihood … both work.
def custom_func(p1,gamma):
n = pt.tensor.shape(p1).value[0]
return (n*pt.tensor.log(gamma-1)-n*pt.tensor.log(p1.min())- pt.tensor.sum(gamma*pt.tensor.log(p1/p1.min())))
Using expanded form.