Roots of log binomial likelihood equation

Question

I have been wondering if there is a way to approximate the roots of the log binomial likelihood equation. To be clear the equation is

$$a \cdot \log\left(x\right) + b \cdot \log\left(1 - x\right) = t$$

where $a$ and $b$ are positive integers, and $x$ is between 0 and 1. When $t < a \cdot \log\left(\frac{a}{a+b}\right) + b \cdot \log\left(\frac{b}{a+b}\right)$, I know there are two real roots between 0 and 1.

I tried a few symbolic solvers such as SymPy, SageMath, Mathematica, Matlab's solver but none of them give me any answers. I can numerically solve it using either Newton's method, secant method or Halley’s method, but there are a few problems. These methods are slow. I need to solve this for many different $a$, $b$, and $t$ values. Another problem is when t is very small, the roots get very closer to the end points 0 and 1, then the algorithms fail to converge.

I believe partly because I am not familiar with the field, and partly because English is not my first language, I couldn't find any literature on this topic. I would appreciate any references to the related literature very much.

Answer 1

I prefer to add another answer, since it is based on a totally different approach.

One of the problem when working with the function $$f(x)=a \, \log\left(x\right) + b \, \log\left(1 - x\right) - t$$ is that $x$ is a limited and bounded interval.

To unbound the variable, consider the change of variable $$x=\frac 1{1+e^{-y}}$$ which allows to vary between $-\infty$ and $+\infty$.

So, replacing and simplifying, now we have to find the zero's of function $$g(y)=a y-(a+b) \log \left(1+e^y\right)-t$$ Its maximum is at $y_*=\log \left(\frac{a}{b}\right)$ and a series expansion around this point gives $$g(y)=k-\frac{a b}{2 (a+b)}\left(y-y_*\right)^2-t+O\left(\left(y-y_*\right)^3\right)$$ $$k=\left(a\log(a)+b\log(b)-(a+b) \log(a+b)\right)$$ and then the approximate solutions $$y_\pm=y_*\pm \sqrt{\frac{2(a+b) (k-t)}{a b}} \implies x_\pm=\frac 1 {1+e^{y_\mp}}$$

Applied to the previous case $(a=1, b=2, t=-5)$, this gives $$y_-=-\log (2)-\sqrt{15-3 \log \left(\frac{27}{4}\right)}\implies x_-=0.023247$$ $$y_+=-\log (2)+\sqrt{15-3 \log \left(\frac{27}{4}\right)}\implies x_+=0.913073$$

Based on the $y$'s, Newton method converge very fast

$$\left( \begin{array}{cc} n & y_n \\ 0 & -3.738040 \\ 1 & -5.018754 \\ 2 & -4.979452 \\ 3 & -4.979437 \end{array} \right)$$

$$\left( \begin{array}{cc} n & y_n \\ 0 & 2.351746 \\ 1 & 2.365367 \\ 2 & 2.365354 \end{array} \right)$$

Update

I think that we could combine the two answers to take advantage of the quality of the estimates of $x$ given in the first answer and the nice conditioning of function $g(y)$ defined here.

For each of the root, the path would be $$x_0 \implies y_0=-\log \left(\frac{1}{x_0}-1\right)\implies y_1=y_0-\frac{g(y_0)}{g'(y_0)}\implies x_1=\frac{1}{1+e^{-y_1}}$$ where $$g'(y)=\frac{a+b}{1+e^y}-b$$

Applied to the worked example, this would give as solutions $$9.14146964097\times 10^{-1} \qquad \text{and} \qquad 6.83095206594\times 10^{-3}$$ to be compared to the exact $$9.14146937658\times 10^{-1} \qquad \text{and} \qquad 6.83095206610\times 10^{-3}$$

Answer 2

Consider that you look for the zero of function $$f(x)=a \, \log\left(x\right) + b \, \log\left(1 - x\right) - t$$ If $a<b$ there is a small root close to $0$ and a large root. Consider the large one and expand using Taylor $$f(x)=b\,\log(1-x)-a(1-x)+O\left((1-x)^2\right)-t$$ and then an approximate solution for the large root $$\color{red}{x_2 \sim 1+\frac b a W\left(-\frac{a }{b}e^{\frac{t}{b}}\right)}$$ where $W(.)$ is Lambert function.

For the small root $$f(x)=a \log (x)-b x+O\left(x^2\right)-t$$ $$\color{red}{x_1 \sim -\frac a b W\left(-\frac{b }{a}e^{\frac{t}{a}}\right)}$$

Let us try with $a=1$, $b=2$ and $t=-5$. This gives $x_2=0.914322$ and $x_1=0.00683063$ while the exact solutions are $x_2=0.914147$ and $x_1=0.00683095$.

With these estimates, Newton method should work like a charm. For the worked example, Newton iterates are.

$$\left( \begin{array}{cc} n & x_n \\ 0 & 0.9143221615 \\ 1 & 0.9141471262 \\ 2 & 0.9141469377 \end{array} \right)$$

$$\left( \begin{array}{cc} n & x_n \\ 0 & 0.006830627434 \\ 1 & 0.006830952058 \\ 2 & 0.006830952066 \end{array} \right)$$