-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathreplicate_figure_06_direction_given_length.R
107 lines (79 loc) · 3.61 KB
/
replicate_figure_06_direction_given_length.R
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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
# ------------------------------------------------------------------------------
# load resources
# ------------------------------------------------------------------------------
source("_packages.R")
source("helpers/_helpers.R")
# ------------------------------------------------------------------------------
# settings
# ------------------------------------------------------------------------------
set.seed(8675309)
M = 25000
x_vals = seq(0, 2 * pi, length.out = 500)
th = runif(1, 0, 2 * pi)
m1 = rnorm(2, sd = sqrt(10))
m2 = c(cos(th), sin(th))
r = rexp(1, rate = 0.5)
# ------------------------------------------------------------------------------
# raw materials
# ------------------------------------------------------------------------------
my_d_vals_m1 = density_angle_given_length(x_vals, m1[1], m1[2], r)
my_d_vals_m2 = density_angle_given_length(x_vals, m2[1], m2[2], 1)
vm_d_vals = dvonmises(x_vals, th, 1)
my_draws_m1 = sample_angle_given_length(M, m1[1], m1[2], r)
my_draws_m2 = sample_angle_given_length(M, m2[1], m2[2], 1)
vm_draws = as.vector( rvonmises(M, th, 1) )
max_d_m1 = maxdensity_angle_given_length(m1[1], m1[2], r)
max_d_m2 = maxdensity_angle_given_length(m2[1], m2[2], 1)
# ------------------------------------------------------------------------------
# check integrating constant
# ------------------------------------------------------------------------------
approx_C = integrate(kernel_angle_given_length, lower = 0, upper = 2 * pi, m1[1], m1[2], r)$value
exact_C = constant_angle_given_length(m1[1], m1[2], r)
c(approx_C, exact_C)
# ------------------------------------------------------------------------------
# compare von Mises density and my density
# ------------------------------------------------------------------------------
par(mfrow = c(1, 1))
plot(x_vals, my_d_vals_m2, type = "l", col = "red")
lines(x_vals, vm_d_vals, col = "blue", lty = 2)
# ------------------------------------------------------------------------------
# compare my draws and my density
# ------------------------------------------------------------------------------
par(mfrow = c(1, 1))
hist(my_draws_m1,
freq = FALSE,
ylim = c(0, max_d_m1),
main = "Sampling angle given length in r[cos a, sin a]' ~ N(m, I)",
breaks = "Scott",
col = "lightblue")
lines(x_vals, my_d_vals_m1, lwd = 2)
legend("right",
c(paste("m = [", round(m1[1], 3), ", ", round(m1[2], 3), "]'", sep = ""),
paste("r = ", round(r, 3), sep = "")), bty = "n")
#main = paste("My draws against my density (r = ", r, ")", sep = ""),
# ------------------------------------------------------------------------------
# compare my draws and von Mises draws
# ------------------------------------------------------------------------------
par(mfrow = c(1, 1))
qqplot(my_draws_m2, vm_draws, cex = 0.1, pch = 19)
abline(a = 0, b = 1, lty = 2, col = "red")
# ------------------------------------------------------------------------------
# compare my draws to von Mises density
# ------------------------------------------------------------------------------
par(mfrow = c(1, 1))
hist(my_draws_m2,
freq = FALSE,
ylim = c(0, max_d_m2),
main = "My draws (r = 1) against von Mises density",
breaks = "Scott")
lines(x_vals, vm_d_vals)
# ------------------------------------------------------------------------------
# compare von Mises draws to my density
# ------------------------------------------------------------------------------
par(mfrow = c(1, 1))
hist(vm_draws,
freq = FALSE,
ylim = c(0, max_d_m2),
main = "My density (r = 1) against von Mises draws",
breaks = "Scott")
lines(x_vals, my_d_vals_m2)