jumps_find <-function(u0,u1,degree,yj,uj,err,sigma2hat,alpha,gamma=0.05){
jumps2 <- c()
n <- length(yj)
k0 <- length(u0)
nu1 <- length(u1)
DRSStotal<- matrix(rep(0),length(u1),1)
# fit1<- lm(yj ~ bs(uj, degree = degree, knots = u0,Boundary.knots = c(0,200),intercept=T)-1)
# fit1<- lm(yj ~ sapply(0:(length(u0)-2+degree),BaSplite,x=uj,degree=degree,u=u0,n=length(u0))-1)
fit1<- lm(yj ~ BaSplite1(uj,degree,u0)-1)
yfit1<- fitted(fit1)
DRSSu<-DRSSu2<-DRSSu3<-DRSSuT<- rep(0,nu1)
for(i in 1:nu1){
us<- u0
m=which(u0==u1[i])
if(length(m)>0){
us<- c(u0[-m],rep(u1[i],degree+1))
us<- sort(us)
DRSStotal[i]<- NA
}else{
us<- c(u0,rep(u1[i],degree+1))
us<- sort(us)
# fit2<- lm(yj ~ bs(uj, degree = degree, knots = us,Boundary.knots = c(0,200),intercept=T)-1)
# fit2<- lm(yj ~ sapply(0:(length(us)-2+degree),BaSplite,x=uj,degree=degree,u=us,n=length(us))-1)
fit2<- lm(yj ~ BaSplite1(uj,degree,us)-1)
yfit2<- fitted(fit2)
DRSStotal[i]<- sum((yfit1-yfit2)^2)
#确定x的区间
it<- order(abs(u0-u1[i]))[1:2]
it1<- min(it)
it2<- max(it)
idM2<- which(uj>=u0[it2]+err|uj<u0[it1])
RSS2<- sum((yfit1[idM2]-yfit2[idM2])^2)*n/(length(idM2))
idM3<- which(uj>=u0[it2]|uj<u0[it1]-err)
RSS3<- sum((yfit1[idM3]-yfit2[idM3])^2)*n/(length(idM3))
DRSSuT[i]<- max(RSS2,RSS3)
}
}
lambda_1<- mean(DRSStotal,na.rm = TRUE)/sigma2hat-degree-1
lambda_1<- max(lambda_1,0)
# lambdaM_1[Nt]<- lambda_1
q2_1<- qchisq(1-alpha,degree+1,lambda_1)*sigma2hat
# lambda_2<- max(mean(DRSSuT)/sigma2hat-degree-1,0)
# lambdaM_2[Nt]<- lambda_2
# q2_2<- qchisq(1-alpha,degree+1,lambda_2)*sigma2hat
#plot(u1,DRSStotal,main = paste('Simulation',Nt))
#abline(h=q2,col='red')
### 找跳点
id_1<- which(DRSStotal>q2_1)
# id_2<- which(DRSStotal>q2_2)
#cat('id_1',id_1,'\n')
#cat('id_2',id_2,'\n')
idtemp_1<-id_1
# idtemp_2<-id_2
DRSStotal_temp<-DRSStotal
# while (length(idtemp_2)>0){
# cat(length(idtemp_1),'\n')
if (length(idtemp_1)>=gamma*nu1){
while (length(idtemp_1)>0){
#从id_1中把包含在nx的(2p+1)hn的邻域去掉
nx <- which(DRSStotal_temp==max(DRSStotal_temp,na.rm = TRUE))
u.jump<- mean(u1[nx])
jumps2<- c(jumps2,u.jump)
#将跳点插入
if(length(u.jump)>0){
du<- abs(u0-u.jump)
m<- which(du==min(du))
m<- m[1]
if(du[m]<err/2){
u0<- c(u0[-m],rep(u.jump,degree+1))
u0<- sort(u0)
}else{
u0<- c(u0,rep(u.jump,degree+1))
u0<- sort(u0)
}
}
idinterval<- c(max(1,ceiling((mean(nx)))-ceiling(nu1*2*(degree+0.5)/(k0+1))):min(ceiling((mean(nx)))+ceiling(nu1*2*(degree+0.5)/(k0+1)),length(u1)))
#cat('idinterval',idinterval,'\n')
idtemp_1 <- setdiff(idtemp_1,idinterval)
# idtemp_2 <- setdiff(idtemp_2,idinterval)
# cat('idtemp_1',idtemp_1,'\n')
# cat('idtemp_2',idtemp_2,'\n')
DRSStotal_temp[max(1,ceiling((mean(nx)))-ceiling(nu1*2*(degree+0.5)/(k0+1))):min(ceiling((mean(nx)))+ceiling(nu1*2*(degree+0.5)/(k0+1)),length(u1))]<-0
# cat('DRSStotal_temp',DRSStotal_temp,'\n')
}
}
# cat('jumps2',jumps2,'\n')
jumps_find_list<-list(u0,jumps2)
return(jumps_find_list)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.