install.packages("car")

In the open source stat program R, you can use the function ellipse. All it needs is the coordinates of the center of the ellipse and the variance-covariance matrix of data XY (variances on the diagonal and covariance on the off-diagonal). Alternatively, the function data.ellipse will plot the data and ellipse together for you. E.g., supposing that "x" and "y" are vectors containing your data in R, the following syntax will get you a plot of the data along with the median and 97.5th percentile:

library(car)

data.ellipse(x, y, levels=c(0.5, 0.975))

An easy way to draw such a plot is to use the freeware statistical program PAST (Hammer, Ø.; Harper, D. & Ryan, P. D. (2001): PAST: Paleontological Statistics Software Package for Education and Data Analysis. Palaeontologia Electronica 4: 9).

If you use x- and y-values and a grouping variable, you can draw 95%-ellipses for each group easily by ticking the respective box in the editor.

http://folk.uio.no/ohammer/past/

You can also use Hotelling's T-squared, for instance http://sites.stat.psu.edu/~ajw13/stat505/fa06/10_pHotel/06_pHotel_CI.html

This way you can decide between different confidence levels (90% or 95%, or other).

Best,

install.packages("car")

In the open source stat program R, you can use the function ellipse. All it needs is the coordinates of the center of the ellipse and the variance-covariance matrix of data XY (variances on the diagonal and covariance on the off-diagonal). Alternatively, the function data.ellipse will plot the data and ellipse together for you. E.g., supposing that "x" and "y" are vectors containing your data in R, the following syntax will get you a plot of the data along with the median and 97.5th percentile:

library(car)

data.ellipse(x, y, levels=c(0.5, 0.975))

Sounds heteroscedastic (WLS) through the origin. I think it would be more useful to try something like the attached. 

Data CRE Prediction 'Bounds' and Graphs Example for Section 4 of ...

Thank you very much to everyone!!

The problem as stated has multiple solutions; there is an infinite number of ellipses that will contain 95% of your data points. Even if you add the constraint that the center of the ellipse has to coincide with the average (x_bar, y_bar) of your data points, you still may have multiple solutions. Using a program in R to draw contours for a given percentile of a bivariate normal distribution does not guarantees that the resulting ellipse will contain 95% of your data points - specially if your data is not normally distrubuted. On the other hand, we usually want confidence regions for multivariate parameters not for individual data points; maybe what you want to construct is a prediction region for a future observation.

From what I can make of the graph, it appears you are not including young children.  Correct?  Perhaps this might appear somewhat nonlinear?  And/Or perhaps some stratification might be useful, depending upon your objective? 

A good explanation of how to do this (with code samples for Matlab and C++): http://www.visiondummy.com/2014/04/draw-error-ellipse-representing-covariance-matrix/

 Have you already solved the problem? If not I think I can help you

This example gives you the best fit of the confidence interval.

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

Here is a Solution, if you prefer R:

If your data is saved as `X` then the following code gives the 95% confidence interval.

Note that, `center = mean vector`, `shape = your covariance matrix` and the `radius ` is the square root of the chisquare value at 0.05 with p degrees of freedom. The following code is for a 2 dimensional data

```{r}

X = matrix(rnorm(100^2, 5, 1), ncol =2 , byrow = F)

plot(X, xlab = 'x', ylab = 'y')

ellipse(center = colMeans(X), shape = cov(X),radius = sqrt(qchisq(.95, df=2)),col = 'darkgreen')

```

See my video on Youtube on how to draw one in Excel:

https://www.youtube.com/watch?v=MZ8JRGy-t1I&t=277s

clear all;

close all;

% Create some random data

s = [2 2];

x = randn(334,1);

y1 = normrnd(s(1).*x,1);

y2 = normrnd(s(2).*x,1);

data = [y1 y2];

% Calculate the eigenvectors and eigenvalues

covariance = cov(data);

[eigenvec, eigenval ] = eig(covariance);

% Get the index of the largest eigenvector

[largest_eigenvec_ind_c, r] = find(eigenval == max(max(eigenval)));

largest_eigenvec = eigenvec(:, largest_eigenvec_ind_c);

% Get the largest eigenvalue

largest_eigenval = max(max(eigenval));

% Get the smallest eigenvector and eigenvalue

if(largest_eigenvec_ind_c == 1)

smallest_eigenval = max(eigenval(:,2))

smallest_eigenvec = eigenvec(:,2);

else

smallest_eigenval = max(eigenval(:,1))

smallest_eigenvec = eigenvec(1,:);

end

% Calculate the angle between the x-axis and the largest eigenvector

angle = atan2(largest_eigenvec(2), largest_eigenvec(1));

% This angle is between -pi and pi.

% Let's shift it such that the angle is between 0 and 2pi

if(angle < 0)

angle = angle + 2*pi;

end

% Get the coordinates of the data mean

avg = mean(data);

% Get the 95% confidence interval error ellipse

chisquare_val = 2.4477;

theta_grid = linspace(0,2*pi);

phi = angle;

X0=avg(1);

Y0=avg(2);

a=chisquare_val*sqrt(largest_eigenval);

b=chisquare_val*sqrt(smallest_eigenval);

% the ellipse in x and y coordinates

ellipse_x_r = a*cos( theta_grid );

ellipse_y_r = b*sin( theta_grid );

%Define a rotation matrix

R = [ cos(phi) sin(phi); -sin(phi) cos(phi) ];

%let's rotate the ellipse to some angle phi

r_ellipse = [ellipse_x_r;ellipse_y_r]' * R;

% Draw the error ellipse

plot(r_ellipse(:,1) + X0,r_ellipse(:,2) + Y0,'-')

hold on;

% Plot the original data

plot(data(:,1), data(:,2), '.');

mindata = min(min(data));

maxdata = max(max(data));

Xlim([mindata-3, maxdata+3]);

Ylim([mindata-3, maxdata+3]);

hold on;

% Plot the eigenvectors

quiver(X0, Y0, largest_eigenvec(1)*sqrt(largest_eigenval), largest_eigenvec(2)*sqrt(largest_eigenval), '-m', 'LineWidth',2);

quiver(X0, Y0, smallest_eigenvec(1)*sqrt(smallest_eigenval), smallest_eigenvec(2)*sqrt(smallest_eigenval), '-g', 'LineWidth',2);

hold on;

% Set the axis labels

hXLabel = xlabel('x');

hYLabel = ylabel('y');

https://www.xarg.org/2018/04/how-to-plot-a-covariance-error-ellipse/