EDIT: Refactored code. Old code can be found in the end of the answer.
Initial settings. We take 20
lattice-free unit steps starting from {10, 0, 0}
. Four such paths are calculated.
steps = 20;
start = {10, 0, 0};
paths = 4;
We do this in three dimensions, and whenever we perform random sampling, we strive for 10000
samples.
dims = 3;
samples = 10000;
An uniformly distributed unit-length lattice-free random walk step in dim
dimensions can be computed using this trick. This function computes the resulting position after steps
.
ClearAll[randomWalk];
randomWalk[steps_, dim_] :=
Total[Normalize /@
RandomVariate[NormalDistribution[0, 1], {steps, dim}]]
The core idea of this answer is to use probability distributions which represent likelihood of reaching origin in certain amount of unit-length random walk steps. These distributions are estimated from sample walks, for which in this fully symmetric case the distance covered is the only interesting parameter.
This function is used to compute parameter estimators for PERTDistribution
, which is really a nicely reparameterized variant of BetaDistribution
. It is left as an exercise to the reader to observe that distributions of distances from origin after N steps can be approximated with beta distribution.
ClearAll[distributionEstimates];
distributionEstimates[steps_, dim_, samples_, estimates_] :=
NonlinearModelFit[#, a Sqrt[x] + b, {a, b}, x]["Function"] & /@
Transpose@
ParallelTable[{{n, a}, {n, b}} /.
FindDistributionParameters[
Table[Norm@randomWalk[n, dim], {samples}],
PERTDistribution[{0, n}, a, b],
MapThread[{#1, #2@n} &, {{a, b}, estimates}]], {n, 2, steps}]
This is slightly tricky, but the goal is to fit parameter estimation functions better and better and produce distribution parameter estimation functions which give robustness to EstimatedDistribution
in the next step.
ClearAll[estimates];
estimates =
Last@NestWhile[
Apply[{2 #1,
distributionEstimates[steps, dims, #1, #2]} &], {100, {Sqrt,
Sqrt}}, First@# < 2 samples &];
Here we compute an array of functions providing PDF of random walk covering certain distance. This information is used to drive the actual source-to-destination random walk later on. Distributions for walks of 0 and 1 steps are not computed, because they're both useless for our method, and also too much for EstimatedDistribution
to handle.
ClearAll[stepWeightFunctions];
stepWeightFunctions[steps_, dim_, samples_, estimates_] :=
ParallelTable[
PDF@EstimatedDistribution[Table[Norm@randomWalk[n, dim], {samples}],
PERTDistribution[{0, n}, a, b],
MapThread[{#1, #2@n} &, {{a, b}, estimates}]], {n, 2, steps}]
We compute weights for our case here. These functions could be reused later on.
ClearAll[weightFunctions];
weightFunctions = stepWeightFunctions[steps, dims, samples, estimates];
This is the core of our solution. It simply performs weighted RandomChoice
walk form source. Weights are calculated on basis of distance of candidate step from the origin, and amount of steps left.
Two last steps are not computed this way, because last step is always origin, and second-last step lies on intersection of two predetermined spheres. This can't be accomplished with dumb random step candidates.
ClearAll[incompleteWalks];
incompleteWalks =
ParallelTable[
Last /@ NestList[
Apply[Function[{level, pt}, {level - 1,
pt + (RandomChoice[
weightFunctions[[level - 1]]@Norm[pt + #] & /@ # -> #] &[
Table[randomWalk[1, dims], {samples}]])}]], {steps - 1,
start}, steps - 2], {paths}];
Walks are completed by FindInstance
which provides a step that satisfies above constraints. FindInstance
isn't really random, but given all the other randomness in this routine, this might be good enough.
ClearAll[completeWalks];
completeWalks =
With[{v = Table[Unique[], {dims}]}, #~Join~{v, Table[0, {dims}]} /.
FindInstance[Norm[v] == 1 && Norm[v - Last@#] == 1, v,
Reals] &] /@ incompleteWalks;
Four resulting walks in 3D:
completeWalks //
MapIndexed[{Hue[(First@#2 - 1)/(paths + 1)],
Tube@#1} &] // Graphics3D

By setting initial dimension and starting point parameters one can also accomplish 2D random walk (and probably other dimensions too):
completeWalks //
MapIndexed[{Hue[(First@#2 - 1)/(paths + 1)], Line@#1} &] // Graphics

Old code:
This, rather messy contraption attempts to drive a random walk process from start
point over paths
to origin in fixed amount of unit-distance steps
. I'm not really convinced; I think it's lazy and spends time near the starting point more than it should, and then hurries to destination when indicators become strong enough it is starting to run out of steps to reach origin. Probability distributions which drive RandomChoice
are estimated from sample random walks.
With[
{steps = 50, paths = 4, start = {20, 0, 0}},
With[
(* estimated probability functions for every amount of steps left *)
{p =
ParallelTable[
Function[x,
Evaluate@
PDF[EstimatedDistribution[#,
PERTDistribution[{0, n}, a, b],
(* these are here just to speed up estimation *)
{{a, 0.3407759964445141` + 0.7140058349745467` Sqrt[n]},
{b, -5.88256361017025` + 5.436462403959963` Sqrt[n]}}],
x]] &[
Norm@Total@# & /@
Map[Normalize,
RandomVariate[
NormalDistribution[0, 1], {10000, n, 3}], {2}]],
{n, 2, steps - 1}]},
ParallelTable[
Last /@ NestList[
Apply[Function[{level, pt},
{level - 1,
pt + (RandomChoice[
p[[level - 1]]@Norm[pt + #] & /@ # -> #] &[
Normalize /@
RandomVariate[
NormalDistribution[0, 1], {10000, 3}]])}]],
{steps - 1, start}, steps - 2], {paths}]] //
(* find a suitable second-to-last point *)
Map[#~Join~{{x, y, z}, {0, 0, 0}} /.
FindInstance[
Norm[{x, y, z}] == 1 && Norm[{x, y, z} - Last@#] == 1,
{x, y, z}, Reals] &]] //
MapIndexed[{Hue[Mod[GoldenRatio #2, 1]], Tube@#1} &] // Graphics3D

BrownianBridgeProcess
seems promising, but the steps are not of the same length. $\endgroup$