1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
|
import numpy as np
import scipy.stats
from scipy.special import i0
import numpy.testing
def von_mises_cdf_series(k,x,p):
x = float(x)
s = np.sin(x)
c = np.cos(x)
sn = np.sin(p*x)
cn = np.cos(p*x)
R = 0
V = 0
for n in range(p-1,0,-1):
sn, cn = sn*c - cn*s, cn*c + sn*s
R = 1./(2*n/k + R)
V = R*(sn/n+V)
return 0.5+x/(2*np.pi) + V/np.pi
def von_mises_cdf_normalapprox(k,x,C1):
b = np.sqrt(2/np.pi)*np.exp(k)/i0(k)
z = b*np.sin(x/2.)
C = 24*k
chi = z - z**3/((C-2*z**2-16)/3.-(z**4+7/4.*z**2+167./2)/(C+C1-z**2+3))**2
return scipy.stats.norm.cdf(z)
def von_mises_cdf(k,x):
ix = 2*np.pi*np.round(x/(2*np.pi))
x = x-ix
k = float(k)
# These values should give 12 decimal digits
CK=50
a = [28., 0.5, 100., 5.0]
C1 = 50.1
if k<CK:
p = int(np.ceil(a[0]+a[1]*k-a[2]/(k+a[3])))
F = np.clip(von_mises_cdf_series(k,x,p),0,1)
else:
F = von_mises_cdf_normalapprox(k,x,C1)
return F+ix
|
