Need help finding out if I wrote this code correctly

N = 365;
n = (1:80);
for M = 1:1000
    j = (1:n);
    pn = 1-(1-((n-1)/N));
end
%%
fig = figure(1); clf
plot(pn)
xlabel('Value')
ylabel('Probability')
ylim(x,0.3)

What I'm trying to do is perform an experiment on (pn) 1000 times, with a value (n) from 1 through 80. I'm supposed to be getting 0.5 at a value greater than 22 but I'm not getting that. First time using matlab so help would be very much appreciated.

[Kevin Nelson]

What I'm trying to do is perform an experiment on (pn) 1000 times, with a value (n) from 1 through 80. I'm supposed to be getting 0.5 at a value greater than 22 but I'm not getting that. First time using matlab so help would be very much appreciated.

Answer 1

[Star Strider]

If you want to save the values of ‘pn’ in the loop, subscript it —

pn(M) = 1-(1-((n-1)/N));

That will create a row vector that you can then plot.

What is the ‘j’ calculation supposed to do? The code calculates it in each iteration and never uses it (at least not that I can see).

[Walter Roberson]

n is a vector, so you would have to assign to either pn(M,:) or pn(:,M)

[Star Strider]

I just now noticed that.

EDIT — (5 Sep 2022 at 13:16)

I believe this is the ‘Birthday Problem’ the solution of which is simply:

N = 365;
V = 1:N;
pn = cumprod((N:-1:1) / N);
Result = [V; pn]

Result = 2×365

1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 16.0000 17.0000 18.0000 19.0000 20.0000 21.0000 22.0000 23.0000 24.0000 25.0000 26.0000 27.0000 28.0000 29.0000 30.0000 1.0000 0.9973 0.9918 0.9836 0.9729 0.9595 0.9438 0.9257 0.9054 0.8831 0.8589 0.8330 0.8056 0.7769 0.7471 0.7164 0.6850 0.6531 0.6209 0.5886 0.5563 0.5243 0.4927 0.4617 0.4313 0.4018 0.3731 0.3455 0.3190 0.2937

Eq = find(pn >= 0.5, 1, 'last')

Eq = 22

I don’t understand the reason that it needs to be simulated 1000 times unless some different sort of calculation is involved.

.

Answer 2

[Walter Roberson]

n = (1:80);

n is a vector.

for M = 1:1000
    j = (1:n);

n is a vector, so j is 1:VECTOR instead of being 1:SCALAR or 1:SCALAR:SCALAR . You have a logic error there. What MATLAB does in this particular case is to define the output as if you had written 1:VECTOR(1) . Since n = 1:80 then n(1) is 1, so you are effectively defining j = 1:1 which is the same as j = 1

    pn = 1-(1-((n-1)/N));

but you never use j in the calculation or afterwards, so there is no point calculating it.

You overwrite all of pn in the calculation. The end result is the same if you had only done the last iteration, when M = 1000.

You did not include any comments, so we cannot tell what you are trying to calculate, so we cannot suggest what correct code might be. (But I get the impression that possibly you are testing out the "Birthday Paradox" ?)

[Kevin Nelson]

I included a comment on my question, just below the code. But I see what you mean. What I was trying to do was have (pn) be calculated for each value of (n), then repeat and do the same calculation for the next value. I realize not (j) wasn't included in the calculation, but what is the correct way to do what I just mentioned? I should be getting an answer for each value of (n), and when I run the code, I do get values for each one, but they aren't correct because I should be getting a 0.5 for values > 22. Do you see what I mean? Will just removing (j) give me that result?

[Walter Roberson]

With your n being a vector, your line pn = 1-(1-((n-1)/N)) is already calculating pn(K) = 1-(1-((n(K)-1)/N)) for all K 1 to number of elements of n .

However, your for M = 1:1000 is causing the exact same calculation to repeat 1000 times, calculating exactly the same thing each time.

You are not using random numbers, so repeating the "experiment" 1000 times is of no benefit. A calculation with a fixed outcome is not typically considered as an "experiment", just a "calculation".

[Kevin Nelson]

Do I set M = rand to have it calculate with random numbers? How do I specifically let it calculate each experiment 1000 times?

[Walter Roberson]

It is not very meaningful to talk about probability without random numbers, or at least without cdf or pdf formulas.

1-(1-((n-1)/N)) 

1-(n-1)/N is 1/N

1 - 1/N is (n-1)/N

Therefore your current formula is calculating the same thing as (n-1)/N

If you were doing something like birthday paradox, you would be doing a calculation similar to

1-(1-((n-1)/N).^POWER)

for appropriate POWER -- if, that is, you were calculating based on the closed-form solution, rather than doing it as an "experiment".