How to fit a boundary to a scatter plot

I am playing around with some diffusion simulations using random walks. For example, if I generate many random walks from the same parent distribution (a Gaussian) as:

 ManyRandomWalks =
Table[
RandomWalkData = RandomVariate[NormalDistribution[0, 1], 100];
RandomWalk = 1 +  Accumulate[RandomWalkData];

RandomWalk,
{i, 1, 200}
]


It will look like this:

One can bound this scatter with the equation: $$f(t) = A + \sqrt{D t} + B t$$

I'd like to make a fit of this to get a more accurate value of $$D$$ -- a diffusion constant -- so far the best method I can think of is to bin the data by index or $$x$$-axis and then perform a statistic or count on the bin and then fit this -- much in a similar how one would fit a histogram.

The other approach might me to do some MLE like FindDistributionParameters and define my function as a PDF, extracting parameter values that way.

Are there any inbuilt features to achieve what I want?

• StandardDeviation[ManyRandomWalks] will give you the sqrt of variance at each of your 100 time steps. You could fit that curve and relate variance as a function of time to your diffusion coefficient. e.g. Show[ ListPlot[ StandardDeviation[ManyRandomWalks] ], ListPlot[ManyRandomWalks], Plot[ Evaluate[ NonlinearModelFit[StandardDeviation[ManyRandomWalks], Sqrt[a*t], {{a, 2}}, t]["BestFit"]] , {t, 0, 100}, PlotStyle -> Red ] , PlotRange -> All ] – N.J.Evans Mar 2 at 16:52
• I was thinking on a similar line of thought except with 2Sqrt[2] StandardDeviation[...] – Q.P. Mar 2 at 17:13
• It depends on what you know about the data generation process and what you want the envelope to represent. If you know all of the parameters of the data generation process AND you want the envelope to contain, say, 99% of all potential (not sample) values, then Plot[{1 + 2.5758293035489004 Sqrt[x], 1 - 2.5758293035489004 Sqrt[x]}, {x, 1, 100}, PlotStyle -> Red]] will do that. – JimB Mar 2 at 17:37
• If you don't know the parameters or the structure of the data generation process but you want the envelope to contain, say, the middle 99% of all sample values, then you could use ListPlot[{Transpose[{Range[100], (Quantile[#, 0.995] &) /@ Transpose[ManyRandomWalks]}], Transpose[{Range[100], (Quantile[#, 0.005] &) /@ Transpose[ManyRandomWalks]}]}, Joined -> True] . – JimB Mar 2 at 17:44
• @JimB Well the point is more that I want to extract a diffusion parameter, $D$, but those are pretty useful commands to plot an envelope. I've added more to my equation in the main text – Q.P. Mar 2 at 17:59

Suppose one had a bunch of time series with no specifics about the data generation mechanism but did know that each time series was generated independently from the others. (Yes, that's using a somewhat loose interpretation of independence.)

Further one wants to estimate an envelope that would contain the central $$100(1-\alpha)$$% of the observations for each time step such that the envelope has the following functional form:

$$f(t)=A\pm \sqrt{D t}+B t$$

One could find the sample estimates of the $$1-\alpha/2$$% and $$\alpha/2$$% quantiles for each time step and then fit a regression to get estimates of the parameters $$A$$, $$B$$, and $$D$$.

Using the OP's simulated data (ManyRandomWalks) we create a dataset that has the time step and a $$-1$$ associated with the lower quantiles and $$+1$$ associated with the upper quantiles and then run NonlinearModelFit.

α = 0.05;
lower = Quantile[#, α/2] & /@ Transpose[ManyRandomWalks];
upper = Quantile[#, 1 - α/2] & /@ Transpose[ManyRandomWalks];
data = Transpose[{Join[Range[n], Range[n]],
Join[ConstantArray[-1, n], ConstantArray[1, n]],
Join[lower, upper]}];
nlm = NonlinearModelFit[data, a + p Sqrt[d t] + b t, {a, b, d}, {t, p}];
nlm["BestFitParameters"]
(* {a -> 0.227842, b -> 0.0196022, d -> 4.00303} *)


Now plot everything:

Show[ListPlot[ManyRandomWalks, Joined -> True, PlotStyle -> Thin],
Plot[{a + Sqrt[d t] + b t, a - Sqrt[d t] + b t} /. nlm["BestFitParameters"], {t, 0, n},
PlotStyle -> {{Thickness[0.01], Red}}],
ListPlot[data[[All, {1, 3}]], PlotStyle -> {{Black, PointSize[0.01]}}]]


• This is a really cool approach. This is what I was after! I have posted my own approach below where I use a multiple of the standard deviation to determine the boundary and then fit that. – Q.P. Mar 3 at 20:00

If the objective is to determine an envelope in the following form

$$f(t)=A\pm \sqrt{D t}+B t$$

that contains a desired expected proportion ($$1-\alpha$$) of the values at each time step, then for the time series mentioned the coefficients are $$A=1$$, $$B=0$$, and $$D=(\Phi^{-1}(1-\alpha/2))^2$$ where $$\Phi^{-1}$$ is the inverse of the standard normal cumulative distribution function. No data is needed.

(* Generate several time series *)
nSim = 200
n = 100
ManyRandomWalks = Table[RandomWalkData = RandomVariate[NormalDistribution[0, 1], n];
RandomWalk = 1 + Accumulate[RandomWalkData];
RandomWalk, {i, 1, nSim}];

(* Set parameters associated with this particular method of generating a time series *)
(* No need to estimate those from the data *)
α = 0.01;  (* Proportion of observations expected to be outside the envelope *)
{a, d, b} = {1,
InverseCDF[NormalDistribution[0, 1], 1 - α/2]^2, 0}

(* Plot the time series and an envelope containing the central 100(1-α)% of the values
for each time step *)
Show[ListPlot[ManyRandomWalks, Joined -> True, PlotStyle -> Thin],
Plot[{a + Sqrt[d t] + b t, a - Sqrt[d t] + b t}, {t, 0, n}, PlotStyle -> {{Thickness[0.01], Red}}]]


I've connected the points as a reminder that there are not 100*200 independent points for which to subject to a regression.

I had a go at my own approach where I run multiple simulations of a random walk, all drawn from the same distribution. I join all the walks and sort them by $$x$$-index, and calculate the standard deviation and mean. Then taking $$a \sigma \pm \mu$$ defines the boundary of the data.

F[A_, B_, C_, D_, t_] := A + B t + C Sqrt[D t];

NSim = 512;
NInd = 128;
\[Sigma]SF = 2Sqrt[2];

\[Mu] = 1;
\[Sigma] = 1;

JoinedManyRandomWalk =
Apply[
Join,
Table[
RandomWalkData = RandomVariate[NormalDistribution[1, 1], NInd];
RandomWalk = Transpose[{Range[NInd], Accumulate[RandomWalkData]}];
RandomWalk,
{i, 1, NSim}
]
];

IndexList = DeleteDuplicates[JoinedManyRandomWalk[[1;;,1]]];

StantisticsManyRandomWalk =
Table[
SortedByIndex = Select[JoinedManyRandomWalk, #[[1]] == IndexList[[i]] &];
\[Mu]SBI = Mean[SortedByIndex[[1;;,2]]];
\[Sigma]SBI = StandardDeviation[SortedByIndex[[1;;,2]]];

{Around[\[Mu]SBI, \[Sigma]SF \[Sigma]SBI], \[Mu]SBI + \[Sigma]SF \[Sigma]SBI, \[Mu]SBI - \[Sigma]SF \[Sigma]SBI, \[Mu]SBI},
{i, 1, Length[IndexList]}
];

UpperFit = NonlinearModelFit[StantisticsManyRandomWalk[[1;;,2]], F[aa, bb, 1, dd, t], {aa, bb, dd}, t];
LowerFit = NonlinearModelFit[StantisticsManyRandomWalk[[1;;,3]], F[aa, bb, -1, dd, t], {aa, bb, dd}, t];
UpperFitParams = {aa, bb, dd} /. UpperFit["BestFitParameters"];
LowerFitParams = {aa, bb, dd} /. LowerFit["BestFitParameters"];

UpperLowerFitsPlot =
Show[
{
Plot[{UpperFit[x], LowerFit[x]},{x, 0, NInd}, PlotStyle->Directive[Red]],
ListPlot[JoinedManyRandomWalk]
},
PlotRange->All, Frame->True
]


I haven't played around with JimB's answer just yet, but the above also works when one has non-zero means for the generating distribution.