python - 如何使用优化功能复制scipy.stats.fit？

我正在尝试使分布适合某些值。这是我的代码

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，但事实并非如此。怎么了？。

最佳答案

怎么了？。



成本函数是错误的。
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的帕累托分布无法容纳numpy生成的数据。

import matplpotlib.pyplot as plt
import scipys as sp
import numpy as np

plt.subplot(2, 1, 1)
plt.hist(np.random.pareto(1.5, 10000), 1000)  # This is a Lomax or Pareto II distribution
plt.xlim(0, 10)

plt.subplot(2, 1, 2)
plt.hist(sp.stats.pareto.rvs(1.5, size=1000), 1000)  # This is a Pareto distribution
plt.xlim(0, 10)

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

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)

根据文档，numpy.random.pareto并非完全来自帕累托分布：


  从具有指定形状的Pareto II或Lomax分布中抽取样本。
  
  Lomax或Pareto II分布是移位的Pareto分布。可以通过将Lomax分布加1并乘以比例参数m来获得经典的Pareto分布（请参见注释）。


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


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