Let $X_1\sim\text{Bernoulli}(p)$ and $X_2\sim\text{Bernoulli}(3p)$ be independent Bernoulli random variables where $p\in[0,1/3]$. Derive the MLE of $p$.
We have that
$$L(p\mid \vec{x})=p^{x_1}(1-p)^{1-x_1}(3p)^{x_2}(1-3p)^{1-x_2}$$
Upon taking the natural log of both sides we get
$$\mathscr{L}(p\mid\vec{x})=x_1\text{log}(p)+(1-x_1)\text{log}(1-p)+x_2\text{log}(3p)+(1-x_2)\text{log}(1-3p)$$
Then
$$\begin{align*} \frac{\partial\mathscr{L}}{\partial p} &=\frac{x_1}{p}-\frac{1-x_1}{1-p}+\frac{x_2}{p}-\frac{3(1-x_2)}{1-3p}\\\\ &=\frac{6p^2-3x_1p-px_2-4p+x_1+x_2}{p(1-p)(1-3p)} \end{align*}$$
which equals zero when
$$p=\frac{3x_1+x_2+4\pm\sqrt{9x_1^2+6x_1x_2+x_2^2-16x_2+16}}{12}$$
From here it's clear that
$$\frac{3x_1+x_2+4+\sqrt{9x_1^2+6x_1x_2+x_2^2-16x_2+16}}{12}\geq\frac{8}{12}\gt\frac{1}{3}$$
and
$$0\leq \frac{3x_1+x_2+4-\sqrt{9x_1^2+6x_1x_2+x_2^2-16x_2+16}}{12}\leq\frac{1}{3}$$
so
$$\hat{p}=\frac{3x_1+x_2+4-\sqrt{9x_1^2+6x_1x_2+x_2^2-16x_2+16}}{12}$$
However, a problem arises when trying to show that this is a global maximum
We have from software that
$$\begin{align*} \frac{\partial^2\mathscr{L}}{\partial p^2} &=\frac{-18p^4+24p^3+18x_1p^3+6p^3x_2-10p^2-21x_1p^2-13p^2x_2+8x_1p+8px_2-x_1-x_2}{p^2\left(-p+1\right)^2\left(-3p+1\right)^2} \end{align*}$$
Attempting to get a negative upper bound on the numerator to show that this is less than zero, we have that
$$-18p^4+24p^3+18x_1p^3+6p^3x_2-10p^2-21x_1p^2-13p^2x_2+8x_1p+8px_2-x_1-x_2$$
equals
$$-18p^4+p^3(24+18x_1+6x_2)+p^2(-10-21x_1-13x_2)+p(8x_1+8x_2)-x_1-x_2$$
which is less than or equal to
$$-18p^4+48p^3-10p^2+16p\leq\frac{48}{27}+\frac{16}{3}$$
so I fail to get a negative upper bound.
How can I show that what I obtained is a global maximum?
