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?


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.

  • $\begingroup$ either site, I think, is fine. So long as you don't cross post (as in post an active version to both sites at once). I would be happy to migrate it to Stats for you if that is what you want, but why not let it hang here for a couple days at least to see if you get a good answer? $\endgroup$ Commented Oct 24, 2013 at 6:47

1 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 $$


$$ \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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.