c1<-seq(0.1,0.5,length.out=1000)^1
#c1<-seq(0.0,100,length.out=1000)^2
#c1<-exp(seq(log(0.1),log(1000),length.out=1000))
#q<-read.csv("data/sikousakugo/siokouji_s.csv")[,3]
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
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*2^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
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)
