Gradient and Hessian of a likelihood function where y is defined implicitly

Question

$\newcommand{\implicit}{\mathrm{implicit}}$

I'm attempting to optimize model parameters $\theta$ by maximizing the likelihood function $$ f(y) = \ln \Bigl(\frac {n!}{k!(n-k)!}y^k(1-y)^{n-k}\Bigl) $$ where $y$ must be calculated iteratively, as it is defined implicitly by $g(\theta)=0$, as follows: $$ g(\theta) = -1 + \frac{a\cdot \theta_2}{\ln \bigl(\frac{-y}{y - 1} \bigl) - \theta_1} + \frac{b\cdot \theta_4}{\ln \bigl(\frac{-y}{y - 1} \bigl) - \theta_3} $$

($n$, $k$, $a$, and $b$ are constants).

As far as I can tell, using partial implicit differentiation, the gradient (first partial derivative with respect to $\theta_i$) would be defined by $$ \frac{\partial \operatorname{f}}{\partial \operatorname{\theta_i}} = \frac{\frac{\partial \operatorname{f}}{\partial \operatorname{y}}\cdot \frac{-\partial \operatorname{g}}{\partial \operatorname{\theta_i}}}{\frac{\partial \operatorname{g}}{\partial \operatorname{y}}} $$ This seems to match numerical approximations computed by R, but when I try to calculate the Hessian by taking the partial derivative of the gradient with respect to $\theta_i$, I get strange results that are nowhere close to numerical approximations.

This is the formula I came up as my attempt at calculating the Hessian (second-order mixed partial derivative of $f$ with respect to $\theta_i$):

$$ \frac{\partial}{\partial \operatorname{\theta_i}} \Bigl(\frac{\partial \operatorname{f}}{\partial \operatorname{\theta_j}}\Bigl) = \frac{\frac{\partial}{\partial \operatorname{\theta_i}} \Bigl(\frac{\partial \operatorname{f}}{\partial \operatorname{y}}\Bigl) \cdot \frac{-\partial \operatorname{g}}{\partial \operatorname{\theta_i}}}{\frac{\partial \operatorname{g}}{\partial \operatorname{y}}}+ \frac{\frac{\partial \operatorname{f}}{\partial \operatorname{y}} \cdot \frac{\partial}{\partial \operatorname{\theta_i}} \Bigl(\frac{-\partial \operatorname{g}}{\partial \operatorname{\theta_i}}\Bigl)}{\frac{\partial \operatorname{g}}{\partial \operatorname{y}}}+ \frac{\frac{\partial \operatorname{f}}{\partial \operatorname{y}} \cdot \frac{-\partial \operatorname{g}}{\partial \operatorname{\theta_i}}\cdot\frac{\partial}{\partial \operatorname{\theta_i}}\Bigl(\frac{\partial \operatorname{g}}{\partial \operatorname{y}}\Bigl)}{\Bigl(\frac{-\partial \operatorname{g}}{\partial \operatorname{y}}\Bigl)^2} $$

Considering that the values I get from these calculations don't match (not even close) the automatically generated numerical approximations returned by R, I suspect that I've done something wrong here. Can anyone spot an error with the Hessian formula? Is there a property of higher-order mixed partial implicit differentiation that requires a different approach?

Thanks!

UPDATE: Alecos Papadopoulos posted a solution that eliminates the need for iteration by solving directly for $g(\theta)$, thus providing an exact value for $y$. This works perfectly for this problem, as the calculation of the gradient and Hessian does not require implicit differentiation in this case!

For proof of concept (and in case I come across any $g(\theta)$ functions that can't be solved directly), I'm still interested in figuring out what went wrong with my attempt at the Hessian. If anyone has any insight into a general solution for the Hessian using implicit differentiation, it is certainly welcome.

Answer 1

This is a digressing answer, but since computers don't do what we want them to do, but only what they tell them to do, I believe that when we can become more specific, we reduce uncertainty (actual or perceived).

The implicit equation that determines $y$ gives a quadratic equation in $y$, with the roots of the quadratic being functions of the parameters (and in fact $y$ is seen to be the logistic cdf).

For comapactness, denote $\ln \bigl(\frac{-y}{y - 1} \bigl)\equiv h$. Then

$$g(\theta) = 0 \Rightarrow -(h-\theta_1)(h-\theta_3) + a\theta_2(h-\theta_3) + (h-\theta_1)b\theta_4 =0$$ $$\Rightarrow -h^2+\theta_3h+\theta_1h-\theta_1\theta_3+a\theta_2h - a\theta_2\theta_3+hb\theta_4-b\theta_1\theta_4=0$$ $$-h^2+(\theta_1+\theta_3+a\theta_2+b\theta_4)h-(\theta_1\theta_3+a\theta_2\theta_3+b\theta_1\theta_4)=0 $$

Set

$$ \phi_1 \equiv\theta_1+\theta_3+a\theta_2+b\theta_4, \;\; \phi_2= \theta_1\theta_3+a\theta_2\theta_3+b\theta_1\theta_4$$

Then the roots of the polynomial $-h^2+\phi_1h-\phi_2=0$ are

$$h^*_A,h^*_B = \frac {-\phi_1 \pm \sqrt {\phi_1^2 -4\phi_2}}{-2}=\frac {\phi_1}{2}\pm\sqrt {\left(\frac {\phi_1}{2}\right)^2-\phi_2}$$

Then we obtain two equations for $y$

$$\ln \bigl(\frac{-y}{y - 1} \bigl) = h^*_A \Rightarrow y_A = \frac 1{1+e^{-h^*_A}}$$ and

$$\ln \bigl(\frac{-y}{y - 1} \bigl) = h^*_B \Rightarrow y_B = \frac 1{1+e^{-h^*_B}}$$

which is the cdf of the logistic distribution, call it $\Lambda_j,\; j=A,B$, and denote $\lambda_j$ the derivative w.r.t its argument, $\lambda_j = \Lambda_j(1-\Lambda_j)$.

Your log-likelihood becomes

$$\ln L = \ln \Bigl(\frac {n!}{k!(n-k)!}\Bigl)+k\ln\Lambda_j + (n-k)\ln(1-\Lambda_j) $$

Now even the Hessian can be calculated by hand (with patience), let alone the gradient. It is also feasible to check concavity of the log-likelihood. Since you have two equations for $y$, you maximize separately and pick the solution that gives the higher value for the log-likelihood.