# estimateLRdistn: Estimate "the" limiting distribution of the likelihood ratio... In logcondens.mode: Compute MLE of Log-Concave Density on R with Fixed Mode, and Perform Inference for the Mode.

## Description

Sampling from a given distribution, we estimate via Monte Carlo the limiting distribution of 2-log-likelihood-ratio of the modally-constrained log-concave MLE to the (unconstrained) log-concave MLE.

## Usage

 1 2 estimateLRdistn(rdist = rnorm, mode = 0, N.MC = 1e2, n.SS = 10000, xgrid=NULL, prec = 10^-10, seedVal = NULL, debugging = NULL) 

## Arguments

 rdist A function taking an integer argument n and returning n values simulated from a distribution. The distribution is generally log-concave (otherwise we are in a misspecified setting). mode fixed/known location of mode for constrained estimator. N.MC Number of Monte Carlo simulations to do for the limiting distribution. n.SS Sample Size used for each Monte Carlo. (Each MC simulates n.SS values from rdist and computes constrained and unconstrained MLE). xgrid Governs the generation of weights for observations. NULL then data are used as they are. Otherwise can be a single numeric of a numeric vector of length n.SS. Please see preProcess for details. prec Precision variable seedVal An optional seed value debugging Turns off/on debugging. Any non-character value turns debugging off. If debugging is a character string, then this string gives the name of an output file to which myxx (the simulated data from rdist), myxx.uniq (the corresponding unique values), and rdist are saved. If the code crashes, this can be examined. If debugging is on (i.e., is a character) then if TLLRs[i] is less than 0, the value of myxx will be saved to a file with name given by paste(debugging, "tmpxxs", i, ".rsav", sep=""), along with corresponding weights myww and the mode passed in.

## Details

Computes an estimate of the asymptotic distribution of the likelihood ratio statistic 2 (\mbox{log} \hat{f}_n - \mbox{log} \hat{f}_n^0) under the assumption that the true log-concave density f_0 satisfies f_0''(m)<0 where m is the true mode of f_0. The estimate is computed based on a sample of size n.SS from rdist via N.MC Monte Carlo iterations.

Note: the object LCTLLRdistn was created by output from this function with n.SS set to 1.2e3 and N.MC set to 1e4. Thus, estimateLRdistn is _NOT_ needed to simply compute fairly accurate quantiles of the limit distribution of the likelihood ratio statistic. estimateLRdistn is more useful for research purposes. For instance, by passing to mode values that are not the true mode of myr, the statistic can be studied under the alternative hypothesis.

## Value

A list(LRs,TLLRs), i.e., "likelihood ratio" and "two log likelihood ratios". Both are numeric vectors of length N.MC.

Note that theoretically all elements of LRs should be nonnegative, but in practice some rounding errors can occur when n.SS is very large.

## References

Duembgen, L, Huesler, A. and Rufibach, K. (2010) Active set and EM algorithms for log-concave densities based on complete and censored data. Technical report 61, IMSV, Univ. of Bern, available at http://arxiv.org/abs/0707.4643.

Duembgen, L. and Rufibach, K. (2009) Maximum likelihood estimation of a log-concave density and its distribution function: basic properties and uniform consistency. Bernoulli, 15(1), 40–68.

Duembgen, L. and Rufibach, K. (2011) logcondens: Computations Related to Univariate Log-Concave Density Estimation. Journal of Statistical Software, 39(6), 1–28. http://www.jstatsoft.org/v39/i06

Doss, C. R. (2013). Shape-Constrained Inference for Concave-Transformed Densities and their Modes. PhD thesis, Department of Statistics, University of Washington, in preparation.

Doss, C. R. and Wellner, J. A. (2013). Inference for the mode of a log-concave density. Technical Report, University of Washington, in preparation.

See activeSetLogCon and activeSetLogCon.mode, which compute the unconstrained and constrained MLEs, which form the likelihood ratio. The object LCTLLRdistn was created by output from this function.
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 myseed <- 561 {if(require(distr)){ mydistn <- Norm() ##demonstrate use of distr package myr <- mydistn@r } else { myr <- rnorm }} hypothesis.mode <- 0 N.MC <- 100 ## should increase these values for better estimate n.SS <- 50 LRres <- estimateLRdistn(rdist=myr, mode=hypothesis.mode, N.MC=N.MC, prec=10^-10, n.SS=n.SS, seedVal=myseed, debugging=FALSE) TLLRs <- sort(LRres\$TLLRs) ##sort is unnecessary, just for examining data negIdcs <- TLLRs<=0; ## rounding errors Nneg <- sum(negIdcs) print(Nneg) TLLRs[negIdcs] <- 0 cdf.empirical.f <- ecdf(TLLRs) xlims <- c(min(TLLRs), max(TLLRs)) xpts <- seq(from=xlims[1], to=xlims[2], by=.001) plot(xpts, cdf.empirical.f(xpts), type="l", xlab="TLLRs", ylab="Probability") #### LCTLLRdistn used 1e4 Monte Carlos with 1.2e3 samples each Monte ####Carlo. ##lines(xpts, LCTLLRdistn@p(xpts), col="blue") ## "object ##'C_R_approxfun' not found" error on winbuilder