1
$\begingroup$

We have some set of data which has been already processed elsewhere and now contains bins with unequal size and counts of occurrence within each bin (around 8-10 bins, with approximately log-like stretching). We are going to analyze the distribution in R. Most of statistical packages rather deal with raw data, building histograms, deriving probability density distributions etc. However we do not have raw, but already pre-counted data. How it is possible

  1. to construct in R the probability density distribution function as histogram (bar chart)?
  2. To perform kernel smoothing to infer the smoothed distribution function?
  3. to visualize the cumulative function of total volume - again as histograms and/or as smoothed curve. In Mathematica there is a HistogramSmooth and allied to it functions.

The problem to 1) is that this is not just bare histogram with height proportional to counts (since bins are unequal). Particular application is grain size analysis (or equivalent objects, for example, distribution of oil fields worldwide over sizes). Here is the example data:

bins: 60-100, 100-200, 200-300, 300-500, 500-1000, 1000-2000, 2000-3000, >3000 
counts: 275, 320, 112, 65, 53, 44, 16, 15

Bins are ranges of size (volume), and counts are numbers of objects of each range. A special question refers to plotting (and programming in R) such data, if the last bin is stated just "above 3000" - how to infer its reasonable size and then height for visualization on histograms? Then, do I need to manually compute bins for some nonlinear scale on x, or some package in R effectively does this job including handling nonlinear scales?

$\endgroup$
2
  • 3
    $\begingroup$ Not sure there's much R does here to help. These are really statistical problems, not programming problems. You should probably ask such questions about binned data over at Cross Validated where statistical questions are on topic. This doesn't really seem to be a programming problem to me specifically. $\endgroup$
    – MrFlick
    Aug 18, 2017 at 18:15
  • $\begingroup$ @MrFlick Telling people to ask elsewhere (instead of explaining to flag for migration there) leads to them reposting, so that when the original is migrated, we get two copies of the same post and have to clean up. Please consider that when commenting; new users don't know the way things work, so we should explain such issues to them. $\endgroup$
    – Glen_b
    Aug 18, 2017 at 23:47

1 Answer 1

2
$\begingroup$

I do not know if you are still interested in this issue (this is an old question). Anyway my attempt to answer might be useful to somebody in the future.

If I read well the original problem, you would like to 'dis-aggregate a coarse histogram'. To make the problem more challenging, the original data come with an 'open group' (> 3000). As pointed in the question, the situation is somehow 'non-standard' (or at least not a text-book example) but this kind of situation is actually really common (in actuarial applications for example).

The way I would approach the problem is by using the penalized composite link model proposed by Eilers (2007). The idea presented in the paper follows, more or less, these lines

  • We assume that there is an unobserved (latent) vector of counts say $\gamma$ wich is defined for a finer set of bins
  • A matrix $C$ creates the aggregated counts as $\mu = C\gamma$
  • What we observe the $m$-dim vector $y$ (i.e. the counts of element 'falling' in a set of bins) are realizaitons of a Poisson process with $E(y) = \mu$ and hence $y \sim \mbox{Pois}(\mu)$
  • The matrix $C$ (dim $m \times n$ with $n >> m$) describes how the latent distribution was mixed before generating the data (the Rcode below should clarify how to build it).

In order to estimate the elements of the latent vector of counts $\gamma$ we need make some extra assumption

  • The vector $\gamma$ is smooth and $\gamma = \exp(B \beta)$: so we just need to impose that $\beta$, a vector of coefficients, is smooth
  • The matrix $B$ is a B-spline matrix built on a set of equally spaced knots (see e.g. Eilers and Marx, 1991)
  • The likelihood is $$ \ell = \sum_{i} (y_{i} \log \mu_{i} - \mu_{i}) - \lambda P $$
  • $P$ is a penalty matrix obtained as $P = D^{\top}D$
  • $D$ is a finite difference matrix of, say, order 2 (see e.g. Eiles and Marx, 1991)
  • $\lambda$ is a regularization parameter to be chosen/selected (in Eilers, 2007 the author uses the AIC criterion to select this parameter if I remember well)

These are the general setting of the solution I would adopt for the original problem. Fortunately, we do not need to do the computations by ourselves. We can use the Rpackage JOPS instead which contains the function pclm (you can download this package here: https://psplines.bitbucket.io/).

Below you will find a code that reproduces your example and plots the inferred latent density for a fixed value of $\lambda$ (you can follow the reference I mentioned before to estimate the optimal value of this parameter). Please also notice that, in the code below, I make the assumption that the class $> 3000$ goes from $3000$ to $4500$. You can change 4500 with the value you prefer (as long as it is a reasonable one for your problem).

    library(JOPS)
    library(colorout)
    library(ggplot2)
    
    # Data
    xl  = c(60, 100, 200, 300, 500, 1000, 2000, 3000) 
    xr  = c(100, 200, 300, 500, 1000, 2000, 3000, 4500)
    y   = c(275, 320, 112, 65, 53, 44, 14, 15)
    
    # Set-ups
    m      = length(y)
    n      =  xr[m]
    widths = xr - xl + 1
    dens   = y/widths/sum(y)
    
    # Composition matrix C  1 if finer bin is included in the 
    # coarse one, 0 otherwise
    C = matrix(0, m, n)
    for(i in 1:m) C[i, xl[i]:xr[i]] = 1
    C = C[, -c(1:(min(xl)-1))] 
    
    # Prepare B-splines bases
    x = (n - ncol(C)+1):n
    B = JOPS::bbase(x)
   
    # Fit
    mod = JOPS::pclm(y, C, B, lambda = 1, pord = 2, show = T)
    
    # Plot
    fit = data.frame(x = x, fits = mod$gamma / sum(mod$gamma))
    dat = data.frame(xl = xl, xr = xr, y = dens)
    plt = ggplot(dat, aes(ymin = 0)) + 
      geom_rect(aes(xmin = xl, xmax = xr, ymax = y, fill = 
        I('white'), color = I('black'))) +
      geom_line(aes(x = x, y = fits), data = fit, col = I('red'))
      +  ylab('Density') + theme_bw()
    
    plt

enter image description here

To conclude, I hope my attempt replies, at least partially, to the original question (and is clear enough ^_^).

$\endgroup$

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.