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
extern crate rand;
pub mod event;
use rand::prelude::*;
use rand::distributions::Poisson;
use event::Event;
pub fn poisson_process(tmax: f64, lambda: f64) -> Vec<Event> {
let mut rng = thread_rng();
assert!(lambda >= 0.0);
let num_events = Poisson::new(tmax*lambda).sample(&mut rng);
let mut result = vec![];
for _ in 0..num_events {
let timestamp= tmax*random::<f64>();
result.push(Event::new(timestamp));
};
result
}
pub fn variable_poisson(tmax: f64, lambda: fn(f64) -> f64, max_lambda: f64) -> Vec<Event> {
let mut rng = thread_rng();
let num_events = Poisson::new(tmax*max_lambda).sample(&mut rng);
let mut result = vec![];
for _ in 0..num_events {
let timestamp = random::<f64>()*tmax;
let lambda_val = random::<f64>()*max_lambda;
if lambda_val < lambda(timestamp) {
let mut event = Event::new(timestamp);
event.add_intensity(lambda_val);
result.push(event);
}
}
result
}
pub fn hawkes_exponential(tmax: f64, alpha: f64, beta: f64, lambda0: f64) -> Vec<Event> {
let mut t = 0.0;
let mut previous_t: f64;
let mut last_lambda = lambda0;
let mut result: Vec<Event> = vec![];
while t < tmax {
let u0 = random::<f64>();
let s0 = -1.0/lambda0*u0.ln();
let d = if last_lambda > lambda0 {
1.0 + beta*u0.ln()/(last_lambda - lambda0)
} else { std::f64::NEG_INFINITY };
let s1: f64;
previous_t = t;
t += if d > 0.0 {
s1 = -1.0/lambda0*random::<f64>().ln();
if s0 <= s1 {
s0
} else {
s1
}
} else {
s0
};
if t > tmax {
break;
}
last_lambda = lambda0 + alpha + (last_lambda - lambda0)*(-beta*(t-previous_t)).exp();
let mut new_event = Event::new(t);
new_event.add_intensity(last_lambda);
result.push(new_event);
}
result
}