Let $R$ be a convex region of the plane with unit area, and choose $n\ge 3$ points at random from the region. We will derive a general expression for the probability $P_{n}$ that the convex hull of these points contains a particular point $\tau \in R$ (which we will take as the origin of our coordinate system). We can then evaluate that expression for the particular question posed (in which $n=3$, $R$ is an equilateral triangle, and $\tau$ is its center).
Suppose the hull does not contain $\tau$, and assume that $\tau$ is not collinear with any two of the points (this is true with probability 1). Then we can draw a ray from the center point through one of the $n$ points (the "leftmost" point) such that the remaining $n-1$ points are to the right of the ray. (Note that this condition is both necessary and sufficient. Note also that if the convex hull does contain $\tau$, then there is no "leftmost" point.) We can calculate the desired probability by integrating over all such configurations. In particular, take $f(\varphi)$ to be the area of the subregion that lies to the right of the radial ray at angle $\varphi$, and take $r(\varphi)$ to be the distance from $\tau$ to the boundary of $R$ along that ray. Then the probability that a particular point will lie in the angular interval $[\varphi, \varphi+d\varphi]$ is $da = \frac{1}{2}r(\varphi)^2 d\varphi$, and the probability that each of the other $n-1$ points will lie to its right is $f(\varphi)$. Since there are $n$ points to choose as the leftmost, the overall probability of not capturing $\tau$ in the convex hull is
$$
1-P_{n} = n \int f(\varphi)^{n-1} da = n \int_{-\pi}^{\pi} \frac{1}{2} f(\varphi)^{n-1} r(\varphi)^2 d\varphi,
$$
We can make a few observations at this point. First, if the region $R$ is symmetric under reflection through $\tau$, then $f(\varphi)$ is identically $1/2$: each radial line splits the region in half. The result then is
$$
P_{n} = 1 - \frac{n}{2^{n-1}},
$$
which gives the known result $P_{3} = 1/4$ for polygons with an even number of sides. Second, if the region $R$ instead has bilateral symmetry across the $x$-axis, then $r(\varphi)$ is an even function, and $f(\varphi) = 1/2 + f^{-}(\varphi)$, where $f^{-}$ is an odd function. Then terms in the integral with odd powers of $f^{-}$ must vanish: in particular, we have
$$
P_{3} = \frac{1}{4} - 3 \int{{f^{-}(\varphi)}^2 da} = \frac{1}{4} - 3\int_{-\pi}^{\pi}\frac{1}{2}{f^{-}(\varphi)}^2 r(\varphi)^2 d\phi \le \frac{1}{4},
$$
and interestingly the relation $P_{4} = 2P_{3}$ continues to hold.
It remains to evaluate $f^{-}(\varphi)$, $r(\varphi)$, and the resulting integral in the case of an equilateral triangle. The figure above shows an equilateral triangle with vertices at $(-x,0)$ and $(x/2, \pm x\sqrt{3}/2)$. For $\varphi \in [0,\pi/3]$, the function $f^{-}(\varphi)$ is the area of the blue triangle minus the area of the red triangle,
$$
f^{-}(\varphi) = \frac{x^2}{8}\tan\varphi - \frac{x^2}{2}\frac{1}{\cot\varphi + \cot\frac{\pi}{6}} = \frac{x^2}{8}\tan\varphi \left(1 - \frac{4}{1 + \sqrt{3}\tan\varphi}\right);
$$
and the radius $r(\varphi) = \frac{1}{2}x\sec\varphi$ over the same domain. By symmetry, the full integral is just six times its value over $[0,\pi/3]$:
$$
\begin{eqnarray}
P_3 &=& \frac{1}{4} - 18\int_{0}^{\pi/3}\frac{x^6}{512}\tan^2\varphi\left(1 - \frac{4}{1 + \sqrt{3}\tan\varphi}\right)^2 \sec^2\varphi d\varphi \\
&=& \frac{1}{4} - \frac{1}{36\sqrt{3}}\int_{0}^{\pi/3}\tan^2\varphi\left(1 - \frac{4}{1 + \sqrt{3}\tan\varphi}\right)^2 \sec^2\varphi d\varphi \\
&=& \frac{1}{4} - \frac{1}{12}\int_{0}^{1} u^2\left(1 - \frac{4}{1+3u}\right)^2 du,
\end{eqnarray}
$$
where we have used $x = 2/\sqrt{3\sqrt{3}}$ (in order for the triangle to have unit area) and introduced the transformed variable $u = \tan\varphi / \sqrt{3}$ (so $du = (\sec^2\varphi / \sqrt{3}) d\varphi$). The final integral is a straightforward exercise for the reader; the result is
$$
P_3 = \frac{1}{4} - \frac{1}{324}\left(57 - 80\ln{2}\right) = \frac{2}{81}\left(3 + 10\ln{2}\right) = 0.2452215...,
$$
in agreement with the brute-force and numerical results already given. This result can be generalized easily to the regular $m$-gon for any odd $m$ (changing only the values of some constants), and to the case where $n>3$ (complicating the final integral).
Place the equilateral triangle s.t its vertices are in the points $A=(0,0), B=(1, \sqrt3), C = (2,0)$ (I'm assuming sidelength $2$ to have nicer constants). Denote by $x_1, x_2$ and $x_3$ the $x$-coordinates of the random points on the sides $AC, AB$ and $BC$ respectively. So
$$x_1 \in [0,2], x_2\in[0,1] x\in[1,2].$$
Then the $y$-coordinates are $y_1 = 0, y_2 = \sqrt3 x_2$ and $y_3 = \sqrt3(2-x_3)$. Now use the formula for the area of a triangle given its vertices to get (see for example here)
$$A(x_1, x_2, x_3) = \frac{\sqrt3}{2} \left|x_1x_2 + x_1x_3-2x_1-2x_2x_3+2_2\right|$$
The area of the original equlateral triangle is $\sqrt3$ so we are lead to the requirement that
$$f(x_1, x_2, x_3) \geq 1$$
where $f(x_1, x_2, x_3) = |x_1x_2 + x_1x_3-2x_1-2x_2x_3+2_2|$. So the probability is given by the volume of the region of $X = [0,2]\times[0,1]\times[1,2]$ where this is satisfied, divided by $2$ (which is the total volume of $X$.
I haven't managed to calculate the volume (yet), but maybe we must split it in two parts $f\geq1$ and $f\leq-1$.
EDIT:
I believe it can be calculated in $4$ parts depending on whether $f\geq1$ or $f\leq-1$ and whether $x_1 > 2x_2$ or $x_1 < 2x_2$. Let's denote by $\alpha$ the function that cuts the value to the inteval, i.e
$\alpha(t) = 1$ if $t<1$, $\alpha(t) = t$ if $t\in[1,2]$ and $\alpha(t) = 2$ if $t>2$.
I got the following:
$$I_1 = \int_0^2 { \int_0^{\frac{1}{2}x_1} {\alpha\left( \frac{-1+x_1x_2-2x_2}{x_1-2x_2}\right)} }dx_2dx_1$$
$$I_2 = \int_0^2 { \int_{\frac{1}{2}x_1}^1 {\alpha\left( \frac{1-x_1x_2+2x_1-2x_2}{x_1-2x_2}\right)} }dx_2dx_1$$
$$I_3 = \int_0^2 { \int_0^{\frac{1}{2}x_1} {\alpha\left( \frac{-1-x_1x_2+2x_1-2x_2}{x_1-2x_2}\right)} }dx_2dx_1$$
$$I_4 = \int_0^2 { \int_{\frac{1}{2}x_1}^1 {\alpha\left( \frac{1+x_1x_2-2x_2}{x_1-2x_2}\right)} }dx_2dx_1$$
If my calculations are correct, then the probability is $\frac{1}{2}(I_1+I_2+I_3+I_4)$. Still, I don't really know how to calculate these integrals. Plotting the the constraint equations $f=1$ and $f=-1$ in Wolfram alpha seems to indicate that they are hyperboloids but they are slanted and I dodn't know, maybe some coordinate transformation would be necessary(?), but then again our box is according to the coordinate axis.
EDIT2: Here's my python code I used in the simulation:
import math, random
#----Functions----------------
def area(a, b, c):
return abs(a[0]*(b[1]-c[1])+b[0]*(c[1]-a[1])+c[0]*(a[1]-b[1]))/2.0
def random_between(a, b):
t = random.random()
return (a[0]+t*(b[0]-a[0]), a[1] + t*(b[1]-a[1]))
def mean(a):
return sum(a)/float(len(a))
def var(a):
m = mean(a)
v = 0
for x in a:
v += (x-m)**2
return v/float(len(a))
#-----------------------------
#The verteces of the original triangle
A = (0,0)
B = (2,0)
C = (1,math.sqrt(3))
AREA = area(A,B,C)
HALF_AREA = AREA/2.0
#----- the simulation---------
n = 1000000
count = 0
area_percents = []
for i in xrange(n):
x = random_between(A, B)
y = random_between(B, C)
z = random_between(C, A)
area_of_triangle = area(x,y,z)
area_percents.append(area_of_triangle/AREA)
if area_of_triangle >= HALF_AREA:
count+=1
print "P(area over half) = "+str(float(count)/n)
print "Mean of the area percentages: " + str(mean(area_percents))
print "Variance of the area percentages: " + str(var(area_percents))
## I got the numbers for the distribution like this
## but there must be a better way to make a graph (with some libraries)
"""
graph = {}
categories = 100
d = 1.0/categories
for i in xrange(categories):
# i ~ [i*d, (i+1)*d)
graph[i] = 0
for x in area_percents:
i = 0
while x > (i+1)*d:
i += 1
graph[i] += 1
for x in graph.keys():
print graph[x]
"""
Best Answer
The answer is in fact always $\frac14$, independent of $n$; you forgot to multiply by $2$ for symmetry and divide by $\tau$ for normalization in the circle result.
Wherever the first point is chosen, the "diameter" on which it lies divides the $n$-gon into two symmetric halves, and the second and third points must be in opposite halves. The line connecting them must also lie above the centre (as seen from the first point), and if the second point is at a distance $x$ from the first point along the perimeter (in units of the length of the perimeter), there's an admissible range of length $\frac12-x$ for the third point. Thus the probability is
$$2\int_0^{1/2}\left(\frac12-x\right)\mathrm dx=2\int_0^{1/2}x\,\mathrm dx=\frac14\;,$$
where the factor of $2$ is for symmetry because the second and third point can be interchanged.