Writing a new distribution (2)

Hi everybody,

I’m going to write a power-law distribution ( y = x^-2) in which x comes from a uniform distribution in [0,100] interval. According to this topic I could write the following code:

import theano.tensor as tt
def power_law(alpha):
def logp(x):
log_density = -alpha*tt.log(x)
return log_density
return logp

(alpha is 2). but I can’t incorporate the limitation on x!. How can I do this?

Thank you

I guess the range of x will be defined outside the function, eg a prior. You can raise an error if x falls outside the range.

1 Like

This is what I wrote with this approach:


but It didn’t lead to the right answer…

@ricardoV94 you are my last hope to solve this issue! :slight_smile:

The way you wrote, the model is not making use of the model variable x anywhere. In the logp function, x is the value the sampler suggests for s1 at each iteration and it has nothing to do with the x variable in the model. You could have called it anything else (the standard is value). I am not sure how you wanted to incorporate x in your model, so I cannot give you a solution.

Can you show how you are generating your data and how the 2 unobserved variables x, s1 are used?

I have the feeling you don’t really want/ need x (s1 follows a power law distribution with your given alpha already, except it also accepts negative values as mentioned above) but I need you to clarify it.

To elaborate, I express my codes in terms of Numpy library. This is my target function: y = x^-2/3., and I want to draw samples from it, so we have:

np.random.power(a = 1/3)

note that in Numpy, the power law function is: y = ax^(a-1) (we can ignore the coefficient a!)
The above code implements what I want except that the default range for x is [0,1]. I want to extend it to [0,100].

x in my problem is a physical quantity which can vary between 0 ,100 and I have to consider this condition

Does this work?

def power(alpha):
  def logp(value):
    return tt.switch(
      tt.and_(value > 0, value < 100), 
      tt.log(alpha) + (a-1)*tt.log(value),
      -np.inf,
    )
  
  return logp

with pm.Model() as m:
  y = pm.DensityDist(`y`, power(1/3)). # an interval transform here would be more efficient than the switch statement
  ...

Alternatively you might manage with just modeling a uniform x and applying the power law function deterministically statistics - python scipy.stats.powerlaw negative exponent - Stack Overflow

Unfortunately, I couldn’t sample from your code!! It comes with big divergences (despite my target_accept which was 0.99!). Could you complete your code and sample from it?

Thank you

That’s expected when not defining the interval transform so that the sampler can work on an unconstrained space.

This seems to sample from the prior without divergences so it might work for you:

import pymc3.distributions.transforms as tr

def power(alpha):
  def logp(value):
    return tt.switch(
      tt.and_(value > 0, value < 100), 
      tt.log(alpha) + (alpha-1)*tt.log(value),
      -np.inf,
    )
  
  return logp

with pm.Model() as m:
  y = pm.DensityDist('y', power(1/3), transform=tr.interval(0, 100), testval=1)
  trace = pm.sample()

With the transform you don’t even need the switch statement in the logp, since invalid values will not ever be proposed by the sampler.

1 Like

Excellent!.. I think you can write a user-friendly documentation about Pymc3 library!.

But I’ve come up with another idea which it doesn’t work and I don’t know why!. These are my codes:


What is wrong with this approach!?

Because applying a transformation to a random variable Y = f(X) does not transform the distribution in the same way (P(Y) \neq f(P(X))). In your case, as you’ve raised a uniform distribution to a power, you will obtain something like a Beta distribution (and not a power-law distribution).

1 Like

Is there a way to obtain a power-law dist with this approach?