我想求解带有约束条件的极大似然问题，用的是constoptim( )包，帮忙看看这个程序是否

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

                                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()调用是否正确？万分感谢
                       

扫码加我 拉你入群

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

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

分享0 收藏1 回帖

关键词：Optim 约束条件 Stop 极大似然 ONS 程序

求帮忙啊