Generating correlated random variables with normal distribution

Question

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?

Answer 1

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}]

Answer 2

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.

Answer 3

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.