使用 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 创建它们的直方图并将直方图值绘制到图中。您可能决定将 bins 的中间作为点的位置（这个假设当然是错误的，但使用的 bins 越多越有效）。

拟合与问题代码中的一样。然后您会发现我们的拟合大致找到了用于数据生成的参数（在本例中低于 ~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

我同意@ImportanceOfBeingErnes 的解决方案，但我想为分发添加一个（众所周知的？）通用解决方案。如果您的分布函数 f 具有积分 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

我认为您实际上询问的是回归问题，这正是 Praveen 的建议。

你有一个沼泽标准指数衰减，它到达 y 轴大约 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)

如果您想进一步 'test that my function is indeed going through the data points' 那么我建议您在残差中寻找模式。但像这样的讨论可能超出了 Whosebug 的欢迎范围：Q-Q 和 P-P 图、残差图与 y 或 x，等等。

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

Generate random numbers from exponential distribution and model using python

python

random

numpy

curve-fitting

pandas