Применение ммп для матриц

@Aprill · Регистрация: 17.10.2015

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

Здравствуйте, помогите, пожалуйста, разобраться с задачей:
На вход подается бинарная матрица, каждый элемент которой Хij, зависит от элемента θi вектора θ и элемента βj вектора β. Требуется по выборке {Хij} определить вектора θ и β. (рис.1)

Определение векторов происходит на основе метода максимального правдоподобия. (β1 полагается равным нулю, т.е остальсные β будут рассчитываться относительно β1) Функция правдоподобия (рис.2), (рис.3)
С целью максимизации используется квазиньютоновский алгоритм Бройдена-Флетчера-Голдфарба-Шенно. Стартовой точкой метода полагается нулевой вектор u0 (рис.4)
Результатом выполнения алгоритма должны являться вектор δ и вектор θ

При реализации алгоритма для примера (с известным ответом) матрицы 10*10 никак не могу получить нужные вектора.

f.m

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
function L = f(u)%f(u_theta, u_delta)
L = 0;
theta = 1:10;
delta = 1:10;
for i = 1:10
    theta(i) = u(i);
end;
 
for i = 11:20
    delta(i-10) = u(i);
end;
x = [1, 0, 1, 1, 0, 1, 1, 1, 0, 1 ; ...
      0, 0, 0, 1, 1, 1, 0, 1, 1, 1 ; ...
      1, 1, 1, 0, 1, 1, 1, 1, 1, 1 ; ...
      1, 1, 0, 0, 0, 0, 1, 0, 1, 1 ; ...
      0, 1, 1, 0, 0, 1, 0, 0, 1, 1 ; ...
      1, 0, 1, 1, 1, 1, 0, 1, 0, 1 ; ...
      0, 1, 1, 1, 1, 1, 1, 1, 1, 1 ; ...
      0, 0, 1, 1, 1, 1, 0, 0, 0, 0 ; ...
      1, 1, 1, 0, 1, 1, 1, 1, 1, 1 ; ...
      1, 1, 1, 1, 1, 1, 0, 1, 1, 0 ; ];
for i = 1:10
   for j = 2:10
    L = L + (theta(i))*x(i,1) - log(1 + exp(theta(i))) + (theta(i)-delta(j))*x(i,j) - log(1 + exp(theta(i) - delta(j)));
   end
end
L = -L;
end %f(u_theta(k2,k1), u_delta(k2,k1) f(u_theta(k2,k1+1),u_delta(k2,k1+1)zz

Un.m

Matlab M
1
2
3
4
clear all;
x0 = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ];
 
fminunc(@f, x0)

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
>> Un
Warning: Gradient must be provided for trust-region algorithm;
  using line-search algorithm instead. 
> In fminunc at 383
  In Un at 4 
 
Local minimum found.
 
Optimization completed because the size of the gradient is less than
the default value of the function tolerance.
 
<stopping criteria details>
 
ans =
  Columns 1 through 18
    1.3808   -1.1211    2.6343    0.6940   -1.4191    1.3808   -0.3430   -1.7571    2.6343    1.8660         0    0.0000   -1.4667   -0.0000   -0.6934   -2.4800    0.6507   -0.6934
 
  Columns 19 through 20
   -0.6934   -1.4667
>>

Очень долго не могу получить требуемые вектора для конкретной матрицы. (рис5)
Пожалуйста, помогите разобраться. Очень надеюсь на вашу помощь!!!

Зосима · 21.10.2015, 16:59

Aprill, а функция невязки $https://www.cyberforum.ru/cgi-bin/latex.cgi?L\left(\theta ,\delta \right)$ правильно записана? и матрица х?

Функция невязки:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
function L = f(u)
theta = u(:,1); % первый столбец
delta = u(:,2); % второй столбец
x = [1, 0, 1, 1, 0, 1, 1, 1, 0, 1 ; ...
    0, 0, 0, 1, 1, 1, 0, 1, 1, 1 ; ...
    1, 1, 1, 0, 1, 1, 1, 1, 1, 1 ; ...
    1, 1, 0, 0, 0, 0, 1, 0, 1, 1 ; ...
    0, 1, 1, 0, 0, 1, 0, 0, 1, 1 ; ...
    1, 0, 1, 1, 1, 1, 0, 1, 0, 1 ; ...
    0, 1, 1, 1, 1, 1, 1, 1, 1, 1 ; ...
    0, 0, 1, 1, 1, 1, 0, 0, 0, 0 ; ...
    1, 1, 1, 0, 1, 1, 1, 1, 1, 1 ; ...
    1, 1, 1, 1, 1, 1, 0, 1, 1, 0 ];
L = 0;
for i = 1:size(x,1)
    for j = 2:size(x,2)
        L = L + (theta(i))*x(i,1) - log(1 + exp(theta(i))) + (theta(i)-delta(j))*x(i,j) - log(1 + exp(theta(i) - delta(j)));
    end
end
L = -L;
end

Скрипт расчета:

Matlab M
1
2
3
4
5
clear all;
clc
x0 = zeros(10,2);
options = optimset('MaxFunEvals',1e5, 'MaxIter', 1e4);
u0 = fminsearch(@f, x0,options)

Результат:

Matlab M
u0 =
 
    0.8430    1.3167
   -0.2372    0.8284
    1.3652   -0.3261
    0.5706   -1.0717
   -0.6487   -1.7048
    1.1067   -0.6679
   -0.1122   -0.2792
   -1.4290   -2.4606
    2.1424    0.5644
    0.7780   -2.0036

@Aprill · 26.10.2015, 17:23 **[ТС]**

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

Матрица и функция записаны вроде бы правильно.
Возможно вы правы, но результат, к сожалению, не совпадает с необходимыми значениями векторов(( (рисунок во вложении)

Зосима · 26.10.2015, 17:43

О!

покурил формулу для $https://www.cyberforum.ru/cgi-bin/latex.cgi?L\left(\theta ,\delta \right)$ и вот что записал:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
function L = f(u)
theta = u(:,1);
delta = u(:,2);
x = [1, 0, 1, 1, 0, 1, 1, 1, 0, 1 ; ...
    0, 0, 0, 1, 1, 1, 0, 1, 1, 1 ; ...
    1, 1, 1, 0, 1, 1, 1, 1, 1, 1 ; ...
    1, 1, 0, 0, 0, 0, 1, 0, 1, 1 ; ...
    0, 1, 1, 0, 0, 1, 0, 0, 1, 1 ; ...
    1, 0, 1, 1, 1, 1, 0, 1, 0, 1 ; ...
    0, 1, 1, 1, 1, 1, 1, 1, 1, 1 ; ...
    0, 0, 1, 1, 1, 1, 0, 0, 0, 0 ; ...
    1, 1, 1, 0, 1, 1, 1, 1, 1, 1 ; ...
    1, 1, 1, 1, 1, 1, 0, 1, 1, 0 ];
L = 0;
for i = 1:size(x,1)
    j = 2:size(x,2);
    L = L + theta(i)*x(i,1) - log(1 + exp(theta(i))) + ...
         sum( (theta(i)-delta(j)).*x(i,j)' - log(1 + exp(theta(i)-delta(j))) );
end
L = -L;
end

и скрипт расчета:

Matlab M
1
2
3
4
5
clear all;
clc
x0 = zeros(10,2);
options = optimset('MaxFunEvals',1e5, 'MaxIter', 1e4,'TolFun',1e-9);
u0 = fminunc(@f, x0,options)

но самое интересное это результат

u0 =

0.4006 0
-0.0827 -0.0000
1.8274 -1.1622
-0.5323 -0.0000
-0.5323 -0.5352
0.4006 -2.0554
1.8274 0.5036
-0.9855 -0.5352
1.8274 -0.5352
0.9788 -1.1622

Левый столбец это тета, правый - дельта.
Тебе эти числа ничего не напоминают?

@Aprill · 26.10.2015, 18:04 **[ТС]**

Ааа, спасибо вам огромноеее!!! Просто не знаю, как вас благодарить!!! Я очень признательна вам за то, что есть такие отзывчивые и добрые люди! Спасибо!!!

Зосима · 26.10.2015, 18:07

Не по теме:

на здоровье! :yes:

@Aprill · 27.10.2015, 12:09 **[ТС]**

Еще один маленький вопрос.
В первом файле я генерирую n-ое количество матриц T.

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
clear;
bet = [7 9 -8 9 3 -8 -5 1 5 -3];
tet = [-7 10 1 0 6 -8 -2 9 -4 -5];
for i = 1:10
    for j = 1:10
        P(i,j) = 1/(1+exp(bet(i)-tet(j)));
    end
end
 
for k = 1:50          %Кол-во матриц T(i,j)
    r = rand;          
    for i = 1:10
        for j = 1:10
            if r <= P(i,j)
                T(i,j,k) = 1;
            else 
                T(i,j,k) = 0;
            end
        end
    end
end
T

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

Мои попытки увы не увенчались успехом(

Зосима · 27.10.2015, 12:39

Смотри

я чуток переделал функцию: ввел глобальную переменную T с матрицей

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
function L = f(u)
global T
theta = u(1,:);
delta = u(2,:);
 
L = 0;
for i = 1:size(T,1)
    j = 2:size(T,2);
    L = L + theta(i)*T(i,1) - log(1 + exp(theta(i))) +...
         sum( (theta(i)-delta(j)).*T(i,j) - log(1 + exp(theta(i)-delta(j))) );
end
L = -L;
end

И программа расчетов:

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
clear, clc
global T
 
bet = [7 9 -8 9 3 -8 -5 1 5 -3];
tet = [-7 10 1 0 6 -8 -2 9 -4 -5];
for i = 1:length(bet)
    for j = 1:length(tet)
        P(i,j) = 1/(1+exp(bet(i)-tet(j)));
    end
end
 
x0 = zeros(2,length(bet));
options = optimset('MaxFunEvals',1e5, 'MaxIter', 1e4,'TolFun',1e-9);
 
K = 8;  %Кол-во матриц T(i,j)
 
Tet = zeros(K,length(tet)); % заготовка
Bet = zeros(K,length(tet)); % заготовка
h = waitbar(0,'Ша! Я считаю...');
for k = 1:K 
    r = rand(size(T));
    T = zeros(size(P)); % заготовка
    for i = 1:size(T,1)
        for j = 1:size(T,2)
            if r <= P(i,j)
                T(i,j) = 1;
            end
        end
    end 
    u0 = fminunc(@f, x0, options);
    Bet(k,:) = u0(1,:);
    Tet(k,:) = u0(2,:);
    waitbar(k/K)
end
close(h)

НО! тут есть один момент:
если в строке 21 оставить как у тебя, т.е. r = rand - одно значение для всей матрицы, то оптимизация не хочет проходить - ругается на что-то

@Aprill · 27.10.2015, 13:02 **[ТС]**

Очень странно( Ведь r = rand это просто число, которое никак на оптимизацию не влияет...

Зосима · 27.10.2015, 13:07

Aprill, почему не влияет

мне кажется, что в таком варианте в матрице становится больше нулей или часть строк/столбцов становятся линейно-зависимыми и система становится не решаемой

@Aprill · 27.10.2015, 13:20 **[ТС]**

Не, все равно что-то с оптимизацией, даже если число поставить(

Matlab M
1
2
fminunc stopped because it cannot decrease the objective function
along the current search direction.

я без идей даже, как это можно устранить(

Зосима · 27.10.2015, 13:27

у тебя красным ошибки пишет?
если нет - значит не беда

@Aprill · 27.10.2015, 13:36 **[ТС]**

Это я в курсе) Я про random)

А как сделать, чтобы в результате все пары векторов отображались корректно?
Прошу прошения уже за глупые вопросы)

Зосима · 27.10.2015, 13:45

я этот момент прошляпил, так что исправляюсь

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
clear, clc
global T
 
bet = [7 9 -8 9 3 -8 -5 1 5 -3];
tet = [-7 10 1 0 6 -8 -2 9 -4 -5];
for i = 1:length(bet)
    for j = 1:length(tet)
        P(i,j) = 1/(1+exp(bet(i)-tet(j)));
    end
end
 
x0 = zeros(2,length(bet));
options = optimset('MaxFunEvals',1e5, 'MaxIter', 1e4,'TolFun',1e-9);
 
K = 8;  %Кол-во матриц T(i,j)
 
U = zeros(2,length(tet),K); % заготовка
h = waitbar(0,'Ша! Я считаю...');
for k = 1:K 
    r = rand(size(T));
    T = zeros(size(P)); % заготовка
    for i = 1:size(T,1)
        for j = 1:size(T,2)
            if r <= P(i,j)
                T(i,j) = 1;
            end
        end
    end 
    u0 = fminunc(@f, x0, options);
    U(:,:,k) = u0;
    waitbar(k/K)
end
close(h)
U(:,:,3) % отображаем третью пару веткоров

@Aprill · 27.10.2015, 13:56 **[ТС]**

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

Зосима · 27.10.2015, 14:04

Aprill, учти, я тоже могу ошибаться

@Aprill · 28.10.2015, 17:05 **[ТС]**

Немного странные результаты все же выходят(
Странно, что оптимизация выдает цифры порядка 937.95 если на входе много матриц, хотя если взять ту же сгенерированную случайно матрицу и запустить тот же процесс, получаются вполне нормальные цифры порядка 2.77 и т. д.

Зосима · 28.10.2015, 17:21

Aprill, может быть такое, что оптимизация выдает ошибочный результат

@Aprill · 28.10.2015, 20:32 **[ТС]**

Да, но ошибочный результат она выдает именно при работе с большим количеством матриц(
Когда на входе одна матрица, все в порядке...

Зосима · 28.10.2015, 22:19

т.е. когда К=1 считает правильно, когда K>1 оптимизация дает ошибочный результат? Весьма странно

Новые блоги и статьи Все статьи Все блоги /
Переходник USB-CAN-GPIO Eddy_Em 20.03.2026 Достаточно давно на работе возникла необходимость в переходнике CAN-USB с гальваноразвязкой, оный и был разработан. Однако, все меня терзала совесть, что аж 48-ногий МК используется так тупо: просто. . .	Оттенки серого Argus19 18.03.2026 Оттенки серого Нашёл в интернете 3 прекрасных модуля: Модуль класса открытия диалога открытия/ сохранения файла на Win32 API; Модуль класса быстрого перекодирования цветного изображения в оттенки. . .	SDL3 для Desktop (MinGW): Рисуем цветные прямоугольники с помощью рисовальщика SDL3 на Си и C++ 8Observer8 17.03.2026 Содержание блога Финальные проекты на Си и на C++: finish-rectangles-sdl3-c. zip finish-rectangles-sdl3-cpp. zip	Символические и жёсткие ссылки в Linux. algri14 15.03.2026 Существует два типа ссылок — символические и жёсткие. Ссылка в Linux — это запись в каталоге, которая может указывать либо на inode «файла-ИСТОЧНИКА», тогда это будет «жёсткая ссылка» (hard link),. . .
[Owen Logic] Поддержание уровня воды в резервуаре количеством включённых насосов: моделирование и выбор регулятора ФедосеевПавел 14.03.2026 Поддержание уровня воды в резервуаре количеством включённых насосов: моделирование и выбор регулятора ВВЕДЕНИЕ Выполняя задание на управление насосной группой заполнения резервуара,. . .	делаю науч статью по влиянию грибов на сукцессию anaschu 13.03.2026 прикрепляю статью	SDL3 для Desktop (MinGW): Создаём пустое окно с нуля для 2D-графики на SDL3, Си и C++ 8Observer8 10.03.2026 Содержание блога Финальные проекты на Си и на C++: hello-sdl3-c. zip hello-sdl3-cpp. zip Результат:	Установка CMake и MinGW 13.1 для сборки С и C++ приложений из консоли и из Qt Creator в EXE 8Observer8 10.03.2026 Содержание блога MinGW - это коллекция инструментов для сборки приложений в EXE. CMake - это система сборки приложений. Здесь описаны базовые шаги для старта программирования с помощью CMake и. . .

@Aprill 0 / 0 / 0 Регистрация: 17.10.2015 Сообщений: 26
	26.10.2015, 17:23 [ТС]
	Большое вам спасибо! К сожалению, не видела что вы ответили раньше. Не думала, что найдется кто-либо, готовый помочь... Матрица и функция записаны вроде бы правильно. Возможно вы правы, но результат, к сожалению, не совпадает с необходимыми значениями векторов(( (рисунок во вложении) Миниатюры 0

@Aprill 0 / 0 / 0 Регистрация: 17.10.2015 Сообщений: 26
	26.10.2015, 18:04 [ТС]
	Ааа, спасибо вам огромноеее!!! Просто не знаю, как вас благодарить!!! Я очень признательна вам за то, что есть такие отзывчивые и добрые люди! Спасибо!!! 0

Зосима
	26.10.2015, 18:07
	Не по теме: на здоровье! :yes: 0

@Aprill 0 / 0 / 0 Регистрация: 17.10.2015 Сообщений: 26
	27.10.2015, 13:02 [ТС]
	Очень странно( Ведь r = rand это просто число, которое никак на оптимизацию не влияет... 0

Зосима 5246 / 3574 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	27.10.2015, 13:07
	Aprill, почему не влияет мне кажется, что в таком варианте в матрице становится больше нулей или часть строк/столбцов становятся линейно-зависимыми и система становится не решаемой 0

Зосима 5246 / 3574 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	27.10.2015, 13:27
	у тебя красным ошибки пишет? если нет - значит не беда 0

@Aprill 0 / 0 / 0 Регистрация: 17.10.2015 Сообщений: 26
	27.10.2015, 13:36 [ТС]
	Это я в курсе) Я про random) А как сделать, чтобы в результате все пары векторов отображались корректно? Прошу прошения уже за глупые вопросы) 0

@Aprill 0 / 0 / 0 Регистрация: 17.10.2015 Сообщений: 26
	27.10.2015, 13:56 [ТС]
	Странно, что появляются одинаковые вектора. Буду думать, как быть с оптимизацией и rand. Надеюсь, не доставляю вам много хлопот. Спасибо! 0

Зосима 5246 / 3574 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	27.10.2015, 14:04
	Aprill, учти, я тоже могу ошибаться 0

@Aprill 0 / 0 / 0 Регистрация: 17.10.2015 Сообщений: 26
	28.10.2015, 17:05 [ТС]
	Немного странные результаты все же выходят( Странно, что оптимизация выдает цифры порядка 937.95 если на входе много матриц, хотя если взять ту же сгенерированную случайно матрицу и запустить тот же процесс, получаются вполне нормальные цифры порядка 2.77 и т. д. 0

Зосима 5246 / 3574 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	28.10.2015, 17:21
	Aprill, может быть такое, что оптимизация выдает ошибочный результат 0

@Aprill 0 / 0 / 0 Регистрация: 17.10.2015 Сообщений: 26
	28.10.2015, 20:32 [ТС]
	Да, но ошибочный результат она выдает именно при работе с большим количеством матриц( Когда на входе одна матрица, все в порядке... 0

Зосима 5246 / 3574 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	28.10.2015, 22:19
	т.е. когда К=1 считает правильно, когда K>1 оптимизация дает ошибочный результат? Весьма странно 0

Применение ммп для матриц

Решение