1
$\begingroup$

First, set the data and prepare the conditional distribution.


pxy is joint distribution of Random Variable x andy.
py is marginal distribution for variable y of pxy
pcxy is conditional distribution(distribution when x==0)

data = RandomFunction[WienerProcess[0, 1], {0, 1, 0.01}, 1000];
data // Normal // #[[All, 2]] & /@ # & // Set[sample, #] &;

Table[Transpose@{sample[[All, i]], sample[[All, i + 1]]}, {i, 1, 
      99}] // Flatten[#, 1] & // KernelMixtureDistribution // 
  Set[pxy, #] &;
px = MarginalDistribution[pxy, 1];
pcxy = PDF[pxy, {0, y}]/PDF[px, 0];

What I want to do is sampling from pcxy quickly,or fast as long as possible.

Plot[pcxy, {y, -1, 1}]

What I've tried is the following(this didn't end at all):

ProbabilityDistribution[pcxy, {y, -Infinity, Infinity}]
RandomVariate[%, 3] // AbsoluteTiming
$\endgroup$
2
$\begingroup$

The reason it's slow is that the PDF of KernelMixture distributions is given symbolically. This can lead to very large expressions if you're using a lot of points, which would then need to be integrated for RandomVariate to work.

Using an interpolated PDF is a decent alternative

data = RandomFunction[WienerProcess[0, 1], {0, 1, 0.01}, 1000];

sample = data["ValueList"];

pxy = KernelMixtureDistribution[
    Flatten[Transpose /@ Partition[Transpose[sample[[;;, ;; 100]]], 2, 1], 1]
];

px = MarginalDistribution[pxy, 1];

pcxypdf = Interpolation[
    Table[{y, PDF[pxy, {0.0, y}] / PDF[px, 0.0]}, {y, -1, 1, 0.01}]
];

pcxy = ProbabilityDistribution[pcxypdf[y], {y, -1, 1}];

RandomVariate[pcxy, 1000] // AbsoluteTiming // First
(* 0.062171 *)

Show[
    Plot[PDF[pxy, {0, y}] / PDF[px, 0], {y, -1, 1}, PlotStyle -> Red],
    SmoothHistogram[RandomVariate[pcxy, 1000], Automatic, "PDF"],
    ImageSize -> Large
]

enter image description here

$\endgroup$
2
$\begingroup$

One expects the conditional distribution to look very much like a normal distribution and indeed your plot of the estimated density shows that. Therefore, you could sample from a normal distribution with the same mean and standard deviation as the estimated conditional distribution.

data = RandomFunction[WienerProcess[0, 1], {0, 1, 0.01}, 1000];
data // Normal // #[[All, 2]] & /@ # & // Set[sample, #] &;

Table[Transpose@{sample[[All, i]], sample[[All, i + 1]]}, {i, 1, 
      99}] //
    Flatten[#, 1] & // KernelMixtureDistribution // Set[pxy, #] &;
px = MarginalDistribution[pxy, 1];
pcxy = PDF[pxy, {0, y}]/PDF[px, 0];

ymean = NIntegrate[y pcxy, {y, -1, 1}]
ysd = NIntegrate[(y - ymean)^2 pcxy, {y, -1, 1}]^0.5
Plot[{pcxy, PDF[NormalDistribution[ymean, ysd], y]}, {y, -1, 1},
 PlotStyle -> {{Thickness[0.02], LightGray}, {Thin, Red}}]

Conditional distribution and normal approximation

Future samples could then be found with

yFuture = RandomVariate[NormalDistribution[ymean, ysd], 100];
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.