绘制日志（n 超过 k）

Question

我以前从未使用过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)

Answer 1

代码有几个问题：

nchoosek 不对第二个输入进行向量化，即不接受数组作为输入。 fplot 对向量化函数工作得更快。否则可以使用，但会发出警告。
对于如此大的第一个输入值，nchoosek 的结果接近溢出。例如，nchoosek(1000,500)给出2.702882409454366e+299，并发出警告。
nchoosek 需要整数输入。 fplot 通常使用指定限制内的 non-integer 值，因此 nchoosek 会发出错误。

您可以利用阶乘与 gamma function and the fact that Matlab has gammaln 之间的关系解决这三个问题，它直接计算伽马函数的对数：

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

请注意，对于指定范围内的所有 x，您会得到一个带有 y 值的图，但实际上二项式系数仅定义对于 non-negative 个整数。

Answer 2

此解决方案使用 arrayfun to deal with the fact that nchoosek(n,k) 要求 k 为标量。这种方法 不需要工具箱 。

此外，这使用 plot instead of fplot since 已经解决了如何使用 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)

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

Warning: Result may not be exact. Coefficient is greater than 9.007199e+15 and is only accurate to 15 digits

因为nchoosek(1000,500) = 2.7029e+299。正如 @Luis Mendo, this is due to realmax = 1.7977e+308 which is the largest real floating-point supported. See here 所指出的以获取更多信息。

Answer 3

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

multiplicative formula for the binomial coefficient 表示

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

（抱歉，无法在 SO 上编写正确的公式，如果不清楚，请参阅维基百科 link）。

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

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)

另请注意 ./，element-wise 划分。 / 在 MATLAB 中执行矩阵除法，这里不是您需要的。

Answer 4

syms 函数类型完全重现你想要的

syms x

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

绘制日志（n 超过 k）

Plot log(n over k)

matlab

plot

logarithm

binomial-coefficients