有人有解决这个R问题的想法吗?它是要在关系中找到一个变化点,比如下面数据中的x=5。
SAS表现出色;
作为一个简单的案例,信号很清晰,R应该可以工作。我想知道我在R中做错了什么(我使用了一些来自https://stats.stackexchange.com/questions/7527/change-point-analysis-using-rs-nls的代码)。有人有建议/解决方案吗?
谢谢!
Willem
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的代码)。有人有建议/解决方案吗?
谢谢!
Willem
set.seed
,以便我们都可以使用相同的输入数据集进行工作。解释一下您如何确定 SAS“有效”,而 R 则“无效”。 - Carl Witthoft