1
$\begingroup$

Question in short: How can I make two histograms equal by explicitly downsampling the buckets?

More detailed version: I have realizations of sequences of two discrete random variables both giving values in some finite space $X = \{1, ..., n\}$, i.e. I have two data vectors $(x_1, x_2, ..., x_A)$ and $(y_1, y_2, ..., y_B)$ with $x_i, y_j \in X$. I consider their distributions in terms of histograms, i.e. for each value $i \in X$ I count $a_i = |\{u \in \{1,...,A\} : x_u = i\}|$ and $b_i = |\{v \in \{1,...,B\} : y_v = i\}|$. Now I draw two histograms in percentage, i.e. the height of bar $i$ is $A_i = a_i/(a_1 + ... + a_n)$, respectively $B_i = b_i/(b_1 + ... + b_n)$. In general, these histograms look different.

How can I make the $A_i$ roughly equal to the $B_i$ by explicitly taking out some (as less as possible) of the data points $x_i$, respectively $y_j$?

More precisely, I want a solution that gives a tradeoff between 'how equal the histograms look alike' and how much data points I took out. It does not make sense for me to take out neither 90% of the $x_i$ nor 90% of the $y_j$.

Results so far:

For one single variable and a given fixed 'desired' histogram I succeeded more or less by doing the following (working R example below). Note: If one does the same thing for two variables (explicitly writing out what one wants as a system of equations), then the system becomes an optimization problem with quadratic constraints which is utterly difficult to solve in general.

CAUTION: even if the desired probabilities (i.e. heights of the bars in the desired histogram) come from another data vector, the answer is not as simple as 'take out as much until you arrive at the size of the desired bucket' as the bucket size of the given distribution is lower than the desired one (in order to demonstrate that I take out some elements from the first bucket in the R code below).

Given target percentages $p_i \in (0,1)$ with $p_1 + ... + p_n = 1$, find the minimal numbers $s_i \in \mathbb{N}_0$ such that $$\frac{a_i - s_i}{a_i - s_1 + ... + a_n - s_n} = p_i$$

rewriting $w_i = a_i - s_i$, this yields a linear system of equations $$ w_i = p_1w_1 + ... + p_nw_n$$ i.e. $$ p_1w_1 + p_{i-1}w_{i-1} +(p_i-1)w_i + p_{i+1}w_{i+1} + ... + p_nw_n = 0$$ which I solve [meaning that I obtain one single solution $w = (w_1, ..., w_n)$, recompute the $s_i = a_i - w_i$ and then just round. The tradeoff is done by searching for the minimal $s_i$, i.e. the system of equations is not uniquely solveable but is one dimensional, so we can adjust a factor $\lambda$ in front of the $w_i$, i.e. put the conditions $$ 0 \leq s_i = a_i - \lambda w_i \leq a_i$$ on $\lambda$ and then take the $\lambda$ that makes the $s_i$ as small as possible. Although this is not 100% correct, it yields good results:

Before downsampling:

enter image description here

After downsampling:

enter image description here

R code:

set.seed(1234)
options(scipen=99)
desired = rexp(1000, rate=1)*3000
desired = desired[desired >= 0 & desired <= 15000]

dTest = hist(desired, plot = F, breaks = seq(from=0, to=15000, length.out = 8))
buckets = 1:(length(dTest$breaks)-1)
probVec = c(0, 0, 0.1, 0.1+0.15, 0.1+0.1, 0.3)
probVec = c(probVec, 1-sum(probVec))
addBuckets = sample(buckets, size = 1000, replace=T, prob = probVec)
actual = c(desired, dTest$mids[addBuckets])
indicesFirstBucket = which(actual >= dTest$breaks[1] & actual <= dTest$breaks[2])
removeIndices = sample(indicesFirstBucket, size=100)
keepIndices = 1:length(actual)
keepIndices = keepIndices[!keepIndices %in% removeIndices]
actual = actual[keepIndices]
dControl = hist(actual, plot = F, breaks = dTest$breaks)


mfRowBefore = par("mfrow")
dTest$density = dTest$counts/sum(dTest$counts)
dControl$density = dControl$counts/sum(dControl$counts)
par(mfrow=c(1,2))
plot(dTest, freq=F, ylim=c(0, 0.6))
plot(dControl, freq=F, ylim=c(0, 0.6))
par(mfrow=mfRowBefore)


mfRowBefore = par("mfrow")
dTest$density = dTest$counts/sum(dTest$counts)
dControl$density = dControl$counts/sum(dControl$counts)
par(mfrow=c(1,2))
plot(dTest, ylim=c(0, 800))
plot(dControl, ylim=c(0, 800))
par(mfrow=mfRowBefore)

# sets up and solves the linear system of equations
helper_makeDistributionsEqual_getW = function(q) {
  n = length(q)
  m = matrix(0, n, n)  
  for (i in 1:n) {
    for (j in 1:n) {
      m[i,j] = if (i==j) (q[i]-1) else q[i]
    }
  }

  w<-eigen(m)
  vals = w$values
  # columns are eigenvectors:
  # m%*%vecs ~= diag(vals)%*%m
  vecs = w$vectors

  kernelIndices = which(abs(vals) < 0.000001)
  kernelBasis = t(t(vecs[,kernelIndices]))
  if (ncol(kernelBasis) > 1) {
    print("Uh oh... we have a >2 dim kernel...")
    browser()
  }
  kernelVector = kernelBasis[,1]
  return(kernelVector)
}

# given two binned densities 'upper' u and 'lower' l,
# try to find minimal integral positive a such that
# u-a has the same (binned!) distribution as l
#
# put q_i = l_i/sum(l_i)
# then we want that
#
# u_i - a_i
# --------------   ==   q_i
# sum(u_i - a_i)
#
# This is a linear system of equations in the variables w_i = u_i - a_i
# giving a 1-D solution space parametrized by some nonzero vector w
# Now 0 <= a_i = u_i - lambda*w_i iff. u_i >= lambda*w_i
#  iff. 
#         for all i, ... lambda <= u_i/w_i   if w_i > 0
#                        lambda >= u_i/w_i   if w_i < 0   (*)
#                        -no condition-      if w_i = 0
#
# Further, it surely does not make sense if a_i > u_i (as we run into negative counts then!)
# --> we want
#       u_i >= a_i = u_i - lambda w_i   <=>  -lambda w_i <= 0  
#                                       <=> lambda and w_i have the same sign
#
# --> we have a problem if w has different signs! but if it has the same signs in
#     every entry then we are fine (and I guess that one can show that 
#     for this particular linear problem, the kernel is 1D and all
#     elements in the kernel have either positive signin each component or 
#     negative sign in every component)
#
Global.makeDistributionsEqual = function(lower, upper) {
  # sometimes there are buckets that are extremely small but put a lot of
  # restrictions on the process (i.e. then it degenerates
  # and subtracts really a LOT of data points from the huge buckets 
  # which is not what we want: subtract only a few from the large buckets
  # even if you 'miss' [i.e. do not correctly reduce] extremely small buckets!)
  # --> only take those buckets covering 95% of the data


  q = lower / sum(lower)


  w = helper_makeDistributionsEqual_getW(q)
  whichAreZero = which(w == 0)
  whichAreNonZero = which(w != 0)
  res = rep(0, length(upper))

  lower = lower[w != 0]
  upper = upper[w != 0]
  w = w[w != 0]


  if (length(unique(sign(w))) > 1) {
    print("Uh oh: sign(w) was not unique!")
    # This really is an error:
    # we cannot satisfy (**) then!!!
    browser()
  }
  if (unique(sign(w)) < 0) {
    w = -w
  }

  # compute upper and lower bounds for lambda w.r.t. (*)
  testUpper = upper[w > 0]/w[w > 0]

  # we dont need this anymore as w >= 0!
  # testLower = upper[w < 0]/w[w < 0]



  #lowerBound = -Inf
  #if (length(testLower) > 0) {
  #  lowerBound = max(testLower)
  #}

  # (**) now just means that -lambda w_i <= 0
  #     iff. lambda w_i >= 0
  #     iff lambda >= 0 (as w_i >= 0!!!)
  # --> lower bound is just 0
  lowerBound = 0

  #upperBound = Inf
  #if (length(testUpper) > 0) {
  # # as w != 0 and w >= 0 we dont need to ask this here anymore!
  upperBound = min(testUpper)
  if (upperBound == 0) {
    # this means that we empty all the buckets
    print("Uh oh... there is something wrong here!!!")
    browser()
  }
  #}

  # lambda must be in the interval [lowerBound, upperBound]
  # (with round instead of '[' brackets if upper resp lowerBound = +-Inf)
  # but caution: clearly, the whole thing does not make sense if
  # a_i > u_i (because then we end up with negative counts!)
  # --> 0 <= a_i <= u_i
  # (*) ensures that 0 <= a_i is satisfied
  # the condition for a_i <= u_i is
  # u_i >= a_i = u_i - lambda w_i
  #
  # 0 >= -lambda w_i
  #
  # 


  # now u_i >= a_i >= 0
  # for all lambda but how to ensure integrality?
  # a_i = u_i - lambda w_i is supposed to be an integral number
  # iff. lambda satisfies the relation
  #   upperBound >= lambda = (u_i - a_i)/w_i >= 0
  # iff.
  #   w_i * upperBound >= (u_i - a_i) >= 0  [w_i >= 0]
  # iff.
  #   w_i*upperBound - u_i >= -a_i >= -u_i
  # iff.
  #   u_i - w_i*upperBound <= a_i <= u_i
  #
  #
  a = ceiling(upper - w*upperBound)
  a[a < 0] = 0
  a[a > upper] = upper[a > upper]

  res = rep(0, length(q))
  res[whichAreNonZero] = a
  removedPerc = a/upper
  return(res)
}

a = Global.makeDistributionsEqual(dTest$counts, dControl$counts)

dControl$counts = dControl$counts - a


mfRowBefore = par("mfrow")
dTest$density = dTest$counts/sum(dTest$counts)
dControl$density = dControl$counts/sum(dControl$counts)
par(mfrow=c(1,2))
plot(dTest, freq=F, ylim=c(0, 0.6))
plot(dControl, freq=F, ylim=c(0, 0.6))
par(mfrow=mfRowBefore)


mfRowBefore = par("mfrow")
dTest$density = dTest$counts/sum(dTest$counts)
dControl$density = dControl$counts/sum(dControl$counts)
par(mfrow=c(1,2))
plot(dTest, ylim=c(0, 800))
plot(dControl, ylim=c(0, 800))
par(mfrow=mfRowBefore)
$\endgroup$

Your Answer

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

Browse other questions tagged or ask your own question.