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)');