Log-Concave Distribution Estimation

Introduction

Let \(n = 1\) and suppose \(x_i\) are i.i.d. samples from a log-concave discrete distribution on \(\{0,\ldots,K\}\) for some \(K \in {\mathbf Z}_+\). Define \(p_k := {\mathbf {Prob}}(X = k)\) to be the probability mass function. One method for estimating \(\{p_0,\ldots,p_K\}\) is to maximize the log-likelihood function subject to a log-concavity constraint , i.e.,

\[ \begin{array}{ll} \underset{p}{\mbox{maximize}} & \sum_{k=0}^K M_k\log p_k \\ \mbox{subject to} & p \geq 0, \quad \sum_{k=0}^K p_k = 1, \\ & p_k \geq \sqrt{p_{k-1}p_{k+1}}, \quad k = 1,\ldots,K-1, \end{array} \]

where \(p \in {\mathbf R}^{K+1}\) is our variable of interest and \(M_k\) represents the number of observations equal to \(k\), so that \(\sum_{k=0}^K M_k = m\). The problem as posed above is not convex. However, we can transform it into a convex optimization problem by defining new variables \(u_k = \log p_k\) and relaxing the equality constraint to \(\sum_{k=0}^K p_k \leq 1\), since the latter always holds tightly at an optimal solution. The result is

\[ \begin{array}{ll} \underset{u}{\mbox{maximize}} & \sum_{k=0}^K M_k u_k \\ \mbox{subject to} & \sum_{k=0}^K e^{u_k} \leq 1, \\ & u_k - u_{k-1} \geq u_{k+1} - u_k, \quad k = 1,\ldots,K-1. \end{array} \]

Example

We draw \(m = 25\) observations from a log-concave distribution on \(\{0,\ldots,100\}\). We then estimate the probability mass function using the above method and compare it with the empirical distribution.

suppressMessages(suppressWarnings(library(CVXR)))
set.seed(1)
## Calculate a piecewise linear function
pwl_fun <- function(x, knots) {
    n <- nrow(knots)
    x0 <- sort(knots$x, decreasing = FALSE)
    y0 <- knots$y[order(knots$x, decreasing = FALSE)]
    slope <- diff(y0)/diff(x0)
    
    sapply(x, function(xs) {
        if(xs <= x0[1])
            y0[1] + slope[1]*(xs -x0[1])
        else if(xs >= x0[n])
            y0[n] + slope[n-1]*(xs - x0[n])
        else {
            idx <- which(xs <= x0)[1]
            y0[idx-1] + slope[idx-1]*(xs - x0[idx-1])
        }
    })
}
## Problem data
m <- 25
xrange <- 0:100
knots <- data.frame(x = c(0, 25, 65, 100), y = c(10, 30, 40, 15))
xprobs <- pwl_fun(xrange, knots)/15
xprobs <- exp(xprobs)/sum(exp(xprobs))
x <- sample(xrange, size = m, replace = TRUE, prob = xprobs)

K <- max(xrange)
counts <- hist(x, breaks = -1:K, right = TRUE, include.lowest = FALSE,
               plot = FALSE)$counts
library(ggplot2)
ggplot() +
    geom_histogram(mapping = aes(x = x), breaks = -1:K, color = "blue", fill = "orange")

We now solve problem with log-concave constraint.

u <- Variable(K+1)
obj <- t(counts) %*% u
constraints <- list(sum(exp(u)) <= 1, diff(u[1:K]) >= diff(u[2:(K+1)]))
prob <- Problem(Maximize(obj), constraints)
result <- solve(prob)
pmf <- result$getValue(exp(u))

The above lines transform the variables \(u_k\) to \(e^{u_k}\) before calculating their resulting values. This is possible because exp is a member of the CVXR library of atoms, so it can operate directly on a Variable object such as u.

Below are the comparison plots of pmf and cdf.

library(ggplot2)
library(tidyr)
dens <- density(x, bw = "sj")
d <- data.frame(x = xrange, True = xprobs, Optimal = pmf,
                Empirical = approx(x = dens$x, y = dens$y, xout = xrange)$y)
plot.data <- gather(data = d, key = "Type", value = "Estimate", True, Empirical, Optimal,
                    factor_key = TRUE)
ggplot(plot.data) +
    geom_line(mapping = aes(x = x, y = Estimate, color = Type)) +
    theme(legend.position = "top")

d <- data.frame(x = xrange, True = cumsum(xprobs),
                Empirical = cumsum(counts) / sum(counts),
                Optimal = cumsum(pmf))
plot.data <- gather(data = d, key = "Type", value = "Estimate", True, Empirical, Optimal,
                    factor_key = TRUE)
ggplot(plot.data) +
    geom_line(mapping = aes(x = x, y = Estimate, color = Type)) +
    theme(legend.position = "top")

From the figures we see that the estimated curve is much closer to the true distribution, exhibiting a similar shape and number of peaks. In contrast, the empirical probability mass function oscillates, failing to be log-concave on parts of its domain. These differences are reflected in the cumulative distribution functions as well.

Session Info

sessionInfo()
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin17.7.0 (64-bit)
## Running under: macOS  10.14.1
## 
## Matrix products: default
## BLAS/LAPACK: /usr/local/Cellar/openblas/0.3.3/lib/libopenblasp-r0.3.3.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices datasets  utils     methods   base     
## 
## other attached packages:
## [1] tidyr_0.8.2   ggplot2_3.1.0 CVXR_0.99-2  
## 
## loaded via a namespace (and not attached):
##  [1] gmp_0.5-13.2      Rcpp_0.12.19      bindr_0.1.1      
##  [4] pillar_1.3.0      compiler_3.5.1    plyr_1.8.4       
##  [7] R.methodsS3_1.7.1 R.utils_2.7.0     tools_3.5.1      
## [10] digest_0.6.18     bit_1.1-14        evaluate_0.12    
## [13] tibble_1.4.2      gtable_0.2.0      lattice_0.20-35  
## [16] pkgconfig_2.0.2   rlang_0.3.0.1     Matrix_1.2-15    
## [19] yaml_2.2.0        blogdown_0.9.2    xfun_0.4         
## [22] bindrcpp_0.2.2    withr_2.1.2       dplyr_0.7.7      
## [25] Rmpfr_0.7-1       ECOSolveR_0.4     stringr_1.3.1    
## [28] knitr_1.20        tidyselect_0.2.5  rprojroot_1.3-2  
## [31] bit64_0.9-7       grid_3.5.1        glue_1.3.0       
## [34] R6_2.3.0          rmarkdown_1.10    bookdown_0.7     
## [37] purrr_0.2.5       magrittr_1.5      backports_1.1.2  
## [40] scales_1.0.0      htmltools_0.3.6   scs_1.1-1        
## [43] assertthat_0.2.0  colorspace_1.3-2  labeling_0.3     
## [46] stringi_1.2.4     lazyeval_0.2.1    munsell_0.5.0    
## [49] crayon_1.3.4      R.oo_1.22.0

Source

R Markdown

References