Re: Length distribution of random secants on a unit
- To: mathgroup at smc.vnet.net
- Subject: [mg95740] Re: [mg95712] Length distribution of random secants on a unit
- From: Bob Hanlon <hanlonr at cox.net>
- Date: Sun, 25 Jan 2009 21:51:02 -0500 (EST)
- Reply-to: hanlonr at cox.net
I did not copy the definition of distPDF
distPDF[d_] = FullSimplify[D[distCDF[d], d]]
Bob Hanlon
---- Bob Hanlon <hanlonr at cox.net> wrote:
=============
When the two points lie on the same side (1/4 of the time)
With[{d = RandomReal[]},
Show[
RegionPlot[Tooltip[Abs[x1 - x2] <= d],
{x1, 0, 1}, {x2, 0, 1}],
Plot[{x1 - d, d + x1}, {x1, 0, 1},
PlotStyle -> {Darker[Blue], Red},
PlotRange -> {0, 1}]]]
Integrating over the valid region to obtain the CDF
Assuming[{0 <= d <= 1},
FullSimplify[Integrate[1, {x, 0, 1}, {y, Max[0, x - d], Min[1, d + x]}]]]
Piecewise[{{3/4, 2*d == 1}, {1, d == 1}, {(-(-2 + d))*d,
0 < d < 1/2 || 1/2 < d < 1}}, 0]
The CDF for this case is then
dist1CDF[d_] = Piecewise[{{1, d >= 1}, {((2 - d))*d,
0 < d < 1}}, 0]
Piecewise[{{1, d >= 1}, {(2 - d)*d, 0 < d < 1}}, 0]
Similarly, when the two points lie on adjacent sides (1/2 of the time)
With[{d = Sqrt[2]*RandomReal[]},
Show[
RegionPlot[Tooltip[Sqrt[x^2 + y^2] < d],
{x, 0, 1}, {y, 0, 1}],
Plot[Sqrt[d^2 - x^2], {x, 0, 1},
PlotStyle -> Red,
PlotRange -> {0, 1}]
]]
dist2CDF[d_] = FullSimplify[Piecewise[{
{1, d >= Sqrt[2]},
{Integrate[1, {x, 0, d}, {y, 0, Sqrt[d^2 - x^2]},
Assumptions -> {0 < d <= 1}], 0 < d <= 1},
{Integrate[1, {x, 0, 1}, {y, 0, Min[1, Sqrt[d^2 - x^2]]},
Assumptions -> {1 < d < Sqrt[2]}], 1 < d < Sqrt[2]}}]]
Piecewise[{{1, d >= Sqrt[2]}, {(d^2*Pi)/4,
Inequality[0, Less, d, LessEqual, 1]},
{Sqrt[-1 + d^2] + (1/2)*d^2*(ArcCsc[d] - ArcCsc[d/Sqrt[-1 + d^2]]),
1 < d < Sqrt[2]}}, 0]
When the two points lie on opposite sides (1/4 of the time)
With[{d = (Sqrt[2] - 1)*RandomReal[] + 1},
Show[
RegionPlot[Tooltip[Sqrt[(x2 - x1)^2 + 1] < d],
{x1, 0, 1}, {x2, 0, 1}],
Plot[{x1 - Sqrt[d^2 - 1], Sqrt[d^2 - 1] + x1}, {x1, 0, 1},
PlotStyle -> {Darker[Blue], Red},
PlotRange -> {0, 1}]]]
Assuming[{1 <= d <= Sqrt[2]}, FullSimplify[Integrate[1, {x1, 0, 1},
{x2, Max[0, x1 - Sqrt[d^2 - 1]], Min[1, Sqrt[d^2 - 1] + x1]}]]]
Piecewise[{{3/4, 2*d == Sqrt[5]}, {1,
d == Sqrt[2]}, {1 - d^2 + 2*Sqrt[-1 + d^2],
1 < d < Sqrt[5]/2 || Sqrt[5]/2 < d < Sqrt[2]}}, 0]
dist3CDF[d_] = Piecewise[{
{1, d >= Sqrt[2]},
{1 - d^2 + 2*Sqrt[-1 + d^2], 1 < d < Sqrt[2]}}]
Piecewise[{{1, d >= Sqrt[2]}, {1 - d^2 + 2*Sqrt[-1 + d^2],
1 < d < Sqrt[2]}}, 0]
The CDF for the distribution is then
distCDF[d_] = FullSimplify[PiecewiseExpand[
dist1CDF[d]/4 + dist2CDF[d]/2 + dist3CDF[d]/4]]
Piecewise[{{1, d >= Sqrt[2]}, {(2 + Pi)/8,
d == 1}, {(1/8)*d*(4 + d*(-2 + Pi)), 0 < d < 1},
{1/2 +
Sqrt[-1 + d^2] + (1/4)*d^2*(-1 + ArcCsc[d] - ArcCsc[d/Sqrt[-1 + d^2]]),
1 < d < Sqrt[2]}}, 0]
Comparing with the experimental data
len = Norm[(First[#] - Last[#])] &;
corners = {{0, 0}, {1, 0}, {1, 1}, {0, 1}};
dir = {{1, 0}, {0, 1}, {-1, 0}, {0, -1}};
p[t_] := Block[{n, r}, n = Mod[IntegerPart[t], 4];
r = FractionalPart[t];
corners[[n + 1]] + r dir[[n + 1]]]
Show[Histogram[Table[len[{p[RandomReal[{0, 4}]],
p[RandomReal[{0, 4}]]}], {100000}],
Automatic, "ProbabilityDensity"],
Plot[distPDF[d], {d, 0, Sqrt[2]},
PlotStyle -> {AbsoluteThickness[2], Red},
PlotRange -> {0, 2.4}]]
Bob Hanlon
---- andreas.kohlmajer at gmx.de wrote:
=============
I need to work with the length distribution of random secants (of two
random points on the perimeter) on a unit square. It's easy to
generate some random data and a histogram. I used the following code
(Mathematica 7.0):
len = Norm[(First[#] - Last[#])] &;
corners = {{0, 0}, {1, 0}, {1, 1}, {0, 1}};
dir = {{1, 0}, {0, 1}, {-1, 0}, {0, -1}};
p[t_] := Block[{n, r},
n = Mod[IntegerPart[t], 4];
r = FractionalPart[t];
corners[[n + 1]] + r dir[[n + 1]]
]
Histogram[
Table[len[{p[RandomReal[{0, 4}]], p[RandomReal[{0, 4}]]}], {100000}],
PlotRange -> All]
The histogram shows a small increase close to 1, a big peak at 1 and
some kind of exponential decay to Sqrt[2] (= maximum).
Does anybody know how to calculate this distribution exactly? What
about moving from a unit square to a random rectangle or a random
polygon? Thanks!
--
Bob Hanlon