clearvars,clc

tic

hit=0; convCt=0;

for trial=1:2000000

  for i=1:4

    x(i)=2*rand-1; y(i)=2*rand-1;

    angle(i)=atan2d(y(i),x(i));

    if angle(i)<0 angle(i)=angle(i)+360;

    end

  end

  if isConvex(x,y)

    convCt=convCt+1;

    [angle,ix]=sort(angle);

    adiff=angle(2:4)-angle(1:3);

    adiff(4)=angle(1)-angle(4)+360;

    if max(adiff)<180

      hit=hit+1;

    end

  end

end

disp([hit convCt  hit/convCt])

toc

function conv=isConvex(x0,y0)

  for pt=1:4

     x=x0;y=y0;

     x(pt)=[]; y(pt)=[];

     x=x-x0(pt); y=y-y0(pt);

     conv=true;

     for i=1:3

       angle(i)=atan2d(y(i),x(i));

       if angle(i)<0 angle(i)=angle(i)+360;

       end

     end

     angle=sort(angle);

     if hasReflex(angle)==false

       conv=false;

       break

     end

  end

end

function t=hasReflex(angle) % requires angle be sorted

  adiff=angle(2:3)-angle(1:2);

  adiff(3)=angle(1)-angle(3)+360;

  if max(adiff)>180

    t=true;

  else

    t=false;

  end

end

The program tests for convexity by finding, for each point, the coordinates of the other points and computing their relative angles. Those angles are sorted, and each must include a reflex angle as the outermost angle, to be considered convex.

The results:

so it was able to do 2 million trials in one minute, of which 1,388,507 were convex, and among those 721,783 included the center of the square (found by a similar method of angles relative to the center). This gave a probability of 52.0% (fairly safe choice of precision).