jeudi 24 février 2022

Quantile function error and high mean square error (MSE)

I want to generate random sample from the following discrete distribution, so I used the CDF technique i.e F(x)=p So, I

#The CDF of continuous distribution is as follows
cdf.con<- function(q, t, a, lower.tail = TRUE, log.p = FALSE)
{
  stopifnot(t > 0, a > 0)
  p <- t / (t + 1)
  if(lower.tail)
  {
    cdf <- p * pexp(q, rate = t, lower.tail = TRUE) + (1 - p) * pgamma(q, shape = a, rate = t, lower.tail = TRUE)
  }
  else
  {
    cdf <- p * pexp(q,  rate = t, lower.tail = FALSE) + (1 - p) * pgamma(q, shape = a, rate = t, lower.tail = FALSE)
  }
  if(log.p) return(log(cdf)) else return(cdf)
}

#PMF 
#************************
pmf<- function(x, t, a, log = FALSE)
{
  stopifnot(t > 0, a > 0, x >= 0)
  if(!isTRUE(all(x == floor(x)))) stop("'x' must only contain nonnegative integers")
  if(log)
  {
    log(cdf.con(q = x, t, a, lower.tail = FALSE) - cdf.con(q = x + 1, t, a, lower.tail = FALSE))
  }
  else
  {
    cdf.con(q = x, t, a, lower.tail = FALSE) - cdf.con(q = x + 1, t, a, lower.tail = FALSE)
  }
}

#CDF $ (1-CDF) of discrete random variable
#************************************************
cdf.dic<- function(q, t, a, lower.tail = TRUE, log.p = FALSE)
{
  if(lower.tail)
  {
    cdf <- sapply(1 : length(q), function(i) sum(pmf(x = 0 : q[i], t, a, log = FALSE)))
  }
  else
  {
    cdf <- 0.1e1 - sapply(1 : length(q), function(i) sum(x = pmf(0 : q[i], t, a, log = FALSE)))
  }
  if(log.p) return(log(cdf)) else return(cdf)
}



# Generate random sample
#******************************
#Quantile of Cont.
qntf.con<- function(p, t, a, lower.tail = TRUE, log.p = FALSE, L = 1e-4, U = 50)
{
  stopifnot(t > 0, a > 0)
  if(lower.tail)
  {
    fx <- function(p)
    {
      tryCatch(uniroot(function(q) cdf.con(q, t, a, lower.tail = TRUE, log.p = FALSE) - p, lower = L, upper = U)$root, error = function(e) NaN)
    }
    qtf <- sapply(p, fx)
    if(log.p) return(log(qtf)) else return(qtf)
  }
  else
  {
    fx <- function(p)
    {
      tryCatch(uniroot(function(q) cdf.con(q, t, a, lower.tail = FALSE, log.p = FALSE) - p, lower = L, upper = U)$root, error = function(e) NaN)
    }
    qtf <- sapply(p, fx)
    if(log.p) return(log(qtf)) else return(qtf)
  }
}

I generate random sample from discrete distribution by making (floor ) to quantile function of continuous distribution as follows

# Quantiles and random sample generator from discrete r.v
#******************************************************
qntf.dis<- function(p, t, a, lower.tail = TRUE, log.p = FALSE)
{
  stopifnot(t > 0, a > 0)
  if(lower.tail)
  {
    qtf <- floor(qntf.con(p, t, a, lower.tail = TRUE, log.p = FALSE))
  }
  else
  {
    qtf <- floor(qntf.con(p, t, a, lower.tail = FALSE, log.p = FALSE))
  }
  if(log.p) return(log(qtf)) else return(qtf)
}

#generate r.s from discrete distribution as follows
r.s<- function(n, t, a) # For discrete r.v
{
  stopifnot(t > 0, a > 0)
  
  qntf.dis(p = runif(n), t, a, lower.tail = TRUE, log.p = FALSE)
}

#generate random sample from discrete r.v, I got NAN in some values
set.seed(110)
r.s(n=200, t=2, a=0.42) #produce NAN
#For the following set of parameters, I did not get NAN 
set.seed(110)

r.s(n=2000, t=0.5, a=2) # do not produce NAN

can someone tell me why NAN appears ? and how can we fix this?

My Second question is, I got high Mean square error

 myx <- r.s(n, t=theta1[1], a=theta1[2]) # generate sample from discrete distribution
    myx<- myx[!is.na(myx)]
# the log LHF
####################
loglik <- function(param) {
      t <- param[1]
      if(t<=0){return(NaN)}
      a <- param[2]
      if(a<=0){return(NaN)}
      sum(pmf(myx, t=t, a= a, log=TRUE))
      # log likelihood function is sum(ln(probability mass function))
    }# logLikFun
    library(maxLik)
    mle <-  maxLik(loglik, start=c(t= 0.5, a= 2))
    summary( mle )



   #Compute MSE
###########################
        #to compute it first compute the bias
        bias_theta= coef(mle)[1]- 0.5
        bias_alpha= coef(mle)[2]- 2
        MSE= c(vcov(mle)[1] + bias_theta^2, vcov(mle)[4] + bias_alpha^2)
        
    MSE
          t      a 
    1.824092 2.845212 



Aucun commentaire:

Enregistrer un commentaire