使用python从指数分布和模型生成随机数

Question

我的目标是创建一个随机点数据集,其直方图看起来像一个指数衰减函数,然后通过这些点绘制指数衰减函数.

首先,我试图从指数分布中创建一系列随机数(但是没有成功,因为这些应该是点,而不是数字).

from pylab import *
from scipy.optimize import curve_fit
import random
import numpy as np
import pandas as pd

testx = pd.DataFrame(range(10)).astype(float)
testx = testx[0]

for i in range(1,11):
   x = random.expovariate(15) # rate = 15 arrivals per second
   data[i] = [x]

testy = pd.DataFrame(data).T.astype(float)
testy = testy[0]; testy

plot(testx, testy, 'ko')

结果看起来像这样.

然后我定义了一个函数来绘制一条直线:

def func(x, a, e):
return a*np.exp(-a*x)+e

popt, pcov = curve_fit(f=func, xdata=testx, ydata=testy, p0 = None, sigma = None) 

print popt # parameters
print pcov # covariance

plot(testx, testy, 'ko')

xx = np.linspace(0, 15, 1000)
plot(xx, func(xx,*popt))

plt.show()

我正在寻找的是:(1)从指数(衰变)分布创建随机数组的更优雅方式,以及(2)如何测试我的函数确实通过数据点.

Answer 1

我猜想以下内容接近你想要的.您可以使用numpy生成从指数分布中提取的一些随机数,

data = numpy.random.exponential(5, size=1000)

然后,您可以使用它们创建直方图,并将numpy.hist直方图值绘制到绘图中.你可以决定把箱子的中间作为点的位置(这个假设当然是错误的,但是你使用的垃圾箱越多越有效).

拟合工作与问题的代码一样.然后,您将发现我们的拟合粗略地找到用于数据生成的参数(在这种情况下低于~5).

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

data = np.random.exponential(5, size=1000)

hist,edges = np.histogram(data,bins="auto",density=True )

x = edges[:-1]+np.diff(edges)/2.
plt.scatter(x,hist)

func = lambda x,beta: 1./beta*np.exp(-x/beta)

popt, pcov = curve_fit(f=func, xdata=x, ydata=hist) 

print(popt)

xx = np.linspace(0, x.max(), 101)
plt.plot(xx, func(xx,*popt), ls="--", color="k", 
         label="fit, $beta = ${}".format(popt))
plt.legend()
plt.show()

Answer 2

我认为你实际上是在问回归问题，这就是 Praveen 所建议的。

\n\n

您有一个沼泽标准指数衰减，在 y 轴大约 y=0.27 处到达。因此它的方程是y = 0.27*exp(-0.27*x). 我可以围绕该函数的值建立高斯误差模型，并使用以下代码绘制结果。

\n\n

import matplotlib.pyplot as plt\nfrom math import exp\nfrom scipy.stats import norm\n\n\nx = range(0, 16)\nY = [0.27*exp(-0.27*_) for _ in x]\nerror = norm.rvs(0, scale=0.05, size=9)\nsimulated_data = [max(0, y+e) for (y,e) in zip(Y[:9],error)]\n\nplt.plot(x, Y, \'b-\')\nplt.plot(x[:9], simulated_data, \'r.\')\nplt.show()\n\nprint (x[:9])\nprint (simulated_data)\n

\n\n

这是情节。请注意，我保存了输出值以供后续使用。

\n\n

\n\n

现在我可以计算自变量上受噪声污染的指数衰减值的非线性回归，这就是curve_fit所做的事情。

\n\n

from math import exp\nfrom scipy.optimize import curve_fit\nimport numpy as np\n\ndef model(x, p):\n    return p*np.exp(-p*x)\n\nx = list(range(9))\nY = [0.22219001972988275, 0.15537454187341937, 0.15864069451825827, 0.056411162886672819, 0.037398831058143338, 0.10278251869912845, 0.03984605649260467, 0.0035360087611421981, 0.075855255999424692]\n\npopt, pcov = curve_fit(model, x, Y)\nprint (popt[0])\nprint (pcov)\n

\n\n

额外的好处是，不仅curve_fit计算参数 \xe2\x80\x94 0.207962159793 \xe2\x80\x94 的估计，还提供了该估计方差的估计 \xe2\x80\x94 0.00086071 \xe2\x80 \x94 作为 的元素pcov。鉴于样本量较小，这似乎是一个相当小的值。

\n\n

以下是计算残差的方法。x请注意，每个残差是数据值与使用参数估计值估计的值之间的差值。

\n\n

residuals = [y-model(_, popt[0]) for (y, _) in zip(Y, x)]\nprint (residuals)\n

\n\n

如果您想进一步“测试我的函数确实正在通过数据点”，那么我建议在残差中寻找模式。但像这样的讨论可能超出了 stackoverflow 上受欢迎的范围：QQ 和 PP 图、残差 vsy或 的图x等等。

\n