number_of_senders = 4;
number_of_timeslots = 1000;
access_probability = 0:0.01:1;

access_probability_length = length(access_probability);

sender_access = zeros(number_of_senders);

simulated_throughput = zeros(1, access_probability_length);
estimated_throughput = zeros(1, access_probability_length);

% simulate all time slots on each probability
for probability_idx = 1 : access_probability_length
    
    current_result = zeros(1, number_of_timeslots);
    probability = access_probability(probability_idx);

    for time_slot = 1 : number_of_timeslots
        % generate random sender accesses by probability
        sender_access = zeros(1, number_of_senders);
        for sender = 1 : number_of_senders
            if rand() < access_probability(probability_idx)
                sender_access(1, sender) = 1;
            end
        end
    
        current_result(1, time_slot) = sum(sender_access(1, :)) == 1;
    end
    
    simulated_throughput(probability_idx) = mean(current_result);
    estimated_throughput(probability_idx) = number_of_senders * probability * (1 - probability)^(number_of_senders - 1);
end

plot(access_probability, estimated_throughput)
hold on
plot(access_probability, simulated_throughput)
legend("theoretical", "simulation")
xlabel("access probability")
ylabel("throughput")