clear all;

% Parameters
N = 100000;
S = 25;

% Generate 200 numbers grid from -1 to 1
x_values = linspace(-1, 1, 200);

% Generate N random numbers using the given function
x_k = rand(1, N) + rand(1, N) - 1;

% Initialize lists
ai_Mean_List = zeros(1, N);
f_dash_List = zeros(1, length(x_values));

% Generate ai_Mean_List
for s = 1:N
    ai_List = fi(s, x_k);
    ai_Mean_List(s) = mean(ai_List);
end

% Calculate f_dash_List
for x = 1:length(x_values)
    for i = 1:S
        f_dash_List(x) = f_dash_List(x) + ai_Mean_List(i) * fi(i, x_values(x));
    end
end

% Calculate theoretical PDF values
pdf_values = max(1 - abs(x_values), 0);
pdf_values = pdf_values / trapz(x_values, pdf_values);

% Plot the theoretical and estimated PDFs
figure;
plot(x_values, pdf_values, 'g', 'LineWidth', 4, 'DisplayName', 'Theoretical PDF');
hold on;
plot(x_values, f_dash_List, 'DisplayName', 'Estimated PDF');
grid on;
legend('show');
title(['Estimated and Theoretical PDFs for S=', num2str(S)]);
hold off;

% Function Definition
function fi_result = fi(i, x)
    if i == 1
        fi_result = 1 / sqrt(2);
    elseif mod(i, 2) == 0
        fi_result = sin(i/2 * pi * x);
    else
        fi_result = cos((i-1)/2 * pi * x);
    end
end