3
I’m doing a study with the bivariate Birnbaum Saunders distribution. To solve my problem I need to create a code that estimates the parameters of a mixed model. Would you like to insert the gradient of my verisimilitude function in the code below via the Numderiv package or otherwise, but I’m not getting it, could someone help me create this gradient via some package or R code? Follow the code:
#Parâmetros iniciais para geração das variáveis aleatórias T1 e T2
###################################
#Tamanho da Amostra
###############################
n<-100
#####################################
#Número de amostras
########################################
N=100
#####################################################################
#Matriz para receber os valores estimados em cada passo do Monte Carlo
#######################################################################
m=matrix(ncol=4,nrow=N)
#########################################################
#Inicio do loop Monte carlo
for (i in 1:N){
mu1<-1.5
phi1<-1
mu2<-2
phi2<-8
rho<-0.3
u1 <-rnorm(n)
u2 <-rnorm(n)
z1 <-(((sqrt(1+rho))+(sqrt(1-rho)))/2)*u1+(((sqrt(1+rho))-(sqrt(1- rho)))/2)*u2
z2 <-(((sqrt(1+rho))-(sqrt(1-rho)))/2)*u1+(((sqrt(1+rho))+(sqrt(1- rho)))/2)*u2
#################################################################
#Variáveis (T1,T2)~BSB(mu1,phi1,mu2,phi2,rho)
##########################################################
T1<-(mu1/(1+(1/phi1)))*((1/2)*(sqrt(2/phi1))*z1+sqrt(1+((1/2)* (sqrt(2/phi1))*z1)^2))^2
T2<-(mu2/(1+(1/phi2)))*((1/2)*(sqrt(2/phi2))*z2+sqrt(1+((1/2)* (sqrt(2/phi2))*z2)^2))^2
################################################################
#Variável indicadora
######################################
u<-1*(T2>1)
#######################################################
#Variável de interesse cuja densidade é obtida apartir da
#distribuição Birnbaum Saunders bivariada, observe que esta
#é uma variável com censura e de acordo com a minha especificação
#dos parâmetros a censura é cerca de 10% a 20%.
############################################################
y<-T1*u
sum(1*(y>0))
rm(mu1,phi1,mu2,rho)
#######################################################
#Função Log de Verossimilhança
#Observe que o phi2 foi considerado fixo para garantir a
#identificabilidade do modelo!
############################################################
log.lik<-function(par){
mu1 <-abs(par[1])
phi1 <-abs(par[2])
mu2 <-abs(par[3])
arho <-par[4]
f1<-(1/(2*sqrt(2*pi)))*exp((-phi1/4)*(sqrt(y[y>0]*(phi1+1)/(phi1*mu1))-sqrt(phi1*mu1/(y[y>0]* (phi1+1))))^2)* (((sqrt(phi1+1))/(sqrt(phi1*mu1*y[y>0])))+ ((sqrt(phi1*mu1))/((y[y>0]^3)*(phi1+1))))* (sqrt(phi1/2))*pnorm(((sqrt(phi2*(phi2+1))* (phi2*mu2-phi2-1))/(sqrt(2)* (phi2+1)*sqrt(phi2*mu2*(1-tanh(arho)^2))))+ (((sqrt(phi1)*tanh(arho))/(sqrt(2*(1-tanh(arho)^2))))*(sqrt(y[y>0]* (phi1+1)/(phi1*mu1))-sqrt((phi1*mu1)/(y[y>0]* (phi1+1))))))
f2<-pnorm(sqrt(phi2/2)*(sqrt((phi2+1)/(phi2*mu2))-sqrt((phi2*mu2)/(phi2+1))))
lv<--(sum(log(f1))+sum(1*(y==0)*log(f2)))
lv
}
#######################################################
#Tentativa de se criar o gradiente (Não deu certo!)
############################################################
require(numDeriv)
# grad<-function(par){
# mu1 <-par[1]
# phi1<-par[2]
# mu2 <-par[3]
# rho <-par[4]
# f1<-(1/(2*sqrt(2*pi)))*exp((-phi1/4)*(sqrt(y[y>0]* (phi1+1)/(phi1*mu1))- sqrt(phi1*mu1/(y[y>0]*(phi1+1))))^2)*(((sqrt(phi1+1))/(sqrt(phi1*mu1*y[y>0])))+ ((sqrt(phi1*mu1))/((y[y>0]^3)*(phi1+1))))*(sqrt(phi1/2))*pnorm(((sqrt(phi2* (phi2+1))*(phi2*mu2-phi2-1))/(sqrt(2)*(phi2+1)*sqrt(phi2*mu2*(1-rho^2))))+ (((sqrt(phi1)*rho)/(sqrt(2*(1-rho^2))))*(sqrt(y[y>0]*(phi1+1)/(phi1*mu1))- sqrt((phi1*mu1)/(y[y>0]*(phi1+1))))))
# f2<-pnorm(sqrt(phi2/2)*(sqrt((phi2+1)/(phi2*mu2))- sqrt((phi2*mu2)/(phi2+1))))
# lv<--(sum(log(f1))+sum(1*(y==0)*log(f2)))
# g1<-(sum(1*(y>0)*D(log(f1),"mu1"))+sum(1*(y==0)*D(log(f2),"mu1")))
# g2<-(sum(1*(y>0)*D(log(f1),"phi1"))+sum(1*(y==0)*D(log(f2),"phi1")))
# g3<-(sum(1*(y>0)*D(log(f1),"mu2"))+sum(1*(y==0)*D(log(f2),"mu2")))
# g4<-(sum(1*(y>0)*D(log(f1),"rho"))+sum(1*(y==0)*D(log(f2),"rho")))
# gd=cbind(eval(g1),eval(g2),eval(g3),eval(g4))
# colSums(gd, na.rm = TRUE)
# }
#######################################################
#Chutes iniciais
############################################################
mu1_0 <-1.5
phi1_0<-1
mu2_0 <-2
arho_0<-atanh(0.3)
start<-c(mu1_0,phi1_0,mu2_0,arho_0)
#######################################################
#Estimação via método BFGS da função optim
############################################################
op<-optim(start,log.lik,method = "BFGS")
#######################################################
#Como fiz uma reparametrização do parâmetro rho afim de não
#ter problemas devido a simulação dos valores a entrada da
#matriz referente a rho recebe tangente hiperbólico da estimativa de rho
#via função optim!
############################################################
m[i,]<-c(op$par[1:3],tanh(op$par[4]))
}
colMeans(m)
Thank you @Erikson K.!
– fsbmat