dist_l1_case_1.jl 1.06 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

using MarkovProcesses
import QuadGK: quadgk
load_model("SIR")

# Case 1 : 6.4

x_obs = [2, 4, 3, 3]
t_x = [0.0, 0.5, 1.2, 2.2]
t_y_bis = [0.0, 0.5, 1.2, 2.2]
values = zeros(length(x_obs), 1)
values[:,1] = x_obs
l_tr = Vector{Nothing}(nothing, length(x_obs))
σ1 = Trajectory(SIR, values, t_x, l_tr)

y_obs = [6, 6]
t_y = [0.0, 2.2]
t_x_bis = [0.0, 2.2]
values = zeros(length(y_obs), 1)
values[:,1] = y_obs
l_tr = Vector{Nothing}(nothing, length(y_obs))
σ2 = Trajectory(SIR, values, t_y, l_tr)

24
test_1 = dist_lp(@view(x_obs[:]), @view(t_x[:]), @view(y_obs[:]), @view(t_y[:]); p=1) == 6.4
25
26
27
28
29
30
31
32
33
34
35
test_1 = test_1 && dist_lp(σ1, σ2; p=1) == 6.4

f_x(t::Real) = MarkovProcesses._f_step(x_obs, t_x, t)
f_y(t::Real) = MarkovProcesses._f_step(y_obs, t_y, t)
diff_f(t) = abs(f_x(t) - f_y(t))
int, err = quadgk(diff_f, 0.0, 2.2)

test_2 = isapprox(6.4, int; atol=err)

# Case 1 bis : inverse of case 1 

36
test_1_bis = dist_lp(@view(y_obs[:]), @view(t_y[:]), @view(x_obs[:]), @view(t_x[:]); p=1) == 6.4 
37
38
39
40
test_1_bis = test_1_bis && dist_lp(σ2, σ1; p=1) == 6.4

return test_1 && test_1_bis && test_2