python - 给定 pdf 的 scipy 中的自定义分发

Question

我尝试使用通过 scipy.stats 给出的 pdf 定义自定义分布

import numpy as np
from scipy.stats import rv_continuous

class CustomDistribution(rv_continuous):
    def __init__(self, pdf=None):
        super(CustomDistribution, self).__init__()
        self.custom_pdf = pdf
        print "Initialized!"

    def _pdf(self, x, *args):
        if self.custom_pdf is None:
            # print 'PDF is not overridden'
            return super(CustomDistribution, self)._pdf(x, *args)
        else:
            # print 'PDF is overridden'
            return self.custom_pdf(x)

def g(x, mu):
    if x < 0:
        return 0
    else:
        return mu * np.exp(- mu * x)

my_exp_dist = CustomDistribution(pdf=lambda x: g(x, .5))
print my_exp_dist.mean()

如您所见，我尝试定义参数 mu=0.5 的指数分布，但输出如下。

Initialized!

D:\Anaconda2\lib\site-packages\scipy\integrate\quadpack.py:357:

IntegrationWarning: The algorithm does not converge. Roundoff error is detected in the extrapolation table. It is assumed that the requested tolerance cannot be achieved, and that the returned result (if full_output = 1) is the best which can be obtained.
warnings.warn(msg, IntegrationWarning)

D:\Anaconda2\lib\site-packages\scipy\integrate\quadpack.py:357:

IntegrationWarning: The maximum number of subdivisions (50) has been achieved.

2.0576933609

If increasing the limit yields no improvement it is advised to analyze the integrand in order to determine the difficulties. If the position of a local difficulty can be determined (singularity, discontinuity) one will probably gain from splitting up the interval and calling the integrator on the subranges. Perhaps a special-purpose integrator should be used. warnings.warn(msg, IntegrationWarning)

我应该怎么做才能改善这一点？

注意:计算精度问题在this GitHub issue中讨论。 .

Answer 1

这似乎可以满足您的要求。每次创建实例时，都必须为类的实例赋予 lambda 参数值。 rv_continuous 足够聪明，可以推断出您没有提供的项目，但您当然可以提供我在这里提供的更多定义。

from scipy.stats import rv_continuous
import numpy

class Neg_exp(rv_continuous): 
    "negative exponential"
    def _pdf(self, x, lambda):
        self.lambda=lambda
        return lambda*numpy.exp(-lambda*x)
    def _cdf(self, x, lambda):
        return 1-numpy.exp(-lambda*x)
    def _stats(self,lambda):
        return [1/self.lambda,0,0,0]

neg_exp=Neg_exp(name="negative exponential",a=0)

print (neg_exp.pdf(0,.5))
print (neg_exp.pdf(5,.5))

print (neg_exp.stats(0.5))

print (neg_exp.rvs(0.5))