# How to estimate discrete probability distribution from a dataset of pairwise frequencies?

My problem is best explained with an example..

Suppose every day I count the cars driving down my street, but because I'm a slow counter, I only count two different types of cars per day. On the first day, I count 30 trucks and 65 cars. On the second day, I count 25 trucks and 12 motorcycles. On the third day I count 19 motorcycles and 72 cars. On the fourth day I count 11 vans and 14 motorcycles. And on the fifth day I count 12 vans and 44 trucks.

While the total volume of traffic changes day to day, I expect the relative proportion of vehicle types to remain near constant.

How do I use this information to develop a discrete probability distribution for the vehicles on my road? Put another way, how can I determine the probability that the next car I observe will be type X?

I suspect there's a nice conditional probability tree solution to this, but I can't quite put my finger on it.

## UPDATE

Upon more thought, I think this may be a good candidate for maximum likelihood estimation. If I assume $$P(\text{car})= p_1$$, $$P(\text{truck})= p_2$$, $$P(\text{motorcycle})= p_3$$, and $$P(\text{van})= p_4$$, then I can calculate the likelihood of my observed pairs of data given some values for $$p_1$$, $$p_2$$, $$p_3$$, and $$p_4$$. Now I just need to define and optimize my likelihood function.. Hmm..

• Each day you make a Binomial observation. Go on from there...
– whuber
Jul 1, 2021 at 14:24
• @whuber, thanks - that's the approach I've been working on. Still trying to get the math to sing to me..
– Ben
Jul 1, 2021 at 14:32

Since you have a good idea, let's formalize it.

Begin with notation. We need to suppose you know all the possible types of vehicle. Let's index them with whole numbers $$1,2,\ldots,k.$$ You assume a constant chance $$p_i$$ that vehicle $$i$$ will appear at any time, so (obviously) these chances are non-negative and sum to unity.

When on a particular day $$t$$ you count only vehicles of types $$i$$ and $$j\ne i,$$ let $$x_{ti}$$ and $$x_{tj}$$ be those counts. I will assume you decide, independently of any of the preceding data or other related information, which vehicles to count. (Thus, for instance, you don't just count the first two types of vehicles you encounter: that would provide more information.)

If the appearances of the vehicles are independent, then the chance of observing $$x_{ti}$$ and $$x_{tj}$$ must be proportional to the powers of the probabilities,

$$\Pr(x_{ti}, x_{tj}) \,\propto\, (1-p_i-p_j)^{-(x_{ti}+x_{tj})}\,p_i^{x_{ti}}\,p_j^{x_{tj}}.$$

Proceeding from day to day, again assuming independence, these probabilities continue to multiply. Thus, for each pair $$\{i,j\},$$ we may once and for all collect our counts $$x_{ti}$$ and $$x_{tj}$$ over all the days $$t$$ when we observed both these vehicle types. Let $$X_{ij}$$ be the sum of all such $$x_{tj}$$ (and, therefore, $$X_{ji}$$ is the sum of all such $$x_{ti}$$). Set $$X_{ii}=0$$ for all $$i.$$

For example, indexing Truck, Car, Motorcycle, and Van with the numbers $$1$$ through $$4,$$ respectively, the data in the question yield the matrix

$$X = \begin{array}{l|rrrr} & \text{Truck} & \text{Car} & \text{Cycle} & \text{Van} \\ \hline\text{Truck} & 0 & 65 & 12 & 12 \\ \text{Car} & 30 & 0 & 19 & 0 \\ \text{Cycle} & 25 & 72 & 0 & 11 \\ \text{Van} & 44 & 0 & 14 & 0 \end{array}$$

For instance, on days when you were observing trucks ($$1$$) and cars ($$2$$), you saw a total of $$X_{12}=65$$ cars and $$X_{21} = 30$$ trucks.

Let $$\mathbf{p} = (p_i)$$ be the vector of probabilities and let $$\mathbf{x} = (x_j) = \left(\sum_{i=1}^k X_{ij}\right)$$ represent the total counts of each type of vehicle (the column sums of $$X$$).

Upon taking logarithms of the total probability, we obtain a log likelihood

$$\Lambda = C + 2 \mathbf{x}\cdot \log \mathbf{p} - \sum_{i\ne j} (X_{ij}+X_{ji})\log(p_i+p_j)$$

where $$C$$ is what became of all the implicit constants of proportionality. Maximize this to obtain the maximum likelihood estimate of $$\mathbf p.$$ I obtain the value $$\hat {\mathbb{p}} = (0.26, 0.53, 0.13, 0.081).$$

Generally, with more than two types of vehicles, the equations for the critical points of the likelihood are polynomials of degree two or higher: you won't obtain closed formulas. This requires numerical optimization.

• Let $N_t$ be the number of vehicles on day $t$, including non counted types. Should not your pair probabilities be proportional to $Pr(x_{ti}, x_{tj}) \propto (1 - p_i - p_j)^{N_t - (x_{ti} + x_{tj})} p_i^{x_{ti}} p_j^{x_{tj}}$ ? The value of $N_t$ might not be known, but why can you drop the term on the equation? Jul 2, 2021 at 16:46
• @Lucas That's a great question. If you could propose what to do about the unknown $N_t,$ it could lead to a good answer. My short answer is, when you go out on day $t$ to count only vehicles $i$ and $j,$ why do any of the other vehicles matter?
– whuber
Jul 2, 2021 at 18:01
• I have trouble maximizing the log-likelihood function. If I understand it correctly, it returns 4 numbers in this example, right? How can you make sure that the four $p_i$s sum to 1? Jul 2, 2021 at 18:26
• @COOL I use the $\log p_i$ for the parameters, noting that it's unnecessary to enforce the sum-to-unity restriction of the $p_i$ because the likelihood is a homogeneous function of them. This creates a function with a one-dimensional set of maxima, but that doesn't cause problems in practice with most numerical algorithms. If by chance there were a problem, you could resort to constrained optimization.
– whuber
Jul 2, 2021 at 18:51
• @COOL Given the array X, here is my code to find an estimate q.hat:p.start <- rep(0, ncol(X)); X.c <- colSums(X); X2 <- X + t(X); Lambda <- function(p.log) { p <- exp(p.log); P.log <- log(outer(p, p, "+")); sum(X2 * P.log) - 2 * sum(X.c * p.log) }; fit <- nlm(Lambda, p.start, gradtol=1e-12, hessian=TRUE); exp(fitestimate) -> q.hat; q.hat <- q.hat / sum(q.hat) – whuber Jul 3, 2021 at 14:02 I have a proposition that is near to what you thought, but with an additional step of working with an exact marginal distribution. This will be a long answer. Skip the mathematics if you are not inclined. The conclusions about the computations is explained after they are made. The model Let $$(N_t, X_{ti}, X_{tj})$$ be the total number of vehicles, the number of vehicles of type $$i$$ and the number of vehicles of type $$j$$, $$i\neq j$$ , on the $$t$$-th day. The random variable $$N_t$$ is latent, meaning it is not observed directly. The idea is to propose a reasonable joint distribution for this triplet, and then compute the marginal of distribution of $$(X_{ti}, X_{tj})$$, which will then be used to compute the proportions. For this idea to work, we have to propose first a distribution for $$N_t$$. Since we are working with a countable variable, I will assume that $$N_t \sim \mbox{Poisson}(\lambda)$$. For the conditional distribution of $$(X_{ti}, X_{tj})$$, we propose $$\mathbb{P}(X_{ti} = x_i, X_{tj} = x_j|N_t = n) = \frac{n!}{x_i!x_j!(n-x_i-x_j)!}p_i^{x_i}p_j^{x_j}(1 - p_i - p_j)^{n-x_i-x_j} \quad.$$ That is, they are a trinomial distribution with $$N_t = n$$ observations, where the $$(1 - p_i - p_j)$$ term correspond to the probability of observing vehicles of type different than $$i$$ and $$j$$. The joint distribution can be written as $$\mathbb{P}(N_t = n, X_{ti} = x_i, X_{tj} = x_j) = \mathbb{P}(N_t = n) \mathbb{P}(X_{ti} = x_i, X_{tj} = x_j|N_t = n) \quad.$$ We have analytical an expression for the RHS. However, we are interest in the marginal $$\mathbb{P}(X_{ti} = x_i, X_{tj} = x_j) = \sum_{n=0}^\infty \mathbb{P}(N_t = n) \mathbb{P}(X_{ti} = x_i, X_{tj} = x_j|N_t = n) \quad.$$ Before evaluating the summation, notice that $$\mathbb{P}(X_{ti} = x_i, X_{tj} = x_j|N_t = n) = 0$$ when $$n < x_i + x_j$$. Now we evaluate the summation \begin{align} \mathbb{P}(X_{ti} = x_i, X_{tj} = x_j) &= \sum_{n=0}^\infty \frac{e^{-\lambda}\lambda^{n}}{n!} \frac{n!}{x_i!x_j!(n-x_i-x_j)!}p_i^{x_i}p_j^{x_j}(1 - p_i - p_j)^{n-x_i-x_j}\\ &= \frac{e^{-\lambda}p_i^{x_i}p_j^{x_j}}{x_i!x_j!}\sum_{n=x_i+x_j}^\infty \frac{\lambda^{n}}{(n-x_i-x_j)!}(1 - p_i - p_j)^{n-x_i-x_j}\\ &= \frac{e^{-\lambda}p_i^{x_i}p_j^{x_j}}{x_i!x_j!}\sum_{k = 0}^\infty \frac{\lambda^{k+x_i+x_j}}{k!}(1 - p_i - p_j)^{k}\\ &= \frac{e^{-\lambda}\lambda^{x_i+x_j}p_i^{x_i}p_j^{x_j}}{x_i!x_j!}\sum_{k = 0}^\infty \frac{(\lambda(1 - p_i - p_j))^{k}}{k!}\\ &= \frac{e^{-\lambda}\lambda^{x_i+x_j}p_i^{x_i}p_j^{x_j}}{x_i!x_j!}e^{\lambda(1-p_i-p_j)}\\ &= \frac{e^{-\lambda p_i}(\lambda p_i)^{x_i}}{x_i!}\frac{e^{-\lambda p_j}(\lambda p_j)^{x_j}}{x_j!}\\ \end{align} After these exausting computations, we have a great result: $$X_{ti}$$ and $$X_{tj}$$ are marginally independent. Moreover, their distributions is $$X_{ti} \sim \mbox{Poisson}(\lambda p_i)$$ for all $$i \, \in \, \{1,\ldots, k\}$$, where $$k$$ is the number of types of vehicles. With this, you can write the likelihood as a product of independent Poissons for each type of vehicle, where the number of observations for each will vary on how you choose the vehicle counting. Estimation Write $$\lambda_i = \lambda p_i$$. Let $$n_i$$ be total number of days you chose to observe the $$i$$-th type of vehicle, and let $$S_i$$ be the sum of the observations of that vehicle type. The maximum likelihood estimator for $$\lambda_i$$ is $$\hat{\lambda}_i = \frac{1}{n_i}S_i \quad,$$ That is, it is just the average of the observations. But we do not want to estimate $$\lambda_i$$, we want to estimate $$p_i$$. Well, we know that $$1 = \sum_{i=1}^k p_i \quad.$$ Multiplying both sides by $$\lambda$$, we have $$\lambda = \sum_{i=1}^k \lambda_i \quad.$$ By the invariance property of the MLE, we have $$\hat{\lambda} = \sum_{i=1}^k \hat{\lambda}_i \quad.$$ But $$p_i = \lambda_i/\lambda$$, hence $$\hat{p}_i = \frac{\hat{\lambda}_i}{\hat{\lambda}} \quad.$$ Therefore, we can estimate the proportions $$p_i$$, and also the expected number of vehicles $$\lambda$$! Your example To show that this approach might work, lets compute each parameter in your data example: \begin{align} &\hat{\lambda}_{truck} = (30+25+44)/3 = 33\\ &\hat{\lambda}_{car} = (65+72)/2 = 68.5\\ &\hat{\lambda}_{cycle} = (12+19+14)/3 = 15\\ &\hat{\lambda}_{van} = (12+11)/2 = 11.5 \end{align} For $$\lambda$$, we have $$\hat{\lambda} = \hat{\lambda}_{truck} + \hat{\lambda}_{car} + \hat{\lambda}_{cycle} + \hat{\lambda}_{van} = 128 \quad.$$ For the probabilities \begin{align} &\hat{p}_{truck} = \hat{\lambda}_{truck}/\hat{\lambda} = 0.2578\\ &\hat{p}_{car} = \hat{\lambda}_{car}/\hat{\lambda} = 0.5351\\ &\hat{p}_{cycle} = \hat{\lambda}_{cycle}/\hat{\lambda} = 0.1171\\ &\hat{p}_{van} = \hat{\lambda}_{van}/\hat{\lambda} = 0.0898 \end{align} The probabilities are strinkingly similar to those provided by @whuber. The difference is that it is very very easy to compute it, no optimization required. Final Analysis Here is a final overall comparison when considering this approach. Advantages: 1. Its very easy to compute the estimators, they are analytical; 2. You probably can perform hypothesis testing, if you wish; 3. You can estimate the total number of vehicle. Disadvantages: 1. We assumed that $$N_t$$ does not vary with the day, which might be false due to weekly or monthly seasonality; 2. I do not know how you could check if the Poisson distribution is adequate for $$N_t$$ with the data at hand. • Nice writeup, but the assumption thatN_t$does not vary with the day is false. The amount of cars on the road changes drastically from day to day depending on factors like weather, day of week, is it a holiday, etc. This is a fundamental part of the challenge. – Ben Jul 5, 2021 at 18:32 • I imagine it is a complicated process. However, you might still be able to apply the analysis above. For instance, if you wished to include day of week and raining, you could suppose$N_t = Poisson(\lambda_{d(t), r(t)})$, where$d(t)$would be the day of the week and$r(t)$a rain boolean. The computations are almost the same, the only difference is that you would have to group the variables for estimation. The only strong hypothesis remaining is independence. In optimizing a "pseudo likelihood" that is decomposed as a product of observations, I think you are implicitly assuming independence. Jul 5, 2021 at 18:49 With some help from @whuber and @Lucas Prates, I've put together this solution (which is pretty close to @whuber's solution). My general approach is to use maximum likelihood estimation in the following manner 1. Index the vehicle types in alphabetical order, car = 1, motorcycle = 2, truck = 3, van = 4. 2. Initialize arbitrary probabilities for $$p_1 = 0.25$$, $$p_2 = 0.25$$, $$p_3$$, and $$p_4 = 0.25$$ where $$p_i$$ gives the true probability of observing vehicle $$i$$ next. 3. Let $$p_{jk}$$ be the probability of observing vehicle $$j$$ next when we restrict ourselves to viewing only vehicles of class $$j$$ and $$k$$. Assuming independence, we get $$p_{jk} = \frac{p_j}{p_j + p_k}$$. 4. Let $$x_{jk}$$ be the number of $$j$$ vehicles we observed when we restricted ourselves to viewing only vehicles of class $$j$$ and $$k$$. 5. $$x_{jk}$$ is a binomial random variable. Thus, if we believe in our discrete probability distribution, the probability of observing $$x_{jk}$$ is $$\binom {n_{jk}}{x_{jk}}p_{jk}^{x_{jk}}(1-p_{jk})^{n_{jk}-x_{jk}}$$ 1. The probability of all our pair-wise observations will just be the product of their individual probabilities. So, our likelihood function looks like this $$L(p,x)=\prod_{jk} \binom {n_{jk}}{x_{jk}}\left( p_{jk}^{x_{jk}}(1-p_{jk})^{n_{jk}-x_{jk}} \right)$$ where $$jk$$ spans all pair-wise indexes that we observed. However, since our goal is merely to find $$\bf p$$ that maximizes the likelihood function, $$n_{jk}$$ and $$x_{jk}$$ terms are constants, so we can drop the $$\binom {n_{jk}}{x_{jk}}$$ terms. 1. Next, we take the log to get $$log(L)= \sum_{jk} x_{jk}log(p_{jk})+ (n_{jk}-x_{jk})log(1-p_{jk})$$ 1. Then we can calculate the partial derivative of $$log(L)$$ w.r.t. $$p_i$$ as $$\frac{\partial log(L)}{\partial p_i} = \sum_{i=j} \left( \frac{x_{jk}}{p_{jk}} - \frac{n_{jk} - x_{jk}}{1-p_{jk}} \right) \frac{p_k}{(p_j + p_k)^2} + \sum_{i=k} \left( \frac{x_{jk}}{p_{jk}} - \frac{n_{jk} - x_{jk}}{1-p_{jk}} \right) \frac{-p_j}{(p_j + p_k)^2}$$ 1. This looks nasty, but we know all these terms so we can plug and chug. I.e. we can use gradient descent to find the $$\bf p$$ that maximizes $$log(L)$$. I've implemented this in R... library(data.table) get_discrete_dist <- function(pairsDT, alpha = 0.0001, iters = 1000){ # Copy input data.table pairsDTCopy <- copy(pairsDT) # Insert total obs pairsDTCopy[, n_jk := qty1 + qty2] # initialize probabilities for every class vehicles <- sort(unique(c(pairsDTCopy$$vehicle1, pairsDTCopy$$vehicle2))) vehicles <- data.table(vehicle = vehicles) vehicles[, p := 1/.N] for(i in seq_len(iters)){ # Insert probabilities into pairsDTCopy pairsDTCopy[vehicles, p1 := i.p, on = c("vehicle1"="vehicle")] pairsDTCopy[vehicles, p2 := i.p, on = c("vehicle2"="vehicle")] # Calculate log likelihood pairsDTCopy[, p_jk := p1/(p1 + p2)] pairsDTCopy[, logl := qty1*log(p_jk) + (n_jk - qty1)*log(1 - p_jk)] logl <- sum(pairsDTCopy$logl)  # -1316.98

# Print
print(paste0("iter: ", i, ", log likelihood: ", logl))

# Calculate the gradient (partial log likelihood / partial p_i for every class)
g1 <- pairsDTCopy[, list(grad = sum((qty1/p_jk - (n_jk - qty1)/(1 - p_jk)) * (p2)/(p1 + p2)^2)), keyby = list(vehicle = vehicle1)]
g2 <- pairsDTCopy[, list(grad = sum((qty1/p_jk - (n_jk - qty1)/(1 - p_jk)) * (-p1)/(p1 + p2)^2)), keyby = list(vehicle = vehicle2)]

vehicles[, p := p + alpha * grad]
vehicles[, p := p/sum(p)]  # normalize
}

# Return
return(vehicles[])
}

pairsDT <- data.table(
vehicle1 = c("truck", "truck", "motorcycle", "van", "van"),
vehicle2 = c("car", "motorcycle", "car", "motorcycle", "truck"),
qty1 = c(30, 25, 19, 11, 12),
qty2 = c(65, 12, 72, 14, 44)
)
get_discrete_dist(pairsDT, alpha = 0.0001, iters = 100)

$$$$
`
• A nice build up from your original idea. However, you are including constant terms such as $log(1/2\pi)$ on the optimization function, this is not necessary. Moreover, why is ${n \choose k}p^{k}(1-p)^{n-k}$ hard to compute? If you are referring to the term ${n \choose k}$, just take the log and it will be a constant as well (the $n$ does not matter, only your parameters). Optimizing this log-likelihood seems easier than using the normal approximation. The approximation might be poor if one of your probabilities is near the extremes of $[0,1]$, which might be the case for your example. Jul 6, 2021 at 12:07