Answered: With what probability do $4$ points placed uniformly randomly in the unit square of $\mathbb{R}^2$ form a convex/concave quadrilateral?

Here is a different approach to the problem. Consider a random triangle illustrated below. Regions are labelled near ($N_i$) and far ($F_i$) from point $i$. $T$ is the triangle.

If the 4th point falls in regions $F_1$, $F_2$, or $F_3$ a convex quadrilateral will result. The probability that 4 points produce a convex shape is equivalent to the probability that the 4th point will be in one of those regions. That is, $p=E[F_1+F_2+F_3]$ where $F_i$ is the random variable expressing the area of that region.

To simplify the calculations, consider the shaded region labelled $A_1$ and the counterparts $A_2$ and $A_3$. $$F_1+N_2+N_3=A_1$$ $$N_1+F_2+F_3+T=1-A_1$$ Combining these one finds, $$F_1+F_2+F_3 = 2 - 2T - A_1 - A_2 - A_3$$ and therefore, $$p=E[F_1+F_2+F_3]=2 - 2E[T] - 3E[A]$$ Fortunately others have worked out $E[T]={11\over144}$

To work out $E[A]$, I considered separately the case where the unit square is bisected across opposite sides (as is the case for $A_1$). You can show that this occurs with probability ${2\over3}$, given two uniform random points. The smallest of the two areas $A_\min$ follows the distribution $$f(a_\min)=4a_\min$$ for $0<a_\min<{1\over2}$. Given the two points $(2,3)$ point $1$ point will fall in the smaller of the two areas with probability $a_\min$, in which case $A_1=1-a_\min$. $$ E[A_{opp}]=16\int_0^{1\over2}a_\min^2(1-a_\min)\,da_\min = {5\over12} $$

For the case where the unit square is cut across a corner (as is the case for $A_2$), $A_\min={1\over2}X_0Y_0$ where $X_0$ and $Y_0$ are independent random variables with the identical distribution, $f(x)=2x$ for $0<x<1$. This gives the distribution for $A_\min$, $$f(a_\min)=-16 a_\min \log 2a_\min$$ for $0<a_\min<{1\over2}$. $$ E[A_{corner}]=-32\int_0^{1\over2}a_\min^2(1-a_\min)\log 2a_\min\,da_\min = {23\over72} $$ And so, overall $$E[A] = P(A_{opp})E[A_{opp}]+P(A_{corner})E[A_{corner}]={2\over3}{5\over12}+{1\over3}{23\over72} = {83\over216}$$ Finally, $$ p = 2 - 2{11\over144} - 3{83\over216}={25\over36}=\left({5\over6}\right)^2=0.69\bar{4} $$ which is close to your simulation result. A very simple form!

$P_{convex}=\int_0^1 da a_{123}p(a_{123})$ where $a_{123}=\frac{1}{2}\left|(\mathbf{x}_{12} \times \mathbf{x}_{13})\right|$ this is the probability that the 4th point falls within the triangle formed by the first three. Normalization included due to unit square, and $$p(a)=\frac{1}{2} \int_0^1 dx_1 \int_0^1 dx_2 \int_0^1 dx_3 \int_0^1 dy_1 \int_0^1 dy_2 \int_0^1 dy_3 \delta \left( a - \sqrt{\left( (x_2 - x_1)^2 (y_3 - y_1)^2 + (y_2 - y_1)^2 (x_3 - x_1)^2 \right) } \right)$$

an integral which might be difficult to evaluate in closed form, where $\delta(x)$ is the Dirac delta function.

The $p(a)$ probability density function for the area of the triangle formed by three uniformly random points in the unit square seems to fit quite well to the beta distribution with $\alpha = 1.18$ and $\beta = 5.1$ and an upper bound of $\frac{1}{2}$

The percentage of fourth points that fall inside the triangle formed by first three is then the expectation of the area $E\left[a\right]=\frac{\alpha / 2}{\alpha+\beta}$ or ~10%, corresponding to concave fraction of 90%. Seems off ...

Also looked at the distribution of the maximum area of the 4 triangles formed by a set of 4 points each picked from $(U[0,1],U[0,1]) Mean moves from 9% to 15%