How to calculate a distribution out of a series of min-max values?

Question

For each minute of an hour, I have a minimum and maximum observed value.

for example:

time   min    max
00:00  12.13  15.10
00:01  14.23  17.02
00:02  11.12  12.10
...
00:59  09.11  09.89

We can assume that, during that minute, the observed values were uniformly distributed between the minimum and the maximum values.

1. How can I calculate the hourly distribution of the observed values out of the 60 min and max values. My goal is to calculate the x-percentile values.

2. Do you know if this type of calculation has a name I can search for or if it is implemented in a statistical software, such as R?

Answer 1

I doubt it's really the case that the data are uniformly distributed within each minute, but if you think that's a good approximation, an efficient solution can be had by means of binary searching.

Analysis

Each time represents a uniform distribution supported on the interval from $\text{min}$ to $\text{max}$. By definition, this has a cumulative distribution function starting at $0$, increasing linearly from $0$ to $1$ on the interval $[\text{min}, \text{max}]$, then remaining at $1$. Notice that it is a piecewise linear function. The average of these CDFs is the CDF of the hourly data. Because each of its components is piecewise linear, the overall CDF is piecewise linear: its graph will be a polyline. We can therefore represent the CDF as an ordered array of vertices of the graph. (In R, this will be a two-column array.) Obtaining a quantile means inverting this function: given a y-value, find the corresponding x-value. That is done by searching through the set of y-values, which must monotonically non-decrease within the graph, and linearly interpolating between the vertices found to surround the target value of $y$.

As an example, consider these data:

    Begin  End
[1,]    0    1
[2,]    1    3
[3,]    3    8
[4,]    6    5
[5,]    5    1
[6,]    1    2
[7,]    4    4
[8,]   10   12

(Some of the minima and maxima are interchanged just to show it will not matter.)

Notice there are no values between $8$ and $10$--the CDF must be horizontal here--and that in the penultimate interval the values are constant at $4$--the CDF will be vertical here. Here is a graph of it:

Solution

The R code to produce this plot is

x <- matrix(c(0,1,3,6,5,1,4,10, 1,3,8,5,1,2,4,12),ncol=2) # Sample data
cuts <- sort(matrix(x, ncol=1))                     # All endpoints in x
e <- apply(x, 1, function(x) edf(x, cuts))          # List of CDFs of rows of x
density <- apply(e, 1, sum) / dim(x)[1]             # Mixture density
cdf <- unique(cbind(cuts, density))                 # CDF array for the mixture

# Plot the CDF *by finding the inverse values*
y <- matrix(0:100/100, ncol=1)                      # Percents 0, 1, ..., 100
q <- apply(y, 1, function(q) quantile(q, cdf))      # The computed quantiles

plot(cdf, xlab="Values")                            # The CDF array
lines(q, y)                                         # Lines *from the quantiles*

The line segments in the plot were computed by the quantile function, not by connecting the points in the CDF array: this visually confirms the correctness of the solution across the full range of percentiles (at least for this dataset).

For example, to obtain the 75th percentile, execute

> quantile(0.75, cdf)

The output is 5.5. Indeed, looking across the graph at a height of $0.75$ shows a crossing around $5.5$. Looking back at the data, we can compute the chance of exceeding $5.5$:

• In the third interval $[3,8]$, there's a 50% chance.

• In the fourth interval $[6,5]$ there is also a 50% chance.

• In the eighth interval $[10,12]$ there is a 100% chance.

Because each interval occurs one-eighth of the time, the total chance of exceeding $5.5$ therefore equals $(.5 + .5 + 1.0)/8 = 0.25 = 1-0.75$, as required.

As another example, quantile(0.875, cdf) returns 10, whereas quantile(0.8749999999, cdf) returns 8. This shows how the quantile function selects the largest possible value it can. (This behavior can be modified by flipping the graph around both axes, making sure to keep the cdf array sorted: -quantile(-0.875, -cbind(rev(cdf[,1]), rev(cdf[,2]))) returns 8, the smallest possible answer.)

Details

This solution exploits two helper functions. The first, edf, computes the uniform CDF for a single interval:

edf <- function(x, cuts) {
    # x is a 2-vector representing an interval
    # cuts is a vector of possible interval endpoints, including those of x
    # Returns the cdf of a uniform density over x.
    e <- seq(0, 0, along.with=cuts)
    x.max = max(x); x.min = min(x)
    delta <- x.max - x.min
    if (delta > 0) {
        y <- (cuts - x.min) / delta
        i <- cuts < x.max & x.min <= cuts
        e[i] <- y[i]
        e[cuts >= x.max] <- 1
    } 
    else {
        s <- binsearch(x.max, cuts)
        e[s[1]:length(e)] <- 1
    }
    e
}

The second, quantile, performs the search and linear interpolation needed to invert a piecewise linear graph:

quantile <- function(q, cdf) {
    # 0 <= q <= 1
    # cdf is a piecewise linear CDF as a two-column matrix cbind(x, CDF(x)).
    # Returns a quantile for q.

    # Search within the range for q.
    y <- cdf[,2]
    s <- binsearch(q, y)

    # Prepare for linear interpolation.
    r <- y[s]
    delta = r[2] - r[1]
    if (delta == 0) {
        delta = 1
        r[2] <- r[1]+1
    }

    # Linearly interpolate the quantile.
    x <- cdf[,1]
    ((q - r[1]) * x[s[2]] + (r[2] - q) * x[s[1]]) / delta
}

Assisting these functions is a search routine, binsearch. I could not find a bug-free binary search (although no doubt some are lurking in the R repositories), so here is one quickly whipped up for the occasion. The tricky part is handling CDFs with jumps: this requires searching within an array that can have repeated values and selecting exactly the correct one (namely, the highest possible index that works).

binsearch <- function(x, s) {
    # x is the target.
    # s is sorted ascending (but can have repeated values).
    # Returns a pair of indexes into s, c(low, high), such that
    #     if min(s) <= x <= max(s), then
    #         s[low] <= x <= s[high] and 
    #             either high == length(s) or x < s[high+1].
    #
    # Output is undefined if x is outside the range of s (but termination
    # is assured).

    low <- 1; high <- length(s)
    while (high - low > 1) {
        mid = floor((high+low)/2)
        if (s[mid] <= x) low <- mid
        else high <- mid
    }
    c(low, high)
}