/
vignettes.Rmd
189 lines (136 loc) · 7.77 KB
/
vignettes.Rmd
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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
---
title: "`synimu` Vignettes"
bibliography: biblio.bib
output:
rmarkdown::html_vignette:
fig_caption: yes
vignette: >
%\VignetteIndexEntry{Simulations from Time Series Models}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
In @zhang2021scalewise, we introduce a method allowing to compute synthetic (or average) signals issued from multiple correlated and/or even non-stationary signals. In a nutshell, this method minimizes a weighted average of the wavelet variance of the synthetic signal whose weights can be chosen depending on the application at hand. This method is illustrated with gyroscopes which often contain non-stationary processes such as random walks. This package implements the method of @zhang2021scalewise, called Scale-wise Variance Optimization or SVO, as well as the approach of @s17020352, which is used as a comparison.
# Real Gyroscope Data Analysis
We consider the real gyroscope data considered in @zhang2021scalewise, where a virtual gyroscope is created by fusing 12 real MEMS from two device families, InvenSense
ICM-20689 and Bosch BMI055. The employed data is included in the package and can be loaded as follows:
```{r}
library(synimu)
data(BMI055)
data(ICM20689)
Xt = cbind(ICM20689, BMI055)
```
We first need to choose which scales we would like to minimize in the wavelet variance plot. In the SVO proposed in @zhang2021scalewise, this is obtained by choosing non-negative weights on the scales that sum to one. For example, if one is interested in minimizing the large scales of the wavelet variance, suitable weights can be obtained as follows:
```{r, fig.width = 7}
n_scales = log2( nrow(Xt) ) - 1
tau = 1:n_scales
weights = exp(tau-15)/(1 + exp(tau-15))
weights = weights/sum(weights)
plot(NA, xlim = range(tau), ylim = range(weights), xlab = "Wavelet variance scales", ylab = "Weight", xaxt="n")
p10s = seq(0,log10(2^(tail(tau,1))))
abline(v = log2(10^p10s), col = "grey", lt = 2 )
abline(h = seq(0,0.15,0.05), col = "grey", lt = 2 )
lines(tau, weights, lwd = 1.3)
# plot x axis labels in samples (powers of 10)
axis(1, at = log2(10^p10s), labels = sapply(p10s, function(i) as.expression(bquote(10^ .(i)))))
```
Once the weights on scales have been specified, one can obtain the coefficients to be used in the linear combination of single gyroscope signals as follows:
```{r, cache = TRUE}
c_hat = find_optimal_coefs(Xt, weights)
c_hat
```
The obtained virtual gyroscope can be compared with what would be obtained with equal coefficients using the "plot_virtual_gyro" function:
```{r, fig.width = 7, fig.height = 10}
plot_virtual_gyro(Xt, c_hat, names = "SVO")
```
# Simulation Studies
In this section we will provide the instructions on how to reproduce the Monte-Carlo simulations presented in @zhang2021scalewise.
## White Noise + Random Walk
In this section we simulate data for six gyroscopes whose noise error model is composed by a white noise and a random walk (with correlated innovations), as in Equation (11) in @zhang2021scalewise.
```{r}
# definition of the noise parameters
sigma_wn = c(2.805556e-07, 1.977778e-07, 1.361111e-07, 1.063889e-07, 1.069444e-07, 2.538889e-07)
sigma_rw = matrix(c(
2.550583e-14, -8.573388e-16, 1.028807e-14, -1.628944e-14, -2.400549e-14, -5.572702e-15,
-8.573388e-16, 4.715364e-14, 1.993313e-14, 1.071674e-15, 3.215021e-15, 2.079047e-14,
1.028807e-14, 1.993313e-14, 3.489369e-13, -1.281722e-13, 5.572702e-15, -3.022119e-14,
-1.628944e-14, 1.071674e-15, -1.281722e-13, 2.091907e-13, 4.093793e-14, 5.444102e-14,
-2.400549e-14, 3.215021e-15, 5.572702e-15, 4.093793e-14, 5.551269e-14, 2.443416e-14,
-5.572702e-15, 2.079047e-14, -3.022119e-14, 5.444102e-14, 2.443416e-14, 2.501286e-13
), nrow = 6, byrow = T)
N = 2^20
set.seed(1)
Xt = simulate_wn(N, sigma_wn) + simulate_corr_rw(N, sigma_rw)
```
We define weights in order to minimize the largest scales of the wavelet variance and we compute the coefficients with the approach of SVO:
```{r, cache = TRUE}
n_scales = log2( N ) - 1
tau = 1:n_scales
weights = exp(tau-13)/(1 + exp(tau-13))
weights = weights/sum(weights)
c_svo = find_optimal_coefs(Xt, weights)
c_svo
```
As a comparison, we also compute the optimal coefficients according to the method presented in @s17020352, which can be calculated as follows:
```{r, cache = TRUE}
# Estimate the Q matrix (the covariance of the Random Walk innovations)
Q_hat = algorithm_2(Xt)$Q
# Compute the coefficients
c_vaccarozaki = find_optimal_coefs_vaccaro(Q_hat)
c_vaccarozaki
```
The results can be seen in the following plot. In this case, since the assumed model for the simulated gyroscopes satisfies the parametric assumptions in @s17020352, the results of both approaches are mostly similar.
```{r, fig.width = 7, fig.height = 10}
plot_virtual_gyro(Xt, c_svo, c_vaccarozaki, names = c("SVO", "VaccaroZaki"))
```
## White Noise + 3 correlated AR(1)
In this example we simulate data for six gyroscopes whose noise error model is composed by a white noise and three first-order auto-regressive processes (with correlated innovations), as in Equation (14) in @zhang2021scalewise.
```{r, cache = TRUE}
phi_ar1_1 = rep(0.9998705, 6)
sigma_ar1_1 = matrix(c(
6.713833e-13, 1.515552e-13, 1.245013e-13, -7.117931e-15, 4.284888e-14, 9.577566e-14,
1.515552e-13, 5.607230e-13, 2.264650e-13, -7.737235e-14, 2.810539e-14, 1.440987e-13,
1.245013e-13, 2.264650e-13, 7.499051e-13, -1.626870e-14, -5.293172e-14, -9.995748e-14,
-7.117931e-15, -7.737235e-14, -1.626870e-14, 1.207920e-13, -1.873054e-14, 6.817015e-14,
4.284888e-14, 2.810539e-14, -5.293172e-14, -1.873054e-14, 7.499327e-13, 2.414762e-13,
9.577566e-14, 1.440987e-13, -9.995748e-14, 6.817015e-14, 2.414762e-13, 6.461646e-13
), nrow = 6, byrow = T)
phi_ar1_2 = rep(0.9975214, 6)
sigma_ar1_2 = matrix(c(
4.008129e-12, -1.806110e-13, -4.477015e-13, -1.048567e-12, 2.546714e-12, -2.349504e-12,
-1.806110e-13, 1.177847e-11, -1.290506e-13, -3.741980e-12, -5.110108e-13, 2.849421e-12,
-4.477015e-13, -1.290506e-13, 1.965775e-11, 2.026288e-12, -1.750890e-12, 1.446290e-12,
-1.048567e-12, -3.741980e-12, 2.026288e-12, 1.309887e-11, 2.317561e-13, 4.127606e-12,
2.546714e-12, -5.110108e-13, -1.750890e-12, 2.317561e-13, 1.928375e-11, 5.714532e-12,
-2.349504e-12, 2.849421e-12, 1.446290e-12, 4.127606e-12, 5.714532e-12, 1.434466e-11
), nrow = 6, byrow = T)
phi_ar1_3 = rep(0.9999933, 6)
sigma_ar1_3 = matrix(c(
5.062006e-14, -6.432804e-15, 3.260125e-15, 6.977995e-15, 1.727949e-14, 7.385838e-15,
-6.432804e-15, 1.612038e-14, -3.184490e-15, -3.836343e-15, -7.682144e-15, -1.000440e-15,
3.260125e-15, -3.184490e-15, 2.102923e-14, 9.145839e-16, 3.985049e-15, -1.377599e-15,
6.977995e-15, -3.836343e-15, 9.145839e-16, 7.881771e-15, -1.291789e-15, 5.191186e-15,
1.727949e-14, -7.682144e-15, 3.985049e-15, -1.291789e-15, 5.442670e-14, 3.786316e-15,
7.385838e-15, -1.000440e-15, -1.377599e-15, 5.191186e-15, 3.786316e-15, 4.990004e-14
), nrow = 6, byrow = T)
set.seed(1)
Xt = simulate_wn(N, sigma_wn) +
simulate_corr_ar1(N, phi_ar1_1, sigma_ar1_1) +
simulate_corr_ar1(N, phi_ar1_2, sigma_ar1_2) +
simulate_corr_ar1(N, phi_ar1_3, sigma_ar1_3)
```
We again estimate the coefficients for the virtual gyroscopes with the SVO approach and the method proposed in @s17020352 and compare the results. In this case, SVO, which does not rely on any parametric assumption on the underlying stochastic processes, achieves better results.
```{r, cache = TRUE}
c_svo = find_optimal_coefs(Xt, weights)
c_svo
```
```{r, cache = TRUE}
# Estimate the Q matrix (the covariance of the Random Walk innovations)
Q_hat = algorithm_2(Xt)$Q
# Compute the coefficients
c_vaccarozaki = find_optimal_coefs_vaccaro(Q_hat)
c_vaccarozaki
```
```{r, fig.width = 7, fig.height = 10}
plot_virtual_gyro(Xt, c_svo, c_vaccarozaki, names = c("SVO", "VaccaroZaki"))
```
# Bibliography