> data<-read.csv("C:/Users/sjr88/Desktop/0421.csv",header=T)
> x1 <- data[,1]
> x2 <- data[,2]
> d <- data[,5]
> mu <- data[,6]
>
> #对数似然函数
> logLikweibullcox <- function(param){
+   beta <- param[1]
+   yita <- param[2]
+   gama1 <- param[3]
+   gama2 <- param[4]
+   ll <-sum(d*(log(abs(beta))+(beta-1)*log(abs(mu))-beta*log(abs(yita))+(x1*gama1+x2*gama2))-((mu/yita)^beta)*exp(x1*gama1+x2*gama2))
+   return(-ll)}
> #一节偏导矩阵
> logLikGradValues <-function(param){
+   beta <- param[1]
+   yita <- param[2]
+   gama1 <- param[3]
+   gama2 <- param[4]
+   logLikGradValues <- numeric(4)
+   logLikGradValues[1]<- sum(d*(1/beta+log(abs(mu))-log(abs(yita)))-(mu/yita)^beta*log(abs(mu/yita))*exp(x1*gama1+x2*gama2))
+   logLikGradValues[2]<- sum(d*(-beta/yita)+beta/yita*(mu/yita)^beta*exp(x1*gama1+x2*gama2))
+   logLikGradValues[3]<- sum(d*(x1)-x1*(mu/yita)^beta*exp(x1*gama1+x2*gama2))
+   logLikGradValues[4]<- sum(d*(x2)-x2*(mu/yita)^beta*exp(x1*gama1+x2*gama2))
+   return(logLikGradValues)}
> #二阶偏导矩阵
> logLikeHessValues <-function(param){
+   beta <- param[1]
+   yita <- param[2]
+   gama1 <- param[3]
+   gama2 <- param[4]
+   logLikeHessValues <- matrix(0,nrow = 4,ncol = 4)
+   logLikeHessValues[1,1] <- sum(d*(-1/beta^2)-(mu/beta)^beta*(log(abs(mu/yita)))^2*exp(x1*gama1+x2*gama2))
+   logLikeHessValues[2,2] <- sum(d*(beta/yita^2)-(mu/yita)^beta*(beta^2+beta)/yita^2*exp(x1*gama1+x2*gama2))
+   logLikeHessValues[3,3] <- sum(-x1^2*(mu/yita)^beta*exp(x1*gama1+x2*gama2))
+   logLikeHessValues[4,4] <- sum(-x2^2*(mu/yita)^beta*exp(x1*gama1+x2*gama2))
+   logLikeHessValues[1,2] <- sum(d*(-1/yita)+1/yita*(mu/yita)^beta*(1+beta*log(abs(mu/yita)))*exp(x1*gama1+x2*gama2))
+   logLikeHessValues[2,1] <- logLikeHessValues[1,2]
+   logLikeHessValues[1,3] <- sum(-x1*(mu/yita)^beta*log(abs(mu/yita))*exp(x1*gama1+x2*gama2))
+   logLikeHessValues[3,1] <- logLikeHessValues[1,3]
+   logLikeHessValues[2,3] <- sum(x1*beta/yita*(mu/yita)^beta*exp(x1*gama1+x2*gama2))
+   logLikeHessValues[3,2] <- logLikeHessValues[2,3]
+   logLikeHessValues[1,4] <- sum(-x2*(mu/yita)^beta*log(abs(mu/yita))*exp(x1*gama1+x2*gama2))
+   logLikeHessValues[4,1] <- logLikeHessValues[1,4]
+   logLikeHessValues[2,4] <- sum(x2*beta/yita*(mu/yita)^beta*exp(x1*gama1+x2*gama2))
+   logLikeHessValues[4,2] <- logLikeHessValues[2,4]
+   logLikeHessValues[3,4] <- sum(-x1*x2*(mu/yita)^beta*exp(x1*gama1+x2*gama2))
+   logLikeHessValues[4,3] <- logLikeHessValues[3,4]
+   return(logLikeHessValues)}
> #?maxLik
>
> mle <- maxLik(logLikweibullcox,grad=logLikGradValues,hess=logLikeHessValues,start = c(beta =3,yita = 500,gama1 =-2, gama2 =3),)
> print(summary(mle))
--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 1 iterations
Return code 3: Last step could not find a value above the current.
Boundary of parameter space?  
Consider switching to a more robust optimisation method temporarily.
Log-Likelihood: 12340.08
4  free parameters
Estimates:
        Estimate Std. error  t value Pr(> t)   
beta   3.000e+00  4.620e-06 649362.6  <2e-16 ***
yita   5.000e+02  2.358e+00    212.0  <2e-16 ***
gama1 -2.000e+00  3.053e-02    -65.5  <2e-16 ***
gama2  3.000e+00  1.814e-02    165.4  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------