- こちらでintegrate()関数を用いて積分をしている
- 以下に示すように回帰すると、回帰式の係数は結構、異なるかもしれない
- 回帰曲線自体は同じくらいなのに
- 推定係数を比較するとずいぶん違うが、積分すると、大差ない
- 雑な作りなので、バグがあるかもしれない
- さらに、続けて、似たようなデータが出たときに回帰係数がどれくらいばらつくかも調べてみよう→こちら
minC<--2
maxC<-2
x<-c(2^(seq(from=minC,to=maxC,length=10)))
a<-5
b<-2
U<-5
L<-2
y<-U/(1+(a/x)^b)+L
s<-0.3
y<-y+rnorm(length(y),sd=sqrt(var(y))*s)
X<-log(x)
fr <- function(ks) {
k1 <- ks[1]
k2 <- ks[2]
k3 <- ks[3]
k4 <- ks[4]
sum((y-(k3/(1+(k1/x)^k2)+k4))^2)
}
fr2 <- function(ks) {
k1 <- ks[1]
k2 <- ks[2]
k3 <- ks[3]
k4 <- ks[4]
sum((y-(k3/(1+exp(-k2*(log(x)-log(k1))))+k4))^2)
}
opt.out<-optim(runif(4)*10, fr)
opt.out2<-optim(runif(4)*10, fr2)
par(mfcol=c(2,2))
plot(x,y)
t<-2^(seq(from=min(X),to=max(X),length=100))
points(t,opt.out$par[3]/(1+(opt.out$par[1]/t)^opt.out$par[2])+opt.out$par[4],col=2,cex=0.01)
plot(log(x),y)
points(log(t),opt.out$par[3]/(1+(opt.out$par[1]/t)^opt.out$par[2])+opt.out$par[4],col=2,cex=0.01)
plot(X,y)
points(log(t),opt.out2$par[3]/(1+exp(-opt.out2$par[2]*(log(t)-log(opt.out2$par[1]))))+opt.out2$par[4],col=2,cex=0.01)
plot(exp(X),y)
points(t,opt.out2$par[3]/(1+exp(-opt.out2$par[2]*(log(t)-log(opt.out2$par[1]))))+opt.out2$par[4],col=2,cex=0.01)
opt.out
opt.out2
hf1<-function(t){opt.out$par[3]/(1+(opt.out$par[1]/t)^opt.out$par[2])+opt.out$par[4]}
hf2<-function(t){opt.out2$par[3]/(1+exp(-opt.out2$par[2]*(log(t)-log(opt.out2$par[1]))))+opt.out2$par[4]}
integrate(hf1,0,4)
integrate(hf2,0,4)
> opt.out
$par
[1] 81.615153 1.145219 62.900702 1.759287
$value
[1] 0.2300158
$counts
function gradient
501 NA
$convergence
[1] 1
$message
NULL
> opt.out2
$par
[1] 7.536972 1.488542 6.623899 1.828171
$value
[1] 0.2518245
$counts
function gradient
343 NA
$convergence
[1] 0
$message
NULL
> hf1<-function(t){opt.out$par[3]/(1+(opt.out$par[1]/t)^opt.out$par[2])+opt.out$par[4]}
> hf2<-function(t){opt.out2$par[3]/(1+exp(-opt.out2$par[2]*(log(t)-log(opt.out2$par[1]))))+opt.out2$par[4]}
>
> integrate(hf1,0,4)
10.6721 with absolute error < 7.8e-05
> integrate(hf2,0,4)
10.66855 with absolute error < 3.8e-05