Please check the Code and let me know the error. When i try to run it in MATLAB. It gives me an error .

%annual probability that the water height exceeds the threshold hcr.
cv=1
nsam=1000
while cv>0.1
  Q=lognrnd(5.1434,0.1416,nsam,1);
  Depth=0.2*Q.^0.4329;
  Depth1=Depth>4.7;
  if any(Depth1)
    Pf=sum(Depth1)/nsam
    cv=sqrt(1/Pf/nsam);
    nsam=round(1/Pf/(0.1^2));
  end
end

[SL: indented the while and if blocks to emphasize them.]

[Walter Roberson]

How long should we let it run? It has been more than a minute of computation now but no error message or output has been produced. What error message are you seeing?

Answer 1

[Walter Roberson]

It is low probability that even a single element of Q will meet the 0.2*Q.^0.4329 > 4.7 test.

But suppose exactly one element does. Then with nsam = 1000 at first, we have Pf=sum(Depth1)/nsam which would be Pf = 1/1000 . Then cv=sqrt(1/Pf/nsam); would be cv = sqrt(1/(1/1000)/1000) which would be cv = 1 (no change from what it was.) Then nsam=round(1/Pf/(0.1^2)) would be nsam = round(1/(1/1000)/(1/10^2)) which would give nsam = 100000 .

Then we spin for a long time waiting for another item to match. When one finally does, we would have Pf=sum(Depth)/nsam which would be Pf=1/10000 . Then cv=sqrt(1/(1/10000)/10000) would be cv = 1 (no change from what it was). THen nsam=round(1/(1/10000)/(0.1^2) which would take nsam to 10000000 . And so on, continually increasing nsam but getting no closer to reducing cv.

Finally after rather some time, having matched single items N times, nsam = 1000*10^(2*N), we get two matches in depth. Pf=sum(Depth1)/nsam which would be Pf = 2/10^(3+2*N), cv = sqrt(1/Pf/nsam) would be 1/(2/10^(3+2*N)/10^(3+2*N) which is going to give 1/2, and nsam would increase by another factor of 100. But does that help? No, because the next time that a single Depth match is made, cv will be reset back to 1 again.

So in order for cv to reach 1/10, there would have to be 10 or more simultaneous matches on depth out of nsam. What is the probability of that? True, it does increase as nsam increases, but what is our initial probability ?

With lognrnd(5.1434,0.1416,nsam,1)

[M,V] = lognstat(5.1434,0.1416)
M =
          173.023129137865
V =
          606.311812331174
sqrt(V)
ans =
        24.6233996907652

and our target is exp(log(4.7/0.2)/0.4329) == 1469.50114573441 so that is (1469.50114573441 - 173.023129137865) / 24.6233996907652 = 52.6522751885792 standard deviations. That is approximately 8E-605 . That is extremely unlikely to happen in your lifetime.

Answer 2

[John BG]

Hi Farooq

this is John BG ( <mailto:[e-mail] [e-mail]> )

1.

increasing Q to 5000

2.

why randomly sampling following a log normal distribution if can enfilade with a straight line?

3.

no need for if any(Depth1) because Depth1 already contains all the values you want to process

cv=1
nsam=1000
while cv>0.1

    Q=[1:1:5e3]

    Depth=0.2*Q.^0.4329;
    Depth1=Depth>4.7;

    Pf=sum(Depth1)/nsam
    cv=sqrt(1/Pf/nsam);
    nsam=round(1/Pf/(0.1^2));

end

Is this of any help?

Pf
= 3.5310

if you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?

To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link

thanks in advance

John BG

<mailto:[e-mail] [e-mail]>

[Walter Roberson]

John, that will not work. Suppose there is no match the first time, so sum(Depth1) is 0. Then Pf is 0, and 1/Pf is infinity. cv would become infinity. nsam would become infinity as well.

The if any() is important there is important to ensure that there are no divisions by 0 and that cv and nsam are only changed when there is at least one match on depth.