Getting NaN greater values in a function
조회 수: 4 (최근 30일)
이전 댓글 표시
Hey there, I am trying to plot a function, that works perfect for smaller values but as I input the larger values I get NaN as output. Below is the code for what coded:
clearvars;
iterations = [10, 100, 1000];
count = 1;
p = 0.5;
for n = iterations
figure();
hold on;
for k = 1:n
binomial = (factorial(n)) / (factorial(k) * factorial(n - k)) * (p^k * (1-p)^(n-k));
plot(k, binomial, '*');
end
end
As shown in the above graphs, I get the output for the iterations 10 &100, but there is no output for for 1000. Even if I replace 10 by 1000 I get the same output so I am sure its nothing related to the indexing.
댓글 수: 0
채택된 답변
Torsten
2022년 10월 7일
편집: Torsten
2022년 10월 7일
warning('off')
iterations = [10, 100, 1000];
count = 1;
p = 0.5;
for n = iterations
figure();
hold on;
for k = 1:n
%binomial = nchoosek(n,k) * p^k * (1-p)^(n-k);
binomial = exp(sum(log(n-(0:k-1)))-sum(log(1:k))) * p^k * (1-p)^(n-k);
plot(k, binomial, '*');
end
end
댓글 수: 2
Stephen23
2022년 10월 7일
편집: Stephen23
2022년 10월 7일
"I believe nchoosek function is same as ... so why does it give NaN?"
Because calculating using factorials will very quickly get values well beyond those which can be represented using binary floating point numbers. Which is why NCHOOSEK uses another algorithm (hint: TYPE).
추가 답변 (1개)
Steven Lord
2022년 10월 7일
What is the factorial of 1000 in double precision?
factorial(1000)
It overflows. If we computed it symbolically, using arbitrary precision arithmetic:
vpa(factorial(sym(1000)))
That's a lot larger than realmax so it is not surprising it overflows.
realmax
Most likely you're dividing Inf by either Inf (if the denominator overflows), NaN, or 0 (if the denominator underflows.)
Inf ./ [Inf, NaN, 0]
The plot function doesn't display NaN values.
If you use nchoosek as @Torsten suggested it avoids computing such large numbers and then dividing a large numerator by a large denominator.
참고 항목
카테고리
Help Center 및 File Exchange에서 Characters and Strings에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!