I think I've got your idea, you simply want to know when all three points are lying uniformly in a unit disk, how to calculate the possibility of the Circumsphere lying fully in the unit sphere, right?
So, there's quite a lot mistakes in your code:
What is If
doing up here?
Probability
can only handle codes that'll generate True
and False
only, so using If
is not a wise choice as it will return nothing if the the first part yield False
. I suggest you to use conditioned possibility.
Code will not do your desired job.
This code will work better I suppose:
f[x : {a_?NumberQ, b_}, y_, z_] := (And @@ Thread[Norm /@ {x, y, z} <= 1])
&& (Norm[#1] + #2 & @@ Circumsphere[{x, y, z}]) < 1;
f[{.5, .5}, {2, 2}, {.3, .3}]
Probability[f[x, y, z],
{x \[Distributed] UniformDistribution[{{-1, 1}, {-1, 1}}],
y \[Distributed] UniformDistribution[{{-1, 1}, {-1, 1}}],
z \[Distributed] UniformDistribution[{{-1, 1}, {-1, 1}}]}]
Possibility
is simply not designed for these questions.
The usage of Possibility
is quite limited, check this:
f[{x1_, x2_}] := x1 > .5;
Probability[f[x], x \[Distributed] UniformDistribution[{{-1, 1}, {-1, 1}}]]
NProbability[f[x], x \[Distributed] UniformDistribution[{{-1, 1}, {-1, 1}}]]
It will return nothing wihile it should return something like 1/4.
I think we can open a new question and ask how to wrestle with this sort of complex definitions.
Edit 1
It seems that Probability will not work when we split those multidimension variables' components. so, simply change the method we use to define those functions can give you a desired result:
f[a : {ax_, ay_}, b : {bx_, by_},
c : {cx_, cy_}] := (And @@
Thread[Norm /@ {a, b, c} <= 1]) && (Norm[#1] + #2 & @@
Circumsphere[{a, b, c}]) < 1;
NProbability[
f[{ax, ay}, {bx, by}, {cx, cy}], {ax \[Distributed]
UniformDistribution[{-1, 1}],
ay \[Distributed] UniformDistribution[{-1, 1}],
bx \[Distributed] UniformDistribution[{-1, 1}],
by \[Distributed] UniformDistribution[{-1, 1}],
cx \[Distributed] UniformDistribution[{-1, 1}],
cy \[Distributed] UniformDistribution[{-1, 1}]}]
(*0.193461*)
But as I've said before, this is just the probability when point are distributed in a unit square, so the result should be multiplied by 4/Pi
which yield the final result of
0.246322
I don't know whether it's correct as some error message was reported(NIntegrate converge slowly), but I think this shall be a correct way.
Edit 2
I used a simple Monte-Carlo method to proof that my code is correct:
n = 1000000;
N[Count[f @@@ RandomReal[{-1, 1}, {n, 3, 2}], True]/n]
(*0.19485*)
It seems that NProbability generate a proper result this time~
Edit 3
Post a new simpler code:
f[{ax_, ay_}, {bx_, by_}, {cx_,
cy_}] = (And @@
Thread[Norm /@ {{ax, ay}, {bx, by}, {cx, cy}} <=
1]) && (Norm[#1] + #2 & @@
Circumsphere[{{ax, ay}, {bx, by}, {cx, cy}}]) < 1;
Probability[
f[{ax, ay}, {bx, by}, {cx, cy}], {ax, ay, bx, by, cx,
cy} \[Distributed] UniformDistribution[Array[{-1, 1} &, 6]]]
Using {ax, ay, bx, by, cx, cy} \[Distributed] UniformDistribution[Array[{-1, 1} &, 6]]
can save you a lot of characters~ :P
Thus, the final result shall be:
0.24632
RegionCentroid[Circumsphere[{{0, 0}, {1, 0}, {0, 1}}]]
returns {1/2,1/2}, which is presumably not what you want. What isNorm[x] <= 1 && Norm[y] <= 1 && Norm[z] <= 1
going to do for you? $\endgroup$