1

I need to draw a bivariate normal distribution ellipse based on this article. It says

In the case of the bivariate normal distribution, both approximate and exact methods are available for calculations of both confidence and tolerance ellipses. We present our modified version of the exact methods using common statistics of the simple linear correlation analysis. Given n pairs of observations $x$ and $y$, with standard deviations $s_x$, and $s_y$, and correlation coefficient $r$, one must fix the $\alpha$ probability level and take the Snedecor's $F_\alpha$ value with 2 and $n-2$ degrees of freedom. The semi-axes $L_1$ and $L_2$, and the slopes $b_1$ and $b_2 = -1/b_1$, of the axes of the $100(1 — \alpha)\%$ confidence and tolerance ellipses can be calculated using equations (1) and (2), respectively.

$$ L_1, L_2 = K\sqrt{(n-1)({s_x}^2 + {s_y}^2) \pm \sqrt{[(n-1)({s_x}^2 + {s_y}^2)]^2 - 4(n-1)^2(1-r^2){s_x}^2 {s_y}^2}} \tag{Eq. 1} $$

where

$K = F/n(n-1)$ for confidence ellipses

$K = F(n+1)/n(n-2)$ for tolerance ellipses

$$ b, -1/b = ({s_y}^2 - {s_x}^2)/2rs_xs_y \pm \sqrt{1 + [({s_y}^2 - {s_x}^2)/2rs_xs_y]^2} \tag{Eq. 2} $$

I won't have access to all the data points, only the means, standard deviations and $r$. But to test the equations I'm using the following data.

x = [19, 25, 22, -1, 4, 14, 21, 22, 23, 27, 29, 25, 29, 15, 29, 24, 0, 2, 26, 17, 19, 9, 20, -6, -13, -13, -11, -4, -4, 11, 23]
y = [28, 28, 26, 19, 16, 24, 26, 24, 24, 29, 29, 27, 31, 26, 38, 23, 13, 14, 28, 19, 19, 17, 22, 2, 4, 5, 7, 8, 14, 14, 23]

From this dataset I get

  • $mean_x = 13$
  • $mean_y = 20.226$
  • $s_x = 13.619$
  • $s_y = 8.808$
  • $r = 0.924$

I'm using this calculator (with the parameters $v_1$=2, $v_2$=29, cumulative prob=0.9) to get

$F_{90\%} = 2.5$

I'm calculating the following results in Clojure (that's very similar to Lisp), and here is the code I'm using to get them

(def F 2.5)

(def K (/ (* F (inc n))
          (* n (- n 2))))
  • $K_{90\%} = 0.089$
(def Sx 13.619)
(def Sx2 (* Sx Sx))
(def Sy 8.808)
(def Sy2 (* Sy Sy))
(def r 0.924)
(def r2 (* r r))
(def n 31)

(def L1 (* k
           (Math/sqrt (+ (* (dec n) (+ Sx2 Sy2))
                         (Math/sqrt (- (Math/pow (* (dec n) (+ Sx2 Sy2)) 2)
                                       (* 4 (Math/pow (dec n) 2) (- 1 r2) Sx2 Sy2)))))))
  • $L_1^{90\%} = 11.002$
  • $L_2^{90\%} = 1.985$
(def b (+ (/ (- Sy2 Sx2)
             (* 2 r Sx Sy))
          (Math/sqrt (inc (Math/pow (/ (- Sy2 Sx2) (* 2 r Sx Sy)) 2)))))
  • $b = 0.625$

This is drawing an ellipse like

My wrong ellipse

where I should get something like this

The ellipse I should get

Jp_
  • 131
  • 5
  • https://stats.stackexchange.com/search?q=ellipse+%5Br%5D – whuber Feb 03 '20 at 19:07
  • I don't see an answer that fits to my case in this link. – Jp_ Feb 03 '20 at 19:15
  • 1
    I see many answers that *look* like they fit your case, so could you provide some additional comments that help us see how this is a different question? On the surface, it looks like you have a bug in your code but you haven't supplied enough information to know any more than that. – whuber Feb 03 '20 at 20:29
  • I said I don't have the data points, I'm suppose to get the semi-axes from the standard deviations and the correlation r. I don't see another answer doing it. Is the value of F correct? – Jp_ Feb 03 '20 at 22:31
  • It would be difficult to draw such an ellipse *without* using the covariance matrix in the first place! *Every* answer goes through that process, even if it begins with a dataset. For instance, the [first hit in the search](https://stats.stackexchange.com/questions/9898) gives code that starts right off with the covariance matrix. – whuber Feb 03 '20 at 23:12
  • Yes, this is what my question is about. Aparently the article that I quoted explains how to draw such ellipse, but I must be doing something wrong in the step about the Snedecor's F value. – Jp_ Feb 03 '20 at 23:26
  • Perhaps--but could you show us how you are doing the calculations rather than just exhibiting the results? BTW, the formula you are applying is just an ugly way to represent a fundamentally simple operation (inversion of a $2\times 2$ matrix) and the equation of a circle, as explained in my post at https://stats.stackexchange.com/a/62147/919. – whuber Feb 04 '20 at 15:11
  • Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/104062/discussion-between-jp-and-whuber). – Jp_ Feb 04 '20 at 17:24

1 Answers1

2

Ellipses drawn the right way

It turned out that I lost a lot of time taking the article's formula too serious, it doesn't explain where the formula comes from. So I had to read a lot to understand the context and I'm providing some references for people who's not from maths or stats.

The interesting thing is that I can find $COV(x, y)$ (the covariance of $x$ and $y$) from the Standard Deviations $s_x$, $s_y$ and $r$ (Pearson Correlation Coeficient) as showed here.

Therefore,

$COV(x, y) = COV(y, x) = s_xs_yr_{xy}$

$COV(x, x) = s_xs_x$

$COV(y, y) = s_ys_y$

$C_{xy} = \begin{bmatrix} COV(x, x) & COV(x, y) \\ COV(y, x) & COV(y, y) \end{bmatrix}$ (as showed here)

I calculate the eigen values and eigen vectors based on this refs.

(let [[[eigen-vec-1-x eigen-vec-2-x]
       [eigen-vec-1-y eigen-vec-2-y]] (eigen-vec c-xy)
      [eigen-val-x eigen-val-y]       (eigen-val c-xy)]))

I get the tolerance 50%, 80%, 90% and 95% based on this table (with df=2).

(def tolerances (map chi-squared [0.5 0.8 0.9 0.95]))

The semi-axes and the slope are calculated using

(def x-axis (map #(Math/sqrt (* % eigen-val-x)) tolerances))
(def y-axis (map #(Math/sqrt (* % eigen-val-y)) tolerances))
(def slope  (Math/atan2 eigen-vec-1-x eigen-vec-1-y))

Some other useful references.

Jp_
  • 131
  • 5