Finding the maximum point of probability density function I'm curious about why we always find mle using the first (partial) derivative without checking the end points or singular point or the second (partial) derivative? Thx a lot!
 A: It's not a stupid question at all. See this post for a case where a likelihood can have two maxima and a minimum. 
When dealing with maximum likelihood in a general theoretical approach, we tend to silently assume that the likelihood is a unimodal function (usually having a maximum). Moreover, many "known" distributions have log-concave densities (in their variable). This, coupled with the fact that the unknown coefficients have in many cases a linear relationship with the variable (or we can make it linear through a one-to-one parametrization, which leaves the MLE unaffected), makes the density log-concave in the unknown coefficients also... which are the arguments with respect to which we maximize the (by now, concave) log-likelihood. Satisfaction of the second-order conditions follows, in such cases.  
But in more specific theoretical works, where novel log-likelihoods arise, the researcher has in my opinion the responsibility to treat specifically the issue of whether the second-order conditions are satisfied or not.  
Finally, in applied work, the software algorithms check on their own whether the Hessian is negative definite at the point that they locate as stationary, (and report on the matter) so at least we know whether we have a local maximum or not.
A: First of all, in response to answer from Alecos Papadopoulos, should the software check for negative definite?  Yes?  Do they? I suspect many don't. But actually, if there are any constraints, including bound constraints such as parameters being nonnegative, and one or more constraints are "active" at a candidate solution (e.g., parameter being estimated is on a bound), then checking for negative definiteness of the Hessian, is NOT what should be done. The correct 2nd order condition is that Z' * Hessian * Z be negative semidefinite , where Z is a basis for the null space of the Jacobian of active constraints.  ( Z' * Hessian * Z is the projection of the Hessian into the null space of the Jacobian of active constraints).  If the only active constraints are bounds, then Z' * Hessian * Z amounts to eliminating the rows and columns of parameters on a bound from the Hessian.  Moreover, the first order conditions require the correct sign Lagrange multiplier for each active bound constraint, which amounts to the requirement that any parameter on a lower bound needs to have its gradient component be nonpositive, and any parameter on an upper bound needs to have its gradient component be nonnegative.  And if all first and second order conditions are satisfied, then that only tells you it's a local maximum, unless you know the likelihood function to be concave (or log-concave).
So let's say you found the GLOBAL maximum, but there are several local maxima with likelihood function values almost as high, do you think that maximum likelihood estimation should provide you great confidence in the solution? The confidence intervals which the software spits out are only "valid" relative to that local (even if global) maximum, and will provide you NO indication that similar, or better values, are far outside any confidence intervals you form. If there are many disparate regions with similar likelihood, the maximum likelihood might not be very high likelihood in absolute terms.
