Solution (source code)

= Solution

The inverse <logit link> is $\mu_i(\beta)=e^{\beta^Tx_i}/(1+e^{\beta^Tx_i})$. <independence> therefore gives the <Likelihood function> and <log-likelihood>
$$
\boxed{L(\beta;y)=\prod_i\mu_i^{y_i}(1-\mu_i)^{1-y_i},\qquad
\ell(\beta)=\sum_i\left[y_i\beta^Tx_i-\log(1+e^{\beta^Tx_i})\right].}
$$
Differentiation yields the score $\nabla_\beta\ell=\sum_i x_i(y_i-\mu_i)$. The matrix whose rows are $x_i^T$ has dimension $n\times p$, not the PDF's printed $p\times n$. With this necessary dimension correction, any finite maximum <Likelihood function> estimate satisfies
$$
\boxed{X^Ty=X^T\widehat\mu,\qquad
\widehat\mu_i=\frac{e^{\widehat\beta^Tx_i}}{1+e^{\widehat\beta^Tx_i}}.}
$$
The <Hessian matrix> is $-X^TWX$, with $W_{ii}=\mu_i(1-\mu_i)$, so the <log-likelihood> is concave. These score equations apply to a finite interior estimate; complete or quasi-complete separation can instead make its supremum occur only as coefficients diverge.

The <deviance> is twice the difference between maximized saturated and fitted <log-likelihoods>. A saturated <Bernoulli distribution> has fitted means $y_i\in\{0,1\}$ and <log-likelihood> zero, interpreted by limits. Hence the <binomial deviance> here is
$$
\boxed{D(y;\widehat\mu)=2\sum_i\left[y_i\log\frac{y_i}{\widehat\mu_i}+(1-y_i)\log\frac{1-y_i}{1-\widehat\mu_i}\right]
=-2\sum_i\left[y_i\log\widehat\mu_i+(1-y_i)\log(1-\widehat\mu_i)\right].}
$$
Zero-multiplier logarithms are interpreted as zero.