# Generating correlated random variables with normal distribution

This site and this site provide some background as to what I am trying to achieve, but with some variations in what is required. First I want three random variables with a normal distribution, however they are correlated as in $d_i= c a_i +(1-c) b_i$, where $c$ is the correlation factor, and $(a_i, b_i, d_i)$ are the variables. The constraint $b_i < d_i < a_i$ is a requirement.

Here I have a Monte Carlo based simulation that yields all three variables which satisfies the criteria mentioned above:

RaN[m_, s_, co_] :=
Module[{me = m, sd = s, c = co},
q1 = RandomVariate[NormalDistribution[me, sd]];
q2 = RandomVariate[NormalDistribution[me, sd]];
p = c*q1 + (1 - c)*q2;
If[TrueQ[q2 < p < q1], {q1, q2, p}, RaN[m, s, co]]]

Manipulate[
ListPointPlot3D[Table[RaN[0, 1, cr], {i, 1000}]], {cr, 0.005, 0.995}]


However this is not an optimized code. It can be expensive time-wise, and also there is the possibility of running into recursion depth problems if total runs is increased.

My question is: Are there better ways of seeking, not just 3 correlated variables, but a sequence of $N$ variables such that the variables $a_i$ is the highest and $b_i$ is the lowest numbers. The site mentioned earlier points to a technique that involves matrices. Any suggestions?

• This seems to have the same output properties than your function and doesn't suffer the "recursive trap", but I'm not sure about the "statistical truth" of your statements RaN[m_, s_, c_] := Append[#, c*#[[1]] + (1 - c)*#[[2]]] &@ Sort@RandomVariate[NormalDistribution[m, s], {2}]; Oct 23, 2015 at 23:51
• If I change "2" to "5", I get {-0.669003, -0.00770918, 0.253854, 0.25824, 0.517137, -0.66239}. How do you get 3 numbers sandwiched between the lowest & highest...... Oct 24, 2015 at 0:05
• Or much faster RaN[m_, s_, c_, n_] := Append[#, c*#[[1]] + (1 - c)*#[[2]]] & /@ (Sort /@ RandomVariate[NormalDistribution[m, s], {n, 2}]); Manipulate[ListPointPlot3D[RaN[0, 1, cr, 1000]], {cr, 0.005, 0.995}] Oct 24, 2015 at 0:07
• Looks gd! How to generalize to sequence of N numbers Oct 24, 2015 at 0:14
• I don't fully understand if the way you're generating the "correlated" rvs is sound. I prefer to leave and see what other more statistically-savvy users have to say. Oct 24, 2015 at 0:23

Analytic approach:

Manipulate[
SeedRandom["five"];
ListPointPlot3D[
RandomVariate[
TransformedDistribution[{a, b, c*a + (1 - c)*b}, {{a, b} \[Distributed]
OrderDistribution[BinormalDistribution[r], {1, 2}]}], 10^3],
PlotLabel -> Row[{"c = ", c, " | ", "r = ", r}]],
{{c, 0.5}, 0, 1}, {{r, 0}, -.99, .99}]


• Neat code. Why "five" in SeedRandom? Nov 5, 2015 at 20:00
• "five" is a random choice whim :) seed is fixed to focus on the dependence on c and r. Nov 5, 2015 at 21:06

You only have to deal with choosing $a$ and $b$ because if $c$ (a weighting factor not a correlation factor) is between 0 and 1, then $d$ has to be between $a$ and $b$. (If $c$ is outside that range, I don't think you can get there from here.)

What you've used above is a type of rejection sampling. If the means are the same for $a$ and $b$, then you're only wasting about half of the samples which doesn't seem too inefficient but you could speed up the process by generating a little over twice the number of needed observations at one time and then select the first $n$ that satisfy the desired ordering.

If $a$ and $b$ are selected from a bivariate normal distribution with the same means (mu) and same standard deviations (sigma) and a correlation coefficient (rho), then the desired triplet can be determined as the following:

$$(\min(a,b), c\max(a,b)+(1-c)\min(a,b), \max(a,b))$$

Here is some Mathematica code that can certainly be made more efficient:

(* Set some parameters for a bivariate normal and take a random sample; *)
mu = 0;    (* Mean *)
sigma=1;   (* Standard deviation *)
rho=-0.5;  (* This correlation coefficient can be between -1 and +1 *)

n=1000;    (* Sample size *)
raw=RandomVariate[BinormalDistribution[{mu,mu},{sigma,sigma},rho],n];

(* Create a table with the minimum (b) and maximum (a) of each pair
of random samples and a placeholder for d *)
data=Table[{Min[raw[[i,All]]],0,Max[raw[[i,All]]]},{i,n}];

(* Determine the weighted average of a and b *)
c = 0.4;
data[[All,2]]=c data[[All,3]]+(1-c)data[[All,1]];

(* Plot the resulting samples *)
ListPointPlot3D[data, BoxRatios->{1, 1, 1}]


Note that none of the elements of the triplet will have a normal distribution but the triplet is generated from a bivariate normal distribution.

• I think you forgot to include c; anyway, data = With[{d = Sort /@ raw}, MapThread[Riffle, {d, d.{1 - c, c}}]]. Oct 24, 2015 at 3:19
• @J.M. Ooops! Still working on cutting-and-pasting properly. Thanks for pointing that out. I've added the definition for c.
– JimB
Oct 24, 2015 at 3:25
• Nimble code! Just wish the elements satisfy normal distribution. Is there a quick way to check if a distribution is normal (even though this is not part of the query) Oct 24, 2015 at 3:36
• @thils. Therapy for removing a feeling of dependence on normality is not as quick but better in the long run. But more seriously, why do you need to have the variables have a normal distribution? That should probably be another question or a discussion in chat.
– JimB
Oct 24, 2015 at 3:43
• @JimBaldwin Because of normal distribution's omnipresence...its almost everywhere Oct 24, 2015 at 4:05

I am not certain what the ultimate aim here (in particular correlation relationship). I post this in the event it may be helpful. In the following a and b are independent (standardized) normal random variables that are correlated with (standardized) normal variable d but in such a way that when a is poorly correlated b is highly correlated.

mn[c_] :=
MultinormalDistribution[{0, 0,
0}, {{1, 0, 1 - c}, {0, 1, c}, {1 - c, c, 1}}]
Manipulate[
Column[{
Show[ListPointPlot3D[RandomVariate[mn[1 - c], 10000],
AxesLabel -> {"a", "b", "d"}, BaseStyle -> 20],
Plot3D[{c, 1 - c}.{x, y}, {x, -3, 3}, {y, -3, 3}, Mesh -> None,
PlotStyle -> {Pink, Opacity[0.3]}], ImageSize -> 400],
Show[ListPointPlot3D[RandomVariate[mn[1 - c], 10000],
AxesLabel -> {"a", "b", "d"}, BaseStyle -> 20,
RegionFunction -> Function[{x, y, z}, y < z < x]],
Plot3D[{c, 1 - c}.{x, y}, {x, -3, 3}, {y, -3, 3}, Mesh -> None,
PlotStyle -> {Pink, Opacity[0.3]}], ImageSize -> 400]
}],
{c, 0.05, 0.995, Appearance -> "Labeled"}]


The lower plot is the truncated sample based on constraint in OP.

Some conditional probabilities illustrating the relationship between variables and "correlation" c.

tab = Table[{j,
NProbability[
z > 0.5 \[Conditioned] x > 0.5, {x, y, z} \[Distributed]
mn[1 - j]],
NProbability[
z > 0.5 \[Conditioned] y > 0.5, {x, y, z} \[Distributed]
mn[1 - j]]}, {j, 0.05, 0.95, 0.05}];
ListPlot[{tab[[All, {1, 2}]], tab[[All, {1, 3}]]}, Frame -> True,
GridLines -> {None, {Probability[x > 0.5,
x \[Distributed] NormalDistribution[]]}},
PlotLegends -> {"P(d>0.5|a>0.5)", "P{d>0.5|b>0.5}"}]


with gridline probability P(Z>0.5), Z is N(0,1) and symmetry around c=0.5

I re-iterate I am not sure what is the aim and whether this post has any relevance. I am happy to delete for lack of relevance.

• Close to what I am after. The aim: correlated random events are part of the essence of quantum entanglement...with specific distributions you can proceed to explore the bizarre quantum world Oct 24, 2015 at 7:03
• @thils I would not accept my answer as it is just exploratory and wait for better answers...I hope you can play but I am sure others will have closer and answers rather than longer comment :) Oct 24, 2015 at 7:08