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
How can I reduce the MSE of estimator especially for alpha "a" parameter?
