R中的变化点检测答案

【问题标题】：changepoint detection in RR中的变化点检测
【发布时间】：2014-05-13 21:52:55
【问题描述】：

有人知道如何解决这个 R 问题吗？就是在关系中找到一个变化点，比如下面数据中的x=5。

fitDat <- data.frame(x=1:10, y=c(rnorm(5, sd=0.2),(1:5)+rnorm(5, sd=0.2)))
plot(fitDat$x, fitDat$y)

SAS 运行良好；

/*simulated example*/
data test;
INPUT id $ x y;
cards;
1   1 -0.22711769
2   2 -0.08712909
3   3  0.06922072
4   4 -0.12940913
5   5 -0.43152927
6   6  1.17685016
7   7  1.83410448
8   8  2.88528795
9   9  4.30078012
10 10  4.84517101
;
run;

proc nlmixed data=test;
    parms B0 = 0 B1=1 t0=5 sigma_r=1;
    mu1 = (x <= t0)*B0;
    mu2 = (x > t0)*(B0 + B1*(x - t0));
    mu=mu1+mu2;
    model y~normal(mu,sigma_r);
run;

R 失败；

changePointOptim <- function(x, int, slope, delta){ # code adapted from https://stats.stackexchange.com/questions/7527/change-point-analysis-using-rs-nls
       int + (pmax(delta-x, 0) * slope)
}

# nls
nls(y ~ changePointOptim(x, b0, b1, delta), 
              data = fitDat, 
              start = c(b0 = 0, b1 = 1, delta = 5))

# optimization
sqerror <- function (par, x, y)
       sum((y - changePointOptim(x, par[1], par[2], par[3]))^2)
minObj <- optim(par = c(0, 1, 3), fn = sqerror,
              x = fitDat$x, y = fitDat$y, method = "BFGS")
minObj$par

# nlme
library(nlme)
nlmeChgpt <- nlme(y ~ changePointOptim(b0, b1, delta,x), 
              data = fitDat, 
              fixed = b0 + b1 + delta, 
              start = c(b0=0, b1=1, delta=5)) 
summary(nlmeChgpt)

因为它是一个简单的案例，有一个清晰的信号，R 应该可以工作。我想知道我在 R 中做错了什么（我使用了来自 https://stats.stackexchange.com/questions/7527/change-point-analysis-using-rs-nls 的一些代码）。有人建议/解决方案吗？

谢谢！

威廉

【问题讨论】：

请提供可重复的输入和结果。使用set.seed，这样我们就可以使用相同的输入数据集。解释您如何确定 SAS“有效”而 R“无效”。

标签： r regression nls

【解决方案1】：

您在 changePointOptim 中的标志有误。那么：

set.seed(0)
fitDat <- data.frame(x=1:10, y=c(rnorm(5, sd=0.2),(1:5)+rnorm(5, sd=0.2)))
plot(fitDat$x, fitDat$y)


changePointOptim <- function(x, int, slope, delta){ # code adapted from http://stats.stackexchange.com/questions/7527/change-point-analysis-using-rs-nls
       int + (pmax(x-delta, 0) * slope)  ####Chamnged sign here
}

# optimization
sqerror <- function (par, x, y)
       sum((y - changePointOptim(x, par[1], par[2], par[3]))^2)
minObj <- optim(par = c(0, 1, 3), fn = sqerror,
              x = fitDat$x, y = fitDat$y, method = "BFGS")

给予：

> minObj$par
[1] 0.1581436 1.1762401 5.5963392

这是关于你开始的。

# nls
nls(y ~ changePointOptim(x, b0, b1, delta),
              data = fitDat,
              start = list(b0 = 0, b1 = 1, delta = 5)) ##Note it is a list

给出相同的结果：

0.1581 1.1762 5.5963

目前不确定 NLME

【讨论】：

确实是代码中的错误。谢谢，这确实节省了时间。