Why can log-linear models for binary data be trained using Poisson regression 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)?
 A: I'm not sure I follow your examples, but I think you are asking about equivalence between (appropriately reparameterized) poisson and multinomial models, of which the bernoulli version is a special case.  The equivalence you might be thinking where some variables are treated as explanatory and therefore conditioned on in the multinomial / logit analysis but treated as random and therefore modelled in the poisson / log-linear analysis.  The equivalence works because of the shared likelihood kernel and the fact that parameters in the log linear model that fix aspects of the data that are conditioned on in the other analysis factor out.  I'd recommend reading Rodriguez's excellent notes, in particular section 6.2.5 on the equivalent log linear model or Agresti 2002 section 8.5 for a slower version.
If you want to see more general math, then google 'Multinomial Poisson Transformation' to get papers by Lang and Baker that I don't have to hand.  This is one of the few cases where having incidental parameters that increase with the data set size does not make ML estimates of the others inconsistent.
