Циклы со сложными математическими выражениям, как правильно записать?!

@EvgenSanich · Регистрация: 17.07.2013

Студворк — интернет-сервис помощи студентам

Привет всем гуру Матлаб! Если у Вас найдется свободная минутка

, подскажите пожалуйста, как правильно сделать ниже указанные записи в функции matlab в соответствии с представленными кусочками алгоритмов на картинках ниже:

[COLOR="rgb(0, 100, 0)"]Картинка 1[/COLOR]
Не могу понять, как присвоить значение P() в соответствии представленным выражением в условии, я сделал так:

Matlab M
1
2
3
4
5
6
7
8
9
10
if jj == 0
   disp('Вывод P[0]:')
   fprintf('%.10f\n', P(1))
   else 
        for kk = 1:1:fix(a) - 1
             P(kk) = ((kk + jj - 1)/kk) * p^kk * (1 - p)^jj * P(kk)* e^-((kk + jj)* r * (tau/n));
             disp('Вывод  P(kk) = ((kk + jj - 1)/kk) * p^k * (1 - p)^jj * P(kk)* e^-((kk + jj)* r * (maturity/n):')
             fprintf('%.10f\n', P(kk))
        end
 end

скорее всего в скобках ((kk + jj - 1)/kk) не знак деления, тогда не пойму как правильно записать?

Картинка 2
Аналогичная ситуация, возможно я не правильно записал, сделал так:

Matlab M
1
2
3
4
5
6
if rem((n - h),2)||((n - h)/2) > 0 || ((n - h)/2) == 0 || ((n - h)/2) < n || ((n - h)/2) == n 
    disp('Удар! барьер достигнут')
    C((n-h)/2) = 0; 
else
    disp('Барьер не достигнут')
end

ругается на эту запись вот так:

Attempted to access C(50.8548); index must be a positive
integer or logical.

Error in downandoutnabezdivid (line 23)
C((n-h)/2) = 0;

изначально в начале функции я задал массив C = zeros(1, n)

Подскажите пожалуйста, в чем возможны ошибки?

@Heidegger · 12.09.2013, 15:54

Сообщение от EvgenSanich

как правильно сделать ниже указанные записи в функции matlab в соответствии с представленными кусочками алгоритмов на картинках ниже

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

Сообщение от EvgenSanich

скорее всего в скобках ((kk + jj - 1)/kk) не знак деления

Да, это не знак деления, а биномиальный коэффициент:

$https://www.cyberforum.ru/cgi-bin/latex.cgi?C_{n}^{k}=\binom{n}{k}=\frac{n!}{k!(n-k)!}$

В Matlab ему соответствует функция nchoosek(n, k).

Во втором случае можно записать так:

Matlab M
1
2
3
4
5
6
7
8
9
if mod(n-h, 2) == 0
    q = (n-h)/2;
    if (q >= 0) && (q <= n)
        disp('Удар! барьер достигнут')
        C(q) = 0; 
    else
        disp('Барьер не достигнут')
    end
end

@EvgenSanich · 19.09.2013, 11:02 **[ТС]**

Доброго Вам дня! Алгоритм с которым я не перестаю бороться, и очень не хочу превращать все в спам, поэтому еще один последний призыв к помощи. Сам алгоритм на рисунке ниже, моя реализация представлена тут:

Данные для ввода:

Matlab M
1
2
3
4
5
6
7
8
S = 50 
X = 50 
n = 100
q = 1
tau = 1/3
r = 8
u = 1.5
d = 0.5

Алгоритм:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
ffunction P = somefunction(S, u, d, tau, X, n, r, q)
qnum = q/100;
P = zeros(1, n);
R = exp(r/100);
p = exp((((r-qnum)*(tau/n)) - d)/(u/d)); %для без дивидендной акции
a = log(X/(S * d^n))/log(u/d);
if r > q || r == q || q > 0
    for ii = 1:1:fix(a) - 1
        P(ii) = -1; %верхняя граница узлов досрочного исполнения
    end
    for jj = n - fix(a):-1:1
        for ii = fix(a) - 1:-1:1
            niz(ii) = P(ii);
            if  niz(ii) < 0
                niz(ii) = X - S * u^ii * d^(jj + 1);
                c(ii) = (p * P(ii + 1) + (1 - p) * niz(ii))/R;
            else
                if c(ii) > (X - S * u^ii * d^jj) && c(ii) == (X - S * u^ii * d^jj) && ii > 0
                   P(ii) = c; %досрочное исполнение отсутствует
                else
                    if c(ii) > (X - S * u^ii * d^jj) && c(ii) == (X - S * u^ii * d^jj) && ii == 0
                       P(ii) = c;
                    else
                        if jj == 0
                           disp('Вывод P[0]:')
                           fprintf('%.10f\n', P(1))
                        else 
                            for kk = 1:1:fix(a) - 1
                                %ранее было так (не правильно)
                                %P(kk) = ((kk + jj - 1)/kk) * p^kk * (1 - p)^jj * P(kk)* e^-((kk + jj)* r * (tau/n));
                                
                                %Биномиальный коэффициент: Функция MATLAB nchoosek(n, k)
                                V(kk) = (kk + jj - 1)* p^kk * (1 - p)^jj * P(kk)* e^-((kk + jj)* r * (tau/n));
                                M(kk) = (kk)*p^kk * (1 - p)^jj * P(kk)* e^-((kk + jj)* r * (tau/n));
                                P = nchoosek(V, M);
                                disp('Вывод  P(kk) = ((kk + jj - 1) kk) * p^k * (1 - p)^jj * P(kk)* e^-((kk + jj)* r * (maturity/n):')
                                fprintf('%.10f\n', P(kk))
                            end
                        end
                    end
                end
            end 
        end
        if P(1) < 0
           P(1) = X - S;
           disp('Вывод P[0]:')
           fprintf('%.10f\n', P(1))
        end
    end
else
   disp('Ошибка условия: r >= q >0') 
end
end

Собственно вопрос, или я не правильно записал указанный алгоритм, или в книжке из которой я его взял что то не так, или???

Помогите разобраться, очень нужно

@Heidegger · 19.09.2013, 11:13

Сообщение от EvgenSanich

или я не правильно записал указанный алгоритм, или в книжке из которой я его взял что то не так, или???

А с чем возникают вопросы? Желательно сразу конкретизировать.

@EvgenSanich · 19.09.2013, 12:23 **[ТС]**

Ну первое это правильно ли интерпретировал запись из книжки в язык матлаб, и второе это ошибка которая выползает по реализации алгоритма:

Последний изм. алгоритма:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
function P = somefunction(S, u, d, tau, X, n, r, q)
qnum = q/100;
P = zeros(1, n);
R = exp(r/100);
p = exp((((r-qnum)*(tau/n)) - d)/(u/d)); %для без дивидендной акции
a = log(X/(S * d^n))/log(u/d);
if r > q || r == q || q > 0
    for ii = 1:1:fix(a) - 1
        P(ii) = -1; %верхняя граница узлов досрочного исполнения
    end
    for jj = n - fix(a):-1:1
        for ii = fix(a) - 1:-1:1
            niz(ii) = P(ii);
            if  niz(ii) < 0
                niz(ii) = X - S * u^ii * d^(jj + 1);
                c(ii) = (p * P(ii + 1) + (1 - p) * niz(ii))/R;
            else
                if c(ii) > (X - S * u^ii * d^jj) && c(ii) == (X - S * u^ii * d^jj) && ii > 0
                   P(ii) = c; %досрочное исполнение отсутствует
                else
                    if c(ii) > (X - S * u^ii * d^jj) && c(ii) == (X - S * u^ii * d^jj) && ii == 0
                       P(ii) = c;
                    else
                        if jj == 0
                           disp('Вывод P[0]:')
                           fprintf('%.10f\n', P(1))
                        else 
                            for kk = 1:1:fix(a) - 1
                                %ранее было так (не правильно)
                                %P(kk) = ((kk + jj - 1)/kk) * p^kk * (1 - p)^jj * P(kk)* e^-((kk + jj)* r * (tau/n));
                                
                                %Биномиальный коэффициент: Функция MATLAB nchoosek(n, k)
                                V(kk) = (kk + jj - 1)* p^kk * (1 - p)^jj * P(kk)* exp(-((kk + jj)* r * (tau/n)));
                                M(kk) = (kk)*p^kk * (1 - p)^jj * P(kk)* exp(-((kk + jj)* r * (tau/n)));
                                P(kk) = nchoosek(V, M);
                                disp('Вывод  P(kk) = ((kk + jj - 1) kk) * p^k * (1 - p)^jj * P(kk)* e^-((kk + jj)* r * (maturity/n):')
                                fprintf('%.10f\n', P(kk))
                            end
                        end
                    end
                end
            end 
        end
        if P(1) < 0
           P(1) = X - S;
           disp('Вывод P[0]:')
           fprintf('%.10f\n', P(1))
        end
    end
else
   disp('Ошибка условия: r >= q >0') 
end

А ошибка такая:

Вывод P[0]:
0.0000000000
Вывод P(kk) = ((kk + jj - 1) kk) * p^k * (1 - p)^jj * P(kk)* e^-((kk + jj)* r * (maturity/n):
1.0000000000

Error using nchoosek (line 25)
The second input has to be a non-negative integer.

Error in somefunction(line 42)
P(kk) = nchoosek(V,
M);

То есть я как то не правильно интерпретирую запись бинома?

Добавлено через 49 минут

Сообщение от EvgenSanich

То есть я как то не правильно интерпретирую запись бинома?

это строчка 36 алгоритма указанного в посте выше

@Heidegger · 19.09.2013, 12:31

Алгоритм сложный, ошибок и неточностей очень много. Прежде всего, нужно переформулировать исходный алгоритм в терминах принятой в Matlab one-based indexing (нумерации элементов массива, начиная с 1). Экспонента исправлена верно, биномиальный коэффициент записан действительно не там и не так. Достаточно подставить nchoosek(k+j-1, k) в нужном месте формулы (с учётом, опять же, one-based indexing).

p.s. Matlab в качестве платформы обязателен? В Mathcad реализация, думаю, будет проще.

@EvgenSanich · 19.09.2013, 12:42 **[ТС]**

Да, Matlab обязателен, требование в институте. Не совсем понял формулировку Matlab one-based indexing, в двух словах как говорят "на пальцах" можете пояснить что это за правило? я попробую разобраться с алгоритмом еще раз. А вот про бином, меня смущает выражение которое по книге идет как (k+j-1, k), а далее умножение на выражение p^kk * (1 - p)^jj * P(kk)* exp(-((kk + jj)* r * (tau/n))). То есть надо и (k+j-1) умножить на это выражение и k?

@Heidegger · 19.09.2013, 13:30

Сообщение от EvgenSanich

Не совсем понял формулировку Matlab one-based indexing

Сообщение от Heidegger

нумерации элементов массива, начиная с 1

В Matlab первый элемент массива всегда имеет индекс, равный 1. В алгоритме используется индексация элементов массива, начиная с 0. Необходимо внести соответствующие изменения.

Сам алгоритм немного путанный. Например, что такое C, R и D?

примерно так

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
function Y = hitf3(S, u, d, tau, X, n, r, q)
if ~((q > 0) && (r >= q))
    disp('Ошибка условия: r >= q >0')
    return
end
P = zeros(1, n+1);
p = (exp((r-q)*(tau/n))-d)/(u-d); %для бездивидендной акции
a = ceil(log(X/(S*d^n))/log(u/d));
P(1:a) = -1; %верхняя граница узлов досрочного исполнения
for j = n-a: -1: 0
for i = a-1: -1: 0
    niz = P(i+1);
    if  niz < 0
        niz = X - S*u^i*d^(j+1);
    end
    c = (p*P(i+2)+(1-p)*niz)/r;
    if (c >= X-S*u^i*d^j) && (i > 0)
        P(i+1) = c; %досрочное исполнение отсутствует
    elseif (c >= X-S*u^i*d^j) && (i == 0)
        P(i+1) = c;
        if j == 0
            Y = P(1);
            return
        else
            Y = 0;
            for k = 0:a-1
                Y = Y + nchoosek(k+j-1, k)*p^k*(1-p)^j* ...
                    P(k+1)*exp(-(k+j)*r*(tau/n));
            end
            return
        end
    else 
        break %узел раннего исполнения
    end
end
end
if P(1) < 0
    P(1) = X - S;
end
Y = P(1);

@EvgenSanich · 19.09.2013, 14:00 **[ТС]**

Это очередной биномиальный алгоритм из финансовой математики по оценке опционов, в котором я пытаюсь разобраться. В данном примере С, D и R это элементы формул (C - колл, D - диагональ, R - cтавка).

Спасибо Вам огромное за подсказки, буду разбираться дальше.

@Heidegger · 19.09.2013, 14:08

Сообщение от EvgenSanich

В данном примере С, D и R это элементы формул (C - колл, D - диагональ, R - cтавка).

Среди задаваемых и вычисляемых параметров нет С, D и R (заглавными буквами), при этом C и R участвуют в вычислениях (заменил на c и r соответственно).

@EvgenSanich · 19.09.2013, 17:22 **[ТС]**

Удалось побороть вышеуказанные алгоритм и опять с толкнулся с проблемой, но на этот раз по второй части вышеуказанного кусочка, а именно попытка реализовать алгоритм представленный на картинке ниже привела к нулевым результатам. Не могу понять в чем хитрость?

Единственное, какие изменения я в нес со своей стороны, это вторая переменная P за исключением небольших изменений в условия в if для нее и расчета. По идеи результаты полученного алгоритма не должны быть нулевыми, это как минимум.

Входные данные произвольные, я брал из примера в учебнике:

S = 95
X = 100
u = 1.5
d = 0.5
r = 10
n = 5
H = 90

Мой алгоритм:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
%Алгоритм БМОО колл и пут типа knock-out   
function [C, P] = knockoutoptions(S, u, d, X, H, n, r)
C = zeros(1, n+1);
P = zeros(1, n+1);
R = exp(r/100);
p = (R-d)/(u-d);
h = ceil(log(H/S)/log(u));
H = S*u^h; %Поскольку уровень H может не совпадать с реальной ценой 
                 %акции, то мы понижаем его до величины S*u^h, т.е. до
                 %максимальной цены акции, не превосходящей H.
 
if ((H<X)&&(H<S)) %условие для колл (вниз-и-выход)
    for ii = 1:1:n
        C(ii+1) = max(0, S*u^(n-ii)*d^ii - X);
    end
    if mod(n - h, 2) == 0
       q = (n-h)/2;
       if (q>=0)||(q<=n)     
           C(q) = 0; 
           disp('Удар! барьер достигнут')
       end
    end
    for jj = n-1:-1:0
        for ii = 1:1:jj
            C(ii+1) = (p*C(ii+1)+(1-p)*C(ii+1))/R;
        end
        if mod(jj-h, 2) == 0
           q = (jj-h)/2;
           if (q>=0)||(q<=jj)
               C(q) = 0;
               disp('Удар! барьер достигнут')
           end
        end
    end
    %fprintf('%.5f\n', C(1))   
 
elseif((H>X)&&(H>S)) %условие для опицона пут (вверх-и-выход)
    for ii = 1:1:n
        P(ii+1) = max(0, X - S*u^(n-ii)*d^ii);
    end
    if mod(n - h, 2) == 0
       q = (n-h)/2;
       if (q>=0)||(q<=n)
           P(q) = 0;
           disp('Удар! барьер достигнут')
       end
    end
    for jj = n-1:-1:0
        for ii = 1:1:jj
            P(ii+1) = (p*P(ii+1)+(1-p)*P(ii+1))/R;
        end
        if mod(jj-h, 2) == 0
           q = (jj-h)/2;
           if (q>=0)||(q<=jj)
               P(q) = 0;
               disp('Удар! барьер достигнут')
           end
        end
    end
    %fprintf('%.5f\n', P(1))  
end
fprintf('%.5f\n', C(1), P(1)) 
end

в чем может быть ошибка?

@Heidegger · 20.09.2013, 19:08

Сообщение от EvgenSanich

в чем может быть ошибка?

А в чём конкретно заключается ошибка? Если результат не соответствует ожиданиям, в проблемных местах ставят breakpoints (F12) и затем пошагово проходят программу (F10), контролируя значения переменных.

Новые блоги и статьи Все статьи Все блоги /
SDL3 для Web (WebAssembly): Реализация движения на Box2D v3 - трение и коллизии с повёрнутыми стенами 8Observer8 20.02.2026 Содержание блога Box2D позволяет легко создать главного героя, который не проходит сквозь стены и перемещается с заданным трением о препятствия, которые можно располагать под углом, как верхнее. . .	Конвертировать закладки radiotray-ng в m3u-плейлист damix 19.02.2026 Это можно сделать скриптом для PowerShell. Использование . \СonvertRadiotrayToM3U. ps1 <path_to_bookmarks. json> Рядом с файлом bookmarks. json появится файл bookmarks. m3u с результатом. # Check if. . .	Семь CDC на одном интерфейсе: 5 U[S]ARTов, 1 CAN и 1 SSI Eddy_Em 18.02.2026 Постепенно допиливаю свою "многоинтерфейсную плату". Выглядит вот так: https:/ / www. cyberforum. ru/ blog_attachment. php?attachmentid=11617&stc=1&d=1771445347 Основана на STM32F303RBT6. На борту пять. . .	Камера Toupcam IUA500KMA Eddy_Em 12.02.2026 Т. к. у всяких "хикроботов" слишком уж мелкий пиксель, для подсмотра в ESPriF они вообще плохо годятся: уже 14 величину можно рассмотреть еле-еле лишь на экспозициях под 3 секунды (а то и больше),. . .
И ясному Солнцу zbw 12.02.2026 И ясному Солнцу, и светлой Луне. В мире покоя нет и люди не могут жить в тишине. А жить им немного лет.	«Знание-Сила» zbw 12.02.2026 «Знание-Сила» «Время-Деньги» «Деньги -Пуля»	SDL3 для Web (WebAssembly): Подключение Box2D v3, физика и отрисовка коллайдеров 8Observer8 12.02.2026 Содержание блога Box2D - это библиотека для 2D физики для анимаций и игр. С её помощью можно определять были ли коллизии между конкретными объектами и вызывать обработчики событий столкновения. . . .	SDL3 для Web (WebAssembly): Загрузка PNG с прозрачным фоном с помощью SDL_LoadPNG (без SDL3_image) 8Observer8 11.02.2026 Содержание блога Библиотека SDL3 содержит встроенные инструменты для базовой работы с изображениями - без использования библиотеки SDL3_image. Пошагово создадим проект для загрузки изображения. . .

@Heidegger 617 / 242 / 16 Регистрация: 31.07.2013 Сообщений: 376
	19.09.2013, 12:31
	Алгоритм сложный, ошибок и неточностей очень много. Прежде всего, нужно переформулировать исходный алгоритм в терминах принятой в Matlab one-based indexing (нумерации элементов массива, начиная с 1). Экспонента исправлена верно, биномиальный коэффициент записан действительно не там и не так. Достаточно подставить nchoosek(k+j-1, k) в нужном месте формулы (с учётом, опять же, one-based indexing). p.s. Matlab в качестве платформы обязателен? В Mathcad реализация, думаю, будет проще. 1

@EvgenSanich 0 / 0 / 0 Регистрация: 17.07.2013 Сообщений: 35
	19.09.2013, 12:42 [ТС]
	Да, Matlab обязателен, требование в институте. Не совсем понял формулировку Matlab one-based indexing, в двух словах как говорят "на пальцах" можете пояснить что это за правило? я попробую разобраться с алгоритмом еще раз. А вот про бином, меня смущает выражение которое по книге идет как (k+j-1, k), а далее умножение на выражение p^kk * (1 - p)^jj * P(kk)* exp(-((kk + jj)* r * (tau/n))). То есть надо и (k+j-1) умножить на это выражение и k? 0

@EvgenSanich 0 / 0 / 0 Регистрация: 17.07.2013 Сообщений: 35
	19.09.2013, 14:00 [ТС]
	Это очередной биномиальный алгоритм из финансовой математики по оценке опционов, в котором я пытаюсь разобраться. В данном примере С, D и R это элементы формул (C - колл, D - диагональ, R - cтавка). Спасибо Вам огромное за подсказки, буду разбираться дальше. 0

Циклы со сложными математическими выражениям, как правильно записать?!

Решение