Why do we need to triangulate a convex polygon in order to sample uniformly from it? Suppose I want to uniformly sample points inside a convex polygon.
One of the most common approaches described here and on the internet in general consists in triangulation of the polygon and generate uniformly random points inside each triangles using different schemes.
The one I find most practical is to generate exponential distributions from uniform ones taking -log(U) for instance and normalizing the sum to one.
Within Matlab, we would have this code to sample uniformly inside a triangle:
vertex=[0 0;1 0;0.5 0.5]; %vertex coordinates in the 2D plane

mix_coeff=rand(10000,size(vertex,1)); 
          %uniform generation of random coefficients
x=-log(x); %make the uniform distribution exponential
x=bsxfun(@rdivide,x,sum(x,2)); 
   %normalize such that sum is equal to one
unif_samples=x*vertex; 
      %calculate the 2D coordinates of each sample inside the triangle

And this works just fine:

However, using the exact same scheme for anything other than a triangle just fails. For instance for a quadrilateral, we get the following result:

Clearly, sampling is not uniform anymore and the more vertices you add, the more difficult it is to "reach" the corners.
If I triangulate the polygon first then uniform sampling in each triangle is easy and obviously gets the job done.
But why? Why is it necessary to triangulate first?
Which specific property have triangle (and simplexes in general since this behaviour seems to extend to n-dimensional constructions) that makes it work for them and not for the other polygons?
I would be grateful if someone could give me an intuitive explanation of the phenomena or just point to some reference that could help me understand what is going on.
 A: The tldr answer is that in the square case, there are multiple ways to express a "deep" interior point as a convex combination of the vertices, but only one way for points that are nearer to the vertices.
To give a bit more detail and explain why this happens for squares and not triangles, let me first re-express your algorithm in more formal language. Given an integer $n$, define the unit simplex as $\Delta_n=\{(x_1,\dots, x_n)\in \mathbb{R}^n: \sum_i x_i=1, x_i\geq 0\}$. Now, let $p_1, \dots, p_n$ denote the vertices of a polygon. Your algorithm first samples a vector $(w_1, \dots, w_n)$ uniformly from the simplex (using a normalized exponential distribution is a standard way to do this), and then transforms this to the point $\sum_i w_i p_i$.
The first question is why does this work for triangles (that is, when $n=3$)? Well, it is quite easy to see geometrically that $\Delta_3$ is simply an equilateral triangle with vertices at the points $e_1=(1,0,0), e_2=(0,1,0), e_3=(0,0,1)$. In fact, $\Delta_3$ lives in the two-dimensional subspace of $\mathbb{R}^3$ defined by the constraint $x_1+x_2+x_3=1$. Let's call this subspace $H$. Now, I claim that there is an invertible affine-linear map $T: H\to\mathbb{R}^2$ such that $T(e_i)=p_i$ (where recall $p_i$ are the vertices of the triangle we wish to sample from). Why does this claim imply that the random variable $\sum_i w_ip_i$ has a uniform distribution over the triangle? Well by construction the random variable $(w_1,w_2,w_3)=\sum_i w_ie_i$ has a uniform distribution on $\Delta_3$; by the standard change-of-variables formula, the transformed variable $T(\sum_i w_ie_i)=\sum_i w_iT(e_i)=\sum_i w_ip_i$ will have a density that is obtained by multiplying the density of $\sum_i w_ie_i$ by a factor involving the Jacobian of T. Since $T$ is affine-linear its Jacobian is constant, therefore the density of the transformed variable is also constant. I leave it as an exercise to (1) derive an explicit expression for $T$ and (2) to justify the step $T(\sum_i w_ie_i)=\sum_i w_iT(e_i)$ (this step would be immediate for a linear map, but recall $T$ is actually affine linear).
The next question is what is the difference for the quadrilateral case ($n=4$)? Well, now the simplex $\Delta_4$ is actually three-dimensional (it is a solid tetrahedron), so it's impossible to invertibly map it to a 2-dimensional quadrilateral. It's fairly clear geometrically that in this case the density of the convex combination at some value $y:=\sum_i w_ip_i$ will be proportional to the area of the set $C_y:=\{w\in \Delta_4: \sum_i w_ip_i=y\}$, that is, the area of the set of all possible convex coefficients that produce the value $y$ (viz coarea formula).
So now we are led to the consider the geometry of the sets $C_y$ for different $y$ in the quadrilateral. For simplicity, let's assume that the vertices are given by $(0,0), (1,0), (0,1), (1,1)$. If we take $y=(0,0)$, then it is obvious that there is only one way to write $y$ as a convex combination of the vertices of the square, namely $y=1*(0,0)+0*(1,0)+0*(0,1)+0*(1,1)$. Thus $C_y=\{(1,0,0,0)\}$ consists of a single point. On the other hand, consider the center of mass of the square: $y=(.5,.5)$. There are now multiple ways to write this as a convex combination of the vertices, such as:
$$y=.25*(0,0)+.25*(1,0)+.25*(0,1)+.25*(1,1)$$
$$y=.5*(0,0)+0*(1,0)+0*(0,1)+.5*(1,1)$$
$$y=({\frac 1 2}-\alpha)*(0,0)+\alpha*(1,0)+\alpha*(0,1)+({\frac 1 2}-\alpha)*(1,1), 0\leq \alpha\leq 1/2$$
So in this case $C_y$ is isomorphic to a line segment $[0,1/2]$, which is clearly much bigger than a single point.As an exercise, you can try other values of $y$ and convince yourself that the closer $y$ is to a corner, the fewer ways there are to express it as a convex combination of the vertices.
