如何复制scipy.stats.fit使用优化函数？

Question

我试图使分布符合某些值。这是我的密码

from __future__ import print_function
import pandas as pd
import numpy as np
import scipy as sp
import scipy.optimize as opt
import scipy.stats
import matplotlib.pyplot as plt

    values = np.random.pareto(1.5, 10000)
    loc = values.min()
    scale = 1


    def cost_function(alpha):
        cost = -sp.stats.pareto(alpha, loc=loc, scale=scale).pdf(values)
        return cost.sum()

    opt_res = opt.fmin(cost_function, 1.5)

    alpha_fit_v = sp.stats.pareto.fit(values, floc=loc, fscale=scale)


    print('opt_res = ', opt_res,
          ' alpha_fit_v = ', alpha_fit_v)

我原以为alpha_fit_v等同于opt_res，但事实并非如此。怎么了？。在

Answer 1

What's wrong?.

1. 成本函数是错误的。在
2. np.random.pareto的分布与sp.stats.pareto不同

1。成本函数错误

求逆概率之和是没有意义的。您需要使用对数：

def cost_function(alpha):
    cost = -sp.stats.pareto(alpha, loc=loc, scale=scale).logpdf(values)
    return cost.sum()

2。np.random.pareto的分布与sp.stats.pareto

不同

这是一个棘手的问题，但是您可能已经注意到即使sp.stats.pareto.fit也没有返回正确的结果。这是因为scipy的Pareto分布无法拟合numpy生成的数据。在

^{pr2}$

也就是说，这将按预期工作：

values = sp.stats.pareto.rvs(1.5, size=1000)
loc = 0
scale = 1

def cost_function(alpha):
    cost = -sp.stats.pareto(alpha, loc=loc, scale=scale).logpdf(values)
    return cost.sum()

opt_res = opt.fmin(cost_function, 1.5)

alpha_fit_v = sp.stats.pareto.fit(values, floc=loc, fscale=scale)

print('opt_res = ', opt_res,
      ' alpha_fit_v = ', alpha_fit_v)

# opt_res =  [ 1.49611816]  alpha_fit_v =  (1.4960937500000013, 0, 1)

根据文献，^{}并不完全来自于帕累托分布：

Draw samples from a Pareto II or Lomax distribution with specified shape.

The Lomax or Pareto II distribution is a shifted Pareto distribution. The classical Pareto distribution can be obtained from the Lomax distribution by adding 1 and multiplying by the scale parameter m (see Notes).

因此，如果使用numpy生成数据，您有两种选择：

• 您可以为scipy分布设置loc=-1。在
• 您可以执行values = np.random.pareto(1.5, 10000) + 1并设置loc=0。在