使用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

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

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

import matplotlib.pyplot as plt
from math import exp
from scipy.stats import norm


x = range(0, 16)
Y = [0.27*exp(-0.27*_) for _ in x]
error = norm.rvs(0, scale=0.05, size=9)
simulated_data = [max(0, y+e) for (y,e) in zip(Y[:9],error)]

plt.plot(x, Y, 'b-')
plt.plot(x[:9], simulated_data, 'r.')
plt.show()

print (x[:9])
print (simulated_data)

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

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

from math import exp
from scipy.optimize import curve_fit
import numpy as np

def model(x, p):
    return p*np.exp(-p*x)

x = list(range(9))
Y = [0.22219001972988275, 0.15537454187341937, 0.15864069451825827, 0.056411162886672819, 0.037398831058143338, 0.10278251869912845, 0.03984605649260467, 0.0035360087611421981, 0.075855255999424692]

popt, pcov = curve_fit(model, x, Y)
print (popt[0])
print (pcov)

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

下面是如何计算残差的。注意，每个残差都是数据值和使用参数估计值从x估计的值之间的差。

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

如果你想进一步“测试我的函数是否确实通过了数据点”，那么我建议在残差中寻找模式。但这样的讨论可能超出了stackoverflow所欢迎的范围：Q-Q和P-P图，残差图vsy或x，等等。

Answer 2

我同意@ImportanceOfBeingErnes的解决方案，但我想添加一个（众所周知的？）分布的一般解。如果有一个带积分的分布函数f（即f = dF / dx），那么通过将随机数与inv F（即积分的反函数）映射，可以得到所需的分布。在指数函数的情况下，积分又是指数，而逆是对数。所以可以这样做：

import matplotlib.pyplot as plt
import numpy as np
from random import random


def gen( a ):
    y=random()
    return( -np.log( y ) / a )


def dist_func( x, a ):
    return( a * np.exp( -a * x) )


data = [ gen(3.14) for x in range(20000) ]
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
ax.hist(data, bins=80, normed=True, histtype="step") 
ax.plot(np.linspace(0,5,150), dist_func( np.linspace(0,5,150), 3.14 ) )
plt.show()

Answer 3

我想下面的内容和你想要的很接近。你可以用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()