New Answer
OK, here is an almost failsafe way. If a SmoothKernelDistribution
gives you a good estimate of the region where your dense triangle points lie, then this way will find the correct answer. I'm showing the behavior on a random set of points without fixing the Seed
.
We use the kernel distribution to decide which points lie inside the triangle.
selectInnerPoints[pts_] := Module[{sk},
sk = SmoothKernelDistribution[pts];
Select[pts, PDF[sk, #] > 5 &]
];
points = Join[triPoint, singular];
ListPlot[{points, selectInnerPoints[points]}]

What we are going to do is simple. We find the smallest triangle that does not contain any of the inside points. Sounds simple and is simple. First, we need a function to check whether a point is inside a given triangle. Bluntly stolen from Eric and even written with the same variable names:
insideQ[{p1_, p2_, p3_}, v_] :=
Module[{v0 = p1, v1 = p2 - p1, v2 = p3 - p1, a, b},
a = (Det[{v, v2}] - Det[{v0, v2}])/Det[{v1, v2}];
b = -(Det[{v, v1}] - Det[{v0, v1}])/Det[{v1, v2}];
a > 0 && b > 0 && a + b < 1
]
Next, we just need the area of a triangle, because this is what we want to minimize. You know this from school, but in case you don't, Eric has the answer again
area[{p1_, p2_, p3_}] := 1/2 Abs[Det[{p2 - p1, p3 - p1}]]
The next would be a one-liner if it wasn't so tedious to write the variables down. We minimize the area of a starting triangle under the condition that no inner point is outside. We start with a triangle large enough to contain all points which can easily be constructed from the bounding box of all your points.
findTriangle[pts_] := Module[{xmin, xmax, ymin, ymax, innerPoints},
innerPoints = selectInnerPoints[pts];
{{xmin, xmax}, {ymin, ymax}} = MinMax /@ Transpose[pts];
FindMinimum[{area[{{p1x, p1y}, {p2x, p2y}, {p3x, p3y}}],
And @@ (insideQ[{{p1x, p1y}, {p2x, p2y}, {p3x, p3y}}, #] & /@
innerPoints)
},
{{p1x, xmin}, {p1y, ymin}, {p2x, xmax + ymax},
{p2y, ymin}, {p3x, xmin}, {p3y, xmax + ymax}}
]
]
Let's try it.
Show[
ListPlot[{triPoint, singular}],
Graphics[{EdgeForm[Black], FaceForm[Opacity[.2, Blue]],
Polygon[{{p1x, p1y}, {p2x, p2y}, {p3x, p3y}} /. res]}],
Frame -> True,
FrameTicks -> False,
Axes -> False
]

Worked for several different random point sets. I hope I deserve the accept with this.
Original Answer
Since you are speaking of the
the most probable triangle
Why not trying a kernel density estimator which gives you an estimation of the density in an instant? Without adjusting any parameters, let's try a quick hack. I create a SmoothKernelDistribution
from your points, and as you can see, it gives a pretty good estimate where your point density is high
sk = SmoothKernelDistribution[Join[singular, triPoint]];
With[{dens = DensityPlot[
Log@(.01 + PDF[sk, {x, y}]), {x, 0, 1}, {y, 0, 1},
MaxRecursion -> 2, PlotPoints -> 50],
lp = ListPlot[Join[singular, triPoint]]
},
Manipulate[
Show[
dens,
ContourPlot[p == PDF[sk, {x, y}], {x, 0, 1}, {y, 0, 1},
ContourStyle -> Red],
lp
],
{p, 1, 10}
]
]

Now, you can use this to find 3 points that are most far apart. You can use the simple Euclidean distance between the 3 points for that.
dist[{{x1_?NumericQ,y1_},{x2_,y2_},{x3_,y3_}}]:=
(x1-x2)^2+(x1-x3)^2+(x2-x3)^2+(y1-y2)^2+(y1-y3)^2+(y2-y3)^2
And now you try to find the 3 points that maximize this distance function
prob = 3;
res = Quiet@
Last@NMaximize[{dist[{{x1, y1}, {x2, y2}, {x3, y3}}],
PDF[sk, #] > prob & /@ {{x1, y1}, {x2, y2}, {x3, y3}}}, {x1, y1,
x2, y2, x3, y3}];
Plotting the result
Show[
DensityPlot[
Log@(.01 + PDF[sk, {x, y}]), {x, 0.2, .8}, {y, 0.1, .6},
MaxRecursion -> 2, PlotPoints -> 50],
ListPlot[{triPoint, singular}],
Graphics[{EdgeForm[Red], FaceForm[],
Polygon[{{x1, y1}, {x2, y2}, {x3, y3}}] /. res}]
]

Not perfect, but not bad either.