2
$\begingroup$

If I have a truncated histogram like this, where the ordinate is the number of occurrences $N_X$ of each value of $X$, and the abscissa is the value of $X$, is there a way to estimate the probability $P_X$ of each value of $X$?

enter image description here

If the histogram were not truncated below $X=10$, I could compute $N=\sum_X N_X$, and then estimate the probability of a value $X$ as $P_X = N_X/N$. Is there an approach to get $P_X$ for $X\geq10$ from the truncated histogram, short of assuming or estimating the scaling of $N_X$ and extrapolating to $X=0$?

$\endgroup$
8
  • 1
    $\begingroup$ Are you saying that you actually have values at $0<X<10$, but for some reason, they can't be measured/observed, so you'd like to extrapolate to values between $0 < X < 10$ (or estimate probabilities in that range)? $\endgroup$ Feb 22, 2019 at 22:37
  • 1
    $\begingroup$ OK. So just to make sure, I'm clear, you don't know the total number of observations at all (or total counts)? $\endgroup$ Feb 22, 2019 at 22:54
  • 1
    $\begingroup$ Correct. I only know the number of observations above the truncation point at $X=10$ in this minimal example $\endgroup$ Feb 22, 2019 at 22:55
  • 1
    $\begingroup$ Do you know if the values of $X$ are discrete or continuous? $\endgroup$ Feb 22, 2019 at 22:57
  • 1
    $\begingroup$ They are actually continuous, but are binned for the sake of computation. So in this example the bins run between multiples of $5$ $\endgroup$ Feb 22, 2019 at 22:58

1 Answer 1

2
$\begingroup$

There are a lot of different ways you could approach this. You could try to fit a polynomial or complex spline to your data, you could fit a a non-linear function, like an exponential function, to your data, or perhaps most simply, you could linearize the relationship between $N_X$ and $X$ through a transformation of your data. There are also some sampling approaches which do not make any distributional assumptions.

I'll highlight a few possibilities, and then explain some pitfalls and offer some suggestions.

"Linearize" the relationship between $N_X$ and $X$

Suppose you are only given the binned data, so you know the number of observations that fall between 10 - 15, 15 - 20, 20 - 25, .... Then, you could model the relationship between the bins and $N_X$ and once this relationship is known, and if we can assume it holds over the range $0<X<10$, then you could extrapolate (this can be dangerous -- see below) to estimate the number of observations in the bins 0 - 5, and 5 - 10. Based on the histogram you've provided it looks like it would be reasonable to assume an exponential relationship between $X$ and $N_X$ (you could try functional relationships too). We could then try to model the following:

\begin{eqnarray*} Y=log(N_{X}) & = & b+mX \end{eqnarray*}

where $b$ and $m$ are to be estimated by a linear regression regression on the transformed data $Y=log(N_X)$ by using the lower boundary of each bin as the input for $X$ (you could use the upper bound, or the midpoint, etc. -- this is just the method I've chosen to demonstrate this method). Then, once you've modeled the relationship based on values where $X>10$, you simply estimate $N_X$ for the $X$ values that are missing. You will then be able to derive the probabilities using a frequentist argument by taking $P_X=N_X/N$.

To make this more concrete, I'll demonstrate with some R code. First, let's just generate some truncated data since you didn't provide data in your question:

set.seed((543))
y <- rexp(5000, .05) + rnorm(5000, 0, 3)
#Force negative values to 0
y[y<0] <- 0
xbin <- y[y>10]
x <- cut(xbin, breaks=seq(10,200, 5), include.lowest=TRUE, labels=FALSE)
mydf <- data.frame(table(x), stringsAsFactors = FALSE)
mydf$x <- as.numeric(mydf$x)*5 + 5
names(mydf) < -c("x", "Nx")
head(mydf)

This results in a dataframe that looks like this:

  x  Nx
1 10 686
2 15 540
3 20 452
4 25 336
5 30 235
6 35 182
    .
    .
    .

As you can see, we are missing the bins that have a left-hand boundary corresponding to $X=0$ and $X=5$, so we must estimate the $N_X$ for these values. So, we can transform the known $N_X$ and then fit a linear regression to the logged values:

print(model1 <- lm(log(Nx) ~ x, data=mydf))

Call:
lm(formula = log(Nx) ~ x, data = mydf)

Coefficients:
(Intercept)            x  
    7.07181     -0.05261

So we have our estimates for $b=7.07181$ and for $m=-0.05261$. All that's left to do is substitute values for $X=0$ and $X=5$ into the estimated equation $Y=log(N_{X}) = 7.07181+-0.05261X$ and take the inverse function to obtain the estimated $N_X$ in the bins:

\begin{eqnarray*} N_{X} & = & e^{7.07181-0.05261X} \end{eqnarray*}

This can be accomplished with R's predict function. We use the predict function below and then sum the $N_X$ for the bin with left boundaries 0 and 5:

sum(exp(predict(model1, list(x=c(0, 5)))))
[1] 2084.038

So, rounding up, our estimate for the total number of observations between $0<X<10$ is given by approximately $2085$. We can check this against the number of observations below $X=10$ from our simulated data:

sum(y<=10)
[1] 1917

Not bad. We overestimated by $2085-1917=168$ observations. The graphic below shows a barplot of the binned frequency counts for the known bins as well as our estimates for the bins defined by 0 - 5 and 5 - 10. The blue line above each bar shows the actual simulated values as well.

enter image description here

We could then proceed by adding this number to the total number of observations above $X=10$ (this is known) and you'd obtain your estimate for $N$. All estimated probabilities can then be derived by taking $P_X=N_X/N$ as you suggested in your question.

Predicting $N_X$ from $X$ using Nonlinear Least Squares

You may be able to get better estimates using nonlinear least squares estimates for an exponential function. We'll use R's nls function to do this for us and fit a model of the following form:

\begin{eqnarray*} N_{X} & = & ae^{bX} \end{eqnarray*}

To fit a nonlinear least squares to the binned data, we must provide nls estimated values for $a$ and $b$. I think it makes sense to use our estimates for $m$ and $b$ from the transformed linear regression analysis to give us starting values for $a$ and $b$ here (you can to to estimate this visually, otherwise). So, we'll use $a=e^{7.07181}$ and $b=-0.05261$ for starting values:

print(model2<-nls(Nx~a*exp(b*x), data=mydf, start=list(a=exp(7.07181),b=-0.05261)))
Nonlinear regression model
  model: Nx ~ a * exp(b * x)
   data: mydf
         a          b 
1161.54763   -0.05093 
 residual sum-of-squares: 2782

Now, we can predict again for the binned values with a left-hand boundaries of 0 and 5 and sum the result:

sum(predict(m, list(x=c(0,5))))
[1] 2061.965

In this analysis, we overestimated the count by only $145$ observations: $2062-1917=145$.

Now, the problem with these methods is that they rely on the assumption that your data trends similarly over $0<X<10$ as it does over $X\ge 10$. You are having to take this on face value as there is no evidence to suggest that something weird doesn't happen to your data at $X=0$ (or some other value of $X<10$. In fact, some models have to be specially created to handle problems with zeros (see inflated poisson and hurdle models on CV for example). There is always some significant risk of extrapolating your model beyond the range of values for your "$X$-values". So, another approach that avoids this assumptions may be best:

Frequentist Sampling Approach

So to get around this extrapolation problem, and assuming you are only trying to estimate probabilities of being in a bin (as opposed to actual counts), you could try sampling methods. I'm not sure what your exact memory/disk limitations are, but one way to estimate the probabilities of being in a bin, for example would be as follows:

  1. Write all your data to disk, including $X>10$ and $X<10$.
  2. Initialize to zero counters for your bins (we'll call these $Z_X$, where $X$ is the lower boundary of the bin.
  3. Randomly sample a large number of observations, say, $n=2000$ observations from the values written to disk
  4. For each sampled value, determine, based on its value, to which bin it belongs and increment the counter by 1. For example, if the sampled value is 4.32, then set $Z_0=Z_0+1$.
  5. When completed, divide the value of each bins counter by the number of times you sampled. These are you estimated probabilities.

Here is some R code to demonstrate this algorithm:

z0 <- 0
z5 <- 0

for (i in 1:2000){
  y_i <- sample(y, 1)
  z0 <- z0 + as.numeric(y_i < 5)   
  z5 <- z5 + as.numeric(y_i >= 5 & y_i < 10)   
}
z0/2000  
z5/2000

sum(y < 5)/5000
sum(y >= 5 & y < 10)/5000

This gives estimated probabilities for bin 0 - 5 as $0.2195$ and for bin 5 - 10 as $0.1745$. This is pretty close to the actual probabilities of $0.2126$ and $0.1708$.

Mark and Recapture Methods

If you need to know the actual counts in these bins, as opposed to probabilities, you may have to rely on something like a mark and recapture technique to estimate the population size of each of your bins. I'm not sure what the exact limitations of your computing environment are, but the Lincoln-Peterson estimator described in the wikipedia page I provided could be adapted to your situation. Here's an example in R of estimating the 0 - 5 and the 5 - 10 bin using this method:

sample1 <- mydf2[sample(nrow(mydf2), 2000), ]
sample2 <- mydf2[sample(nrow(mydf2), 2000), ]


sample1$bin0 <- sample1$y < 5
sample2$bin0 <- sample2$y < 5
sample1$bin5 <- sample1$y >= 5 & sample1$y < 10
sample2$bin5 <- sample2$y >= 5 & sample1$y < 10

n_bin0 <- sum(sample1$bin0)
n_bin5 <- sum(sample1$bin5)
K_bin0 <- sum(sample2$bin0)
K_bin5 <- sum(sample2$bin5)

recaptured <- sample2[sample2$id %in% sample1$id,]
k_bin0 <- sum(recaptured$bin0)
k_bin5 <- sum(recaptured$bin5)

N_hat_bin0 <- ceiling(K_bin0*n_bin0/k_bin0)
N_hat_bin5 <- ceiling(K_bin5*n_bin5/k_bin5)
N_hat_bin0
N_hat_bin5
#Difference between estimate and reality for 0-5 bin:
N_hat_bin0 - sum(y < 5)
#Difference between estimate and reality for 5-10 bin:
N_hat_bin5 - sum(y < 10 & y >= 5)

This yields:

N_hat_bin0-sum(y < 5)
[1] -6
N_hat_bin5-sum(y < 10 & y >= 5)
[1] 78

So with this method, we underestimate the number of observations in the 0 - 5 bin by $6$ and overestimate the number of observations in the 5 - 10 bin by $78$. This leads to a net overestimation for $X<10$ of $72$. Once you have these counts in hand, you could simply add them to the counts where $X>10$ and use this as $N$ in the calculations for $P_X$.

As I said in the introduction, these are just some of the methods you could try, and many others exist which I didn't explain. Different techniques will be better/worse depending on your exact limitations and assumptions you may be able to make about your data.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.