5
$\begingroup$

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

some points

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?

$\endgroup$
  • 1
    $\begingroup$ 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}]; $\endgroup$ – Dr. belisarius Oct 23 '15 at 23:51
  • $\begingroup$ 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...... $\endgroup$ – thils Oct 24 '15 at 0:05
  • 1
    $\begingroup$ 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}] $\endgroup$ – Dr. belisarius Oct 24 '15 at 0:07
  • $\begingroup$ Looks gd! How to generalize to sequence of N numbers $\endgroup$ – thils Oct 24 '15 at 0:14
  • 2
    $\begingroup$ 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. $\endgroup$ – Dr. belisarius Oct 24 '15 at 0:23
4
$\begingroup$

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

Out

$\endgroup$
  • $\begingroup$ Neat code. Why "five" in SeedRandom? $\endgroup$ – thils Nov 5 '15 at 20:00
  • 1
    $\begingroup$ "five" is a random choice whim :) seed is fixed to focus on the dependence on c and r. $\endgroup$ – Gosia Nov 5 '15 at 21:06
5
$\begingroup$

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

Random sample

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

$\endgroup$
  • $\begingroup$ I think you forgot to include c; anyway, data = With[{d = Sort /@ raw}, MapThread[Riffle, {d, d.{1 - c, c}}]]. $\endgroup$ – J. M. is away Oct 24 '15 at 3:19
  • $\begingroup$ @J.M. Ooops! Still working on cutting-and-pasting properly. Thanks for pointing that out. I've added the definition for c. $\endgroup$ – JimB Oct 24 '15 at 3:25
  • $\begingroup$ 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) $\endgroup$ – thils Oct 24 '15 at 3:36
  • 1
    $\begingroup$ @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. $\endgroup$ – JimB Oct 24 '15 at 3:43
  • $\begingroup$ @JimBaldwin Because of normal distribution's omnipresence...its almost everywhere $\endgroup$ – thils Oct 24 '15 at 4:05
3
$\begingroup$

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

enter image description here

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

enter image description here

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.

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – thils Oct 24 '15 at 7:03
  • $\begingroup$ @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 :) $\endgroup$ – ubpdqn Oct 24 '15 at 7:08

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.