Write a scipy function without using a standard library (exponential power)

Question

My question might come across as stupid or so simple, but I could not work towards finding a solution. Here is my question: I want to write an exponential power distribution function which is available in scipy. However, I don't want to use the scipy for this. How do I go about it?

Here are my efforts so far:

import math
import numpy as np

def ExpPowerFun(x,b, size=1000):
        distribution = b*x**(b-1)*math.exp(1+x**b-math.exp(x**b))
        return distribution

I used this equation based on this scipy doc. To be fair, using this equation and writing a function using it doesn't do much. As you can see, it returns only one value. I want to generate a distribution of random numbers based on scipy's exponential power distribution function without using scipy. I have looked at class exponpow_gefrom github code. However, it uses scipy.special(-sc), so it's kind of useless for me, unless there is any workaround and avoids the use of scipy.

I can't figure out how to go about it. Again, this might be a simple task, but I am stuck. Please help.

Answer 1

the simplest way to generate a random number for a given distribution is using the inverse of the CDF of that function, the PPF (Percent point function) will give you the distribution you need when you apply it on uniform distributed numbers.

for you case the PPF (taken directly from scipy source code with some modifications) is:

np.power(np.log(1-np.log(1-x)), 1.0/b)

hence you code should look like this

def ExpPowerFun(b, size=1000):
    x = np.random.rand(size)
    return np.power(np.log(1-np.log(1-x)), 1.0/b)

import matplotlib.pyplot as plt
plt.hist(ExpPowerFun(2.7,10000),20)
plt.show()

Edit: the uniform distribution has to be from 0 to 1 ofc since the probabilities are from 0% to 100%