_ | 覦覈襦 | 豕蠏手 | 殊螳 | 譯殊碁
FrontPage › ChangePoint

Contents

1 覲蟆曙
2 豕蠏 覲蟆曙


slope螳 覲 讌 谿城.

1 覲蟆曙 #

library(segmented)
options(warn = -1)

change_point <- function(y, chart=F, n=10){
  #y <- sin(sort(runif(100,0,5)))*10+100
  #n=5
  #chart=T
  
  x <- 1:length(y)
  #x <- (1:length(y)) * (max(y) - min(y)) / length(y)

  for(cp_cnt in 1:n){
    fit_lm = lm(y ~ x)  # intercept-only model
    fit_tmp = try(segmented(fit_lm, seg.Z = ~x, npsi = cp_cnt),silent=T)
    
    if (class(fit_tmp)[1] == "try-error"){
      break
      cp_cnt <- cp_cnt - 1
    } else if (length(fit_tmp$psi) == 0){
      break
      cp_cnt <- cp_cnt - 1
    }
    
    r2 <- cor(y, fitted(fit_tmp))^2  
    if(r2 >= 0.9){
      fit_segmented <- fit_tmp
      break
    }
    
    fit_segmented <- fit_tmp
  }
  
  if(1==chart){
    plot(fit_segmented)
    points(x,y)
    lines.segmented(fit_segmented)
    points.segmented(fit_segmented)  
  }
  
  cp <- fit_segmented$psi[,2]
  return(cp)
}


set.seed(1024)
y <- sin(sort(runif(100,0,5)))*10+100 + sin(sort(runif(50,0,5)))*10+100
cp <- change_point(y, T)

蟆郁骸
change_point_2.png


2 豕蠏 覲蟆曙 #

library(segmented)
options(warn = -1)

recent_change_point <- function(y, chart=F, n=10){
  #y <- sin(sort(runif(100,0,5)))*10+100
  #n=5
  #chart=T
  
  x <- 1:length(y)
  #x <- (1:length(y)) * (max(y) - min(y)) / length(y)

  for(cp_cnt in 1:n){
    fit_lm = lm(y ~ x)  # intercept-only model
    fit_tmp = try(segmented(fit_lm, seg.Z = ~x, npsi = cp_cnt),silent=T)
    
    if (class(fit_tmp) == "try-error"){
      break
      cp_cnt <- cp_cnt - 1
    } else if (length(fit_tmp$psi) == 0){
      break
      cp_cnt <- cp_cnt - 1
    }
    
    r2 <- cor(y, fitted(fit_tmp))^2  
    if(r2 >= 0.9){
      fit_segmented <- fit_tmp
      break
    }
    
    fit_segmented <- fit_tmp
  }
  
  if(1==chart){
    plot(fit_segmented)
    points(x,y)
    lines.segmented(fit_segmented)
    points.segmented(fit_segmented)  
  }
  
  cp <- fit_segmented$psi[,2]
  
  start_x <- round(tail(cp,1), 2)
  y2 <- y[x >= start_x]
  x2 <- x[x >= start_x]
  return(data.frame(x=x2, y=y2))
}


set.seed(1024)
y <- sin(sort(runif(100,0,5)))*10+100 + sin(sort(runif(50,0,5)))*10+100
x <- 1:length(y)
plot(x, y)
lines(recent_change_point(y))

蟆郁骸
recent_change_point_3.png
蠍 蠍郁鍵..
企: : るジ讓曙 襦螻豺 企Ν 譯殊語. 襦螻豺
EditText : Print : Mobile : FindPage : DeletePage : LikePages : Powered by MoniWiki : Last modified 2022-11-15 10:49:30

伎 願 願 願 伎 伎螳.