Any arbitrary distribution of a multivariate binary variable can be represented as a log-linear model. That is, for $X = (X_1,\dots,X_d)$ a $d$ dimensional binary rv, the distribution can be written as $\log(\text{p}(x)) = \sum\limits_{A\subseteq V} \beta_A (\prod\limits_{i \in A} x_i)$ where $V={1,\dots,d}$. E.g. for $X = (X_1,X_2)$, $\log(\text{p}(x)) = \beta_\emptyset + \beta_1 x_1 + \beta_2 x_2 + \beta_{1,2} x_1 x_2$. Note, $\beta_\emptyset$ is a normalizing term.
I've read that this model can be trained using Poisson regression. For example, in R
m_out = glm(count ~ x1 + x2 + x1*x2, family="poisson")
for the 2d case above.
I'm unclear about why this is so: why does modeling the counts as being distributed as Poisson with mean = $\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_{1,2} x_1 x_2$ work?
In particular, let $\beta(x)=\log(\text{p}(x))-\beta_\emptyset$ if we write out the likelihood for our betas we have that: $L(\beta) = \prod\limits_{i=1}^m e^{\beta(x_i)+\beta_\emptyset}=\exp( \sum\limits_{i=1}^m (\beta(x_i)+\beta_\emptyset))$ so the log likelihood is : $l(\beta)=\sum\limits_{i=1}^m (\beta(x_i)+\beta_\emptyset) = \sum\limits_{A\subseteq V} y_A(\beta(x_A)+\beta_\emptyset) = m\beta_\emptyset+\sum\limits_{A\subseteq V} y_A\beta(x_A)$ where $x_A$ is the vector $x$ such that $x_i=1$ if $i \in A$ $x_i=0$ otherwise and $y_A$ is the empirical count of how many observations had the pattern $x_A$ and $\beta_\emptyset = -\log(\sum\limits_{A\subseteq V} e^{\beta(x_A)})$ is the normalizer.
However, according to Wikipedia what we are doing with glm is modeling our counts $y_A \sim \text{Pois}(\beta(x_A)+\beta_0)$ where $\beta_0$ is an intercept term. Ignoring constants, this will lead to a log likelihood of: $\sum\limits_{A\subseteq V} (y_A(\beta(x_A)+\beta_0)-e^{\beta(x_A)+\beta_0})$.
Is there a reason why maximizing both likelihoods will give the same coefficients (at least up to constants)?
