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

Question

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..

Answer 1

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.

Answer 2

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.

Answer 3

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}} $$

6. 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.

7. Next, we take the log to get

$$ log(L)= \sum_{jk} x_{jk}log(p_{jk})+ (n_{jk}-x_{jk})log(1-p_{jk}) $$

8. 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} $$

9. 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)]
    grads <- rbind(g1, g2)[, list(grad = sum(grad)), keyby = vehicle]
    
    # Update probabilities (gradient step)
    vehicles[grads, grad := i.grad, on = "vehicle"]
    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)

      vehicle          p       grad
1:        car 0.52790741  0.5766159
2: motorcycle 0.12749926 -0.2156412
3:      truck 0.26296780 -0.9608602
4:        van 0.08162554 -0.2951460
```