python:从自定义概率函数中随机抽样 [英] python: random sampling from self-defined probability function

问题描述

我有一个具有概率密度函数的分段四次分布:

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

p(x)= c(x/a)^2 if 0≤x<a; 
      c((b+a-x)^2/b)^2 if a≤x≤b;
      0 otherwise

假设c，a，b是已知的，我试图从分布中抽取100个随机样本.我该如何用numpy/scipy做到?

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?

推荐答案

一种标准方法是为累积分布函数的逆找到一个明确的公式G = F^-1.这是可行的(尽管自然会是分段定义的)，然后使用G(U)，其中U在[0,1]上是统一的，以生成您的样本.

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

和

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

然后您的p(x)由

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

我集成了f来找到累积分布，然后求逆，得到以下结果:

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)

在这种情况下，可以得出一个明确的公式.当您无法计算出这样的逆公式时，请考虑使用@lucianopaz

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

这篇关于python:从自定义概率函数中随机抽样的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！