close all;
clear all;
clc;

%% Parameters

% Dimension
D = 2; 
% Number of samples
N = 1000; 
% Offset
Offset = [2,7];

%% Sample
% Samples, each column is a random variable
X = randn(N,D); 
% Offset
X = X + Offset ;

%% Weight

% Arbitrary weights, e.g. Gaussian, unnormalised
%w = exp(-([1:N]'/N-0.5).^2); 
%w = max(X(:,1),0.1);
mu = 0.7;
w = 0.01+(1/(mu*sqrt(2*pi))) * exp(-0.5*((X(:,1)/mu).*(X(:,1)/mu)))


%% For plotting the gaussian used for calculating the weights
%x=[-2:0.1:6];
%y = 0.01+(1/(mu*sqrt(2*pi))) * exp(-0.5*((x/mu).*(x/mu)))






%% Compute non-weighted covariance

% row containing mean of each column
C = sum(X) / N;
% Covariance matrix
sigma = cov(X-C);
% Singular value decomposition
[U,S,D]=svd(sigma)




%% Compute weighted covariance

% row containing weighted mean of each column
C_w = sum(w.*X) / sum(w) ;

%% Covariance matrix
% Biased weighted sample covariance:
sigma_w =  (w.*(X-C_w))'*(X-C_w) / sum(w);
% Singular value decomposition
[U_w,S_w,D_w]=svd(sigma_w)


%% Plots samples

figure;
hold on;
%yyaxis left
for i=1:N
    % Non-weighted samples
    plot (X(i,1), X(i,2), '.r','MarkerFaceColor', 'r');
    % Weighted samples
    %scatter(X(i,1), X(i,2), 300*w(i), 'ob', 'MarkerEdgeColor', 'w', 'MarkerFaceColor', 'b', 'MarkerFaceAlpha',0.5);
    h=plot (X(i,1), X(i,2), 'ob', 'MarkerSize', 10*w(i));
    
end
alpha(.5)
axis square equal;    

%% Plot non-weighted center of mass
plot (C(1), C(2), 'xr', 'MarkerSize', 50, 'LineWidth', 10);
plot (C(1), C(2), 'or', 'MarkerSize', 50, 'LineWidth', 10);


%% Plot weighted center of mass
plot (C_w(1), C_w(2), 'xb', 'MarkerSize', 50, 'LineWidth', 10);
plot (C_w(1), C_w(2), 'ob', 'MarkerSize', 50, 'LineWidth', 10);



%% Plot non-weighted PCA

a=[0:0.1:2*pi];
Xu=[cos(a);sin(a)];
Xu= (U*sqrt(S)) * Xu;
patch (C(1) + 3*Xu(1,:),C(2) + 3*Xu(2,:),'r','FaceAlpha',0.1);

%% Plot weighted PCA

Xu_w=[cos(a);sin(a)];
Xu_w= (U_w*sqrt(S_w)) * Xu_w;
patch (C_w(1) + 3*Xu_w(1,:),C_w(2) + 3*Xu_w(2,:),'b','FaceAlpha',0.1);


%% Plot the gaussian used for calculating the weights
%yyaxis right
%area (x,y, 'FaceAlpha', 0.2);
%grid on;