Learning poisson parameters from aggregated data
================================================

In this note, we will study the problem of learning parameters of poisson distributions from *aggregated* data. In other words, the data does not directly reflect realizations of these distributions, but rather contains aggregates (i.e. sums) of such values. 

Consider the following scenario as a concrete example: A city consists of five areas: **A, B, C, D, E**. We are interested in the distribution of morning trips. We know that there exist three types of travelers:

1. *Workers*: going from home to work
2. *Shoppers*: going from home to a shopping area
3. *Kids*: going from home to school

Assume we only have the data for trips that originate from area **A**. Our data only contains the number of trips originating from **A** in the morning. We can not say anything about the purpose of trips by merely relying on the data. 

Know assume that we have the following pieces of background knowledge:

1. The trips of each group follows a poisson distribution
2. For each travel type, the rates of poisson distributions for different destinations are equal. 
3. We know the types of trips for each destination (They are presented in the image below).

**_The second assumption requires us to adjust for the population of areas when applying the following approach. More on this in another note._**


```
from IPython.display import SVG
SVG(filename='areas.svg')
```




    

    



Our problem consists of finding three rates: $\lambda_w$ for working trips, $\lambda_s$ for shopping trips, and $\lambda_k$ for school trips. 

## Reminder: Poisson density function
Remember that the poisson distribution has the following density function:
\begin{align}
Pr(X=k) = \frac{\lambda^k e^{-\lambda}}{k!}
\end{align}

#### example
This equation implies that under the assumption $\lambda_w = 10$, the probability of having 7 trips from **A** to **E** in the morning is $$\frac{{10}^{7} * e^{-10}}{7!} = 0.09$$

The value 0.09 is the *likelihood* of observation $X_{A \to E}=7$. 

In the next section we'll see an interesting property that will help us in our analysis.


## Sum of independent poisson variables
*For independent random variables $X_1,\ldots,X_n$ which follow poisson distributions with parameters $\lambda_1,\ldots,\lambda_n$, The variable $Y=\sum_1^n X_i$ follows a poisson distribution with parameter $\lambda = \sum_1^n \lambda_i$*

#### example
According to this property, under the assumption $\lambda_w = 10, \lambda_s = 3$, the trips from **A** to **B** have a poisson distribution with parameter $\lambda = 10 + 3 = 13$, and hence the probability of observing 15 trips is:

\begin{equation}
\frac{{13}^{15}*e^{-13}}{15!} = 0.088
\end{equation}

The value 0.088 is the *likelihood* of the observation $X_{A \to B} = 15$.

# Learning parameters from data

Now assume that we have the aggregated data for the morning trips originating from area **A** for a number of days:

|Destination | Day1 | Day2 | Day3 |
|---|:---:|:---:|:---:|
|*B*|99|92|78|
|*C*|7|6|10|
|*D*|119|94|98|
|*E*|89|85|81|

We have three unknowns: $\lambda_w, \lambda_s, \lambda_k$. Each of them can take a value in the range $[0, \infty)$. We want to pick values that are most consistent with the data. So we select the parameters that maximize the likelihood of data. 

To compute the likelihood, we first compute the probability of each observation (for each destination-day). Assuming that *trips to different destination* and *trips of different days* are independent, the likelihood will be the product of these probabilities. 

So the probability of the observation for destination **B** in day 1 is:
\begin{equation}
\frac{e^{-(\lambda_w + \lambda_s)} * {(\lambda_w + \lambda_s)}^{99}}{99!}
\end{equation}

Remember that the trips from **A** to **B** are taken by *shoppers* and *workers*:


```
from IPython.display import SVG
SVG(filename='only_b.svg')
```




    

    



As another example, consider the probability of observation for destination **C** in day 2:
\begin{equation}
\frac{e^{-\lambda_s} * {\lambda_s}^{6}}{6!}
\end{equation}

Remember that the trips from **A** to **C** are taken only by *shoppers*. 


```
from IPython.display import SVG
SVG(filename='only_c.svg')
```




    

    



By the same principle, the likelihood for the whole data ($\mathcal{D}$) is computed as:

\begin{align}
L(\mathcal{D}|\lambda_w, \lambda_s, \lambda_k) = 
&\frac{e^{-3*(\lambda_w + \lambda_s)} * {(\lambda_w + \lambda_s)}^{99+92+87}}{99!*92!*87!} + \\
&\frac{e^{-3*\lambda_s} * {\lambda_s}^{7+6+10}}{7!*6!*10!} + \\
&\frac{e^{-3*(\lambda_w + \lambda_k)} * {(\lambda_w + \lambda_k)}^{119+94+98}}{119!*94!*98!} + \\
&\frac{e^{-3*\lambda_w} * {\lambda_w}^{89+85+81}}{89!*85!*81!} 
\end{align}

Our problem is to find the values $\lambda_w, \lambda_s, \lambda_k$ such that the likelihood is maximized:
\begin{equation}
\lambda^* = \underset{\lambda}{\operatorname{argmax}} L(\mathcal{D}|\lambda)
\end{equation}

To facilitate the computations, we can formulate the problem as that of maximizing the *logarithm of likelihood*. 

\begin{align}
LL(\mathcal{D}|\lambda_w, \lambda_s, \lambda_k) \equiv &\log(L(\mathcal{D}|\lambda_w, \lambda_s, \lambda_k)) = \\
&-9 \lambda_w - 6 \lambda_s - 3 \lambda_k + \\
&278 \log (\lambda_w + \lambda_s) +
23 \log (\lambda_s) + \\
&311 \log (\lambda_w + \lambda_k) +
255 \log (\lambda_w) + \mathcal{C}
\end{align}

Since the constant $\mathcal{C}$ (emerged from the factorial terms in denominators) does not affect the optimizaiton problem, we safely drop it.

## Parameter learning as an optimization problem

We formulated our problem as an instance of continuous optimization problem:

\begin{align}
\text{maximize} &LL(\mathcal{D}|\lambda) \\
&\lambda \geq 0
\end{align}

For our example dataset, the actual optimization problem will be:

\begin{align}
\text{maximize} & -9 \lambda_w - 6 \lambda_s - 3 \lambda_k + \\
&278 \log (\lambda_w + \lambda_s) +
23 \log (\lambda_s) + \\
&311 \log (\lambda_w + \lambda_k) +
255 \log (\lambda_w) \\
&\lambda_w, \lambda_s, \lambda_k \geq 0
\end{align}

Note that increasing the number of days only changes the constant coefficients in the objective function of this optimization problem. For example, assume that the number of trips for 10 days are reported in the following table:

|Destination | Day1 | Day2 | Day3 |Day4 |Day5 |Day6 |Day7 |Day8 |Day9 |Day10 |
|---|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|
|*B*|99|92|78|87 | 78 | 76 | 74 | 107 | 75 | 97 |
|*C*|7|6|10|13 | 11 | 11 | 13 | 8 | 17 | 10 |
|*D*|119|94|98|126 | 118 | 117 | 125 | 127 | 108 | 120 |
|*E*|89|85|81|74 | 79 | 68 | 78 | 75 | 80 | 68 |

Then the new optimization problem will be:
    
\begin{align}
\text{maximize} & -30 \lambda_w - 20 \lambda_s - 10 \lambda_k + \\
&846 \log (\lambda_w + \lambda_s) +
100 \log (\lambda_s) + \\
&1094 \log (\lambda_w + \lambda_k) +
831 \log (\lambda_w) \\
&\lambda_w, \lambda_s, \lambda_k \geq 0
\end{align}    

Generalizing to the data of $n$ days is straightforward: Denote the total number of trips to **B**, **C**, **D**, and **E** by $N_b$, $N_c$, $N_d$, and $N_e$. Then the optimization problem will be:

\begin{align}
\text{maximize} & -3n \lambda_w - 2n \lambda_s - n \lambda_k + \\
&N_b \log (\lambda_w + \lambda_s) +
N_c \log (\lambda_s) + \\
&N_d \log (\lambda_w + \lambda_k) +
N_e \log (\lambda_w) \\
&\lambda_w, \lambda_s, \lambda_k \geq 0
\end{align} 

Now it is time to solve this problem.

# Solving the optimization problem

In this section we try to solve the optimization problem using the 10-day data. The actual values used for sampling the data were $\lambda_w = 80$, $\lambda_s = 10$, and $\lambda_k = 30$.


```
import numpy as np
def ll_data(n, N):
    def loglikelihood(lambdas):
        lw = lambdas[0]
        ls = lambdas[1]
        lk = lambdas[2]
        res = -3*n*lw - 2*n*ls - n*lk
        res += N[0] * np.log(lw+ls) + N[1] * np.log(ls)
        res += N[2] * np.log(lw+lk) + N[3] * np.log(lw)
        return -res #reversing the sign to maximize
    return loglikelihood
```


```
from scipy.optimize import minimize
def max_loglikelihood(N, n):
    first_guess = [1, 1, 1]
    bounds = [(0, None)]*3 #lambdas should be non-negative
    res = minimize(ll_data(n, N), x0, bounds=bounds, method='L-BFGS-B')
    print 'lw:%.2f ls:%.2f lk:%.2f  --> %.2f'%(res.x[0], res.x[1], res.x[2],-res.fun)

n = 10
N = [846, 100, 1094, 831]
max_loglikelihood(N, n)
```

    lw:79.31 ls:9.54 lk:30.09  --> 9921.12


Solving the optimization problem gives us the values $\lambda_w = 79.31$, $\lambda_s = 9.54$, and $\lambda_k = 30.09$. 

We can also solve the problem with the data for 100 days:


```
n = 100
N = [8979, 992, 11185, 7926]
max_loglikelihood(N, n)
```

    lw:79.53 ls:9.95 lk:32.32  --> 100996.61


The value for $\lambda_k$ (32.32) has diverged from the actual value (30). However, this is not due to the quality of solutions to the optimization problem. Replacing the actual values in the log-likelihood function evaluates to a lower result:


```
print '%.3f' % -ll_data(n, N)([80, 10, 30])
```

    100994.771


We can conclude that the information available in the 100-day data does not allow for a perfect prediction of $\lambda_k$.

The partitioning issue
========================

Let us consider another scenario. This scenario is equivalent to the previous one in all aspects, except for the composition of trips: The trips to **B** and **D** are taken by *workers* and *shoppers*, and the trips to **C** and **E** are taken by *kids*. The figure below shows the composition of trips to each area.


```
from IPython.display import SVG
SVG(filename='areas2.svg')
```




    

    



Denoting the total number of trips to **B**, **C**, **D**, and **E** by $N_b$, $N_c$, $N_d$, and $N_e$, the optimization problem for this scenario will be:

#### (version 1)
\begin{align}
\text{maximize} & -2n (\lambda_w + \lambda_s) - 2n \lambda_k + \\
&N_b \log (\lambda_w + \lambda_s) +
N_c \log (\lambda_k) + \\
&N_d \log (\lambda_w + \lambda_s) +
N_e \log (\lambda_k) \\
&\lambda_w, \lambda_s, \lambda_k \geq 0
\end{align} 

As you see, the variables $\lambda_w$ and $\lambda_s$ always appear in the term $\lambda_w + \lambda_s$. One can substitute $\lambda_{ws}$ for $\lambda_w + \lambda_s$ and re-write the problem as:

#### (version 2)
\begin{align}
\text{maximize} & -2n (\lambda_{ws}) - 2n \lambda_k + \\
&N_b \log (\lambda_{ws}) +
N_c \log (\lambda_k) + \\
&N_d \log (\lambda_{ws}) +
N_e \log (\lambda_k) \\
&\lambda_{ws}, \lambda_k \geq 0
\end{align} 

The problem is that, once $\lambda_{ws}$ is obtained, each partitioning of $\lambda_{ws}$ into $\lambda_w$ and $\lambda_s = 1 - \lambda_w$ gives the same likelihood. This means that the values $\lambda_w$ and $\lambda_s$ obtained by solving the first version of the optimization problem are arbitrary conditioned on their sum. 

Let us check this by solving the first version:


```
def ll_data(n, N):
    def loglikelihood(lambdas):
        lw = lambdas[0]
        ls = lambdas[1]
        lk = lambdas[2]
        res = -2*n*lw - 2*n*ls - 2*n*lk
        res += N[0] * np.log(lw+ls) + N[1] * np.log(lk)
        res += N[2] * np.log(lw+ls) + N[3] * np.log(lk)
        return -res #reversing the sign to maximize
    return loglikelihood
```


```
n = 100
N = [8929, 2904, 9031, 2973]
max_loglikelihood(N, n)
```

    lw:44.90 ls:44.90 lk:29.39  --> 76806.73


Solving the optimization problem gives $\lambda_w = 44.9$, $\lambda_s = 44.9$, and $\lambda_k = 29.39$. However, all we know from the data about $\lambda_w$ and $\lambda_s$ is that $\lambda_w + \lambda_s = 89.8$. Any partitioning of this sum into non-negative pairs of $\lambda_w$ and $\lambda_s$ will give the same log-likelihood:


```
print '%.3f' % -ll_data(n, N)([44.9, 44.9, 29.39])
print '%.3f' % -ll_data(n, N)([89.8, 0, 29.39])
print '%.3f' % -ll_data(n, N)([0, 89.8, 29.39])
```

    76806.733
    76806.733
    76806.733


## Is this really a problem?
Tias made this important observation:
> The partionting problem does not matter in our application, as long as we are going to sample according to the same trip decomposition. 

It is true! To sample from the distribution specified in the current scenario, we only need to know $\lambda_k$ and $\lambda_s + \lambda_w$. Hence the partitioning issue does not pose a practical limitation. 

349d23132817a324d8c96ae333bd900c307424d3
ipynb
Jupyter Notebook
Hanoi/Learning_Poisson_Rates.ipynb
Behrouz-Babaki/notebooks
9e84263879ba17e8581b8bdee60aaf3006fa4a73
Unlicense
Hanoi/Learning_Poisson_Rates.ipynb
Behrouz-Babaki/notebooks
9e84263879ba17e8581b8bdee60aaf3006fa4a73
Unlicense
Hanoi/Learning_Poisson_Rates.ipynb
Behrouz-Babaki/notebooks
9e84263879ba17e8581b8bdee60aaf3006fa4a73
Unlicense
Qwen/Qwen-72B
1. YES
2. YES