Maximum of a distribution not invariant in different coordinate systems?

Question

Let's say we have a probability distribution in x,y space:

$$ p(x, y)=\frac{1}{4 \pi} \sqrt{x^2+y^2} \exp \left(- \sqrt{x^2+y^2}\right) $$

This can be converted into polar coordinates as:

$$ p(r, \theta)=\frac{1}{4\pi} r^2 \exp (-r), r\geq0 $$

Now let's say we want to know the most probable distance from the centre, that is r = argmax p. This can be done in two ways.

1. We differentiate $p(r,\theta)$ with respect to $r$, and solve for $$ \frac{\partial p(r, \theta)}{\partial r} = 0 $$ this gives $r=2$

2. We rewrite $p(x,y)$ as

$$ p(x, y)=\frac{1}{4 \pi} r \exp \left(-r\right) $$

Then solving $ \frac{\partial p(x, y)}{\partial r} = 0 $ instead gives $r=1$.

This yields two different answers -- even though they should refer to the same quantity (the most probable distance from the centre).

Questions are -- which is correct? and why? and if both are correct then how do we face this contradiction?

Related questions seem to be on bayesian inference so I included those tags -- though my question is on generic distributions.

Answer 1

I wrote some R code to help visualize this function, it is a very pretty surface.

library(plotly)

x = seq(from = -3, to = 3, by = .2)
y = x
A = matrix(NA,nrow=length(x),ncol=length(y))
for(i in 1:length(x)){
  for(j in 1:length(y)){
    A[i,j] = sqrt(x[i]^2+y[j]^2)*exp(-sqrt(x[i]^2 + y[j]^2))
  }
}


plot_ly(x=x,y=y,z=A)%>%add_surface()

Here is the picture that you get,

As you can see, the highest value of the PDF occurs on the ring where $r=1$. One should be a little careful and call this the ``most likely'' region, as opposed to "most probable", as the probability is technically zero at all points.

You wrote that, $$ p(r,\theta) = (\text{constant})\times r^2\exp\big( -r \big) $$ Where does the $r^2$ come from? When you convert to polar coordinates you replace $x^2+y^2$ by $r^2$. Perhaps, you multiplied by an extra factor of $r$? You only do this multiplication when you replace $dx ~ dy$ by $r ~ dr ~ d\theta$. However, this is not an integration problem, and so you do the usual substitution.

If you wrote it instead as, $$ p(r,\theta) = (\text{constant})\times r\exp\big( -r \big) $$ Then the maximum is achieved exactly at $r=1$. Which we can check with WolframAlpha,

You can see from the R picture that the maximum value of the PDF is approximately $0.36$, the mathematical perfect answer is $\frac{1}{e}$, which agrees with your calculus calculation.