forked from WinVector/RcppDynProg
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathSegmentationL.Rmd
More file actions
84 lines (69 loc) · 2.19 KB
/
Copy pathSegmentationL.Rmd
File metadata and controls
84 lines (69 loc) · 2.19 KB
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
---
title: "Linear Segmentation"
author: "John Mount"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Linear Segmentation}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
In this example we fit a piecewise linear function to example data.
Please see [here](https://github.com/WinVector/RcppDynProg) for a discussion of the methodology.
```{r r1, fig.height = 6, fig.width = 8, fig.align = "center"}
library("RcppDynProg")
set.seed(2018)
g <- 100
d <- data.frame(
x = 0.05*(1:(3*g))) # ordered in x
d$y_ideal <- sin((0.3*d$x)^2)
d$y_observed <- d$y_ideal + 0.25*rnorm(length(d$y_ideal))
# plot
plot(d$x, d$y_observed,
xlab = "x", ylab = "y",
main = "raw data\ncircles: observed values, dashed line: unobserved true values")
lines(d$x, d$y_ideal,
type = "l",
lty = "dashed")
x_cuts <- solve_for_partition(d$x, d$y_observed, penalty = 1)
print(x_cuts)
d$estimate <- approx(x_cuts$x, x_cuts$pred, xout = d$x, method = "linear", rule = 2)$y
d$group <- as.character(findInterval(d$x, x_cuts[x_cuts$what=="left", "x"]))
```
```{r r2, fig.height = 6, fig.width = 8, fig.align = "center"}
print(sum((d$y_observed - d$y_ideal)^2))
```
```{r r3, fig.height = 6, fig.width = 8, fig.align = "center"}
print(sum((d$estimate - d$y_ideal)^2))
```
```{r r4, fig.height = 6, fig.width = 8, fig.align = "center"}
print(sum((d$estimate - d$y_observed)^2))
```
```{r r5, fig.height = 6, fig.width = 8, fig.align = "center"}
# plot
plot(d$x, d$y_observed,
xlab = "x", ylab = "y",
main = "RcppDynProg piecewise linear estimate\ndots: observed values, segments: estimated shape")
points(d$x, d$y_ideal,
type = "l",
lty = "dashed")
cmap <- c("#a6cee3",
"#1f78b4",
"#b2df8a",
"#33a02c",
"#fb9a99",
"#e31a1c",
"#fdbf6f",
"#ff7f00",
"#cab2d6",
"#6a3d9a",
"#ffff99",
"#b15928")
names(cmap) <- as.character(seq_len(length(cmap)))
points(d$x, d$y_observed, col = cmap[d$group], pch=19)
groups <- sort(unique(d$group))
for(gi in groups) {
di <- d[d$group==gi, , drop = FALSE]
lines(di$x, di$estimate, col = cmap[di$group[[1]]], lwd=2)
}
```