我试了BB和rootSolve效果都不是很好,不知道是不是我的式子的问题。
我要解的是一个三元非线性方程,
trifun<-function(x,theta,ksi,delta){
# target equations of the simulation
n<-length(theta)
F<-rep(NA,n)
F[1]<-(x^theta[2]/(theta[1]+x^theta[2]))*(ksi*delta*(x^theta[2]+theta[1])*(1+exp(theta[3]*x))
-theta[1])/(x^theta[2]*(1+exp(theta[3]*x))+theta[1]*exp(theta[3]*x));
F[2]<-(x^theta[2]*log(x)/(theta[1]+x^theta[2]))*(ksi*delta*(x^theta[2]+theta[1])*(1+exp(theta[3]*x))
-theta[1])/(x^theta[2]*(1+exp(theta[3]*x))+theta[1]*exp(theta[3]*x));
F[3]<-(exp(theta[3]*x)*x/(1+exp(theta[3]*x)))*(ksi*delta*(x^theta[2]+theta[1])*(1+exp(theta[3]*x))
-theta[1])/(x^theta[2]*(1+exp(theta[3]*x))+theta[1]*exp(theta[3]*x));
F
}
先定义一个多元函数,其中theta为参数,其他参数后面都通过模拟给出
alpha1<-2
bet1<-2
alpha2<-2
bet2<-3
theta1<-bet1*alpha1^bet1/(bet2*alpha2^bet2)
theta2<-bet2-bet1
alpha<-3
para.theta<-c(theta1,theta2,alpha)
n<-1000#??
Tvar<-rweibull(n,shape=bet1,scale=1/alpha1)# where shape=a, scale=b in the density function
Cvar<-rweibull(n,shape=bet2,scale=1/alpha2)
delta<-as.numeric(Tvar<=Cvar)
Xvar<-apply(cbind(Tvar,Cvar),1,min)
ksi<-vector()
for(i in 1:n) ksi<-rbinom(1,1,delta/(1+exp(alpha*Xvar))+(1-delta)/(1+exp(2*alpha*Xvar)))
triequation<-function(theta){
Fun<-rep(0,3)
for(i in 1:n) Fun<-Fun+trifun(Xvar,theta,ksi,delta);
Fun
} # 不知道这样定义函数是否合适??请版上大侠给点建议或意见,谢谢!!
para.theta
triequation(para.theta)
library(BB)
theta0<-rep(1,3)
#dfsane(par=theta0,fn=triequation,method=2)
sane(par=theta0,fn=triequation,method=2,control=list(maxit=5000))
[此贴子已经被作者于2008-8-5 21:01:27编辑过]