CODEWITHSUNDEEP
pythonscipydistribution
VincentN
VincentN

Reputation: 111

Creating sample from custom distribution using scipy

I am trying to get some sample from a given distribution, in fact, it is a 3-parameter Pareto distribution. Here are the codes:

from scipy.stats import gamma, rv_continuous

class pareto3_pdf(rv_continuous):
    def _pdf(self,x,alpha,lambd,k):
        return (gamma(alpha + k) * lambd**alpha * x**(k - 1)) / (gamma(alpha) * gamma(k) * (lambd + x)**(alpha + k))
pareto3 = pareto3_pdf(name="pareto")


x = pareto3.rvs(alpha = 3,lambd = 4,k = 2)
print(x)

and the output: TypeError: unsupported operand type(s) for *: 'rv_frozen' and 'int'

I am not quite sure how to fix this. If anyone has any suggestion it would be appreciated.

Thank you in advance.

edit:

I have now changed the code, but it keeps giving negative values.

import scipy.stats as stats
from scipy.stats import rv_continuous
from scipy.special import gamma

class pareto3_pdf(rv_continuous):
    def _pdf(self,x,alpha,lambd,k):
        return (gamma(alpha + k) * lambd**alpha * x**(k - 1)) / (gamma(alpha) * gamma(k) * (lambd + x)**(alpha + k))
pareto3 = pareto3_pdf(name="pareto")
pare3 = pareto3.rvs(alpha = 5,lambd = 4,k = 2)
print(pare3)

and if I try to simplify this into a 2-parameter model, OverflowError: (34, 'Result too large') error popup.

import scipy.stats as stats
from scipy.stats import rv_continuous
from scipy.special import gamma

class pareto2_pdf(rv_continuous):
    def _pdf(self,x,alpha,lambd):
        return (alpha * lambd**alpha / (lambd + x)**(alpha + 1))
pareto2 = pareto2_pdf(name="pareto2")
pare2 = pareto2.rvs(alpha = 2,lambd = 2)
print(pare2)

Upvotes: 1

Views: 1502

Answers (2)

user6655984
user6655984

Reputation:

As I wrote elsewhere your distribution is available in SciPy as betaprime(k, alpha, scale=lamda), so sampling it is built-in. A little test:

from scipy.stats import betaprime
alpha, lamda, k = 5, 4, 2
sample = betaprime.rvs(k, alpha, scale=lamda, size=1000)
print(sample.mean())
print(betaprime.mean(k, alpha, scale=lamda))

prints 

2.0134570579012108
2.0

Close enough. (Of course, the mean of a random sample is random.)

Upvotes: 0

Einar A
Einar A

Reputation: 141

You have to import gamma from scipy.special instead of scipy.stats.
The reason is that scipy.stats.gamma is distribution and scipy.special.gamma is the gamma function.

from scipy.stats import rv_continuous 
from scipy.special import gamma 

class pareto3_pdf(rv_continuous):
    def _pdf(self,x,alpha,lambd,k):
        return (gamma(alpha + k) * lambd**alpha * x**(k - 1)) /(gamma(alpha) * gamma(k) * (lambd + x)**(alpha + k))
pareto3 = pareto3_pdf(name="pareto")
x = pareto3.rvs(alpha = 3,lambd = 4,k = 2)

Upvotes: 5

Related Questions