Rで平滑化スプライン
データ x,z
lambda: 平滑化パラメータ
罰則付き最小二乗法
#c1<-seq(0.1,100,length.out=1000) c1<-exp(seq(log(0.1),log(1000),length.out=1000)) n<-rnorm(length(c1),mean=1,sd=0.015) z<-rpois(length(c1),lambda=n*c1) x<-1:length(z) #x<-seq(1,10,length.out=100) x1<-x yy<-z #yy<-3*(x-5)^2+rnorm(length(x),mean=0,sd=10) #y1<-y #plot(x1,yy) knot<-150 f12<-function(ytmp){ #y<-3*x^2+rnorm(length(x),mean=0,sd=10) # plot(x,y) h<-diff(x_knots) len<-length(h) v<-6*(diff(ytmp)[2:len]/h[2:len]-diff(ytmp)[1:(len-1)]/h[1:(len-1)]) m<-matrix(nrow=len-1,ncol=len-1) for(i in 1:(len-1)){ if(i==1){ m[i,]<-c(rep(0,i-1),2*(h[i]+h[i+1]),h[i+2],rep(0,len-3-i+2-1)) } if(i>1 & i<len-1){ m[i,]<-c(rep(0,i-2),h[i],2*(h[i]+h[i+1]),h[i+2],rep(0,len-3-1-(i-2))) } if(i==(len-1)){ m[i,]<-c(rep(0,i-2),h[i],2*(h[i]+h[i+1])) } } u<-solve((m),v) b<-c(0,u,0)/2 a<-diff(c(0,u,0))/6/diff(x_knots) d<-ytmp c<-diff(ytmp)/diff(x_knots)-1.0/6.0*diff(x_knots)*(2*c(0,u,0)[1:(length(c(0,u,0))-1)]+c(0,u,0)[2:length(c(0,u,0))]) f<-function(z){ sum<-0 for(i in 1:(length(x_knots)-1)){ sum<-sum+dunif(z,min=x_knots[i],max=x_knots[i+1])*(x_knots[i+1]-x_knots[i])*(a[i]*(z-x_knots[i])^3+b[i]*(z-x_knots[i])^2+c[i]*(z-x_knots[i])+d[i]) if(length(which(z==x_knots[i]))>=1 & i!=1){ sum[which(z==x_knots[i])]<-sum[which(z==x_knots[i])]-d[i] } } sum } ff<-function(z){ sum<-0 for(i in 1:(length(x_knots)-1)){ sum<-sum+dunif(z,min=x_knots[i],max=x_knots[i+1])*(x_knots[i+1]-x_knots[i])*(6*a[i]*(z-x_knots[i])^1+2*b[i]) if(length(which(z==x_knots[i]))>=1 & i!=1){ sum[which(z==x_knots[i])]<-sum[which(z==x_knots[i])]-6*a[i]*(z[z==x_knots[i]]-x_knots[i])^1-2*b[i] } } sum } par(mfcol=c(2,1)) #plot(x1,yy) plot(x_knots,ytmp) curve(f(x),add=T,col=2) #xxx<-seq(1,10,length.out=100) xxx<-x_knots+0.00 points(xxx,f(xxx),col=4) #xxx<-seq(1,10,length.out=100) points(xxx,f(xxx),col=3) points(x1,yy,col=4) xxx<-x1 curve(ff(x),xlim=c(1,10)) points(xxx+0.000,ff(xxx+0.000),col=3) #xxx<-x1 #points(xxx,f(xxx),col=4) #lambda<-0.1 lambda<-0.0 qq1<-sum((yy-f(x1))^2) qq2<-sum(ff(x1)^2) #qq1<-sum(rnorm()) qq<-qq1+lambda*qq2 #lambda<-1 #beta<a-0.5 #beta*sum((yy-f(x1))^2)+ cat(qq1,qq2,qq,"\n") #lambda*sum(ff(x1+0.00)^2) #f(xxx) qq } #f12(rep(1,length(yy))) #f12(yy) x_knots<-seq(1,max(x1),length.out=15) #y_knots<-predict(smooth.spline(x1,yy),x_knots)$y y_knots<-rep(1,length(x_knots)) #x1<- #yyy<- #yy_sp<-smooth.spline(x_knots,y_knots) #z2<-optim(yy,f12,method="BFGS") #z2<-optim(yy_sp$y,f12,method="BFGS") ###z2<-optim(y_knots,f12,method="BFGS") z2<-optim(y_knots,f12,method="CG") #z2<-optim(yy_sp$y,f12,method="CG") #"CG" #z2<-optim(yy_sp$y,f12) #z2<-optim(rep(1,length(yy)),f12,method="BFGS") #optim(rep(1,length(yy)),f12) #optim(yy,f12) #f12(rep(1,length(y1))) #curve(f(x),col=2,add=T)
罰則付き最尤法
c1<-seq(0.1,0.5,length.out=1000)^1 #c1<-seq(0.0,100,length.out=1000)^2 n<-rnorm(length(c1),mean=1,sd=0.015) z<-rpois(length(c1),lambda=n*c1) x<-1:length(z) #x<-seq(1,10,length.out=100) x1<-x yy<-z #yy<-3*(x-5)^2+rnorm(length(x),mean=0,sd=10) #y1<-y #plot(x1,yy) knot<-150 f12<-function(ytmp){ #y<-3*x^2+rnorm(length(x),mean=0,sd=10) # plot(x,y) h<-diff(x_knots) len<-length(h) v<-6*(diff(ytmp)[2:len]/h[2:len]-diff(ytmp)[1:(len-1)]/h[1:(len-1)]) m<-matrix(nrow=len-1,ncol=len-1) for(i in 1:(len-1)){ if(i==1){ m[i,]<-c(rep(0,i-1),2*(h[i]+h[i+1]),h[i+2],rep(0,len-3-i+2-1)) } if(i>1 & i<len-1){ m[i,]<-c(rep(0,i-2),h[i],2*(h[i]+h[i+1]),h[i+2],rep(0,len-3-1-(i-2))) } if(i==(len-1)){ m[i,]<-c(rep(0,i-2),h[i],2*(h[i]+h[i+1])) } } u<-solve((m),v) b<-c(0,u,0)/2 a<-diff(c(0,u,0))/6/diff(x_knots) d<-ytmp c<-diff(ytmp)/diff(x_knots)-1.0/6.0*diff(x_knots)*(2*c(0,u,0)[1:(length(c(0,u,0))-1)]+c(0,u,0)[2:length(c(0,u,0))]) f<-function(z){ sum<-0 for(i in 1:(length(x_knots)-1)){ sum<-sum+dunif(z,min=x_knots[i],max=x_knots[i+1])*(x_knots[i+1]-x_knots[i])*(a[i]*(z-x_knots[i])^3+b[i]*(z-x_knots[i])^2+c[i]*(z-x_knots[i])+d[i]) if(length(which(z==x_knots[i]))>=1 & i!=1){ sum[which(z==x_knots[i])]<-sum[which(z==x_knots[i])]-d[i] } } sum } ff<-function(z){ sum<-0 for(i in 1:(length(x_knots)-1)){ sum<-sum+dunif(z,min=x_knots[i],max=x_knots[i+1])*(x_knots[i+1]-x_knots[i])*(6*a[i]*(z-x_knots[i])^1+2*b[i]) if(length(which(z==x_knots[i]))>=1 & i!=1){ sum[which(z==x_knots[i])]<-sum[which(z==x_knots[i])]-6*a[i]*(z[z==x_knots[i]]-x_knots[i])^1-2*b[i] } } sum } par(mfcol=c(2,1)) #plot(x1,yy) plot(x_knots,ytmp) curve(f(x),add=T,col=2) #xxx<-seq(1,10,length.out=100) xxx<-x_knots+0.00 points(xxx,f(xxx),col=4) #xxx<-seq(1,10,length.out=100) points(xxx,f(xxx),col=3) points(x1,yy,col=4) xxx<-x1 curve(ff(x),xlim=c(1,10)) points(xxx+0.000,ff(xxx+0.000),col=3) #xxx<-x1 #points(xxx,f(xxx),col=4) #lambda<-0.1 lambda<-0.5 yyb<-predict(smooth.spline(x_knots,ytmp),x1)$y #qq1<-sum((yy-f(x1))^2) mm1<-f(x1[which(yyb<=16)]) mm1[mm1<=0]<-0 lm1<--log(dpois(yy[which(yyb<=16)],lambda=mm1)) lm1[lm1==Inf]<-735 qq1<-sum(lm1) cat("A1",qq1,sum(mm1),"\n") mm1<-f(x1[which(yyb>16)]) mm1[mm1<=0]<-0 sd1<-sqrt(mm1+0.015^2*mm1^2) lm2<--log(dnorm(yy[which(yyb>16)],mean=mm1,sd=sd1)) lm2[lm2==Inf]<-735 qq1<-qq1+sum(lm2) cat("B1",qq1,sum(mm1),"\n") qq2<-sum(ff(x1)^2) #qq1<-sum(rnorm()) if(any(ytmp<0)){ qq3<-1*exp(10*mean(ytmp[ytmp<0])) }else{ qq3<-0 } lambda2<-1 qq<-qq1+lambda*qq2+lambda2*qq3 #lambda<-1 #beta<a-0.5 #beta*sum((yy-f(x1))^2)+ cat(qq1,qq2,qq,"\n") #lambda*sum(ff(x1+0.00)^2) #f(xxx) qq } #f12(rep(1,length(yy))) #f12(yy) x_knots<-seq(1,max(x1),length.out=15) y_knots<-predict(smooth.spline(x1,yy),x_knots)$y y_knots[y_knots<=0]<-0 #y_knots<-rep(1,length(x_knots)) #x1<- #yyy<- yy_sp<-smooth.spline(x_knots,y_knots) #z2<-optim(yy,f12,method="BFGS") z2<-optim(yy_sp$y,f12,method="BFGS") #z2<-optim(yy_sp$y,f12) ####z2<-optim(y_knots,f12) #z2<-optim(rep(1,length(y_knots)),f12,method="SANN") #z2<-optim(rep(1,length(y_knots)),f12) # #z2<-optim(rep(1,length(y_knots)),f12,method="CG") #z2<-optim(rep(1,length(y_knots)),f12,method="BFGS") #z2<-optim(y_knots,f12,method="BFGS") #z2<-optim(y_knots,f12,method="CG") #z2<-optim(yy_sp$y,f12,method="CG") #"CG" #z2<-optim(yy_sp$y,f12) #z2<-optim(rep(1,length(yy)),f12,method="SANN") #optim(rep(1,length(yy)),f12) #optim(yy,f12) #f12(rep(1,length(y1))) #curve(f(x),col=2,add=T) plot(yy,ylim=c(0,max(yy))) points(x_knots,y_knots,col=0) #points(yy,col=2) points(x_knots,y_knots,col=3) points(x_knots,z2$par,col=4) qq_sp1<-smooth.spline(c1) qq_sp1<-predict(qq_sp1,x_knots)$y points(x_knots,qq_sp1,col=5) mean((qq_sp1-y_knots)^2) mean((qq_sp1-z2$par)^2)
