function X = f_model_matrix(y, cluster)
    unique_clusters = unique(cluster);
    num_clusters = numel(unique_clusters);
    X = zeros(numel(y), num_clusters);

    for i = 2:num_clusters
        X(:, i - 1) = y .* (cluster == unique_clusters(i));
    end
end