1
$\begingroup$

Suppose that $y_1,y_2,.....,y_N$ form a random sample from the CSUG maximum distribution (like order statistic but on a random number of observations, $N$) with unknown parameter $\theta$, where the pdf of each observation has the following form $f(y|\theta)=\frac{\theta}{(1-\theta)(1-y)^2}$ for $0\leq y\leq 1-\theta$ .

I want to estimate $\theta$ using maximum likelihood estimation. (I am working in R)

I begin to make maximum likelihood estimator of parameter in the CSUG maximum distribution

I have a problem: how can I generate a sample belonging to this max distribution?

Can I find the distribution of the maximum?

$\endgroup$
  • $\begingroup$ I think you need to give us more detail and emphasise the statistical part of your question as advice about R programming is off-topic here but available on other sites. $\endgroup$ – mdewey Apr 29 '17 at 14:11
  • $\begingroup$ ok! i have random variable belong to distribution call "the correlated standard uniform geometric " which have pdf contain parameter theta i need to make maximum likelihood method in R to estimate theta $\endgroup$ – Mary Yousef Apr 29 '17 at 14:41
  • 1
    $\begingroup$ This is R coding question, must be closed as such $\endgroup$ – Aksakal Apr 29 '17 at 15:07
  • $\begingroup$ please edit to make your post clearer $\endgroup$ – Glen_b -Reinstate Monica Apr 30 '17 at 15:53
  • 1
    $\begingroup$ I assume this is the paper your question arises from: Hao, Jie and Godbole, Anant (2016) "Distribution of the Maximum and Minimum of a Random Number of Bounded Random Variables" Open Journal of Statistics, Vol.06 No.02, 12 pages Article ID:65869, 10.4236/ojs.2016.62023 ... is that correct? $\endgroup$ – Glen_b -Reinstate Monica Apr 30 '17 at 16:12
1
$\begingroup$

This paper seems to cover what you are interested in including parameter estimation.

If $y_1,\ldots,y_n$ is a sample from the CSUG distribution then the MLE of $\theta$ is $\hat{\theta}= 1 - y_{(n)}$ where $y_{(n)}$ is the max of the sample.

From your comments you want to use inverse transform. The necessary equations are below. From the histograms, it looks like this form of simulation resembles the pdf.

$$ \begin{align*} f(y) &= \dfrac{\theta}{(1-\theta)(1-y)^2}\mathbf{1}\left\{0\leq y \leq 1-\theta\right\}\\ F(y) &= \begin{cases} 1 & y > 1-\theta\\ \dfrac{\theta y}{(1-\theta)(1-y)} & 0 \leq y \leq 1-\theta\\ 0 &y<0 \end{cases}\\ F^{-1}(u) &= \dfrac{(1-\theta)u}{\theta + (1-\theta)u} \quad 0 \leq u \leq 1 \end{align*} $$

## Simulate n values from the CSUG max distribution 
## using inverse transformation
rcsug2 <- function(n, theta){
  u <- runif(n)
  return((1-theta) * u / (theta + (1-theta) * u))
}

par(mfrow=c(1,3))
## theta 0.75
theta <- 0.75
y_seq <- seq(0,1-theta,length.out = 1000)

y <- rcsug2(1000, theta)
hist(y, prob=TRUE, ylim=c(0,10), 
     xlim=c(0,1),
     main=paste("theta=0.75","\n","theta_hat=",1-max(y)))

lines(y_seq,theta/((1-theta) * (1 - y_seq)^2) * (y_seq <= 1 - theta), lwd=2)

## theta 0.5
theta <- 0.5
y_seq <- seq(0,1-theta,length.out = 1000)

y <- rcsug2(1000, theta)
hist(y, prob=TRUE, ylim=c(0,10), 
     xlim=c(0,1),
     main=paste("theta=0.5","\n","theta_hat=",1-max(y)))

lines(y_seq,theta/((1-theta) * (1 - y_seq)^2) * (y_seq <= 1 - theta), lwd=2)

## theta 0.1
theta <- 0.1
y_seq <- seq(0,1-theta,length.out = 1000)

y <- rcsug2(1000, theta)
hist(y, prob=TRUE, ylim=c(0,10), 
     xlim=c(0,1),
     main=paste("theta=0.1","\n","theta_hat=",1-max(y)))

lines(y_seq,theta/((1-theta) * (1 - y_seq)^2) * (y_seq <= 1 - theta), lwd=2)

enter image description here

EDIT:

You can optimize the log-likelihood using optim or optimize by writing out the log-likelihood as a function. Note the function should just take in the theta value

neg.log.lik <- function(theta){
  n <- length(y)
  return(-n*log(theta) + n * log(1-theta) + sum(2*log(1-y)) - log(sapply(theta, function(x){ all(y >= 0 & y <= 1-x)})))
}

Here I specified the negative log-likelihood so I can use optim or optimize directly. The key to the function is checking whether the y values are within their specified range. You HAVE to include this because this is the only part of the log-likelihood that involves both the parameter and the data. Using optim you get

> y <- rcsug2(1000,0.5) 
> optim(0.1, neg.log.lik)
$par
[1] 0.5000271

$value
[1] -756.1681

$counts
function gradient 
      78       NA 

$convergence
[1] 0

$message
NULL

Warning message:
In optim(0.1, neg.log.lik) :
  one-dimensional optimization by Nelder-Mead is unreliable:
use "Brent" or optimize() directly

Heeding the warning message we can also try optimize

> optimize(neg.log.lik, c(0,1))
$minimum
[1] 0.4999542

$objective
[1] -755.8762

But I want to show you that YOU DO NOT NEED TO USE OPTIM in this case. Consider this simulation:

mle1 <- NULL
mle2 <- NULL
for(i in 1:10000){
  y <- rcsug2(1000,0.75) 
  mle1 <- c(mle1, 1-max(y))
  mle2 <- c(mle2,  optimize(neg.log.lik, c(0,1))$minimum)
}

The theoretically derived MLE looks exactly the same as the one obtained using opimize:

> summary(mle1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.7500  0.7501  0.7501  0.7502  0.7503  0.7518 
> summary(mle2)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.7500  0.7500  0.7501  0.7502  0.7503  0.7518 

> par(mfrow=c(1,2))
> hist(mle1)
> hist(mle2)

enter image description here

But calculating the MLE using optimize takes CONSIDERABLY more time.

library(microbenchmark)
f1 <- function(){
  y <- rcsug2(10000,runif(1))
  1-max(y)
}
f2 <- function(){
  y <- rcsug2(10000,runif(1))
  optimize(neg.log.lik, c(0,1))$minimum
}
microbenchmark(f1,f2)
Unit: nanoseconds
 expr min lq  mean median uq  max neval
   f1   0  0  0.62      1  1    2   100
   f2   0  0 29.50      0  1 2550   100

That is 29.5/0.62 = 47.6 times faster

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ suppose that y1,y2,.....,yn from arandom sample from the CSUG maximum distribution (like order statistic but on random number N of obsrvation)with unkown theta .the pdf of each observation has the following form f(y|theta)=theta/(1-theta)(1-y)^2 for 0=<y=<1-theta . I want estimate theta by using maximum likelihood method in R $\endgroup$ – Mary Yousef Apr 29 '17 at 16:58
  • $\begingroup$ ok! i have this paper but i want to apply the estiamation part in R $\endgroup$ – Mary Yousef Apr 30 '17 at 15:31
  • $\begingroup$ the paper says that the MLE is $\hat{\theta} = 1 - y_{(n)}$ where $y_{(n)}$ is the maximum of your sample. $\endgroup$ – bdeonovic Apr 30 '17 at 16:10
  • $\begingroup$ how i can produce this estimator of theta in R program ? first idea i have , I must generating sample belong to the distribution in problem , second I wil use MLE in R programme !! Are you have any idea about R programme? $\endgroup$ – Mary Yousef Apr 30 '17 at 18:07
  • 1
    $\begingroup$ @Mary Note that bdeonovic already gave you the formula for the MLE in comments and then he edited his code to calculate the MLE of theta. It's all done already. Please spend a little effort to try to understand it. $\endgroup$ – Glen_b -Reinstate Monica May 2 '17 at 5:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.