Xi <- seq(-1,by=0.001,0.999)
f1 <- (1/(sqrt(2*pi)*sqrt(0.01)))*exp(-(Xi+0.5)^2/0.02)
f2 <- 1/(sqrt(2*pi)*sqrt(0.05))*exp(-(Xi-0.5)^2/0.1)
t <-(0.3+0.25*Xi+f1+f2)^(-1)
#plot(Xi,t^(-1),col=”blue”,type=”l”)
library(actuar)
Yi <- rburr(2000,1,t,1)
z <- cbind(Xi,Yi)
o <- order(Yi,Xi)
sort_z <- z[o, ]
v <- sort_z[,2]
k <-380
n<-length(Yi)
u <- v[n-k]
z.samp = sort_z[(2000-k+1):2000,]
X <- z.samp[,1]
Y <- z.samp[,2]
a <- 0.1
library(numDeriv)
library(alabama)
gamma <- rep(0,n)
for(j in 1:n)
{
x=Xi[j]
#the objective function to minimize
fn <- function(b)
{
(1/k)*sum( ( log(b[1]+b[2]*(X-x)) + (1+(b[3]+b[4]*(X-x))^(-1)) *
log(1 + Y*( (b[3]+b[4]*(X-x))/(b[1]+b[2]*(X-x)))) )*exp((-1/2)*((X-x)/a)^2) )
}
# 关于四个分量的导函数
gr <- function(b)
{
g <- rep(NA, 4)
g[1] <- (1/k)*sum( (1/(b[1] + b[2] * (X – x)) – (1 + (b[3] + b[4] * (X – x))^(-1)) * (Y *
    ((b[3] + b[4] * (X – x))/(b[1] + b[2] * (X – x))^2)/(1 + Y * ((b[3] +
    b[4] * (X – x))/(b[1] + b[2] * (X – x)))))) * exp((-1/2) * ((X –
    x)/a)^2) )
g[2] <- (1/k)*sum(((X – x)/(b[1] + b[2] * (X – x)) – (1 + (b[3] + b[4] * (X – x))^(-1)) *
    (Y * ((b[3] + b[4] * (X – x)) * (X – x)/(b[1] + b[2] * (X – x))^2)/(1 +
        Y * ((b[3] + b[4] * (X – x))/(b[1] + b[2] * (X – x)))))) * exp((-1/2) * ((X – x)/a)^2))
g[3] <- (1/k)*sum(((b[3] + b[4] * (X – x))^((-1) – 1) * (-1) * log(1 + Y * ((b[3] + b[4] *
    (X – x))/(b[1] + b[2] * (X – x)))) + (1 + (b[3] + b[4] * (X – x))^(-1)) *
    (Y * (1/(b[1] + b[2] * (X – x)))/(1 + Y * ((b[3] + b[4] * (X – x))/(b[1] +
        b[2] * (X – x)))))) * exp((-1/2) * ((X – x)/a)^2))
g[4] <- (1/k)*sum( ((b[3] + b[4] * (X – x))^((-1) – 1) * ((-1) * (X – x)) * log(1 +
    Y * ((b[3] + b[4] * (X – x))/(b[1] + b[2] * (X – x)))) + (1 + (b[3] +
    b[4] * (X – x))^(-1)) * (Y * ((X – x)/(b[1] + b[2] * (X – x)))/(1 +
    Y * ((b[3] + b[4] * (X – x))/(b[1] + b[2] * (X – x)))))) * exp((-1/2) *
    ((X – x)/a)^2))
g
}
hin <- function(b)
{
h <- rep( NA,(2*k)+2 )
h[1:k] <- b[1]+b[2]*(X-x)
h[(k+1):(2*k)] <- b[3]+b[4]*(X-x)
h[2*k+1] <- b[1]
h[2*k+2] <- b[3]
h
}
# initial value
p0 <- c(1,0,0.1,0)
gamma[j] <- ans <- constrOptim.nl(par=p0, fn=fn, gr=gr, hin=hin)$par[3]  # 只需返回b[3]即可
}
运行结果可以出来，就是和理论值误差较大，所以，希望大家看看程序是否正确，重点是constoptim()调用是否正确？万分感谢