Журнал участка (n над k)

Question

Я никогда не использовал Matlab прежде, и я действительно не знаю, как исправить код. Мне нужно построить лог (1000 на 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 обычно использует нецелые значения в указанных пределах, поэтому nchoosek выдает ошибку.

Вы можете решить эти три проблемы, используя взаимосвязь между факториалом и гамма-функцией и тот факт, что Matlab имеет gammaln, который непосредственно вычисляет логарифм гамма-функции :

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

Обратите внимание, что вы получаете график со значениями y для всех x в указанном диапазоне, но фактически биномиальный коэффициент определяется только для неотрицательных целых чисел.

Answer 2

syms тип функции воспроизводит именно то, что вы хотите

syms x

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

Answer 3

ОК, так как вы все равно получили спойлеры для домашнего задания, я опубликую ответ, который, как мне кажется, легче понять.

Мультипликативная формула для биномиального коэффициента говорит, что

n над k = произведение _{i = от 1 до k} ((n + 1-i) / i)

(извинитеНельзя написать правильные формулы для SO, смотрите ссылку на Википедию, если это неясно).

Чтобы вычислить логарифм произведения, мы можем вычислить сумму логарифмов:

log (product (x _i )) = сумма (log (x _i ))

Таким образом, мы можем вычислить значения (n+1-i)/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, и здесь не то, что вам нужно.

Answer 4

Это решение использует arrayfun, чтобы справиться с тем фактом, что для nchoosek(n,k) требуется k, чтобы быть скаляром.Этот подход не требует наборов инструментов .

Кроме того, здесь используется plot вместо fplot, поскольку этот умный ответ уже посвящен тому, как поступить с 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, который является самой большой поддерживаемой реальной плавающей точкой.См. здесь для получения дополнительной информации.