ドースレスポンスカーブ

> opt.out
$par
[1] 5.271906 2.112786 5.655760 1.899040

$value
[1] 0.1945089

$counts
function gradient 
     205       NA 

$convergence
[1] 0

$message
NULL

> opt.out2
$par
[1] 13.855497  1.723488 19.658422  1.871615

$value
[1] 0.1943919

$counts
function gradient 
     399       NA 

$convergence
[1] 0

$message
NULL

> 

# 指数的に濃度を振る

minC<--2
maxC<-2
x<-c(2^(seq(from=minC,to=maxC,length=10)))
# Hill's eq のパラメタ
# a EC50
# b n
# U upper limit
# L lower limit
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)

# optim()関数での推定をするための関数を２つ
# y=U*1/(1+(a/x)^b)+L
fr <- function(ks) {   
    k1 <- ks[1]
    k2 <- ks[2]
    k3 <- ks[3]
    k4 <- ks[4]
    sum((y-(k3/(1+(k1/x)^k2)+k4))^2)
    #sum((y-(x^k2/(k1^k2+x^k2)))^2)
}
# y=U*1/(1+exp(-b*(log(x)-log(a))))+L
fr2 <- function(ks) {   ## Rosenbrock Banana function
    k1 <- ks[1]
    k2 <- ks[2]
    k3 <- ks[3]
    k4 <- ks[4]
    sum((y-(k3/(1+exp(-k2*(log(x)-log(k1))))+k4))^2)
}
# optim()関数による推定
# 推定パラメタの初期値は適当にして推定する
opt.out<-optim(runif(4)*10, fr)
opt.out2<-optim(runif(4)*10, fr2)

# plotする
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)
#t<-seq(from=min(x),to=max(x),length=100)
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