Uniformly distributed points on unit disk Consider a circle of radius $R$. Generate a homogeneous distribution of random points within the circle using the transformation method (inverse sampling).
Hint: If we draw $(r,\varphi)$ uniformly in $[0,R)\times[0,2\pi)$ the points are not homogeneously distributed within the circle.
We apply the transformation method to the problem. Let $z_1$ and $z_2$ be uniformly distributed in $[0,1)$. Then we have to solve the following integral equation for $r$ and $«varphi$
$$z_2z_2 = \frac{1}{\pi R^2}\int_0^r\int_0^\varphi r' dr'd\varphi = \frac{1}{2\pi R^2}r^2\varphi \tag{1}$$
We chose $\varphi=2\pi z_2$ and $r=R\sqrt{z_1}$.
E.g. for $N=1000$ and $R=1$ we get

Now I don't fully understand how we end up with (1). I understand the hint and I know that we have to "preserve area" or rather "preserve density". But I was't able to get (1).
I'd like to write down a simple equality that has to hold and then apply a change of variables which I then expect to result in (1).
How do you get (1)?
I think the starting point would be: The CDF needs to be equal before and after transforming, leading us to:
$\int_0^y P(y')dy' = \int_0^z P_u(z')dz' \tag{2}$
whereas $P_u$ is our uniformly distributed. $P$ could be differently distributed. In our case, $P_u$ and $P$ are uniformly distributed.
$P_u$ is uniformly distributed on the unit rectangle while $P$ is uniformly distributed on a circle with radius $R$. Leading to:
$P_u(z_1, z_2) = \begin{cases} \frac{1}{\text{Vol(Unit Rectangle)=1}}=1, & (z_1, z_2)\in [0,1)^2 \\ 0, & \text{else}\end{cases} \tag{3}$
$P(z_1, z_2) = \begin{cases} \frac{1}{\text{Vol(Circle)}=\pi R^2}=\frac{1}{\pi R^2}, & (z_1, z_2) \in \text{Circle(R)},\\ 0, &\text{else}\end{cases} \tag{4}$
So we get
$$\int_0^{z_1}\int_0^{z_2} 1 dz_1'dz_2' \stackrel{!}{=} \int_0^{z_1}\int_0^{z_2} \frac{1}{\pi R^2} dz_1'dz_2' = \int_0^r\int_0^\varphi \frac{1}{\pi R^2} r' dr' d\varphi \tag{5}$$
Note that we transformed to polar coordinates in the second equal sign. (I'm unsure about this step thought)
and then we end up with:
$$z_1z_2=\frac{1}{\pi R^2}\frac{1}{2}r^2 \varphi \tag{6}$$
Is that correct?
 A: Be aware that all your double integrals must be evaluated according to the standard methods of double integrals, according to which
$$  \int_0^{z_1}\!\!\int_0^{z_2} \frac{1}{\pi R^2}\,\mathrm dz_1'\mathrm dz_2'
 = \frac{1}{\pi R^2} \int_0^{z_1}\!\!\int_0^{z_2} \mathrm dz_1'\mathrm dz_2'
 \neq \int_0^{z_1}\!\!\int_0^{z_2} 1\,\mathrm dz_1'\mathrm dz_2',$$
so your "$\stackrel{!}{=}$" claim is incorrect.
There is another problem with the way you have written the integrals, which is that as written above, the integration variable $z_1'$ ranges from $0$ to $z_2$, not from $0$ to $z_1$ as you probably intended.
This turns out not to affect the result of the integrals over Cartesian coordinates, since the only place $z_1'$ and $z_2'$ occur is to name the variables of integration, and renaming these variables is permitted.
But in the integral over polar coordinates, this does matter.
In fact,
\begin{align} \int_0^r\!\int_0^\varphi \frac{1}{\pi R^2} r'
 \mathrm dr' \mathrm d\varphi'
&= \int_0^r\left(\int_0^\varphi \frac{1}{\pi R^2} r'
 \mathrm dr'\right) \mathrm d\varphi'\\
&= \int_0^r\left(\frac{\varphi^2}{2\pi R}\right) \mathrm d\varphi'\\
&= \frac{r \varphi^2}{2\pi R},
\end{align}
which is not what you wanted at all.
Remember that in a multiple integral the variables of integration named by d match the integral signs from the outside inward (or from the inside outward), not left to right or right to left.
(I have also changed your $\mathrm d\varphi$ to $\mathrm d\varphi'$
in order to be consistent with the way you name your other variables of integration.)
In the rest of the answer I will put the variables of integration in the order in which you presumably intended, not in the order you wrote them, so as to reproduce the correct integrals for this exercise.

So what is actually going on in this exercise?
In this context, each double integral is also computing the probability
of an event, namely, the event that some coordinates will fall within some
region.
For example, the integral
$$ \int_0^{1/4}\!\!\int_0^{1/4} 1\,\mathrm dz_2'\mathrm dz_1' $$
is the probability that the Cartesian coordinates $(z_1',z_2')$ of a point
in the unit square will fall in the region where
$0 < z_1' < \frac14$ and $0 < z_2' < \frac14,$
that is, the probability that a point randomly sampled from the red square below
(with uniform distribution) will fall within the purple region.

Notice that this has some nice properties. In particular, if we increase or decrease the width of the rectangle, the probability increases or decrease in proportion to the width; and likewise, the probability increases or decreases in proportion to the height of the rectangle.
Similarly, if $R = 1,$ the integral
$$ \int_0^{1/2}\!\!\int_0^{\pi/2}
 \frac{1}{\pi R^2} r'\, \mathrm d\varphi' \mathrm dr' $$
is the probability that the polar coordinates $(r',\varphi')$ of a point in the disk of radius $R=1$ (uniformly distributed) will fall in the region where
$0 < r' < \frac12$ and $0 < \varphi' < \frac\pi2.$
That is, given a uniform distribution over the red disk in the figure below, this particular integral gives the probability that a randomly-selected point falls in the blue sector.

Now notice that as we increase or decrease the angle of the sector, the probability increases or decreases in proportion to the angle.
But if we increase or decrease the radius of the sector, the probability does not increase or decrease in proportion to the radius.
Instead, it increases or decreases in proportion to the square of the radius.
This is related to the fact that as you get farther from the origin, the grid cells in the figure above get larger.
This is also related to the factor $r'$ in the integral.
What we're actually trying to do here is to set up an invertible continuous function between $z_1$ and $r$ and another such function between $z_2$ and $\varphi$ with a certain additional requirement.
Such functions will give us a mapping from a rectangle inside the unit square
(such as the purple region above) to a sector of a circle within a circle of radius $R$ (such as the blue region in the previous figure).
The "additional requirement" is that the probability of choosing a point in the circular sector, given sampling with the desired distribution over the circle of radius $R,$ is the same as the probability of choosing a point in the rectangle when we sample with uniform distribution over the unit square.
This requirement is not trivial to enforce, because we're taking a region of a particular shape (a rectangle) over one probability distribution and trying to equate its probability with that of a region of a very different shape (a circular sector) over a different probability distribution.
We are not just lifting the rectangle out of the square and dropping it onto the circle; we are shrinking one edge of the rectangle down to a point, curving the other edge out into a circular arc, and distorting all the angles and distances between points inside the rectangle in order to fit them within this new shape.
In other words, the probability described by
$$ \int_0^{z_1}\!\!\int_0^{z_2} 1\,\mathrm dz_2'\mathrm dz_1' \tag A$$
is the probability measure of the rectangular region
$[0,z_1]\times[0,z_2]$ in Cartesian coordinates on the unit square.
On the other hand, the probability described by
$$ \int_0^r\!\!\int_0^\varphi \frac{1}{\pi R^2} r'\,d\varphi' dr' \tag B$$
is the probability measure of the sector
$[0,r]\times[0,\varphi]$ in polar coordinates on the circle of radius $R.$
The properties we would like to obtain are that the points in the rectangle correspond one-to-one with the points in the sector, and that the probability of choosing a point in the rectangle is the same as the probability of choosing a point in the sector.
And that property is satisfied if the two integrals are always equal.
We get the equation
$$z_2z_2 = \frac{1}{\pi R^2}\int_0^r\int_0^\varphi r'\,d\varphi' dr'$$
simply by setting integral (A) equal to integral (B) and simplifying the integral on each side.
In summary, what I am saying is that if we map the random variables $z_1$ and $z_2$ onto the Cartesian coordinates of the circle,
any rectangle $[0,z_1]\times[0,z_2]$ will map onto some region within the circle,
and if we do the mapping correctly then our randomly selected point (randomly selected by choosing $z_1$ and $z_2$ uniformly from $[0,1]$
and mapping to the circle)
will have the correct probability of landing in that region within the circle.
That is why we want the probabilities to be equal, and therefore we decide what mapping to use by making the probabilities equal and then deducing what mapping we must have used in order to make this true.
The equation between the two integrals is not something deduced by necessity from first principles. Rather, it is something we set up because it is relatively easy to work with and because it gives us an algorithm for distributing points uniformly in the circle, given only the ability to distribute points uniformly within the segment $[0,1]$ as a prerequisite.

Now you might ask why we cannot use some intermediate form of coordinates to relate the two distributions (one on the unit square, the other on the circle). In particular, we can use Cartesian coordinates to describe points in the circle, the probability function on those points is
$1/(\pi R^2),$ and you can compute probabilities of regions within the circle using integrals of the form
$\int_s^t\int_u^v 1/(\pi R^2)\,\mathrm dx\mathrm dy,$ as you might expect.
The difficulty is that if we try to use Cartesian coordinates $x$ and $y$ to describe the points in the circle, let's say by assigning $x$ to the horizontal axis and $y$ to the vertical axis, the range of $x$ is from $-R$ to $R$ rather than from $0$ to $1.$
But worse still, for any given value of $x,$ we must let $y$ range from
$-\sqrt{R^2 - x^2}$ to $\sqrt{R^2 - x^2}.$
So to measure the probability of the entire circle, the integral would be
$$ \int_{-R}^R\int_{-\sqrt{R^2 - x^2}}^{\sqrt{R^2 - x^2}}
  \frac{1}{\pi R^2}\,\mathrm dy\mathrm dx. $$
For the region described by polar coordinates $[0,r]\times[0,2\pi]$,
which is a circle of radius $r,$ the integral would be
$$ \int_{-r}^r\int_{-\sqrt{r^2 - x^2}}^{\sqrt{r^2 - x^2}}
  \frac{1}{\pi R^2}\,\mathrm dy\mathrm dx, $$
whereas for the blue sector in the figure illustrating the circular distribution, which has polar coordinates $[0,\frac12]\times[0,\frac\pi2]$,
the integral would be
$$ \int_0^r\int_0^{\sqrt{(1/4) - x^2}}
  \frac{1}{\pi R^2}\,\mathrm dy\mathrm dx. $$
For sectors where $\varphi$ is not $0,$ $\frac\pi2,$ $\pi,$ or $2\pi$ the integrals become even more complicated, and generally you would need to integrate two parts of the region separately and add the integrals together.
In short, trying to use Cartesian coordinates to represent the same region described by the integral
$$ \frac{1}{\pi R^2}\int_0^r\int_0^\varphi r'\,d\varphi' dr'$$
is an insanely complicated mess compared to anything we actually need to do.
Worse still, there is no obvious way to relate the Cartesian coordinates
$x$ and $y$ to the coordinates $z_1$ and $z_2$ uniformly distributed in the unit square. Merely using the name $z_1$ for $x$ and $z_2$ for $y$
and using integrals of the form $\int_0^{z_1}\!\!\int_0^{z_2}$
absolutely will not work, because you end up only being able to describe rectangular regions in the first quadrant of the circle, and you cannot even exactly describe the entire first quadrant using such a rectangle.
In short, your attempted bridge between the Cartesian coordinates of the square and the polar coordinates of the circle in Equations $(4)$ and $(5)$ is not only arithmetically incorrect (the integrals you say are equal are not equal), it is geometrically incorrect: it does not enable you to express the probabilities of regions within the circle that you need to express.
