python随机抽号,python：从自定义概率函数中随机抽样

I have a piecewise quartic distribution with a probability density function:

p(x)= c(x/a)^2 if 0≤x

c((b+a-x)^2/b)^2 if a≤x≤b;

0 otherwise

Suppose c, a, b are known, I am trying to draw 100 random samples from the distribution. How can I do it with numpy/scipy?

解决方案

One standard way is to find an explicit formula, G = F^-1 for the inverse of the cumulative distribution function. That is doable here (although it will naturally be piecewise defined) and then use G(U) where U is uniform on [0,1] to generate your samples.

In this case, I think that I worked out the details, but you will need to check the Calculus/Algebra.

First of all, to streamline things it helps to introduce a couple of new parameters. Let

f(a,b,c,d,x) = c*x**2 #if 0 <= x <= a

and

f(a,b,c,d,x) = d*(x-e)**4 #if a < x <= b

Then your p(x) is given by

p(x) = f(a,b,c/a**2,c/b**2,a+b)

I integrated f to find the cumulative distribution and then inverted and got the following:

def Finverse(a,b,c,d,e,x):

if x <= (c*a**3)/3:

return (3*x/c)**(1/3)

else:

return e + ((a-e)**5 - (5*c*a**3)/(3*d))**(1/5)

Assuming this is right, then simply:

def randX(a,b,c):

u = random.random()

return Finverse(a,b,c/a**2,c/b**2,a+b,u)

In this case it was possible to work out an explicit formula. When you can't work out such a formula for the inverse, consider using the Monte Carlo methods described by @lucianopaz