python - 在 Python 中使用周期性正态分布(von mises)提取时间特征
问题描述
我试图找到周期性/包裹正态分布(von Mises)的均值、方差和置信区间,但在一个时间间隔内(与传统的 pi 间隔相反)。我在这里查看了关于堆栈溢出的解决方案,它很接近,但我不确定它是否正是我正在寻找的。
我在这里找到了我正在寻找的东西,它使用了 R(见下面的代码摘录)。我希望在 Python 中复制它。
> data(timestamps)
> head(timestamps)
[1] "20:27:28" "21:08:41" "01:30:16" "00:57:04" "23:12:14" "22:54:16"
> library(lubridate)
> ts <- as.numeric(hms(timestamps)) / 3600
> head(ts)
[1] 20.4577778 21.1447222 1.5044444 0.9511111 23.2038889 22.9044444
> library(circular)
> ts <- circular(ts, units = "hours", template = "clock24")
> head(ts)
Circular Data:
[1] 20.457889 21.144607 1.504422 0.950982 23.203917 4.904397
> estimates <- mle.vonmises(ts)
> p_mean <- estimates$mu %% 24
> concentration <- estimates$kappa
> densities <- dvonmises(ts, mu = p_mean, kappa = concentration)
> alpha <- 0.90
> quantile <- qvonmises((1 - alpha)/2, mu = p_mean, kappa = concentration) %% 24
> cutoff <- dvonmises(quantile, mu = p_mean, kappa = concentration)
> time_feature <- densities >= cutoff
与库通告一样,python 有一个包 scipy.stats.vonmises 但位于间隔 pi 而不是时间之内。是否有任何替代软件包可以提供帮助?
解决方案
我构建了一个 python 函数,它可以满足我的需要,从这个pdf中获取公式
希望这对社区有所帮助。如果我错了,请提供更正。
注意:这适用于区间 [0,2pi] 或 360 度内的值。
import pandas as pd
import numpy as np
from scipy.stats import chi2
def random_dates(start, end, n, unit='D', seed=None):
if not seed:
np.random.seed(0)
ndays = (end - start).days + 1
return pd.to_timedelta(np.random.rand(n) * ndays, unit=unit) + start
def vonmises(df, field):
N = len(df[field])
s = np.sum(np.sin(df[field]))
c = np.sum(np.cos(df[field]))
sbar = (1/N)*s
cbar = (1/N)*c
if cbar > 0:
if sbar >= 0:
df['mu_vm'] = np.arctan(sbar/cbar)
else:
df['mu_vm'] = np.arctan(sbar/cbar) + 2*np.pi
elif cbar < 0:
df['mu_vm'] = np.arctan(sbar/cbar) + np.pi
else:
df['mu_vm'] = np.nan
R = np.sqrt(c**2 + s**2)
Rbar = (1/N)*R
if Rbar < 0.53:
kstar = 2*Rbar + Rbar**3 + 5*(Rbar**5)/6
elif Rbar >= 0.85:
kstar = 1/(3*Rbar -4*(Rbar**2) + Rbar**3)
else:
kstar = -0.4 + 1.39*Rbar + 0.43/(1-Rbar)
if N<=15:
if kstar < 2:
df['kappa_vm'] = np.max([kstar - 2/(N*kstar),0])
else:
df['kappa_vm'] = ((N-1)**3)*kstar/(N*(N**2+1))
else:
df['kappa_vm'] = kstar
if Rbar <= 2/3:
df['vm_plus'] = df['mu_vm'] + np.arccos(np.sqrt(2*N*(2*(R**2) -
N*chi2.isf(0.9,1))/((R**2)*(4*N - chi2.isf(0.9,1)))))
df['vm_minus'] = df['mu_vm'] - np.arccos(np.sqrt(2*N*(2*(R**2) -
N*chi2.isf(0.9,1))/((R**2)*(4*N - chi2.isf(0.9,1)))))
else:
df['vm_plus'] = df['mu_vm'] + np.arccos(np.sqrt((N**2) -
((N**2) - (R**2))*np.exp(chi2.isf(0.9,1)/N))/R)
df['vm_minus'] = df['mu_vm'] - np.arccos(np.sqrt((N**2) -
((N**2) - (R**2))*np.exp(chi2.isf(0.9,1)/N))/R)
df['vm_conft'] = np.where((df['vm_plus'] < df[field]) |
(df['vm_minus'] > df[field]), True, False)
return df
df = pd.concat([pd.DataFrame({'A':[1,1,1,1,1,2,2,2,2,2]}), pd.DataFrame({'B':random_dates(pd.to_datetime('2015-01-01'), pd.to_datetime('2018-01-01'), 10)})],axis=1)
df['C'] = (df['B'].dt.hour*60+df['B'].dt.minute)*60 + df['B'].dt.second
df['D'] = df['C']*2*np.pi/(24*60*60)
df = df.groupby('A').apply(lambda x : vonmises(x, 'D'))
例如,要回到小时,只需乘以 24 并除以 2pi
推荐阅读
- python - 使用 for 循环在两个数据集中查找最接近值的值
- apache-spark - livy.rsc.jars 和 livy.repl.jars 有什么区别?
- javascript - 无法使用新的 angular/cli 应用程序运行 e2e 测试
- amazon-web-services - 有没有办法找出每个调用 lambda 的资源?
- swift - 如何在 Mac 上使用带有效果的 appleScript 发送 iMessage(如 slam 效果)
- r - 与 spatstat 的多重比较
- css - 如何使用 twitter-bootstrap 4.2 创建 100% 高度的左侧菜单?
- css - 如何使用 css 在“内联”中显示数据?
- python - AWS Sagemaker SKlearn 入口点允许多个脚本
- excel - Excel 中带有和 + 或条件的 SUMIFS