Метод наименьших квадратов

@q74 · Регистрация: 14.11.2016

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

Хочу реализовать МНК для отыскания неизвестных коэффициентов.

У меня на входе линейная функция и я хочу найти параметры a и b.

Написал следующий код, но почему то матлаб выдает ошибку ('a = Empty sym: 0-by-1')

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
syms x y a b
 
f=(y-a.*x+b).^2;
df1=diff(f,a);
df2=diff(f,b);
 
X=-5:0.1:5;
Y=X.*2;
 
df1=subs(df1,x,X);
df2=subs(df2,y,Y);
 
[a,b] = solve(sum(df1), sum(df2))

Что я делаю не так?

@nuHrBuH · 23.03.2020, 14:35

Какой-то сбой пошел

@SSC · 23.03.2020, 15:48

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

Что я делаю не так?

Решайте не символьно, а численно, и проблем не будет

@Centurio · 23.03.2020, 15:52

q74, а не лучше ли было CFT воспользоваться?

@q74 · 23.03.2020, 18:47 **[ТС]**

Переписал код следующим образом:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
h=histogram(d,'BinWidth',10,'Normalization','pdf');
x=h.BinLimits(1):h.BinWidth:h.BinLimits(2)-h.BinWidth;
 
syms mu1 sigma1
f1=exp(-((x-mu1).^2)./(2.*sigma1.^2))./(sigma1.*sqrt(2.*pi));
f=(h.Values-f1).^2;
 
df1=sum(diff(f,mu1));
df2=sum(diff(f,sigma1));
 
[mu1,sigma1]=solve(df1,df2)

Если функция h простая, например, линейная, то программа выдает результат.
В моем случае, программа виснет на вычислении 'sum'.

Я хочу использовать свою функцию, написанную самостоятельно, а не встроенные функции. Как можно исправить зависание программы в моем случае?

@nuHrBuH · 23.03.2020, 20:22

А значений h много (штук)? Это же символьные вычисления

@q74 · 23.03.2020, 20:56 **[ТС]**

h содержит 38 значений

@nuHrBuH · 23.03.2020, 21:51

Мое предложение

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
clc,clear
x = (1:40)/100;
h = sin(x);
syms mu1 sigma1
f1=exp(-((x-mu1).^2)./(2.*sigma1.^2))./(sigma1.*sqrt(2.*pi));
f=(h-f1).^2;
df1=sum(diff(f,mu1));
df2=sum(diff(f,sigma1));
f1 = matlabFunction(df1);
f2 = matlabFunction(df2);
F = @(x)[ f1(x(1),x(2));
         f2(x(1),x(2))];
res = fsolve(F,[1 1])

@q74 · 24.03.2020, 10:04 **[ТС]**

nuHrBuH, спасибо тебе большое!!!

Добавлено через 22 минуты
nuHrBuH, можешь, пожалуйста, объяснить эту строчку:

Matlab M
1
F = @(x)[ f1(x(1),x(2)); f2(x(1),x(2))];

Почему ты в функцию подставляешь только два первых значения вектора x?

Добавлено через 1 час 27 минут
nuHrBuH, еще один вопрос появился. Я дописал твою функцию с учетом своей задачи, вот что получилось:

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
h1=histogram(d,'BinWidth',1,'Normalization','pdf');
x=h1.BinLimits(1):h1.BinWidth:h1.BinLimits(2)-h1.BinWidth;
h=h1.Values;
 
syms mu1 sigma1 mu2 sigma2 k
f1=exp(-((x-mu1).^2)./(2.*sigma1.^2))./(sigma1.*sqrt(2.*pi));
f2=exp(-((x-mu2).^2)./(2.*sigma2.^2))./(sigma2.*sqrt(2.*pi));
f3=x.^((k./2)-1).*exp(-x./2).*((1./2).^(k./2))./gamma(k./2);
ff=f1+f2+f3;
f=(h-ff).^2;
df1=sum(diff(f,mu1));
df2=sum(diff(f,sigma1));
df3=sum(diff(f,mu2));
df4=sum(diff(f,sigma2));
df5=sum(diff(f,k));
f1 = matlabFunction(df1);
f2 = matlabFunction(df2);
f3 = matlabFunction(df3);
f4 = matlabFunction(df4);
f5 = matlabFunction(df5);
 
F = @(x)[f1(x(1),x(2),x(3),x(4),x(5));...
         f2(x(1),x(2),x(3),x(4),x(5));...
         f3(x(1),x(2),x(3),x(4),x(5));...
         f4(x(1),x(2),x(3),x(4),x(5))];
     
res = fsolve(F,[1 1 1 1 1]);

Матлаб выдает ошибку в fsolve. Насколько я понял, ошибка в начальном приближении параметра 'k'. Как исправить эту ситуацию? Ведь не всегда удается подобрать начальные точки.

@nuHrBuH · 24.03.2020, 11:37

q74, Я на провокации не отвечаю
У меня вопрос (скорее к форуму): как эту часть кода можно записать молее изящно?

Matlab M
1
2
3
4
5
6
7
8
9
10
f1 = matlabFunction(df1);
f2 = matlabFunction(df2);
f3 = matlabFunction(df3);
f4 = matlabFunction(df4);
f5 = matlabFunction(df5);
 
F = @(x)[f1(x(1),x(2),x(3),x(4),x(5));...
         f2(x(1),x(2),x(3),x(4),x(5));...
         f3(x(1),x(2),x(3),x(4),x(5));...
         f4(x(1),x(2),x(3),x(4),x(5))];

@SSC · 24.03.2020, 11:46

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

Мое предложение

Вы хотя бы проверяйте свои рекомендации
Вот мой код (в одном файле)

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 TMP95  
clc,clear
%x = (100:5:200)/100;
x = (1:40)/100;
h = sin(x);
syms mu1 sigma1
f1=exp(-((x-mu1).^2)./(2.*sigma1.^2))./(sigma1.*sqrt(2.*pi));
f=(h-f1).^2;
df1=sum(diff(f,mu1));
df2=sum(diff(f,sigma1));
f1 = matlabFunction(df1);
f2 = matlabFunction(df2);
F = @(x)[ f1(x(1),x(2));
         f2(x(1),x(2))];
res = fsolve(F,[1 1])
%===================================
y=exp(-((x-res(1)).^2)./(2.*res(2).^2))./(res(2).*sqrt(2.*pi));
 
[AA,resnorm]=lsqcurvefit(@FF,[1 1],x,h)
figure
y2=FF(AA,x);
 
plot(x,h,'or',x,y,'b',x,y2,'g')
grid on
end
function y=FF(a,x)
  y=exp(-((x-a(1)).^2)./(2.*a(2).^2))./(a(2).*sqrt(2.*pi));
end

Прилагаю график результатов
Красные исходные точки для подбора аппроксимации
nuHrBuH, Ваш график синий
Мой график зеленый (получен численным методом)

@nuHrBuH · 24.03.2020, 12:35

Ну а при чем тут я. Написано же правильно

Добавлено через 16 минут
Я попробывал увеличить точность

Matlab M
1
opt = optimset('tolfun',1e-10);

Но не помогло

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
No solution found.
 
fsolve stopped because the last step was ineffective. However, the vector of function
values is not near zero, as measured by the selected value of the function tolerance. 
 
<stopping criteria details>
fsolve stopped because the sum of squared function values, r, changed by 8.661077e-19
relative to its initial value; this is less than max(options.TolFun^2,eps) = 2.220446e-16.
However, r = 2.895904e-03, exceeds sqrt(options.TolFun) = 1.000000e-05.
 
Optimization Metric                                                 Options
relative change r =   8.66e-19                  max(TolFun^2,eps) =   2e-16 (selected)
r =   2.90e-03                                      sqrt(TolFun) =  1.0e-05 (selected)

Как побороть проблему?

@SSC · 24.03.2020, 13:08

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

Как побороть проблему?

Надо в место Вашего начального приближения [1 1], поставить [0.91 0.38], которое находится численным методом, и тогда все нормально.
Но правда возникает вопрос - если решение уже найдено, то накой фиг решать еще раз

Новые блоги и статьи Все статьи Все блоги /
SDL3 для Web (WebAssembly): Синхронизация спрайтов SDL3 и тел Box2D 8Observer8 04.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-sync-physics-sprites-sdl3-c. zip На первой гифке отладочные линии отключены, а на второй включены:. . .	SDL3 для Web (WebAssembly): Идентификация объектов на Box2D v3 - использование userData и событий коллизий 8Observer8 02.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-collision-events-sdl3-c. zip https:/ / www. cyberforum. ru/ blog_attachment. php?attachmentid=11680&d=1772460536 Одним из. . .	Реалии Hrethgir 01.03.2026 Нет, я не закончил до сих пор симулятор. Эта задача сложнее. Не получилось уйти в плавсостав, но оно и к лучшему, возможно. Точнее получалось - но сварщиком в палубную команду, а это значит, в моём. . .	Ритм жизни kumehtar 27.02.2026 Иногда приходится жить в ритме, где дел становится всё больше, а вовлечения в происходящее — всё меньше. Плотный график не даёт вниманию закрепиться ни на одном событии. Утро начинается с быстрых,. . .
SDL3 для Web (WebAssembly): Сборка библиотек: SDL3, Box2D, FreeType, SDL3_ttf, SDL3_mixer и SDL3_image из исходников с помощью CMake и Emscripten 8Observer8 27.02.2026 Недавно вышла версия 3. 4. 2 библиотеки SDL3. На странице официальной релиза доступны исходники, готовые DLL (для x86, x64, arm64), а также библиотеки для разработки под Android, MinGW и Visual Studio. . . .	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. На борту пять. . .

@nuHrBuH 483 / 427 / 205 Регистрация: 04.03.2011 Сообщений: 1,259
	23.03.2020, 14:35
	Какой-то сбой пошел Миниатюры 0

@Centurio Модератор 1765 / 1610 / 541 Регистрация: 13.09.2015 Сообщений: 5,629
	23.03.2020, 15:52
	q74, а не лучше ли было CFT воспользоваться? 0

@nuHrBuH 483 / 427 / 205 Регистрация: 04.03.2011 Сообщений: 1,259
	23.03.2020, 20:22
	А значений h много (штук)? Это же символьные вычисления 0

@q74 0 / 0 / 0 Регистрация: 14.11.2016 Сообщений: 35
	23.03.2020, 20:56 [ТС]
	h содержит 38 значений 0

Метод наименьших квадратов

Решение