Calculating the probability of an event being early or late by a certain amount I'm very new to this so I assume there's a better way to ask this that I just don't know about.  Please point me in the right direction.
Suppose I'm running a library and books are due 30 days after they're borrowed.  And there's a library member named Joe, and I want to figure out how likely it is that he'll return his books on time, and early or late by certain amounts.  Joe has returned his books on the following intervals after borrowing them: [10, 24, 27, 27, 28. 29, 29, 29, 30, 30, 30, 31, 31, 32, 34].  This is just an example and I'd probably be looking at more 
How can I figure out the probability that Joe will have returned his books on-or-before each day 25-35.  My naive attempt would be to sum up the number of returns by preceding intervals like this:


*

*25: 02/15

*26: 02/15

*27: 04/15

*28: 05/15

*29: 08/15

*30: 11/15

*31: 13/15

*32: 14/15

*33: 14/15

*34: 15/15

*35: 15/15


Is the a "smarter" way to do this, and possibly one that would work with less-discrete numbers than whole days?  Ideally I think I'd like to end up with some kind of continuous function so that I could track the returns down to the second and it would be able to tell me the probability of Joe returning his books at 33.1 days and the probability at 33.9 days, and 33.9 days would have a higher probability.
In particular if there's a good way to do this with a readily-available python package that would be a huge help.
 A: There are several approaches you may use, depending on what assumptions are you willing to make and how much data do you have at hand. I present three possible methods with explantion and code. Note that I calculated only the "probability of returning book on day 25" and "the probability of returning book after noon on day 25", but you can choose arbitrary limits (even a second). 
Non Parametric Estimation


*

*In general, the function you wrote as example is the empirical distribution function $$ \hat{F}_n\left( x \right) = \frac{1}{n}\sum_{i=1}^{n} 1_{\left\lbrace X_i <x \right \rbrace} $$ Where $1_A = \begin{cases}1 & A \text{ holds} \\ 0 & \text{otherwise}\end{cases}$.  It convergence to the real distribution with rate of convergence of $O\left( \sqrt{n} \right)$, meaning that for large enough $n$, the errors $ \left| \hat{F}_n\left( x \right) - F\left( x \right) \right| $ have standard deviation of $ O\left(\frac{1}{\sqrt{n}}\right) $. More suprisingly, the convergence is uniform (holds for all values of $x$ together). However, this is all asymptotics. I'm not sure how it works on a small sample. In any case, it really easy to calculate this function:
def F(data, x)
    return np.mean(data < x)

# Probability of returning book on day 25
F(data, 25) - F(data, 24)

# Probability of returning book after noon on day 25
F(data, 25) - F(data, 24.5)


*Note however that $ \hat{F}_n\left( x \right) $ is still quite discrete, and it is possible to estimate the probabilities to a finer resolution using kernel density estimation. This approach approximates the derivative of $ \hat{F}_n\left( x \right) $ and is quite similar to producing histograms. In essense, you calculate $$ \hat{f}_h \left(x\right) = \frac{1}{n}\sum_{i=1}^{n} \frac{1_{\left\lbrace \left|x_i - x\right| < h \right\rbrace}}{2h}$$ This makes sense as soon as you note that $ 1_{\left\lbrace \left|x_i - x\right| < h \right\rbrace} = \hat{F}_n\left( x + h \right) - \hat{F}_n\left( x  - h\right)$. You have to choose the parameter $h$ which sort of acts as resolution. In practice, you can replace $1_{\left\lbrace \left|x_i - x\right| < h \right\rbrace}$ with other functions, all called kernel, which affect the smoothness of the resulting function. This method is more complex, and required more data to achieve the same performance as $ \hat{F}_n $. The actual rate actually depends on the underline distribution.
Although more complex, this estimator has an implementation in sklearn:
from sklearn.neighbors import KernelDensity
from scipy import integrate

clf = KernelDensity(bandwidth=0.5, kernel="tophat")
clf.fit(data)

# clf.score returns the log density
density = lambda x: np.exp(clf.score(x))

# Probability of returning book on day 25
probability, error = integrate.quad(density, 24, 25)

# Probability of returning book after noon on day 25
probability, error = integrate.quad(density, 24.5, 25)

Note that you might want to play around with the parameters, and that score returns the log of the probability.
Parametric Estimation


*Now we are ready to assume a model. It's a stronger assumption, but works better for smaller amounts of data (assuming you choose your model wisely). Two possibilities for such simple model are exponential distribution and gamma distribution. Both distribution are modeling continuous waiting times with small differences (exponential models the time until the first event, gamma until the $\alpha$'s th event). Our goal is to find the rate in which the event happens, i.e., how long until the costumer returned the book. (There might be better models, so don't get stuck on this one if it fails to predict properly).
Exponential distribution has density $$ f\left(x ; \lambda \right) = \lambda e^{-\lambda x}$$ where $x$ is non-negative. Note that given data points that were generated independently, $$ f\left(x_1,\ldots, x_n ; \lambda \right) = \prod_{i=1}^{n} \lambda e^{-\lambda x} =  \lambda^n e^{-\lambda \sum_{i=1}^{n} x}$$ A common estimator for $ \lambda $ using the data is the maximum likelihood estimator (MLE), that is, the paramter $\lambda$ that maximizes the probability density for the given data. In our case, $$ \hat{\lambda} = \frac{n}{\sum_{i=1}^{n} x_i }$$ After estimating the rate, we can calculate the required probability using the density function. scipy actually provides it for us:
from scipy.stats import expon

# Note that expon uses scale = 1 / lambda
scale = mean(data)

# Probability of returning book on day 25
expon.cdf(25) - expon.cdf(24)

# Probability of returning book after noon on day 25
expon.cdf(25) - expon.cdf(24.5)

Final Remark
There is no one size fit all. After estimating, you should test to see if the estimation makes sense. A few advised steps:


*

*Plot the histogram to know what to expect

*Plot the estimated density vs. the density of exponential with the estimated parameter

*Check the average error between points and their estimation. Preferably, do it on dataset that wasn't used for the estimation task (test set)

*If high resolutions are important, it might be advisable to convert the data to seconds.

