Estimation of parameters via monte carlo and optimum function

Navigation
3

I'm doing a bivariate Birnbaum Saunders distribution study. To solve my problem I need to create a code that estimates the parameters of a blend model. I would like to insert the gradient of my likelihood function in the code below via the NumDeriv package or otherwise, however I am not getting it, could anyone help me to create this gradient via some package or R code? Here is 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)
    
asked by anonymous 12.10.2016 / 02:00

1 answer

1

Important: I'm not familiar with the Birnbaum-Saunders distribution. The entire code is based on the Kundu, Balakrishnan and Jamalizadeh (2010) article. The article does not address the problem of censoring variable values, so it is possible that the results are not directly applicable to your problem.

Comparing your code to simulate the data with the article quoted, I checked a discrepancy in passing the intermediate values z1 and z2 to T1 and T2 (Eq. 8). Therefore, in the simulations below, I followed the proposal of the article, because the results differ.

In the first function, the data are simulated from the proposed parameters and a sample size N. The second function, for maximum likelihood estimation, has as its argument a set of data X with two columns and a vector of initial values for the phi parameters - therefore, with only two values.

The optimization function only follows the equations proposed in the article and uses the optim function to find the MLEs of phi1 and phi2 . At the end, the values of mu1 , mu2 and rho are calculated from the solutions in closed form. I tried to identify the number of equations in the original article. Caution: The function does not check if the initial values are positive and the optimization algorithm used (default optim , Nelder-Mead) does not restrict the parameter space to positive numbers. 

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))
    
27.05.2017 / 01:41
How to apply the_content filter in function get_the_content WordPress? Accumulated Interest Quarterly using R