I derived this with the following Euler code.
>: assume(s>0)$
>: ah : (integrate(cos(t)-cos(s),t,0,s)+integrate(cos(s)-cos(t),t,s,%pi))/%pi
2 sin(s) - s cos(s) + (%pi - s) cos(s)
--------------------------------------
%pi
>: integrate(ah*sin(s),s,0,%pi)/%pi, ratsimp
3
-----
2 %pi
>n=10000; A=random(3,n); x=cos(A*2pi); y=sin(A*2pi);
>function F(x1,y1,x2,y2,x3,y3) := abs(det([x2-x1,y2-y1;x3-x1,y3-y1]))/2
>mean(Fmap(x[1],y[1],x[2],y[2],x[3],y[3]))
0.4777744668333
>3/(2pi)
0.4774648292757
It is easier to do this in the continuum (though it is possible for
discrete points too). If you fix two points, forming a chord in the
circle, the average height can be computed, then you got to compute
the average area for all chords. I did that user Euler Math Toolbox
below, including a numerical simulation to make sure, I got everything
right. My result was 3/(2pi).
>: assume(t>0)$
>: s1 : (integrate(cos(t)-cos(s),s,t,%pi)+integrate(cos(s)-cos(t),s,0,t))/%pi
2 sin(t) - t cos(t) + (%pi - t) cos(t)
--------------------------------------
%pi
>: integrate(s1*2*sin(t),t,0,%pi)/2/%pi
3
-----
2 %pi
>function F(x1,y1,x2,y2,x3,y3) := abs(det([x1-x2,y1-y2;x1-x3,y1-y3]))/2
>n=10000; A=random(3,n);
>x=sin(A*2*pi); y=cos(A*2*pi);
>S=Fmap(x[1],y[1],x[2],y[2],x[3],y[3]);
>mean(S)
0.4746595823366
>3/(2*pi)
0.4774648292757
>