from scipy.optimize import curve_fit
import numpy as np
from matplotlib import pyplot as plt
# Define the von mises kernel density estimator
def circular_von_mises_kde(x,mu,sigma):
# Adjust data to take it to range of 2pi
x = [(hr)*2*np.pi/24 for hr in x]
mu*=2*np.pi/24
sigma*=2*np.pi/24
# Compute kappa for vm kde
kappa = 1/sigma**2
return np.exp((kappa)*np.cos((x-mu)))/(2*np.pi*i0(kappa))
# Assuming your data is occurences of some event at the given hour of the day
frequencies= np.zeros((24))
frequencies[data]=1
hr_data = np.linspace(1,24, 24)
fit_params, cov = curve_fit(circular_von_mises_kde, hr_data, data_to_fit, bounds=(0,24))
plt.plot(hr_data, frequencies, 'k.',label='Raw data')
plt.plot(np.linspace(1,25, 1000), circular_von_mises_kde(np.linspace(1,25, 1000), *fit_params), 'r-',label='Von Mises Fit')
plt.legend()
plt.xlabel('Hours (0-24)')
plt.show()
print('The predicted mean is {mu} and the standard deviation is {sigma}'.format( mu=round(fit_params[0],3), sigma=round(fit_params[1], 3)))
很好地找到了VM分布。这就是战斗的一半。 但是,除非我被scipy.stats.vonmises docs中的公式弄错了,否则公式假设数据的中心是0,而事实可能并非如此。所以我们应该建立自己的VM发行版。对于我们的虚拟机分发,我们将确保它是周期性的,在24小时的范围内,而不是传统的2pi范围。请参阅下面的代码和注释。另外,我假设你的数据是你看到一些事件发生的时间,如果不是这样,你需要重新调整。在
Click to see the result of the above code *作为一个快速警告,你可能需要一个更大的数据集来做一些适当的拟合,并真正建立一个人口趋势。在
相关问题 更多 >
编程相关推荐