MLE for Independent Exponentials Here is the problem statement:

Let $(X_i,Y_i)$ be a random sample from a distribution with pdf
  $$f(x,y;\theta)=
> \frac{1}{\theta^3}\exp\left(\frac{-x}{\theta}-\frac{y}{\theta^2}\right)\qquad
> 0<x, 0<y.$$ (a) Find the MLE for $\theta$.
(b) Give the asymptotic distribution of
  $\sqrt{n}(\hat{\theta}-\theta).$

Here are my thoughts:
We basically have two random samples, $X\sim \exp(\theta)$ and $Y\sim \exp(\theta^2)$. This pdf is the joint pdf for $X$ and $Y$. The parameter has no impact on the support of the variables $X$ and $Y$, so calculus should suffice for this problem.
$$L(x,y;\theta) = \prod_{i=1}^n f(x,y;\theta) = \frac{1}{\theta^{3n}}\exp\left(\frac{-1}{\theta^2}\sum_{i=1}^n \theta x_i+y_i\right),$$
so 
$$l(x,y;\theta) = -3n\ln(\theta)-\frac{1}{\theta^2}\sum_{i=1}^n \theta x_i + y_i,$$ which has partial $\theta$ derivative
$$ \frac{-3n\theta^2+\sum_{i=1}^n \theta x_i +2y_i}{\theta^3}.$$
Equating to zero and solving for $\theta$ requires solving a quadratic in $\theta$. 
The second $\theta$ derivative is
$$\frac{3n\theta^2 - \left(\sum_{i=1}^n 2\theta x_i + 6y_i\right)}{\theta^4}.$$
The solutions to the quadratic are 
$$\theta = \frac{\sum_{i=1}^n x_i\pm \sqrt{\left(\sum_{i=1}^nx_i\right)^2+24n\sum_{i=1}^n y_i}}{6n}.$$
Is there no easier way to do the problem? Am I forced to make this substitution to find the maximum?
 A: I hesitated to 'answer' the question before, because it was not clear if this was a homework assignment, ... but as you have now answered it yourself (and neatly so), I thought I might briefly reply to the question posed:  
$$\text{"Is there no easier way to do the problem?"}$$
In exactly this vein, my co-author, Murray Smith, and I sought to automate these sorts of MLE problems, symbolically and exactly, in precisely the way you have proceeded - step by step - but getting the computer to do the grunt work. See, for instance:


*

*Rose, C and Smith, M.D. (2000), Symbolic maximum likelihood estimation with Mathematica, Journal of the Royal Statistical Society, Series D: The Statistician, 49(2), 2000, 229-240. Download available: here
A more sophisticated version of same was then built into the mathStatica software package which we later developed.
To illustrate, for your example:
GIVEN: $(X,Y)$ have a bivarate Exponential distribution, with parameter $\theta  >0$, with joint pdf $f(x,y)$:

(source: tri.org.au) 
Activate the special SuperLog function:

(source: tri.org.au) 
For a random sample of size $n$ drawn on $(X,Y)$, the log-likelihood function for parameter $\theta$ is:

(source: tri.org.au) 
The score function is the gradient of the log-likelihood with respect to $\theta$:

(source: tri.org.au) 
Setting the score to zero and solving for $\theta$ corresponds to the first-order condition:

(source: tri.org.au) 
We require the second (positive) solution. As for second-order conditions, the Hessian, evaluated at the optimal solution, is:

(source: tri.org.au) 
... which is strictly negative, given that $X>0$ and $Y>0$.
