AR1seg_func: Segmentation of an AR(1) Gaussian process
In AR1seg: Segmentation of an autoregressive Gaussian process of order 1
Description
Usage
Arguments
Value
Note
Author(s)
References
Examples
Description
This function consists in an implementation of a robust approach to solve the problem
of multiple change-point estimation in the mean of a Gaussian AR(1) process.
A robust estimator of the autoregression parameter is proposed and used to build a decorrelated series on which a classical penalized least-square approach is applied.
Usage
1
AR1seg_func(y, Kmax = 15, rho = TRUE)
Arguments
y
Vector of observations
Kmax
Maximal number of segments
rho
It corresponds to the autoregression parameter. If it is equal to TRUE then it is estimated using a robust approach, otherwise the user has to provide a numerical value.
By default, the value of rho is TRUE.
Value
Contains the following attributes:
data
Vector of observations
rho
The estimator of rho if the argument rho=TRUE, otherwise the value provided by the user
decorrelated
The decorrelated series using rho
breaks
Matrix of size Kmax*Kmax. The line K=1,...,Kmax corresponds to the optimal segmentation of the series with K segments. By convention, the last break of each line is the length of the series.
selected
Selected number of segments using the modified BIC criterion proposed by Zhang and Siegmund (2007)
SelectedBreaks
Optimal segmentation with a number of segments equal to the value selected
PPbreaks
Matrix of breaks obtained after the post-processing step
PPSelectedBreaks
Result of the post-processing step applied to SelectedBreaks: it is the resulting segmentation of our approach
PPselected
Length of the resulting segmentation (PPSelectedBreaks)
PPmean
Empirical mean of the series on each segment of the resulting segmentation
Note
This package depends on the package Segmentor3IsBack
Author(s)
S. Chakar, E. Lebarbier, C. Levy-Leduc, S. Robin
References
This function corresponds to the implementation of the robust approach for estimating
change-points in the mean of an AR(1) Gaussian process by using the methodology described in the paper arXiv:1403.1958
Examples
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
##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function (y, Kmax = 15, rho = TRUE)
{
l = length(y)
if (rho)
rho = median((diff(y, lag = 2))^2)/median(diff(y)^2) -
1
x = y[2:l] - rho * y[1:(l - 1)]
S = Segmentor(x, model = 2, Kmax = Kmax)
breaks = S@breaks
for (i in 1:Kmax) {
for (j in 1:i) breaks[i, j] = breaks[i, j] + 1
}
rm(i, j)
parameters = S@parameters
PP = function(t) {
x = t
l = length(x)
i = 2
while (l > 2 && i < l) {
if (x[i] == x[i - 1] + 1 && x[i] != x[i + 1] - 1) {
x = c(x[1:(i - 1)], x[(i + 1):l])
l = l - 1
}
else i = i + 1
}
if (l > 1 && x[l - 1] == x[l] - 1)
x = x[1:(l - 1)]
x
}
PPbreaks = matrix(0, nrow = Kmax, ncol = Kmax, dimnames = dimnames(breaks))
PPbreaks[1, ] = breaks[1, ]
for (i in 2:Kmax) {
t = PP(breaks[i, 1:(i - 1)])
PPbreaks[i, ] = c(t, l, rep(0, Kmax - length(t) - 1))
}
rm(i, t)
fMa = function(t, mu) {
M = c()
t = c(0, t)
for (i in 2:length(t)) {
M = c(M, rep(mu[i - 1], t[i] - t[i - 1]))
}
M
}
sswg = function(br, param, series) {
sum((series - fMa(br, param))^2)
}
sswgseg = function(seg, seri) {
res = c()
for (i in 1:(Kmax)) {
res = c(res, sswg(seg@breaks[i, 1:i], seg@parameters[i,
1:i], seri))
}
res
}
minushalflogB = function(t, u) {
t = t[t != 0]
l = length(t)
b = log(t[1]/u)
if (l > 1) {
for (i in 2:l) {
b = b + log(t[i] - t[i - 1])
}
}
b = -b/2
}
ZS = function(seg, seri) {
u = length(seg@data)
Kmax = seg@Kmax
f = function(t) minushalflogB(t, u)
wg = sswgseg(seg, seri)
criterion = -(((u + 1):(u - Kmax + 2))/2) * log(wg) +
lgamma(((u + 1):(u - Kmax + 2))/2) - (0:(Kmax - 1)) *
log(u) + apply(seg@breaks, 1, f)
selected = which.max(criterion)
selected
}
selected = ZS(S, x)
SelectedBreaks = breaks[selected, 1:selected]
PPSelectedBreaks = PPbreaks[selected, ]
PPSelectedBreaks = PPSelectedBreaks[PPSelectedBreaks != 0]
PPselected = length(PPSelectedBreaks)
vec1 = c(1, PPSelectedBreaks[1:(PPselected - 1)] + 1)
vec2 = PPSelectedBreaks[1:(PPselected)]
m = c()
for (i in 1:PPselected) {
m[i] = mean(y[vec1[i]:vec2[i]])
}
list(data = y, rho = rho, decorrelated = x, breaks = breaks,
PPbreaks = PPbreaks, selected = selected, SelectedBreaks = SelectedBreaks,
PPSelectedBreaks = PPSelectedBreaks, PPselected = PPselected, PPmean = m)
}
Example output
Loading required package: Segmentor3IsBack
Segmentor3IsBack v2.0 Loaded
function (y, Kmax = 15, rho = TRUE)
{
l = length(y)
if (rho)
rho = median((diff(y, lag = 2))^2)/median(diff(y)^2) -
1
x = y[2:l] - rho * y[1:(l - 1)]
S = Segmentor(x, model = 2, Kmax = Kmax)
breaks = S@breaks
for (i in 1:Kmax) {
for (j in 1:i) breaks[i, j] = breaks[i, j] + 1
}
rm(i, j)
parameters = S@parameters
PP = function(t) {
x = t
l = length(x)
i = 2
while (l > 2 && i < l) {
if (x[i] == x[i - 1] + 1 && x[i] != x[i + 1] - 1) {
x = c(x[1:(i - 1)], x[(i + 1):l])
l = l - 1
}
else i = i + 1
}
if (l > 1 && x[l - 1] == x[l] - 1)
x = x[1:(l - 1)]
x
}
PPbreaks = matrix(0, nrow = Kmax, ncol = Kmax, dimnames = dimnames(breaks))
PPbreaks[1, ] = breaks[1, ]
for (i in 2:Kmax) {
t = PP(breaks[i, 1:(i - 1)])
PPbreaks[i, ] = c(t, l, rep(0, Kmax - length(t) - 1))
}
rm(i, t)
fMa = function(t, mu) {
M = c()
t = c(0, t)
for (i in 2:length(t)) {
M = c(M, rep(mu[i - 1], t[i] - t[i - 1]))
}
M
}
sswg = function(br, param, series) {
sum((series - fMa(br, param))^2)
}
sswgseg = function(seg, seri) {
res = c()
for (i in 1:(Kmax)) {
res = c(res, sswg(seg@breaks[i, 1:i], seg@parameters[i,
1:i], seri))
}
res
}
minushalflogB = function(t, u) {
t = t[t != 0]
l = length(t)
b = log(t[1]/u)
if (l > 1) {
for (i in 2:l) {
b = b + log(t[i] - t[i - 1])
}
}
b = -b/2
}
ZS = function(seg, seri) {
u = length(seg@data)
Kmax = seg@Kmax
f = function(t) minushalflogB(t, u)
wg = sswgseg(seg, seri)
criterion = -(((u + 1):(u - Kmax + 2))/2) * log(wg) +
lgamma(((u + 1):(u - Kmax + 2))/2) - (0:(Kmax - 1)) *
log(u) + apply(seg@breaks, 1, f)
selected = which.max(criterion)
selected
}
selected = ZS(S, x)
SelectedBreaks = breaks[selected, 1:selected]
PPSelectedBreaks = PPbreaks[selected, ]
PPSelectedBreaks = PPSelectedBreaks[PPSelectedBreaks != 0]
PPselected = length(PPSelectedBreaks)
vec1 = c(1, PPSelectedBreaks[1:(PPselected - 1)] + 1)
vec2 = PPSelectedBreaks[1:(PPselected)]
m = c()
for (i in 1:PPselected) {
m[i] = mean(y[vec1[i]:vec2[i]])
}
list(data = y, rho = rho, decorrelated = x, breaks = breaks,
PPbreaks = PPbreaks, selected = selected, SelectedBreaks = SelectedBreaks,
PPSelectedBreaks = PPSelectedBreaks, PPselected = PPselected,
PPmean = m)
}
AR1seg documentation built on May 2, 2019, 9:51 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Value Note Author(s) References Examples
Description
This function consists in an implementation of a robust approach to solve the problem of multiple change-point estimation in the mean of a Gaussian AR(1) process. A robust estimator of the autoregression parameter is proposed and used to build a decorrelated series on which a classical penalized least-square approach is applied.
Usage
1 | AR1seg_func(y, Kmax = 15, rho = TRUE)
|
Arguments
y |
Vector of observations |
Kmax |
Maximal number of segments |
rho |
It corresponds to the autoregression parameter. If it is equal to TRUE then it is estimated using a robust approach, otherwise the user has to provide a numerical value. By default, the value of rho is TRUE. |
Value
Contains the following attributes:
data |
Vector of observations |
rho |
The estimator of rho if the argument rho=TRUE, otherwise the value provided by the user |
decorrelated |
The decorrelated series using rho |
breaks |
Matrix of size Kmax*Kmax. The line K=1,...,Kmax corresponds to the optimal segmentation of the series with K segments. By convention, the last break of each line is the length of the series. |
selected |
Selected number of segments using the modified BIC criterion proposed by Zhang and Siegmund (2007) |
SelectedBreaks |
Optimal segmentation with a number of segments equal to the value selected |
PPbreaks |
Matrix of breaks obtained after the post-processing step |
PPSelectedBreaks |
Result of the post-processing step applied to SelectedBreaks: it is the resulting segmentation of our approach |
PPselected |
Length of the resulting segmentation (PPSelectedBreaks) |
PPmean |
Empirical mean of the series on each segment of the resulting segmentation |
Note
This package depends on the package Segmentor3IsBack
Author(s)
S. Chakar, E. Lebarbier, C. Levy-Leduc, S. Robin
References
This function corresponds to the implementation of the robust approach for estimating change-points in the mean of an AR(1) Gaussian process by using the methodology described in the paper arXiv:1403.1958
Examples
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 | ##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function (y, Kmax = 15, rho = TRUE)
{
l = length(y)
if (rho)
rho = median((diff(y, lag = 2))^2)/median(diff(y)^2) -
1
x = y[2:l] - rho * y[1:(l - 1)]
S = Segmentor(x, model = 2, Kmax = Kmax)
breaks = S@breaks
for (i in 1:Kmax) {
for (j in 1:i) breaks[i, j] = breaks[i, j] + 1
}
rm(i, j)
parameters = S@parameters
PP = function(t) {
x = t
l = length(x)
i = 2
while (l > 2 && i < l) {
if (x[i] == x[i - 1] + 1 && x[i] != x[i + 1] - 1) {
x = c(x[1:(i - 1)], x[(i + 1):l])
l = l - 1
}
else i = i + 1
}
if (l > 1 && x[l - 1] == x[l] - 1)
x = x[1:(l - 1)]
x
}
PPbreaks = matrix(0, nrow = Kmax, ncol = Kmax, dimnames = dimnames(breaks))
PPbreaks[1, ] = breaks[1, ]
for (i in 2:Kmax) {
t = PP(breaks[i, 1:(i - 1)])
PPbreaks[i, ] = c(t, l, rep(0, Kmax - length(t) - 1))
}
rm(i, t)
fMa = function(t, mu) {
M = c()
t = c(0, t)
for (i in 2:length(t)) {
M = c(M, rep(mu[i - 1], t[i] - t[i - 1]))
}
M
}
sswg = function(br, param, series) {
sum((series - fMa(br, param))^2)
}
sswgseg = function(seg, seri) {
res = c()
for (i in 1:(Kmax)) {
res = c(res, sswg(seg@breaks[i, 1:i], seg@parameters[i,
1:i], seri))
}
res
}
minushalflogB = function(t, u) {
t = t[t != 0]
l = length(t)
b = log(t[1]/u)
if (l > 1) {
for (i in 2:l) {
b = b + log(t[i] - t[i - 1])
}
}
b = -b/2
}
ZS = function(seg, seri) {
u = length(seg@data)
Kmax = seg@Kmax
f = function(t) minushalflogB(t, u)
wg = sswgseg(seg, seri)
criterion = -(((u + 1):(u - Kmax + 2))/2) * log(wg) +
lgamma(((u + 1):(u - Kmax + 2))/2) - (0:(Kmax - 1)) *
log(u) + apply(seg@breaks, 1, f)
selected = which.max(criterion)
selected
}
selected = ZS(S, x)
SelectedBreaks = breaks[selected, 1:selected]
PPSelectedBreaks = PPbreaks[selected, ]
PPSelectedBreaks = PPSelectedBreaks[PPSelectedBreaks != 0]
PPselected = length(PPSelectedBreaks)
vec1 = c(1, PPSelectedBreaks[1:(PPselected - 1)] + 1)
vec2 = PPSelectedBreaks[1:(PPselected)]
m = c()
for (i in 1:PPselected) {
m[i] = mean(y[vec1[i]:vec2[i]])
}
list(data = y, rho = rho, decorrelated = x, breaks = breaks,
PPbreaks = PPbreaks, selected = selected, SelectedBreaks = SelectedBreaks,
PPSelectedBreaks = PPSelectedBreaks, PPselected = PPselected, PPmean = m)
}
|
Example output
Loading required package: Segmentor3IsBack
Segmentor3IsBack v2.0 Loaded
function (y, Kmax = 15, rho = TRUE)
{
l = length(y)
if (rho)
rho = median((diff(y, lag = 2))^2)/median(diff(y)^2) -
1
x = y[2:l] - rho * y[1:(l - 1)]
S = Segmentor(x, model = 2, Kmax = Kmax)
breaks = S@breaks
for (i in 1:Kmax) {
for (j in 1:i) breaks[i, j] = breaks[i, j] + 1
}
rm(i, j)
parameters = S@parameters
PP = function(t) {
x = t
l = length(x)
i = 2
while (l > 2 && i < l) {
if (x[i] == x[i - 1] + 1 && x[i] != x[i + 1] - 1) {
x = c(x[1:(i - 1)], x[(i + 1):l])
l = l - 1
}
else i = i + 1
}
if (l > 1 && x[l - 1] == x[l] - 1)
x = x[1:(l - 1)]
x
}
PPbreaks = matrix(0, nrow = Kmax, ncol = Kmax, dimnames = dimnames(breaks))
PPbreaks[1, ] = breaks[1, ]
for (i in 2:Kmax) {
t = PP(breaks[i, 1:(i - 1)])
PPbreaks[i, ] = c(t, l, rep(0, Kmax - length(t) - 1))
}
rm(i, t)
fMa = function(t, mu) {
M = c()
t = c(0, t)
for (i in 2:length(t)) {
M = c(M, rep(mu[i - 1], t[i] - t[i - 1]))
}
M
}
sswg = function(br, param, series) {
sum((series - fMa(br, param))^2)
}
sswgseg = function(seg, seri) {
res = c()
for (i in 1:(Kmax)) {
res = c(res, sswg(seg@breaks[i, 1:i], seg@parameters[i,
1:i], seri))
}
res
}
minushalflogB = function(t, u) {
t = t[t != 0]
l = length(t)
b = log(t[1]/u)
if (l > 1) {
for (i in 2:l) {
b = b + log(t[i] - t[i - 1])
}
}
b = -b/2
}
ZS = function(seg, seri) {
u = length(seg@data)
Kmax = seg@Kmax
f = function(t) minushalflogB(t, u)
wg = sswgseg(seg, seri)
criterion = -(((u + 1):(u - Kmax + 2))/2) * log(wg) +
lgamma(((u + 1):(u - Kmax + 2))/2) - (0:(Kmax - 1)) *
log(u) + apply(seg@breaks, 1, f)
selected = which.max(criterion)
selected
}
selected = ZS(S, x)
SelectedBreaks = breaks[selected, 1:selected]
PPSelectedBreaks = PPbreaks[selected, ]
PPSelectedBreaks = PPSelectedBreaks[PPSelectedBreaks != 0]
PPselected = length(PPSelectedBreaks)
vec1 = c(1, PPSelectedBreaks[1:(PPselected - 1)] + 1)
vec2 = PPSelectedBreaks[1:(PPselected)]
m = c()
for (i in 1:PPselected) {
m[i] = mean(y[vec1[i]:vec2[i]])
}
list(data = y, rho = rho, decorrelated = x, breaks = breaks,
PPbreaks = PPbreaks, selected = selected, SelectedBreaks = SelectedBreaks,
PPSelectedBreaks = PPSelectedBreaks, PPselected = PPselected,
PPmean = m)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.