全部版块 我的主页
论坛 数据科学与人工智能 数据分析与数据科学 R语言论坛
3070 1
2008-08-05

我试了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编辑过]

二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

全部回复
2014-12-24 17:59:10
multiroot可以做,具体函数是——
multiroot.1D(f, start, maxiter = 100,
             rtol = 1e-6, atol = 1e-8, ctol = 1e-8,
             nspec = NULL, dimens = NULL,
             verbose = FALSE, positive = FALSE, names = NULL, ...)

   
nleqslv()也可以做,具体函数是——
nleqslv(x, fn, jac=NULL,...,
               method = c("Broyden", "Newton"),
               global = c("dbldog", "pwldog", "qline", "gline", "none"),
               xscalm = c("fixed","auto"),
               control = list()
        )

二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

相关推荐
栏目导航
热门文章
推荐文章

分享

扫码加好友,拉您进群