colors = lines(length(unique_z2) * length(unique_z3));
for i = 1:length(unique_z2)
for j = 1:length(unique_z3)
condition_idx = (z2 == unique_z2(i)) & (z3 == unique_z3(j));
Z4_condition = z4(condition_idx);
if ~isempty(Z4_condition)
[m, s] = normfit(Z4_condition);
x_values = min(Z4_condition):0.1:max(Z4_condition);
y_values = normpdf(x_values, m, s);
plot(x_values, y_values, 'Color', colors(color_idx,:), 'DisplayName', ...
['z2 = ', num2str(unique_z2(i)), ', z3 = ', num2str(unique_z3(j))], 'LineWidth', 1.5);
plot(Z4_condition, zeros(size(Z4_condition)), 'o', 'MarkerFaceColor', colors(color_idx,:), 'MarkerEdgeColor', 'k');
disp(['Condition z2 = ', num2str(unique_z2(i)), ', z3 = ', num2str(unique_z3(j))]);
disp(['Mean: ', num2str(m), ', Std Dev: ', num2str(s)]);
disp(['Z4 values: ', num2str(Z4_condition')]);
color_idx = color_idx + 1;
ylabel('Probability Density');
title('Gaussian Fit for Different Conditions (z2, z3)');