% my Data
Ith_MinT = rand(1,100) ;
Can_MinT = rand(1,100) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Calculate moments of the data sets %%%
% Observed
Ith_mean = mean(Ith_MinT); % x mean
Can_mean= mean(Can_MinT); % y mean
CV = cov(Ith_MinT,Can_MinT); % the angle of ellipse is determined by the cov of data
[Evec, Eval]=eig(CV); % Eigen values and vectors of covariance matrix
%%% Plot observed multivariate contours %%
% Observed data
xCenter = Ith_mean; % ellipses centered at sample averages, center at x,y data mean
yCenter = Can_mean;
theta = 0 : 0.01 : 2*pi; % angles used for plotting ellipses 0:360 degree
% compute angle for rotation of ellipse
% rotation angle will be angle between x axis and first eigenvector
x_vec= [1 ; 0]; % vector along x-axis
cosrotation =dot(x_vec,Evec(:,1))/(norm(x_vec)*norm(Evec(:,1)));
rotation =pi/2-acos(cosrotation); % rotation angle
R = [sin(rotation) cos(rotation); ...
-cos(rotation) sin(rotation)]; % create a rotation matrix
% create chi squared vector
chisq = [1.368 3.2188 4.605 5.991 9.210]; % percentiles of chi^2 dist df=2% the scale of the ellipsoid depends on confidence interval chosen
% chisq = [1.368 3.2188 4.605 5.991];
% size ellipses for each quantile
for i = 1:length(chisq)
% calculate the radius of the ellipse
xRadius(i)=(chisq(i)*Eval(1,1))^.5; % primary
yRadius(i)=(chisq(i)*Eval(2,2))^.5; % secondary
% lines for plotting ellipse
x{i} = xRadius(i)* cos(theta);
y{i} = yRadius(i) * sin(theta);
% rotate ellipse
rotated_Coords{i} = R*[x{i} ; y{i}];
% center ellipse
x_plot{i}=rotated_Coords{i}(1,:)'+xCenter;
y_plot{i}=rotated_Coords{i}(2,:)'+yCenter;
end
% Set up plot
figure
set(gca,'Title',text('String','Probability Ellipses of M1 full weighted of patient plan',...
'FontAngle', 'italic', 'FontWeight', 'bold'),...
'xlabel',text('String', '$\mathbf{M1 full weighted bladder}$', 'Interpreter', 'latex'),...
'ylabel',text('String', '$\mathbf{M1 full weighted rectum}$', 'Interpreter', 'latex'))
hold on
% Plot data points
plot(Ith_MinT,Can_MinT,'o');% scattered data plot
% Plot contours
for j = 1:length(chisq)
P = [x_plot{j} y_plot{j}] ;
plot(x_plot{j},y_plot{j})
% Get center
C = mean(P) ;
plot(C(1),C(2),'*')
% Get distances
d = sqrt(sum((C-[P(:,1) P(:,2)]).^2,2)) ;
M1 = [C(1)-min(d),C(2) ; C(1)+min(d),C(2)] ;
plot(M1(:,1),M1(:,2),'r')
M2 = [C(1),C(2)-max(d) ; C(1) ,C(2)+max(d)] ;
plot(M2(:,1),M2(:,2),'r')
end
legend('Data points','50%', '80%', '90%', '95%', '99%')
hold on
