Simulating a bounded continuous random walk in $n$-dimensions

Question

You're given three $n$-dimensional real vectors, $\mathbf{x_0}$, $\mathbf{a}$ and $\mathbf{b}$, with $\mathbf{a} \le \mathbf{x_0} \le \mathbf{b}$ (vector inequality means component-wise inequality). Starting at $\mathbf{x_0}$, simulate a bounded continuous random walk, such that the end point $\mathbf{x}$ of the walk is always inside the bounds $\mathbf{a} \le \mathbf{x} \le \mathbf{b}$.

I want to program this with Mathematica. The parameters of the simulation are the bounds $\mathbf{a}, \mathbf{b}$, the starting point $\mathbf{x_0}$, and a parameter $\sigma$ that controls the standard deviation of the step size. What's a simple and efficient way to do this?

Here's what I've tried:

The following function Walk takes a vector $\mathbf{x}$ and returns the next point on the random walk.

Walk[x_List, a_List, b_List, sigma_] := 
 x + RandomVariate[
   TruncatedDistribution[Transpose[{a - x, b - x}],
    MultinormalDistribution[ConstantArray[0, Length[x]], 
     sigma*IdentityMatrix[Length[x]]]]]

For some reason this is exceedingly slow. You can try it yourself. After executing:

a = ConstantArray[0, 10];
b = ConstantArray[1, 10];
x = RandomReal[{0, 1}, 10];

The following line of code takes forever on my PC:

Walk[x, a, b, 4]

Why?

Answer 1

I do not know why, but the following are much faster than Walk:

ClearAll[walkF,walkF2];
walkF[x_List, a_List, b_List, sigma_] := 
     With[{f =  TruncatedDistribution[Transpose[{a - x, b - x}], 
      ProductDistribution @@  ConstantArray[NormalDistribution[0, sigma], Length[x]]]},
      x + RandomVariate[f]]

walkF2[x_List, a_List, b_List, sigma_] := 
  With[{f = ProductDistribution @@ 
    Thread[TruncatedDistribution[Transpose[{a - x, b - x}], NormalDistribution[0, sigma]]]},
      x + RandomVariate[f]]

a = ConstantArray[0, 10];
b = ConstantArray[1, 10];
x = RandomReal[{0, 1}, 10];

walkF[x, a, b, 4] // Timing
(* {0.015625, {0.836782, 0.333754, 0.402554, 0.953688, 0.680472, 
               0.578463, 0.0481464, 0.728613, 0.0854795, 0.0142773}} *)

3D

x = RandomReal[{0, 1}, 3];
nl = NestList[walkF[#, {0, 0, 0}, {1, 1, 1}, 2] &, x, 100]; // Timing
(* {0.171875,Null} *)
nl2 = NestList[walkF[#, {0, 0, 0}, {1, 1, 1}, 2] &, x, 1000]; // Timing
(* {1.703125,Null} *)

{g1, g2} = Graphics3D[{Directive[Orange, Opacity[0.8], Specularity[White, 30]],
            CapForm[None], JoinForm[None], Tube[#]}, Boxed -> False, 
            Background -> Black, ImageSize -> 400] & /@ {nl, nl2};

Row[{g1, g2}, Spacer[5]]