Probability density function while using spherical coordinates I'm following this tutorial where at somepoint the derived PDF for spherical coordinates for a Lambertian surface is

\begin{array}{l}
p(\theta, \phi) = \dfrac{\sin \theta}{2 \pi}.
\end{array}

But as soon as they compute a sample, the result is instead divided by $ \dfrac{1}{2\pi} $, which as they say is the "pdf of the integral"

Why isn't it divided by $ \frac{\sin \theta}{2 \pi} $ instead?

If we were using differential steradians over the unit hemisphere, the only possible probability density function integrating to 1 is infact $ \frac{1}{2\pi} $

But if we separate the integral over the hemisphere traced by spherical coordinates

$$ 
\int_{\phi = 0}^{2\pi} \int_{\theta = 0}^{\frac{\pi}{2}}  \,  \frac{\sin{\theta}}{2\pi}d\theta d\phi= 1. 
$$

The PDF now becomes $ \frac{\sin \theta}{2 \pi} $ yet they still divide by $ \frac{1}{2\pi} $

EDIT: after carefully reviewing the concept of PDFs and integration over the hemisphere I'm starting to think the article I've linked is making a substantial error, mixing the idea of importance sampling with the pdf of choosing a direction of reflectance from a lambertian surface

Radiance is defined as
$$ L_{(x,\omega)} = \frac{\mathrm{d}^2\Phi}{\mathrm{d}\omega\ cos\theta\ \mathrm{d}x}  $$
Since it's defined over differential solid angles, we can interpret the result of one sample as if it was the flux density "over a unit steradian"

If we use monte carlo estimation and find "the average flux density over a single steradian" and multiply the result by $2\pi$, we get irradiance:

$$ (2\pi) \frac{1}{n}\sum^n L_{(x,\omega)} $$. 
But in this particular case, $2\pi$ has nothing to do with pdfs! since it's the integral domain used for the monte carlo estimation!

Instead, the real pdf is computed for $\theta$ and $\phi$ because that is the probability density function of choosing one direction over the other, according to the particular properties of a lambertian surface. A mirror-like surface has a different probability of choosing one direction over the other, but this has nothing to do us with dividing the sample with $\frac{1}{2\pi}$. It would be different if we were using importance sampling, but in this case it seems like we're not

Is my reasoning correct? If not, what am I missing?
Normally $sin(\theta)$ isn't written as part of the pdf, because it's an artifact from using double integral for spherical integration. The double integral does the integration over a rectangular domain (in you case of size $2\pi$ x $\pi/2$), while you want integration over a sphere/hemisphere, so $sin(\theta)$ is added there as a weight function to transform rectangular to spherical integration.

A common notation is to use integration over solid angle with single integral and differential solid angle $d\omega$ denoting 2D integration over sphere, and you won't see the $sin(\theta)$ term in there.

To clarify, the notation with single integral has constant infinitesimally small solid angle $d\omega$ over the sphere. However, for the double integral over $\theta$ and $\phi$ the solid angle isn't constant (think of sphere tessellation) but a function of $sin(\theta)$. In case of Monte Carlo integration and even distribution of samples over the hemisphere, each sample represents solid angle of $2\pi/N$ steradians and there's no $sin(\theta)$ solid angle weighting required.