gfpop Vignette
Vincent Runge
LaMME, Evry University
Quick Start
we present a basic use of the main functions of the
gfpop package. More details about optional arguments are
given in later sections.
We install the package from Github:
#devtools::install_github("vrunge/gfpop")
library(gfpop)We simulate some univariate gaussian data (n = 1000
points) with relative change-point positions
0.1, 0.3, 0.5, 0.8, 1 and means 1, 2, 1, 3, 1
with a variance equal to 1.
n <- 1000
myData <- dataGenerator(n, c(0.1,0.3,0.5,0.8,1), c(1,2,1,3,1), sigma = 1)We define the graph of constraints to use for the dynamic programming
algorithm. A simple case is the up-down constraint with a penalty here
equal to a classic 2 log(n).
myGraph <- graph(penalty = 2*log(n), type = "updown")The gfpop function gives the result of the segmentation using
myData and myGraph as parameters. We choose a
gaussian cost.
gfpop(data = myData, mygraph = myGraph, type = "mean")## $changepoints
## [1] 101 298 500 802 1000
## $states
## [1] "Dw" "Up" "Dw" "Up" "Dw"
## $forced
## [1] FALSE FALSE FALSE FALSE
## $parameters
## [1] 0.9342741 2.0140290 1.0729112 2.8885869 0.9897591
## $globalCost
## [1] 989.4949The vector changepoints gives the last index of each
segment. It always ends with the length of the vector
vectData.
The vector states contains the states in which lies each
mean. The length of this vector is the same as the length of
changepoint.
The vector forced is a boolean vector. A forced element
means that two consecutive means have been forced to satisfy the
constraint. For example, the “up” edge with parameter c is forced if
m(i+1) - m(i) = c.
The vector parameters contains the inferred
means/parameters of the successive segments.
The number globalCost is equal to the non-penalized
cost, that is the value of the fit to the data ignoring the penalties
for adding changes.
Some examples
Isotonic regression
The isotonic regression infers a sequence of nondecreasing means.
n <- 1000
mydata <- dataGenerator(n, c(0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1), c(0, 0.5, 1, 1.5, 2, 2.5, 3), sigma = 1)
myGraphIso <- graph(penalty = 2*log(n), type = "isotonic")
gfpop(data = mydata, mygraph = myGraphIso, type = "mean")## $changepoints
## [1] 97 209 305 589 806 1000
## $states
## [1] "Iso" "Iso" "Iso" "Iso" "Iso" "Iso"
## $forced
## [1] FALSE FALSE FALSE FALSE FALSE
## $parameters
## [1] -0.1763525 0.6378637 1.1915420 1.8010801 2.4938585 2.9264333
## $globalCost
## [1] 1022.422In this example, we use in gfpop function a robust
biweight gaussian cost with K = 1 and the min
parameter in order to infer means greater than 0.5.
Fixed number of change-points
This algorithm is called segment neighborhood in the change-point litterature. In this example, we fixed the number of segments at 3 with an isotonic constraint. The graph contains two “up” edges with no cycling.
n <- 1000
mydata <- dataGenerator(n, c(0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1), c(0, 0.5, 1, 1.5, 2, 2.5, 3), sigma = 1)
beta <- 0
myGraph <- graph(
Edge(0, 1,"up", beta),
Edge(1, 2, "up", beta),
Edge(0, 0, "null"),
Edge(1, 1, "null"),
Edge(2, 2, "null"),
StartEnd(start = 0, end = 2))
gfpop(data = mydata, mygraph = myGraph, type = "mean")## $changepoints
## [1] 299 600 1000
## $states
## [1] "0" "1" "2"
## $forced
## [1] FALSE FALSE
## $parameters
## [1] 0.4612321 1.7358436 2.7450942
## $globalCost
## [1] 1028.186Robust up-down with constrained starting and ending states
In presence of outliers we need a robust loss (biweight). We can also
force the starting and ending state and a minimal gap between the means
(here equal to 1)
n <- 1000
chgtpt <- c(0.1, 0.3, 0.5, 0.8, 1)
myData <- dataGenerator(n, chgtpt, c(0, 1, 0, 1, 0), sigma = 1)
myData <- myData + 5 * rbinom(n, 1, 0.05) - 5 * rbinom(n, 1, 0.05)
beta <- 2 * log(n)
myGraph <- graph(
Edge("Dw", "Up", type = "up", penalty = beta, gap = 1, K = 3),
Edge("Up", "Dw", type = "down", penalty = beta, gap = 1, K = 3),
Edge("Dw", "Dw", type = "null", K = 3),
Edge("Up", "Up", type = "null", K = 3),
StartEnd(start = "Dw", end = "Dw"))
gfpop(data = myData, mygraph = myGraph, type = "mean")## $changepoints
## [1] 113 306 503 791 1000
## $states
## [1] "Dw" "Up" "Dw" "Up" "Dw"
## $forced
## [1] TRUE FALSE FALSE FALSE
## $parameters
## [1] -0.03001725 0.96998275 -0.09252292 0.99072932 -0.04796669
## $globalCost
## [1] 1016.062If we skip all these constraints and use a standard fpop algorithm, the result is the following
myGraphStd <- graph(penalty = 2*log(n), type = "std")
gfpop(data = myData, mygraph = myGraphStd, type = "mean")## $changepoints
## [1] 5 6 42 46 47 48 49 62 63 108 109 161 162 163 184
## [16] 185 187 188 197 198 200 201 228 229 278 279 298 299 317 318
## [31] 326 327 418 419 443 444 453 454 458 459 485 486 489 490 498
## [46] 499 500 595 603 604 619 620 743 744 752 753 771 772 796 799
## [61] 800 806 807 847 848 850 851 900 901 911 916 920 921 945 946
## [76] 962 963 976 977 996 998 1000
## $states
## [1] "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std"
## [13] "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std"
## [25] "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std"
## [37] "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std"
## [49] "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std"
## [61] "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std"
## [73] "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std" "Std"
## $forced
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [37] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [49] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## $parameters
## [1] -0.24503906 -5.86832380 -0.04982604 -2.81924277 -7.07086760 0.36683911
## [7] 6.95296775 0.81881151 5.12154107 0.04565478 5.66905393 1.03151283
## [13] 5.69380520 -3.39387927 0.80573467 -4.77933344 0.91332395 6.78488480
## [19] 0.59452694 6.13004046 -0.50462682 6.05301009 1.21195155 7.20070280
## [25] 1.11291690 6.79539005 0.99488775 6.83262582 0.62526582 5.39814597
## [31] -0.73813191 -6.34028241 -0.15513629 5.69003521 0.17941235 -6.45128179
## [37] -0.19033582 5.48282573 -0.15618268 6.63888341 0.19416319 5.70552418
## [43] -0.07552765 5.68190887 -0.53834613 4.00510048 -4.79009728 0.91481021
## [49] 2.74483747 -3.54019118 1.27241946 -5.56196498 0.85582942 7.89726278
## [55] 0.80639359 -4.29399043 1.39367563 -4.18297880 0.80303480 -2.92049606
## [61] 3.15154055 -0.86770333 5.91856050 -0.14788897 -5.94575811 -0.26558236
## [67] -5.67045121 -0.04160088 -5.99451647 1.43908894 -2.49351589 1.18003684
## [73] -5.90374958 0.14758700 -7.13576991 -0.04901286 5.91494910 -0.35799947
## [79] 5.10114009 0.07146912 -4.22583843 0.01397304
## $globalCost
## [1] 1374.186abs edge
With a unique "abs" edge, we impose a difference between
the means of size at least 1.
n <- 10000
myData <- dataGenerator(n, c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1), c(0, 1, 0, 2, 1, 2, 0, 1, 0, 1), sigma = 0.5)
beta <- 2*log(n)
myGraph <- graph(
Edge(0, 0,"abs", penalty = beta, gap = 1),
Edge(0, 0,"null"))
gfpop(data = myData, mygraph = myGraph, type = "mean")## $changepoints
## [1] 1002 2000 3000 4000 5000 6000 7001 7999 8998 10000
## $states
## [1] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
## $forced
## [1] TRUE FALSE FALSE TRUE FALSE FALSE FALSE TRUE TRUE
## $parameters
## [1] -0.003358477 0.996641523 -0.022999674 1.988818878 0.988818878
## [6] 2.044028631 -0.032402659 1.007058356 0.007058356 1.007058356
## $globalCost
## [1] 2468.456Notice that some of the edges are forced, the vector
forced contains non-zero values.
Exponential decay
The null edge corresponds to an exponential decay state if its parameter is not equal to 1.
n <- 1000
mydata <- dataGenerator(n, c(0.2, 0.5, 0.8, 1), c(5, 10, 15, 20), sigma = 1, gamma = 0.966)
beta <- 2*log(n)
myGraphDecay <- graph(
Edge(0, 0, "up", penalty = beta),
Edge(0, 0, "null", 0, decay = 0.966)
)
g <- gfpop(data = mydata, mygraph = myGraphDecay, type = "mean")
g## $changepoints
## [1] 200 500 800 1000
## $states
## [1] "0" "0" "0" "0"
## $forced
## [1] FALSE FALSE FALSE
## $parameters
## [1] 0.0052304010 0.0003113861 0.0004739095 0.0195570079
## $globalCost
## [1] 971.721and we plot the result
gamma <- 0.966
len <- diff(c(0, g$changepoints))
signal <- NULL
for(i in length(len):1)
{signal <- c(signal, g$parameters[i]*c(1, cumprod(rep(1/gamma,len[i]-1))))}
signal <- rev(signal)
ylimits <- c(min(mydata), max(mydata))
plot(mydata, type ='p', pch ='+', ylim = ylimits)
par(new = TRUE)
plot(signal, type ='l', col = 4, ylim = ylimits, lwd = 3)
Graph construction
In the gfpop package, graphs are represented by a
dataframe with 9 features and build with the R functions
Edge, Node, StartEnd and
graph.
emptyGraph <- graph()
emptyGraph## [1] state1 state2 type parameter penalty K a
## [8] min max
## <0 lignes> (ou 'row.names' de longueur nulle)state1 is the starting node of an edge,
state2 its ending node. type is one of the
available edge type ("null", "std",
"up", "down", "abs").
penalty is a nonnegative parameter: the additional cost
\(\beta_i\) to consider when we move
within the graph using a edge (or stay on the same node).
parameter is annother nonnegative parameter, a
characteristics of the edge, depending of its type (it is a decay if
type is “null” and a gap otherwise). K and a
are robust parameters. min and max are used to
constrain the rang of value for the node parameter.
We add edges into a graph as follows
myGraph <- graph(
Edge("E1", "E1", "null"),
Edge("E1", "E2", "down", 3.1415, gap = 1.5)
)
myGraph## state1 state2 type parameter penalty K a min max
## 1 E1 E1 null 1.0 0 Inf 0 NA NA
## 2 E1 E2 down 1.5 0 Inf 0 NA NAwe can only add edges to this dataframe using the object
Edge.
The graph can contain information on the starting and/or ending edge
to use with the StartEnd function.
beta <- 2 * log(1000)
myGraph <- graph(
Edge("Dw", "Dw", "null"),
Edge("Up", "Up", "null"),
Edge("Dw", "Up", "up", penalty = beta, gap = 1),
Edge("Dw", "Dw", "down", penalty = beta),
Edge("Up", "Dw", "down", penalty = beta),
StartEnd(start = "Dw", end = "Dw"))
myGraph## state1 state2 type parameter penalty K a min max
## 1 Dw Dw null 1 0.00000 Inf 0 NA NA
## 2 Up Up null 1 0.00000 Inf 0 NA NA
## 3 Dw Up up 1 13.81551 Inf 0 NA NA
## 4 Dw Dw down 0 13.81551 Inf 0 NA NA
## 5 Up Dw down 0 13.81551 Inf 0 NA NA
## 6 Dw <NA> start NA NA NA NA NA NA
## 7 Dw <NA> end NA NA NA NA NA NASome graphs are often used: they are defined by default in the
graph function. To use these graphs, we specify a string
type equal to "std", "isotonic",
"updown" or "relevant". For example,
myGraphIso <- graph(penalty = 12, type = "isotonic")
myGraphIso## state1 state2 type parameter penalty K a min max
## 1 Iso Iso null 1 0 Inf 0 NA NA
## 2 Iso Iso up 0 12 Inf 0 NA NAThe function Node can be used to restrict the range of
value for parameter associated to a node (called also a vertex). For
example the following graph is an isotonic graph with inferred
parameters between 0 et 1 only.
myGraph <- graph(
Edge("Up", "Up", "up", penalty = 3.1415),
Edge("Up", "Up"),
Node("Up", min = 0, max = 1)
)
myGraph## state1 state2 type parameter penalty K a min max
## 1 Up Up up 0 3.1415 Inf 0 NA NA
## 2 Up Up null 1 0.0000 Inf 0 NA NA
## 3 Up Up node NA NA NA NA 0 1Supplementary R functions
Data generator function
the dataGenerator function is used to simulate
n data-points from a distribution of type
equal to "mean", "poisson",
"exp", "variance" or "negbin".
Standard deviation parameter sigma and decay
gamma are specific to the Gaussian mean model.
size is linked to the R rnbinom function from
R stats package.
Standard deviation estimation
We often need to estimate the standard deviation from the observed
data to normalize the data or choose the edge penalties. The
sdDiff returns such an estimation with the default HALL
method [Hall et al., 1990] well suited for time series with
change-points.