#分段回归的代码

x <- c(1:10, 13:22)

y <- numeric(20)

## Create first segment

y[1:10] <- 20:11 + rnorm(10, 0, 1.5)

## Create second segment

y[11:20] <- seq(11, 15, len=10) + rnorm(10, 0, 1.5)

## Plot it

library(segmented)

lin.mod <- lm(y~x)

#psi starting values for the breakpoints to be estimated

segmented.mod <- segmented(lin.mod, seg.Z = ~x, psi=14)

plot(x,y, pch=16, ylim=c(5,20))

plot(segmented.mod, add=T)