2016-05-26 60 views
2

我得到100个经验值:比较直方图泊松分布和高斯曲线

N = np.array([ 6080, 5604, 5972, 5566, 5683, 5619, 5582, 5912, 5614,5837, 5961, 5972, 5807, 5829, 5881, 5757, 5725, 5809, 5626, 5995, 5793, 5608, 5880, 5982, 5748, 6071, 6181, 6034, 6117, 5903, 6190, 5735, 6109, 6126, 6012, 5948, 6139, 6103, 6108, 6031, 6200, 6091, 6199, 6165, 6591, 5803, 6093, 5921, 6194, 5799, 6020, 6156, 6129, 6344, 6243, 6122, 5926, 5904, 5579, 5881, 6157, 5925, 5835, 5778, 6125, 5737, 5703, 5809, 6109, 5978, 5881, 6250, 6143, 5658, 5815, 5633, 5780, 5620, 6180, 5770, 6058, 5688, 5792, 6170, 5915, 6147, 5727, 6300, 6049, 6263, 6168, 6156, 6071, 6196, 6078, 5848, 5847, 6248, 6243, 6084]) 

现在我想用泊松分布和高斯分布进行比较。 我想把曲线上顶我的直方图:

enter image description here]

但它不工作。 例如与泊松分布:我要让泊松分布,奠定了我的直方图

b = np.sqrt(len(N)) 
w = np.ones_like(N)/float(len(N)) 

entries, bin_edges, patches = plt.hist(N, bins=b, weights=w) 


def poisson(x, lamb): 
    return (lamb**x/factorial(x)) * np.exp(-lamb) 

parameters, cov_matrx = curve_fit(poisson, bin_middles, entries) 

的飞度但是当我绘制它,它甚至没有接近

plt.plot(x_plot, poisson(x_plot, *parameters), 'r-', lw=2) 

Where is Poisson?

你会如何解决这个问题?我认为这是因为我的价值观太高。泊松分布以k = 0,1,2,...运行... 矿井就像k'= 6000,6200,...

+0

我一直想在自己的环境中测试这一点,但我理解你如何将泊松传递给curve_fit?我得到一个警告'不能估计参数的协方差'。你能展示一个更完整的代码摘录吗? – Grr

+0

是的,对我也一样。我用这个指南工作..我不知道为什么它不与我的价值http://stackoverflow.com/questions/25828184/fitting-to-poisson-histogram – rfaenger

回答

1

在你的图中有一条红线......它是泊松曲线。你在用什么拉姆达?这很可能太小,你只能得到零。 无论如何,我建议你使用:从SciPy的

def poisson(x, lamb): 
    from scipy.stats import poisson as _poisson  
    return _poisson.pmf(x,lamb) 

parameters, cov_matrx = curve_fit(poisson, bin_middles, entries, p0=6000) 

poisson功能,您还应该添加normed=True你的第一个直方图

entries, bin_edges, patches = plt.hist(N, bins=b, weights=w, normed=True) 
+0

我试着类似于本指南工作:http:// stackoverflow.com/questions/25828184/fitting-to-poisson-histogram ...我以为'curve_fit'计算我的lambda。不是吗?我为每个bin计算了bin_middle。 bin_middles和计数或可能性是我的x,y我从中制作curve_fit ... – rfaenger

+0

我链接了错误的函数...我编辑了我的文章...尝试使用它而不是您的poisson实现...反正估计的lambda是多少? – elabard

+0

是不是:'lambda = np.mean(N)'?我的数组的意思?无论如何:我想curve_fit为我计算最好的lambda ....我现在很困惑。感谢您的链接。 – rfaenger