Извлечение временных характеристик с использованием периодического нормального распределения (фон мизес) в Python

Question

Я пытаюсь найти среднее, дисперсию и доверительный интервал периодического / свернутого нормального распределения (фон Мизеса), но в пределах временного интервала (в отличие от традиционного интервала числа пи). Я посмотрел на решение по переполнению стека здесь , оно закрыто, но я не уверен, что это именно то, что я ищу.

Я нашел именно то, что искал здесь , который использует 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 вместо времени. Есть ли альтернативные пакеты, которые могут помочь?

Answer 1

Я создал функцию 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