绘制日志（n 超过 k）答案

【问题标题】：Plot log(n over k)绘制日志（n 超过 k）
【发布时间】：2019-05-23 14:56:52
【问题描述】：

我以前从未使用过 Matlab，我真的不知道如何修复代码。我需要绘制 log(1000 over k)，k 从 1 到 1000。

y = @(x) log(nchoosek(1000,x));

fplot(y,[1 1000]);

错误：

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly
vectorize your function to return an output with the same size and shape as the input
arguments. 
In matlab.graphics.function.FunctionLine>getFunction
In matlab.graphics.function.FunctionLine/updateFunction
In matlab.graphics.function.FunctionLine/set.Function_I
In matlab.graphics.function.FunctionLine/set.Function
In matlab.graphics.function.FunctionLine
In fplot>singleFplot (line 241)
In fplot>@(f)singleFplot(cax,{f},limits,extraOpts,args) (line 196)
In fplot>vectorizeFplot (line 196)
In fplot (line 166)
In P1 (line 5)

【问题讨论】：

nchoosek 需要一个标量输入。您应该使用循环来分别计算每个k 的结果。但话又说回来，我认为你的练习是为了让你不要使用nchoosek，因为中间值是巨大的，可能不适合双精度浮点数（或者至少给出错误的结果）。 log(n over k) 的方程是什么？
这会帮助你tutorialspoint.com/matlab/matlab_plotting.htm
fitting and plotting a parabola, matlab的可能重复
根本不是重复的。这个问题使用fplot（它接受一个匿名函数作为输入），而另一个问题使用plot（它接受数字向量）。此外，由于 OP 的方法涉及大量数字，这个问题存在潜在的数值精度问题
@LuisMendo 考虑到他的数据都是离散的，对于他来说，在这个应用程序中使用 plot 可能比 fplot 更有意义。

标签： matlab plot logarithm binomial-coefficients

【解决方案1】：

代码有几个问题：

nchoosek 不对第二个输入进行矢量化，也就是说，它不接受数组作为输入。 fplot 对矢量化函数的工作速度更快。否则可以使用，但会发出警告。
对于如此大的第一个输入值，nchoosek 的结果接近溢出。例如，nchoosek(1000,500) 给出2.702882409454366e+299，并发出警告。
nchoosek 需要整数输入。 fplot 通常使用指定范围内的非整数值，因此 nchoosek 会发出错误。

您可以利用阶乘和gamma function 之间的关系以及Matlab 具有直接计算伽马函数对数的gammaln 的事实来解决这三个问题：

n = 1000;
y = @(x) gammaln(n+1)-gammaln(x+1)-gammaln(n-x+1);
fplot(y,[1 1000]);

请注意，您会得到一个图，其中包含指定范围内所有 x 的 y 值，但实际上二项式系数仅针对非负整数定义。

【讨论】：

对阶乘和伽马函数的出色而巧妙的利用 (+1)。另外，我不知道gammaln 函数。谢谢！

【解决方案2】：

好的，既然你现在已经为你的家庭作业做剧透了，我会发布一个我认为更容易理解的答案。

multiplicative formula for the binomial coefficient 这么说

n over k = product_{i=1 to k}( (n+1-i)/i )

（抱歉，无法在 SO 上编写正确的公式，如果不清楚，请参阅 Wikipedia 链接）。

要计算产品的对数，我们可以计算对数的和：

log(product(x_i)) = sum(log(x_i))

因此，我们可以计算所有i 的(n+1-i)/i 的值，取对数，然后将第一个k 值相加得到给定k 的结果。

这段代码使用cumsum，即累积和来实现这一点。它在数组元素k 的输出是从1 到k 的所有输入数组元素的总和。

n = 1000;
i = 1:1000;
f = (n+1-i)./i;
f = cumsum(log(f));
plot(i,f)

还要注意./，即按元素划分。 / 在 MATLAB 中执行矩阵除法，这里不需要。

【讨论】：

这很聪明。我没有考虑过 对数总和 方法，但自从你把它拼出来后就完美了。非常清晰，干净的方法。 (+1)
其他 SE 站点支持使用 MathJax 渲染 LaTex。这似乎也是 SO 上非常有用的功能。
@John：该建议已在 Meta 上多次讨论并被否决。不幸的是，他们不想这样做。
@CrisLuengo 是的，我刚刚进入了一个兔子洞调查这个问题。似乎出现了很多，我们只是生活在丑陋的排版和低效的沟通中，因为原因......

【解决方案3】：

syms 函数类型准确再现您想要的内容

syms x

y = log(nchoosek(1000,x));
fplot(y,[1 1000]);

【讨论】：

【解决方案4】：

这个解决方案使用arrayfun 来处理nchoosek(n,k) 要求k 是一个标量的事实。这种方法不需要工具箱。

此外，这使用plot 而不是fplot，因为this clever answer 已经解决了如何处理fplot。

% MATLAB R2017a
n = 1000;
fh=@(k) log(nchoosek(n,k));  
K = 1:1000;
V = arrayfun(fh,K);    % calls fh on each element of K and return all results in vector V

plot(K,V)

请注意，k 的某些值大于或等于 500，您将收到警告

警告：结果可能不准确。系数大于9.007199e+15，只精确到15位

因为nchoosek(1000,500) = 2.7029e+299。正如@Luis Mendo 所指出的，这是由于realmax = 1.7977e+308 是支持的最大实浮点数。请参阅here 了解更多信息。

【讨论】：