作者askuse (问问题)
看板Statistics
标题[问题] 在R中使用二分法求根与if else回圈的问题
时间Fri Oct 19 23:12:15 2012
[软体程式类别]:R
[程式问题]:使用二分法求根与if else回圈时出现"需要 TRUE/FALSE 值的地方有缺值"
的错误讯息。
[软体熟悉度]:不到3个月
(Q1/K1)*ln(1+B1*(P*y1/x1)^K1)+(q1/k1)*ln(1+b1*(P*y1/x1)^k1)
-(Q2/K2)*ln(1+B2*(P*(1-y1)/(1-x1))^K2)-(q2/k2)*ln(1+b2*(P*(1-y1)/(1-x1))^k2)=0
主要问题是想藉由二分法求出不同P值时的x1。
Q, q, B, b, K, k都是在拟合曲线後所得到的参数，会在跑程式之前输入进去。
y1的部分也会先输入进去。而且已经知道x1的范围都在0~1之间。
程式码的档案在此
http://www.sendspace.com/file/23078m
f <- function(P,Q,q,B,b,K,k,y,x1) {
f=Q[1]/K[1]*log(1+B[1]*(P*y[1]/x1)^K[1])
+q[1]/k[1]*log(1+b[1]*(P*y[1]/x1)^k[1])
-Q[2]/K[2]*log(1+B[2]*(P*(1-y[1])/(1-x1))^K[2])
-q[2]/k[2]*log(1+b[2]*(P*(1-y[1])/(1-x1))^k[2])
}
bisection <- function(P,Q,q,B,b,K,k,y) {
lt=0.0
rt=1.0
while (abs(rt-lt)>1E-6) {
mi=(lt+rt)/2
if (f(P,Q,q,B,b,K,k,y,lt)*f(P,Q,q,B,b,K,k,y,mi)<0)
rt=mi
else
lt=mi
}
bisection=mi
}
IAST <- function(n,y,PP,Q,q,B,b,K,k) {
x1=c()
x2=c()
P1=c()
P2=c()
N1=c()
N2=c()
P=seq(PP[1],PP[2],length=n)
S=c()
y1=y[1]
y2=y[2]
for ( i in 1:n) {
x1[i]=bisection(P[i],Q,q,B,b,K,k,y)
x2[i]=1-x1[i]
}
for ( i in 1:n) {
P1[i]=P[i]*y[1]/x1[i]
P2[i]=P[i]*y[2]/x2[i]
N1[i]=Q[1]*B[1]*P1[i]^K[1]/(1+B[1]*P1[i]^K[1])
+q[1]*b[1]*P1[i]^k[1]/(1+b[1]*P1[i]^k[1])
N2[i]=Q[2]*B[2]*P2[i]^K[2]/(1+B[2]*P2[i]^K[2])
+q[2]*b[2]*P2[i]^k[2]/(1+b[2]*P2[i]^k[2])
N1[i]=N1[i]*N2[i]/(N2[i]+N1[i]*(1/x1[i]-1))
N2[i]=N1[i]*(1/x1[i]-1)
S[i]=(x1[i]/y[1])/(x2[i]/y[2])
}
IAST=data.frame(P=P,x1=x1,y1=y1,x2=x2,y2=y2,S=S,N1=N1,N2=N2)
}
input则在这里
http://www.sendspace.com/file/a8hv66
n=10
PP=c(10,100)
y=c(0.15,1.85)
Q=matrix(c(20,5),1,2,TRUE)
q=matrix(c(10,2.5),1,2,TRUE)
B=matrix(c(0.009,0.00045),1,2,TRUE)
b=matrix(c(0.0009,0.000045),1,2,TRUE)
K=matrix(c(0.8,0.4),1,2,TRUE)
k=matrix(c(0.2,0.1),1,2,TRUE)
iso=IAST(n,PP,y,Q,q,B,b,K,k)
想请教大家程式码与input的错误在哪？
先谢谢大家了！
--
※ 发信站: 批踢踢实业坊(ptt.cc)
※ 编辑: askuse 来自: 140.109.90.178 (10/21 15:59)
1F：→ gsuper:while (abs(rt-lt)>1E-6) 这附近检查看看 10/21 21:14