1
$\begingroup$

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

$\endgroup$
8
  • $\begingroup$ stats.stackexchange.com/search?q=ellipse+%5Br%5D $\endgroup$ – whuber Feb 3 '20 at 19:07
  • $\begingroup$ I don't see an answer that fits to my case in this link. $\endgroup$ – Jp_ Feb 3 '20 at 19:15
  • 1
    $\begingroup$ 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. $\endgroup$ – whuber Feb 3 '20 at 20:29
  • $\begingroup$ 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? $\endgroup$ – Jp_ Feb 3 '20 at 22:31
  • $\begingroup$ 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 gives code that starts right off with the covariance matrix. $\endgroup$ – whuber Feb 3 '20 at 23:12
2
$\begingroup$

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.

$\endgroup$

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.