3

estou fazendo um estudo com a distribuição Birnbaum Saunders bivariada. Para resolver meu problema preciso criar um código que estima os parâmetros de um modelo mistura. Gostaria de inserir o gradiente de minha função de verossimilhança no código abaixo via o pacote NumDeriv ou de outra forma, porém não estou conseguindo, alguém poderia me ajudar a criar esse gradiente via algum pacote ou código do R? Segue o código:

#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) 

1 Resposta 1

1

Importante: eu não tenho familiaridade com o a distribuição Birnbaum-Saunders. Todo o código é baseado no artigo de Kundu, Balakrishnan e Jamalizadeh (2010). O artigo não aborda o problema da censura dos valores da variável, então é possível que os resultados não sejam diretamente aplicáveis ao teu problema.

Comparando teu código para simulação dos dados com o artigo citado, verifiquei uma discrepância na passagem dos valores intermediários z1 e z2 para T1 e T2 (Eq. 8). Por isso, nas simulações abaixo, segui a proposta do artigo, pois os resultado diferem.

Na primeira função, os dados são simulados a partir dos parâmetros propostos e um tamanho de amostra N. A segunda função, para estimativa de máxima verossimilhança, tem como argumento um conjunto de dados X com duas colunas e um vetor de valores iniciais para os parâmetros phi -- portanto, com apenas dois valores.

A função de otimização apenas segue as equações propostas no artigo e utiliza a função optim para encontrar os MLEs de phi1 e phi2. No fim, os valores de mu1, mu2 e rho são calculados a partir das soluções em forma fechada. Procurei identificar o número das equações no artigo original. Cuidado: a função não verifica se os valores iniciais são positivos e o algoritmo de otimização empregado (padrão do optim, Nelder-Mead) não restringe o espaço dos parâmetros para números positivos.

simData <- function(mu1, phi1, mu2, phi2, rho, N){
  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.1<-(mu1/(1+(1/phi1)))*((1/2)*(sqrt(2/phi1))*z1+sqrt(1+((1/2)*    (sqrt(2/phi1))*z1)^2))^2
  #T2.1<-(mu2/(1+(1/phi2)))*((1/2)*(sqrt(2/phi2))*z2+sqrt(1+((1/2)*    (sqrt(2/phi2))*z2)^2))^2
  T1 <- phi1 * ((1/2) * mu1 * z1 + sqrt(((1/2) * mu1 * z1)^2 + 1))^2
  T2 <- phi2 * ((1/2) * mu2 * z2 + sqrt(((1/2) * mu2 * z2)^2 + 1))^2

  data.frame(T1, T2)
}

X <- simData(1.5, 1, 3, 8, 0.5, 1000)

mleBVBS <- function(X, phiInits) {
  #Funções s e r, definidas em (11)
  s <- function(data){
    apply(data, 2, mean)
  }

  r <- function(data){
    apply(data, 2, function(x) 1/mean(1/x))
  }

  # Função em forma fechada para computar mu, dado phi (9)
  muHat <- function(phi, s, r){
    sqrt(s/phi + phi/r - 2)
  }

  # Função em forma fechada para computar rho, dado phi (10)
  rhoHat <- function(phi){
    phi1 <- phi[1]
    phi2 <- phi[2]

    num <- sum(
      (sqrt(X[, 1]/phi1) - sqrt(phi1/X[, 1])) *
        (sqrt(X[, 2]/phi2) - sqrt(phi2/X[, 2]))
    )
    den <- sqrt(
      sum((sqrt(X[, 1]/phi1) - sqrt(phi1/X[, 1]))^2)
    ) * sqrt(
      sum((sqrt(X[, 2]/phi2) - sqrt(phi2/X[, 2]))^2)
    )

    num/den
  }


  sData <- s(X)
  s1 <- sData[1]
  s2 <- sData[2]

  rData <- r(X)
  r1 <- rData[1]
  r2 <- rData[2]
  n <- nrow(X)

  #Função a ser otimizada para obter phi1 e phi2 (12)

  profileLL <- function(phi, s1, s2, r1, r2, n, X){
    phi1 <- phi[1]
    phi2 <- phi[2]

    (-n)*log(muHat(phi1, s1, r1)) - n*log(phi1) -
      n*log(muHat(phi2, s2, r2)) - n*log(phi2) -
      (n/2)*log(1 - rhoHat(phi)) +
      sum(log(
        sqrt(phi1/X[, 1]) + (phi1/X[, 1])^(3/2)
      )) + sum(log(
        sqrt(phi2/X[, 2]) + (phi2/X[, 2])^(3/2)
      ))
  }

  fit <- optim(phiInits, profileLL, control=list(fnscale=-1),
               s1=s1, s2=s2, r1=r1, r2=r2, n=n, X=X)
  phi <- fit$par

  c(mu1=muHat(phi[1], s1, r1), phi1=phi[1], mu2=muHat(phi[2], s2, r2), phi2=phi[2], rho=rhoHat(phi))

}

mleBVBS(X, c(1, 5))
1
  • Obrigado @Erikson K.!
    – fsbmat
    Commented 27/05/2017 às 0:02

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .