55 lines
1.4 KiB
Python
55 lines
1.4 KiB
Python
import numpy as np
|
|
|
|
|
|
def randn(mu, sigma):
|
|
"""Normal random number generator
|
|
mean mu, standard deviation sigma"""
|
|
return sigma * np.random.randn() + mu
|
|
|
|
|
|
def rand_von_mises(mu, kappa):
|
|
"""Generate random samples from the Von Mises distribution"""
|
|
if kappa < 1e-6:
|
|
# kappa is small: sample uniformly on circle
|
|
theta = 2.0 * np.pi * np.random.ranf() - np.pi
|
|
else:
|
|
a = 1.0 + np.sqrt(1.0 + 4.0 * kappa * kappa)
|
|
b = (a - np.sqrt(2.0 * a)) / (2.0 * kappa)
|
|
r = (1.0 + b*b) / (2.0 * b)
|
|
|
|
not_done = True
|
|
while not_done:
|
|
u1 = np.random.ranf()
|
|
u2 = np.random.ranf()
|
|
|
|
z = np.cos(np.pi * u1)
|
|
f = (1.0 + r * z) / (r + z)
|
|
c = kappa * (r - f)
|
|
|
|
not_done = (u2 >= c * (2.0 - c)) and (np.log(c) - np.log(u2) + 1.0 - c < 0.0)
|
|
|
|
u3 = np.random.ranf()
|
|
theta = mu + np.sign(u3 - 0.5) * np.arccos(f)
|
|
|
|
return theta
|
|
|
|
|
|
if __name__ == '__main__':
|
|
# Tests
|
|
|
|
print("Gaussian distribution:")
|
|
r = np.zeros(1000)
|
|
for i in range(1000):
|
|
r[i] = randn(1.0, 2.0)
|
|
|
|
print("True mean 1.0 == Estimated mean ", np.mean(r))
|
|
print("True std 2.0 == Estimated std ", np.std(r))
|
|
|
|
print("Von Mises distribution:")
|
|
|
|
t = np.zeros(1000)
|
|
for i in range(1000):
|
|
t[i] = rand_von_mises(1.0, 8)
|
|
|
|
print("True mean 1.0 == Estimated mean ", np.mean(t))
|