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

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

*

*Index the vehicle types in alphabetical order, car = 1, motorcycle = 2, truck = 3, van = 4.

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

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

*Let $x_{jk}$ be the number of $j$ vehicles we observed when we restricted ourselves to viewing only vehicles of class $j$ and $k$.

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


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


*Next, we take the log to get

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


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


*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
```

