Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a list of $\{x_i,y_i\}$ pairs and I want to show that $x_i$'s and $y_i$'s have positive correlation. The result of Correlation[Transpose[data][[1]], Transpose[data][[2]]] is 0.856906 and ListPlot[data, PlotRange -> All, AxesLabel -> {X, Y}] results in the following plot:

ListPlot of the simulation data

I just wonder if there is any other (better) type of plot in Mathematica to show the correlation between the variables?

share|improve this question
2  
You might get better answers asking at stats.stackexchange.com Then you can come back with the plot suggestions to this site and ask about implementing them in Mathematica. –  Szabolcs Feb 28 '12 at 13:52

4 Answers 4

up vote 20 down vote accepted

You can also calculate the Coefficient of Determination, R Squared.

This is the same as the correlation squared, but by making use of LinearModelFit you can create some additional graphics.

To make a sample distribution you can use this:

CreateDistribution[] := DynamicModule[{savepts = {{-1, -1}}},
  Dynamic[
   EventHandler[
    ListPlot[pts, AxesOrigin -> {0, 0}, 
     PlotRange -> {{0, 7}, {0, 5}}], 
    "MouseDown" :> (savepts = 
       pts = DeleteCases[
         Append[pts, MousePosition["Graphics"]], {-1, -1}])],
   Initialization :> (pts = savepts)]]

CreateDistribution[]

Just click to add some points. The data is collected in the variable pts.

enter image description here

Then calculate R Squared:

lm = LinearModelFit[Sort@pts, a, a]; r2 = lm["RSquared"];
Show[Plot[lm[x], {x, 0, 7}], ListPlot[pts], AxesOrigin -> {0, 0}, 
 PlotRange -> {{0, 7}, {0, 5}}, 
 PlotLabel -> 
  "The correlation is " <> 
   If[D[lm["BestFit"], a] < 0, "negative", "positive"], 
 Epilog -> 
  Inset[Style[
    "\!\(\*SuperscriptBox[\"R\", \"2\"]\) = " <> ToString[r2], 
    11], {1.5, 4.5}]]

enter image description here

Whether the correlation is positive or negative is obtained here from the derivative of the BestFit.

You can add standard deviation bands for sigma = 1, 2 & 3 like this.

lm = LinearModelFit[Sort@pts, {1, x (*, x^2 *)}, x];
{bands68[x_], bands95[x_], bands99[x_]} = 
  Table[lm["SinglePredictionBands", 
    ConfidenceLevel -> cl], {cl, {0.6827, 0.9545, 0.9973}}];
Show[ListPlot[Sort@pts], 
 Plot[{lm[x], bands68[x], bands95[x], bands99[x]}, {x, -0.15, 7.2}, 
  Filling -> {2 -> {1}, 3 -> {2}, 4 -> {3}}], AxesOrigin -> {0, 0}, 
 PlotRange -> {{0, 7}, {0, 5}}, ImageSize -> 480, Frame -> True]

enter image description here

Uncomment x^2 for a quadratic fit.

enter image description here

share|improve this answer

Basic, but useful options

Something that certainly helps in representing your data is to set the range of x and y to be equal:

ListPlot[data, PlotRange -> {{0,3},{0,3}}, AxesLabel -> {X, Y}]

Moreover, if the slope is also important you should set the aspect ratio to 1:

ListPlot[data, PlotRange -> {{0,3},{0,3}}, AxesLabel -> {X, Y}, AspectRatio -> 1]

Confidence ellipses

In some cases it is convenient to add a confidence ellipse.

Let's generate some data:

data = Table[
        RandomVariate[BinormalDistribution[{50, 50}, {5, 10}, .8]], {1000}];

Here we get the estimates of the generated data:

estDist = EstimatedDistribution[data, 
            BinormalDistribution[
               {\[Mu]1, \[Mu]2}, {\[Sigma]1, \[Sigma]2}, \[Rho]]]

And, here there are two useful functions that I found somewhere some time ago. Unfortunately, I'm not able to refer to the original author.

Ellipse[{{r_, M_}, m_}, {x_, y_}] := 
 Show[ContourPlot[({x, y} - r).M.({x, y} - r) == m, 
       {x, 0, 100}, {y, 0, 100}, ContourStyle -> {Red, Thick}], 
      ListPlot[{r}]
 ];

showEllipse[r0_, s1_, s2_, rho_, perc_] := 
  Module[{dchi2, fi, fic, V, InvV}, 
    Clear[\[Sigma], \[Rho]];
    dchi2 = Quantile[ChiSquareDistribution[2], perc]; 
    fi = Which[s1 == s2, \[Pi]/4 Sign[rho], 
               s1 > s2, 1/2 ArcTan[(2 rho s1 s2)/(s1^2 - s2^2)], 
               s1 < s2, 1/2 ArcTan[(2 rho s1 s2)/(s1^2 - s2^2)] - \[Pi]/2]; 
    fic = If[fi < 0, \[Pi] + fi, fi]; 
    V = {{s1^2, rho s1 s2}, {rho s1 s2, s2^2}}; 
    InvV = Inverse[V]; 
    Ellipse[{{r0, InvV}, dchi2}, {x, y}]
  ]

Now we can plot the data and the 95% confidence ellipse:

Show[{
 ListPlot[data, PlotRange -> {{0, 100}, {0, 100}}, AspectRatio ->1],
 showEllipse[estDist[[1]], estDist[[2, 1]], estDist[[2, 2]], estDist[[3]], .95]
}]

enter image description here

Or, even multiple confidence ellipses (95% and 99%):

Show[{
 ListPlot[data, PlotRange -> {{0, 100}, {0, 100}}, AspectRatio ->1],
 showEllipse[estDist[[1]], estDist[[2, 1]], estDist[[2, 2]], estDist[[3]], .99],
 showEllipse[estDist[[1]], estDist[[2, 1]], estDist[[2, 2]], estDist[[3]], .95]
}]

enter image description here

share|improve this answer
5  
Very interesting. I managed to find some built-in functionality for this. It gives smoother ellipses for a small data set. Needs["MultivariateStatistics`"]; Show[{ListPlot[data, PlotRange -> {{0, 100}, {0, 100}}, AspectRatio -> 1], Graphics[EllipsoidQuantile[data, {0.95, 0.99}]]}] –  Chris Degnen Feb 28 '12 at 21:20
    
@ChrisDegnen, the example you're providing is even more interesting, I've never came across that function before. Thanks for pointing it out. –  vikkor Feb 29 '12 at 7:16
    
This method was subsequently used here and here. –  Chris Degnen Sep 3 '12 at 12:42

Another approach could be also based on the smooth kernel distribution.

If this are our data:

data = Table[RandomVariate[BinormalDistribution[{50, 50}, {5, 7}, .8]], {1000}];

We can obtain a representation of a multivariate smooth kernel distribution based on the data values:

skd = SmoothKernelDistribution[data];

And here we get the plot:

ContourPlot[
    Evaluate@PDF[skd, {x, y}], {x, 30, 70}, {y, 30, 70},
    PlotRange -> All, PlotPoints -> 100]

Mathematica graphics

This approach can be helpful to get a glimpse on how the data are distributed. If you consider a case in which the data originate from two different distributions, plotting them in this way gives you a clear picture.

data1 = Table[RandomVariate[BinormalDistribution[{50, 50}, {5, 7}, .8]], {500}];
data2 = Table[RandomVariate[BinormalDistribution[{45, 55}, {7, 5}, .8]], {500}];
data = Partition[Flatten[Join[{data1, data2}]], 2];
skd2 = SmoothKernelDistribution[data];
ContourPlot[Evaluate@PDF[skd2, {x, y}], {x, 30, 70}, {y, 30, 70},
    PlotRange -> All, PlotPoints -> 100]

Mathematica graphics

share|improve this answer

I tried to use the above mentioned codes for plotting error ellipses for 2D-Data. However, I did not get the anticipated results, because an error occured when Mathematica tried to solve the equality in the function Counterplot for my data.

I found another solution based on the explicit calculation of the ellipse by means of covariance analysis. The ellipse is then rotated to coincide with the principal axis of the correlated data. The original explanation can be found at:

http://www.visiondummy.com/2014/04/draw-error-ellipse-representing-covariance-matrix/

Since I did not find an equivalent solution at stackexchange, I adapted the provided matlab code given by Vincent Spruyt to mathematica and answered to this 'historic' post.

define some random correlated data:

(*Random correlated data generation*)
s = 3;
var = 0.3;
mus= Table[RandomReal[], {i, 10000}];

x = RandomVariate[NormalDistribution[#, var]] & /@ (+s mus);
y = RandomVariate[NormalDistribution[#, var]] & /@ (-s mus);
data = {x, y}\[Transpose];

Perform a Covariance analysis and create a 2D-list, that contains the ellipse as points (percLevel is the percentage covering a given confidence interval)

(*define error Ellipse*)
ErrorEllipse[data_, percLevel_] :=
  Module[
   {eVa, eVec, coors, dchi, thetaGrid, ellipse, rEllipse},
   {eVa, eVec} = Eigensystem@Covariance[data];

   (* Get the coordinates of the data mean*)
   coors = Mean[data];

   (*Get the perc [%] confidence interval error ellipse*)
   dchi = \[Sqrt]Quantile[ChiSquareDistribution[2], percLevel/100];

   (* define error ellipse in x and y coordinates*)
  thetaGrid = Table[i, {i, 0, 2 \[Pi], 2 \[Pi]/99}];
   ellipse = {dchi \[Sqrt]eVa[[1]] Cos[thetaGrid], 
     dchi \[Sqrt]eVa[[2]] Sin[thetaGrid]};

   (* rotate the ellipse and center ellipse at coors*)
   rEllipse = coors + # & /@ (ellipse\[Transpose].eVec)
   ];

Use LineListPlot to show the error ellipses:

(* visualize results*)
percentOneSigma = 66.3;
percentTwoSigma = 95.4;

Show[
 ListPlot[data],
 ListLinePlot[ErrorEllipse[data, percentOneSigma],PlotStyle -> {Thick, Red}],
 ListLinePlot[ErrorEllipse[data, percentTwoSigma],PlotStyle -> {Thick, Blue}],
 ListPlot[{Mean[data]}, PlotStyle -> {Thick, Red}]
 ]

enter image description here

Hope, this might be useful.

Cheers

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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