python - 如何通过Von Mises分布求周期间隔和周期均值？

Question

我有一些时间数据(一天中的几个小时)。我想对这些数据进行冯·米塞斯分布拟合，并找到周期平均值。如何在 python 中使用 scipy 执行此操作？

例如:

from scipy.stats import vonmises
data = [1, 2, 22, 23]
A = vonmises.fit(data)

我不确定如何使用拟合或均值或区间方法获得该数据的分布(可能是区间)和周期平均值。

Answer 1

在查找 VM 发行版方面做得很好。这就是战斗的一半。但除非我误会了 scipy.stats.vonmises docs 中的公式，该公式假设数据以 0 为中心，但情况可能并非如此。所以我们应该构建自己的虚拟机发行版。对于我们的 Vm 分布，我们将确保它在 24 小时范围内具有周期性，而不是传统的 2pi 范围。请参阅下面的代码和注释。另外，我假设您的数据是您看到某些事件发生的时间，如果情况并非如此，您将需要重新调整。

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)))

Click to see the result of the above code* 需要注意的是，您可能需要更大的数据集来进行适当的拟合并真正建立人口趋势。