knitr::opts_chunk$set(
collapse = TRUE, comment = "#>", eval = TRUE, cache = FALSE,
warning = FALSE, fig.width = 8, fig.height = 5
)
# Univariate mean change
set.seed(1)
p <- 1
mean_data_1 <- rbind(
mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
mvtnorm::rmvnorm(400, mean = rep(50, p), sigma = diag(100, p)),
mvtnorm::rmvnorm(300, mean = rep(2, p), sigma = diag(100, p))
)
plot.ts(mean_data_1)
# Univariate mean and/or variance change
set.seed(1)
p <- 1
mv_data_1 <- rbind(
mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
mvtnorm::rmvnorm(300, mean = rep(10, p), sigma = diag(100, p))
)
plot.ts(mv_data_1)
# Multivariate mean change
set.seed(1)
p <- 3
mean_data_3 <- rbind(
mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
mvtnorm::rmvnorm(400, mean = rep(50, p), sigma = diag(100, p)),
mvtnorm::rmvnorm(300, mean = rep(2, p), sigma = diag(100, p))
)
plot.ts(mean_data_3)
# Multivariate mean and/or variance change
set.seed(1)
p <- 4
mv_data_3 <- rbind(
mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
mvtnorm::rmvnorm(300, mean = rep(10, p), sigma = diag(100, p))
)
plot.ts(mv_data_3)
# Linear regression
set.seed(1)
n <- 300
p <- 4
x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
theta_0 <- rbind(c(1, 3.2, -1, 0), c(-1, -0.5, 2.5, -2), c(0.8, 0, 1, 2))
y <- c(
x[1:100, ] %*% theta_0[1, ] + rnorm(100, 0, 3),
x[101:200, ] %*% theta_0[2, ] + rnorm(100, 0, 3),
x[201:n, ] %*% theta_0[3, ] + rnorm(100, 0, 3)
)
lm_data <- data.frame(y = y, x = x)
plot.ts(lm_data)
# Logistic regression
set.seed(1)
n <- 500
p <- 4
x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
theta <- rbind(rnorm(p, 0, 1), rnorm(p, 2, 1))
y <- c(
rbinom(300, 1, 1 / (1 + exp(-x[1:300, ] %*% theta[1, ]))),
rbinom(200, 1, 1 / (1 + exp(-x[301:n, ] %*% theta[2, ])))
)
binomial_data <- data.frame(y = y, x = x)
plot.ts(binomial_data)
# Poisson regression
set.seed(1)
n <- 1100
p <- 3
x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
delta <- rnorm(p)
theta_0 <- c(1, 0.3, -1)
y <- c(
rpois(500, exp(x[1:500, ] %*% theta_0)),
rpois(300, exp(x[501:800, ] %*% (theta_0 + delta))),
rpois(200, exp(x[801:1000, ] %*% theta_0)),
rpois(100, exp(x[1001:1100, ] %*% (theta_0 - delta)))
)
poisson_data <- data.frame(y = y, x = x)
plot.ts(log(poisson_data$y))
plot.ts(poisson_data[, -1])
# Lasso
set.seed(1)
n <- 480
p_true <- 6
p <- 50
x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
theta_0 <- rbind(
runif(p_true, -5, -2),
runif(p_true, -3, 3),
runif(p_true, 2, 5),
runif(p_true, -5, 5)
)
theta_0 <- cbind(theta_0, matrix(0, ncol = p - p_true, nrow = 4))
y <- c(
x[1:80, ] %*% theta_0[1, ] + rnorm(80, 0, 1),
x[81:200, ] %*% theta_0[2, ] + rnorm(120, 0, 1),
x[201:320, ] %*% theta_0[3, ] + rnorm(120, 0, 1),
x[321:n, ] %*% theta_0[4, ] + rnorm(160, 0, 1)
)
lasso_data <- data.frame(y = y, x = x)
plot.ts(lasso_data[, seq_len(p_true + 1)])
# AR(3)
set.seed(1)
n <- 1000
x <- rep(0, n + 3)
for (i in 1:600) {
x[i + 3] <- 0.6 * x[i + 2] - 0.2 * x[i + 1] + 0.1 * x[i] + rnorm(1, 0, 3)
}
for (i in 601:1000) {
x[i + 3] <- 0.3 * x[i + 2] + 0.4 * x[i + 1] + 0.2 * x[i] + rnorm(1, 0, 3)
}
ar_data <- x[-seq_len(3)]
plot.ts(ar_data)
# GARCH(1, 1)
set.seed(1)
n <- 400
sigma_2 <- rep(1, n + 1)
x <- rep(0, n + 1)
for (i in seq_len(200)) {
sigma_2[i + 1] <- 20 + 0.5 * x[i]^2 + 0.1 * sigma_2[i]
x[i + 1] <- rnorm(1, 0, sqrt(sigma_2[i + 1]))
}
for (i in 201:400) {
sigma_2[i + 1] <- 1 + 0.1 * x[i]^2 + 0.5 * sigma_2[i]
x[i + 1] <- rnorm(1, 0, sqrt(sigma_2[i + 1]))
}
garch_data <- x[-1]
plot.ts(garch_data)
# VAR(2)
set.seed(1)
n <- 800
p <- 2
theta_1 <- matrix(c(-0.3, 0.6, -0.5, 0.4, 0.2, 0.2, 0.2, -0.2), nrow = p)
theta_2 <- matrix(c(0.3, -0.4, 0.1, -0.5, -0.5, -0.2, -0.5, 0.2), nrow = p)
x <- matrix(0, n + 2, p)
for (i in 1:500) {
x[i + 2, ] <- theta_1 %*% c(x[i + 1, ], x[i, ]) + rnorm(p, 0, 1)
}
for (i in 501:n) {
x[i + 2, ] <- theta_2 %*% c(x[i + 1, ], x[i, ]) + rnorm(p, 0, 1)
}
var_data <- x[-seq_len(2), ]
plot.ts(var_data)
(fastcpd_result <- fastcpd::fastcpd.mean(mean_data_1, r.progress = FALSE)@cp_set)
testthat::expect_equal(fastcpd_result, c(300, 700), tolerance = 0.2)
(CptNonPar_result <- CptNonPar::np.mojo(mean_data_1, G = floor(length(mean_data_1) / 6))$cpts)
testthat::expect_equal(CptNonPar_result, c(300, 700), tolerance = 0.2)
strucchange::breakpoints(y ~ 1, data = data.frame(y = mean_data_1))$breakpoints
#> [1] 300 700
ecp::e.divisive(mean_data_1)$estimates
#> [1] 1 301 701 1001
(changepoint_result <- changepoint::cpt.mean(c(mean_data_1))@cpts)
testthat::expect_equal(changepoint_result, c(300, 1000), tolerance = 0.2)
(breakfast_result <- breakfast::breakfast(mean_data_1)$cptmodel.list[[6]]$cpts)
testthat::expect_equal(breakfast_result, c(300, 700), tolerance = 0.2)
(wbs_result <- wbs::wbs(mean_data_1)$cpt$cpt.ic$mbic.penalty)
testthat::expect_equal(wbs_result, c(300, 700), tolerance = 0.2)
mosum::mosum(c(mean_data_1), G = 40)$cpts.info$cpts
#> [1] 300 700
(fpop_result <- fpop::Fpop(mean_data_1, nrow(mean_data_1))$t.est)
testthat::expect_equal(fpop_result, c(300, 700, 1000), tolerance = 0.2)
gfpop::gfpop(
data = mean_data_1,
mygraph = gfpop::graph(
penalty = 2 * log(nrow(mean_data_1)) * gfpop::sdDiff(mean_data_1) ^ 2,
type = "updown"
),
type = "mean"
)$changepoints
#> [1] 300 700 1000
invisible(
suppressMessages(
capture.output(
result_InspectChangepoint <- InspectChangepoint::inspect(
t(mean_data_1),
threshold = InspectChangepoint::compute.threshold(
nrow(mean_data_1), ncol(mean_data_1)
)
)
)
)
)
result_InspectChangepoint$changepoints[, "location"]
#> [1] 300 700
(jointseg_result <- jointseg::jointSeg(mean_data_1, K = 2)$bestBkp)
testthat::expect_equal(jointseg_result, c(300, 700), tolerance = 0.2)
Rbeast::beast(
mean_data_1, season = "none", print.progress = FALSE, quiet = TRUE
)$trend$cp
#> [1] 701 301 NaN NaN NaN NaN NaN NaN NaN NaN
(stepR_result <- stepR::stepFit(mean_data_1, alpha = 0.5)$rightEnd)
testthat::expect_equal(stepR_result, c(300, 700, 1000), tolerance = 0.2)
(cpm_result <- cpm::processStream(mean_data_1, cpmType = "Student")$changePoints)
testthat::expect_equal(cpm_result, c(299, 699), tolerance = 0.2)
(segmented_result <- segmented::stepmented(
as.numeric(mean_data_1), npsi = 2
)$psi[, "Est."])
testthat::expect_equivalent(segmented_result, c(298, 699), tolerance = 0.2)
plot(
mcp::mcp(
list(y ~ 1, ~ 1, ~ 1),
data = data.frame(y = mean_data_1, x = seq_len(nrow(mean_data_1))),
par_x = "x"
)
)
plot(not::not(mean_data_1, contrast = "pcwsConstMean"))
plot(bcp::bcp(mean_data_1))
(fastcpd_result <- fastcpd::fastcpd.mv(mv_data_1, r.progress = FALSE)@cp_set)
testthat::expect_equal(fastcpd_result, c(300, 700, 1001, 1300, 1700), tolerance = 0.2)
ecp::e.divisive(mv_data_1)$estimates
#> [1] 1 301 701 1001 1301 1701 2001
(changepoint_result <- changepoint::cpt.meanvar(c(mv_data_1))@cpts)
testthat::expect_equal(changepoint_result, c(300, 2000), tolerance = 0.2)
(CptNonPar_result <- CptNonPar::np.mojo(mv_data_1, G = floor(length(mv_data_1) / 6))$cpts)
testthat::expect_equal(CptNonPar_result, c(333, 700, 1300), tolerance = 0.2)
(cpm_result <- cpm::processStream(mv_data_1, cpmType = "GLR")$changePoints)
testthat::expect_equal(cpm_result, c(293, 300, 403, 408, 618, 621, 696, 1000, 1021, 1024, 1293, 1300, 1417, 1693, 1700, 1981), tolerance = 0.2)
invisible(
suppressMessages(
capture.output(
result_InspectChangepoint <- InspectChangepoint::inspect(
t(mv_data_1),
threshold = InspectChangepoint::compute.threshold(
nrow(mv_data_1), ncol(mv_data_1)
)
)
)
)
)
result_InspectChangepoint$changepoints[, "location"]
#> [1] 300 700 701 702 704 707 708 712 715 716 717 718 721 722
#> [15] 723 726 727 729 731 732 734 736 740 742 744 746 748 750
#> [29] 753 755 756 757 759 760 762 764 765 768 769 771 772 774
#> [43] 776 777 784 785 786 789 791 792 794 797 798 799 801 802
#> [57] 803 807 809 810 813 815 817 819 826 827 828 829 831 833
#> [71] 835 836 837 838 840 841 842 843 845 848 849 852 854 860
#> [85] 862 864 866 868 870 872 875 879 881 884 886 887 888 889
#> [99] 896 897 898 899 901 903 904 905 906 909 912 913 915 917
#> [113] 919 921 922 923 927 928 932 934 936 937 940 944 945 947
#> [127] 948 949 951 956 958 959 961 962 963 964 966 967 968 972
#> [141] 974 976 978 979 986 988 990 992 995 1000 1300 1700 1702 1703
#> [155] 1704 1705 1708 1710 1712 1714 1716 1717 1718 1720 1721 1723 1725 1726
#> [169] 1727 1729 1731 1733 1735 1736 1737 1739 1742 1745 1747 1748 1752 1754
#> [183] 1756 1758 1759 1760 1766 1768 1770 1771 1773 1775 1778 1782 1784 1785
#> [197] 1790 1792 1793 1795 1796 1797 1799 1800 1802 1803 1804 1805 1806 1807
#> [211] 1808 1809 1821 1824 1825 1827 1828 1829 1833 1835 1837 1840 1841 1842
#> [225] 1848 1849 1851 1852 1854 1855 1857 1859 1860 1862 1863 1865 1867 1868
#> [239] 1876 1878 1879 1880 1882 1883 1884 1886 1887 1889 1894 1898 1899 1905
#> [253] 1906 1907 1908 1909 1912 1919 1926 1927 1928 1930 1933 1934 1935 1936
#> [267] 1938 1940 1944 1947 1950 1952 1954 1955 1956 1960 1962 1963 1965 1966
#> [281] 1967 1969 1970 1974 1976 1977 1978 1980 1985 1987 1988 1990 1997 1998
Rbeast::beast(
mv_data_1, season = "none", print.progress = FALSE, quiet = TRUE
)$trend$cp
#> [1] 1794 1855 1986 1301 301 703 1981 1769 1860 1834
plot(
mcp::mcp(
list(y ~ 1, ~ 1, ~ 1, ~ 1, ~ 1, ~ 1),
data = data.frame(y = mv_data_1, x = seq_len(nrow(mv_data_1))),
par_x = "x"
)
)
plot(not::not(mv_data_1, contrast = "pcwsConstMeanVar"))
(fastcpd_result <- fastcpd::fastcpd.mean(mean_data_3, r.progress = FALSE)@cp_set)
testthat::expect_equal(fastcpd_result, c(300, 700), tolerance = 0.2)
(CptNonPar_result <- CptNonPar::np.mojo(mean_data_3, G = floor(nrow(mean_data_3) / 6))$cpts)
testthat::expect_equal(CptNonPar_result, c(300, 700), tolerance = 0.2)
invisible(
suppressMessages(
capture.output(
result_InspectChangepoint <- InspectChangepoint::inspect(
t(mean_data_3),
threshold = InspectChangepoint::compute.threshold(
nrow(mean_data_3), ncol(mean_data_3)
)
)
)
)
)
result_InspectChangepoint$changepoints[, "location"]
#> [1] 300 700
(jointseg_result <- jointseg::jointSeg(mean_data_3, K = 2)$bestBkp)
testthat::expect_equal(jointseg_result, c(300, 700), tolerance = 0.2)
invisible(
capture.output(
result_Rbeast <- Rbeast::beast123(
mean_data_3,
metadata = list(whichDimIsTime = 1),
season = "none"
)
)
)
result_Rbeast$trend$cp
#> Warning message:
#> In Rbeast::beast123(mean_data_3, metadata = list(whichDimIsTime = 1), :
#> WARNING: If the input data is regular and ordered in time,the times of individual datapoints are determined fully by 'metadata$startTime' and 'metadata$deltaTime'. But startTime and deltaTime are missing and a default value 1 is used for both!
#> [,1] [,2] [,3]
#> [1,] 301 301 701
#> [2,] 701 701 301
#> [3,] NaN 225 NaN
#> [4,] NaN 684 NaN
#> [5,] NaN NaN NaN
#> [6,] NaN NaN NaN
#> [7,] NaN NaN NaN
#> [8,] NaN NaN NaN
#> [9,] NaN NaN NaN
#> [10,] NaN NaN NaN
strucchange::breakpoints(
cbind(y.1, y.2, y.3) ~ 1, data = data.frame(y = mean_data_3)
)$breakpoints
#> [1] 300 700
ecp::e.divisive(mean_data_3)$estimates
#> [1] 1 301 701 1001
plot(bcp::bcp(mean_data_3))
(fastcpd_result <- fastcpd::fastcpd.mv(mv_data_3, r.progress = FALSE)@cp_set)
testthat::expect_equal(fastcpd_result, c(300, 700, 1000, 1300, 1700), tolerance = 0.2)
ecp::e.divisive(mv_data_3)$estimates
#> [1] 1 301 701 1001 1301 1701 2001
invisible(
suppressMessages(
capture.output(
result_InspectChangepoint <- InspectChangepoint::inspect(
t(mv_data_3),
threshold = InspectChangepoint::compute.threshold(
nrow(mv_data_3), ncol(mv_data_3)
)
)
)
)
)
result_InspectChangepoint$changepoints[, "location"]
#> [1] 300 700 701 703 705 707 708 709 711 712 714 715 717 718
#> [15] 720 721 723 724 726 727 729 731 733 734 736 737 739 740
#> [29] 742 743 744 746 747 749 750 752 753 754 755 756 758 760
#> [43] 762 763 765 766 767 769 770 772 773 774 775 777 779 780
#> [57] 782 784 786 788 790 791 793 795 797 799 801 803 804 806
#> [71] 808 809 810 811 813 814 816 817 818 820 821 823 825 827
#> [85] 828 830 831 833 835 836 837 838 840 842 843 845 846 848
#> [99] 849 850 852 853 854 855 856 858 859 860 862 863 865 866
#> [113] 868 869 871 872 874 876 877 878 879 881 883 885 887 888
#> [127] 889 891 893 894 895 897 898 900 901 903 904 906 908 909
#> [141] 911 913 914 916 917 918 920 921 923 924 925 927 928 929
#> [155] 931 932 934 936 937 938 939 941 942 943 945 946 947 949
#> [169] 950 952 954 955 956 957 958 959 961 962 964 965 967 968
#> [183] 970 972 973 974 975 977 979 981 982 984 985 986 987 988
#> [197] 990 991 992 994 995 997 999 1000 1300 1700 1702 1703 1704 1705
#> [211] 1706 1708 1709 1710 1712 1713 1714 1715 1717 1719 1721 1722 1723 1725
#> [225] 1727 1729 1730 1732 1734 1735 1737 1738 1739 1741 1742 1744 1746 1748
#> [239] 1750 1752 1753 1754 1755 1757 1758 1759 1761 1762 1763 1764 1766 1767
#> [253] 1769 1770 1771 1773 1774 1775 1777 1779 1781 1782 1783 1785 1786 1788
#> [267] 1789 1791 1793 1794 1796 1798 1800 1803 1804 1805 1806 1808 1809 1811
#> [281] 1812 1814 1815 1817 1818 1819 1821 1822 1824 1825 1827 1828 1829 1831
#> [295] 1833 1835 1836 1838 1839 1841 1843 1844 1846 1847 1848 1850 1851 1853
#> [309] 1854 1856 1857 1858 1859 1860 1862 1863 1864 1865 1867 1869 1870 1872
#> [323] 1873 1874 1876 1878 1879 1881 1882 1884 1885 1887 1889 1891 1893 1894
#> [337] 1896 1898 1899 1900 1901 1902 1904 1906 1907 1909 1911 1913 1914 1916
#> [351] 1917 1918 1919 1921 1923 1924 1925 1927 1928 1930 1932 1933 1935 1936
#> [365] 1938 1939 1941 1942 1944 1946 1948 1950 1951 1952 1954 1956 1957 1959
#> [379] 1961 1963 1965 1967 1968 1970 1972 1973 1974 1976 1977 1979 1981 1982
#> [393] 1984 1985 1987 1989 1990 1992 1993 1995 1996 1998
invisible(
capture.output(
result_Rbeast <- Rbeast::beast123(
mv_data_3,
metadata = list(whichDimIsTime = 1),
season = "none"
)
)
)
result_Rbeast$trend$cp
#> Warning message:
#> In Rbeast::beast123(mv_data_3, metadata = list(whichDimIsTime = 1), :
#> WARNING: If the input data is regular and ordered in time,the times of individual datapoints are determined fully by 'metadata$startTime' and 'metadata$deltaTime'. But startTime and deltaTime are missing and a default value 1 is used for both!
#> [,1] [,2] [,3] [,4]
#> [1,] 701 301 301 710
#> [2,] 1301 1301 1301 301
#> [3,] 301 701 702 1301
#> [4,] 814 1993 1829 1986
#> [5,] 1968 767 1822 785
#> [6,] 1994 781 810 774
#> [7,] 761 884 845 1912
#> [8,] 1873 755 798 792
#> [9,] 1865 747 771 879
#> [10,] 1962 924 1700 1919
(fastcpd_result <- fastcpd::fastcpd.lm(lm_data, r.progress = FALSE)@cp_set)
testthat::expect_equal(fastcpd_result, c(97, 201), tolerance = 0.2)
strucchange::breakpoints(y ~ . - 1, data = lm_data)$breakpoints
#> [1] 100 201
(segmented_result <- segmented::segmented(
lm(
y ~ . - 1,
data.frame(y = lm_data$y, x = lm_data[, -1], index = seq_len(nrow(lm_data)))
),
seg.Z = ~ index
)$psi[, "Est."])
testthat::expect_equivalent(segmented_result, c(233), tolerance = 0.2)
(fastcpd_result <- fastcpd::fastcpd.binomial(binomial_data, r.progress = FALSE)@cp_set)
testthat::expect_equal(fastcpd_result, 302, tolerance = 0.2)
strucchange::breakpoints(y ~ . - 1, data = binomial_data)$breakpoints
#> [1] 297
(fastcpd_result <- fastcpd::fastcpd.poisson(poisson_data, r.progress = FALSE)@cp_set)
testthat::expect_equal(fastcpd_result, c(498, 805, 1003), tolerance = 0.2)
strucchange::breakpoints(y ~ . - 1, data = poisson_data)$breakpoints
#> [1] 935
(fastcpd_result <- fastcpd::fastcpd.lasso(lasso_data, r.progress = FALSE)@cp_set)
testthat::expect_true(sum(fastcpd_result - c(79, 199, 320)) <= 1)
strucchange::breakpoints(y ~ . - 1, data = lasso_data)$breakpoints,
#> [1] 80 200 321
(fastcpd_result <- fastcpd::fastcpd.ar(ar_data, 3, r.progress = FALSE)@cp_set)
testthat::expect_equal(fastcpd_result, c(614), tolerance = 0.2)
(CptNonPar_result <- CptNonPar::np.mojo(ar_data, G = floor(length(ar_data) / 6))$cpts)
testthat::expect_equal(CptNonPar_result, numeric(0), tolerance = 0.2)
(segmented_result <- segmented::segmented(
lm(
y ~ x + 1, data.frame(y = ar_data, x = seq_along(ar_data))
),
seg.Z = ~ x
)$psi[, "Est."])
testthat::expect_equivalent(segmented_result, c(690), tolerance = 0.2)
plot(
mcp::mcp(
list(y ~ 1 + ar(3), ~ 0 + ar(3)),
data = data.frame(y = ar_data, x = seq_along(ar_data)),
par_x = "x"
)
)
(fastcpd_result <- fastcpd::fastcpd.garch(garch_data, c(1, 1), r.progress = FALSE)@cp_set)
testthat::expect_equal(fastcpd_result, c(205), tolerance = 0.2)
(CptNonPar_result <- CptNonPar::np.mojo(garch_data, G = floor(length(garch_data) / 6))$cpts)
testthat::expect_equal(CptNonPar_result, c(206), tolerance = 0.2)
strucchange::breakpoints(x ~ 1, data = data.frame(x = garch_data))$breakpoints
#> [1] NA
(fastcpd_result <- fastcpd::fastcpd.var(
var_data, 2, cost_adjustment = NULL, r.progress = FALSE
)@cp_set)
testthat::expect_equal(fastcpd_result, c(500), tolerance = 0.2)
(VARDetect_result <- VARDetect::tbss(var_data)$cp)
testthat::expect_equal(VARDetect_result, c(501), tolerance = 0.2)
well_log <- well_log[well_log > 1e5]
result <- list(
fastcpd = fastcpd.mean(well_log, trim = 0.003)@cp_set,
changepoint = changepoint::cpt.mean(well_log)@cpts,
CptNonPar =
CptNonPar::np.mojo(well_log, G = floor(length(well_log) / 6))$cpts,
strucchange = strucchange::breakpoints(
y ~ 1, data = data.frame(y = well_log)
)$breakpoints,
ecp = ecp::e.divisive(matrix(well_log))$estimates,
breakfast = breakfast::breakfast(well_log)$cptmodel.list[[6]]$cpts,
wbs = wbs::wbs(well_log)$cpt$cpt.ic$mbic.penalty,
mosum = mosum::mosum(c(well_log), G = 40)$cpts.info$cpts,
# fpop = fpop::Fpop(well_log, length(well_log))$t.est, # meaningless
gfpop = gfpop::gfpop(
data = well_log,
mygraph = gfpop::graph(
penalty = 2 * log(length(well_log)) * gfpop::sdDiff(well_log) ^ 2,
type = "updown"
),
type = "mean"
)$changepoints,
InspectChangepoint = InspectChangepoint::inspect(
well_log,
threshold = InspectChangepoint::compute.threshold(length(well_log), 1)
)$changepoints[, "location"],
jointseg = jointseg::jointSeg(well_log, K = 12)$bestBkp,
Rbeast = Rbeast::beast(
well_log, season = "none", print.progress = FALSE, quiet = TRUE
)$trend$cp,
stepR = stepR::stepFit(well_log, alpha = 0.5)$rightEnd
)
package_list <- sort(names(result), decreasing = TRUE)
comparison_table <- NULL
for (package_index in seq_along(package_list)) {
package <- package_list[[package_index]]
comparison_table <- rbind(
comparison_table,
data.frame(
change_point = result[[package]],
package = package,
y_offset = (package_index - 1) * 1000
)
)
}
most_selected <- sort(table(comparison_table$change_point), decreasing = TRUE)
most_selected <- sort(as.numeric(names(most_selected[most_selected >= 4])))
for (i in seq_len(length(most_selected) - 1)) {
if (most_selected[i + 1] - most_selected[i] < 2) {
most_selected[i] <- NA
most_selected[i + 1] <- most_selected[i + 1] - 0.5
}
}
(most_selected <- most_selected[!is.na(most_selected)])
#> [1] 6.0 314.0 434.0 704.0 776.0 1021.0 1057.0 1347.0 1405.0 1502.0 1661.0 1842.0 2023.0 2202.0
#> [15] 2384.5 2445.0 2507.0 2567.5 2738.0 2921.0 3094.0 3251.0 3464.0 3499.0 3622.0 3709.0 3820.0 3976.0
ggplot2::ggplot() +
ggplot2::geom_point(
data = data.frame(x = seq_along(well_log), y = c(well_log)),
ggplot2::aes(x = x, y = y)
) +
ggplot2::geom_vline(
xintercept = most_selected,
color = "black",
linetype = "dashed",
alpha = 0.2
) +
ggplot2::geom_point(
data = comparison_table,
ggplot2::aes(x = change_point, y = 50000 + y_offset, color = package),
shape = 17,
size = 1.9
) +
ggplot2::geom_hline(
data = comparison_table,
ggplot2::aes(yintercept = 50000 + y_offset, color = package),
linetype = "dashed",
alpha = 0.1
) +
ggplot2::coord_cartesian(
ylim = c(50000 - 500, max(well_log) + 1000),
xlim = c(-200, length(well_log) + 200),
expand = FALSE
) +
ggplot2::theme(
panel.background = ggplot2::element_blank(),
panel.border = ggplot2::element_rect(colour = "black", fill = NA),
panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank()
) +
ggplot2::xlab(NULL) + ggplot2::ylab(NULL)
microbenchmark_result <- microbenchmark::microbenchmark(
fastcpd = fastcpd::fastcpd.mean(well_log, trim = 0.003, r.progress = FALSE),
changepoint = changepoint::cpt.mean(well_log, method = "PELT"),
CptNonPar = CptNonPar::np.mojo(well_log, G = floor(length(well_log) / 6)),
strucchange =
strucchange::breakpoints(y ~ 1, data = data.frame(y = well_log)),
ecp = ecp::e.divisive(matrix(well_log)),
breakfast = breakfast::breakfast(well_log),
wbs = wbs::wbs(well_log),
mosum = mosum::mosum(c(well_log), G = 40),
fpop = fpop::Fpop(well_log, nrow(well_log)),
gfpop = gfpop::gfpop(
data = well_log,
mygraph = gfpop::graph(
penalty = 2 * log(length(well_log)) * gfpop::sdDiff(well_log) ^ 2,
type = "updown"
),
type = "mean"
),
InspectChangepoint = InspectChangepoint::inspect(
well_log,
threshold = InspectChangepoint::compute.threshold(length(well_log), 1)
),
jointseg = jointseg::jointSeg(well_log, K = 12),
Rbeast = Rbeast::beast(
well_log, season = "none", print.progress = FALSE, quiet = TRUE
),
stepR = stepR::stepFit(well_log, alpha = 0.5),
not = not::not(well_log, contrast = "pcwsConstMean"),
times = 10
)
#> Unit: milliseconds
#> expr min lq mean median uq max neval
#> fastcpd 8.411284e+01 8.692226e+01 1.011440e+02 1.044509e+02 1.089672e+02 118.05842 10
#> changepoint 3.206773e+01 3.377081e+01 6.843465e+01 3.857181e+01 5.243834e+01 244.76672 10
#> CptNonPar 1.827381e+04 1.901094e+04 2.002752e+04 1.985180e+04 2.076803e+04 22511.59316 10
#> strucchange 5.955079e+04 6.059315e+04 6.185727e+04 6.131291e+04 6.312073e+04 64638.93090 10
#> ecp 7.590543e+05 7.707573e+05 7.859080e+05 7.830752e+05 8.093015e+05 810339.52140 10
#> breakfast 9.170992e+03 9.344041e+03 9.628236e+03 9.382078e+03 9.628663e+03 11073.79318 10
#> wbs 1.139078e+02 1.145472e+02 1.178167e+02 1.166746e+02 1.201676e+02 127.27064 10
#> mosum 1.172847e+00 1.231747e+00 1.740727e+01 1.416854e+00 1.919586e+00 160.76997 10
#> fpop 2.588228e+00 2.630407e+00 4.587742e+00 2.832556e+00 3.312986e+00 18.52067 10
#> gfpop 5.971245e+01 6.072684e+01 6.533492e+01 6.172578e+01 6.839653e+01 87.89407 10
#> InspectChangepoint 1.698673e+02 1.909034e+02 2.392539e+02 2.117010e+02 3.004474e+02 329.87724 10
#> jointseg 2.000894e+01 2.136878e+01 2.551210e+01 2.167757e+01 2.403593e+01 47.98397 10
#> Rbeast 6.533998e+02 6.625203e+02 6.783586e+02 6.792646e+02 6.875840e+02 723.45376 10
#> stepR 2.793709e+01 2.902084e+01 4.380857e+01 3.068416e+01 3.227125e+01 164.81082 10
#> not 9.599763e+01 9.701856e+01 1.028601e+02 1.012292e+02 1.049974e+02 120.73529 10
ggplot2::autoplot(microbenchmark_result)
