谢谢楼上的！

如果用牛顿迭代法解方程，r软件显示结果是：错误在drop(.Call("La_dgesv", a, as.matrix(b), tol, PACKAGE = "base")) :
系统计算上是奇异的: 倒条件数=0
此外: Warning messages:
1: 产生了NaNs in: log(x)
2: 产生了NaNs in: log(x)
3: 产生了NaNs in: log(x)
4: 产生了NaNs in: log(x)
5: 产生了NaNs in: log(x)

是不是初值取得不好？1——5个错误是不是因为x是负数呢？

1.程序看来似乎有误:

function kekaoxingfx并没用上.

or只是将y sort,然后取1~r,不需放在kekaoxing内.

而且kekaoxing也没用上or.

fun看来不像二元二次方程.

2.要解决near singular matrix,可参考e-views user guide:

1)The objective function must be defined at the starting values. For example if the objective function contains LOG(C(1)), then C(1) must be greater than zero.

2)If the starting values are such that the derivatives are all zero, you will immediately see an error message indicating that encountered a “Near

Singular Matrix＂, and the estimation procedure will stop.

谢谢。我把程序改了一下

y<-rweibull(150,1,2)

n<-length(y)

allorder<-sort(y)

r<-120

or<-allorder[1:r]

kekaoxingfx<-function(y,n,ind,opt,opp)

{

if(r>n)

list(fail="输入长度大于样本量，重新输入")

else{

B1<-a-0.789*qnorm(1-(1+opp)/2,0,1) *a/sqrt(n)

B2<-a+0.789*qnorm(1-(1+opp)/2,0,1)*a/sqrt(n)

S1<-b+1.053* qnorm(1-(1+opp)/2,0,1)*b/(sqrt(n)*a)

S2<-b-1.053* qnorm(1-(1+opp)/2,0,1)*b/(sqrt(n)*a)

if(ind==0){

data.frame(Beta=a,Sigma=b, Beta_left=B1, Beta_right=B2, Sigma_left=S1, Sigma_right=S2)

}

else {

R<-exp(-(opt/b)^a)

R1<-exp(-(opt/S1)^B2)

R2<-exp(-(opt/S2)^B1)

t<-S1*(log(1/opp)^B1)

if(opp==0.95) data.frame(R=R,R_left=R1, R_right=R2,t=t)

else data.frame(R=R,R_left =R1)

}

}

}

fun<-function(x) #解二元方程组

{

f<-c((sum(or^x[1]*(log(or)-log(x[2])))+(n-r)*allorder[r]^x[1]*(log(allorder[r]-log(x[2]))))/(r*x[2]^x[1])-1/r*sum(log(or))+log(x[2])-1/x[1],((sum(or^x[1])+ (n-r)*allorder[r]^x[1])/r)^(1/x[1])-x[2]) #二元方程组

j<-matrix(c((sum(or^x[1]*(log(or)-log(x[2]))*log(or))+ (n-r)*allorder[r]^x[1]*(log(allorder[r]-log(x[2])))*log(allorder[r])-( sum(or^x[1]*(log(or)-log(x[2])))+(n-r)*allorder[r]^x[1]*(log(allorder[r]-log(x[2])))*log(x[2])))/x[2]^x[1]*r+1/x[1]^2,(-sum(or^x[1])-(n-r) *allorder[r]^x[1]- (sum(or^x[1]*(log(or)-log(x[2])))+(n-r)*allorder[r]^x[1]*(log(allorder[r]-log(x[2]))))*x[1])/(r*x[2]^(x[1]+1))+1/x[2], ((sum(or^x[1])+ (n-r)*allorder[r]^x[1])/r)^(1/x[1])*log((sum(or^x[1])+ (n-r)*allorder[r]^x[1])/r)*(-1/x[1]^2)*(sum(or^x[1])*log(or)+(n-r)* allorder[r]^x[1]*log(or)) /r,-1),nrow=2,byrow=T) #雅克比行列式

list(f=f,j=j)

}

newton<-function(fun,x,ep=1e-5,it_max=100) #牛顿迭代法解上述二元方程组

{

index<-0;k<-1

while(k<=it_max)

{

x1<-x;obj<-fun(x)

x<-x-solve(obj$j,obj$f);

norm<-sqrt((x-x1)%*%(x-x1))

if(norm<ep)

{

Index<-1;break

}

k<-k+1

}

obj<-fun(x);

list(Beta=x[1], Sigma=x[2],it=k,index=index,funval=obj$f)

}

z<-newton(fun,c(1.32305,1.342913))

a<-as.numeric(z[1]) #把向量强制转化成数

b<-as.numeric(z[2]) #把向量强制转化成数

kekaoxingfx(y,150,1,100,0.95)

结果报错说：

错误在drop(.Call("La_dgesv", a, as.matrix(b), tol, PACKAGE = "base")) :
'a' (2 x 62)必需是正文形的
此外: Warning message:
数据长度[123] 不是矩阵行数[2] 的整倍数

请大家帮我看看，谢谢