Аппроксимировать экспериментальные данные степенной функцией a*x^b методом наименьших квадратов

@Maksim45554 · Регистрация: 08.02.2016

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

Имеется файл с сотней значений x и y (прикреплен). Найти , используя метод наименьших квадратов(!), коэффиценты a и b аппроксимирующей a*x^b функции.
Запрещено использовать какие-либо встроенные функции ( например, datafit очень легко и верно справляется ). Дело в том , что в файле , x в самом начале принимает значение 0, а при попытке линеаризировать степенную функцию используется натуральный логарифм (этот ноль выводится из логарифма как -Inf) попытка как-то "подогнать" этот "0", заменив его на какое-то конкретное число ни к чему не приводит,аппроксимирующая функция получается неверной ( к тому же, почему-то начинает визуально очень сильно походить на прямую) ,подскажите как решить проблему, возможно, идея , предложенная в коде, изначально не верна.

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
Fd=read('C:\Users\Username\Desktop\probe.txt',-1,2);
x=[];z=[];
for i = 1:100
    z(1,i)=(Fd(i,1)); //Считываются координаты для x и y 
    z(2,i)=(Fd(i,2));
end
for i = 1:100
    x(i,1)=(log(z(1,i))); //Считываются прологарифмированные координаты для x и y 
    x(i,2)=(log(z(2,i)));
end
disp(x)
SUMMx2=0;SUMMx1=0;SUMMy1=0;SUMMX1Y1=0;SUMMX2Y1=0;
for i = 1:100
    SUMMx2=SUMMx2+x(i,1)^2;
    SUMMx1=SUMMx1+x(i,1);
    SUMMy1=SUMMy1+x(i,2);
    SUMMX1Y1=SUMMX1Y1+x(i,1)*x(i,2);
end
Kx = [100,SUMMx1
      SUMMx1,SUMMx2];
Ky = [SUMMy1
      SUMMX1Y1];
disp(Kx,Ky);
a=linsolve(Kx,Ky); //Временное использование встроенной , будет заменена на метод Крамера 
disp(a)
for v = 1:100
    xX(v)=z(1,v);
    yY(v)=z(2,v);
end
plot2d(xX,yY,-4);
k=0.00:0.05:5;
M=exp(a(1));
PTC=-(a(2))*k^(M);
plot2d(k,PTC);

Код с использованием datafit (все работает идеально), однако ее нельзя использовать

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Fd=read('C:\Users\Username\Desktop\probe.txt',-1,2);
x=[];
for i = 1:100
    x(1,i)=Fd(i,1); 
    x(2,i)=Fd(i,2);
    plot2d(x(1,i),x(2,i),-4);
end
c=[0;0];
function [ZR]=G(c,z)
    ZR = z(2)-c(1)*(z(1)^c(2))
endfunction
[a,err]=datafit(G,x,c);
disp(a);
disp(err);
PTC1=a(1)*(k)^a(2)
plot2d(k,PTC1);

@Krasme · 11.12.2018, 01:47

есть вариант - фильтровать "неудобные" точки

Кликните здесь для просмотра всего текста

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
clc; clear; clf;
Fd=read('probe.txt',-1,2);
k=0
for i=1:size(Fd,1)
    if (Fd(i,1)>0 & Fd(i,2)>0) then
        k=k+1
        z(k,:)=Fd(i,:)
    end
end
plot(Fd(:,1),Fd(:,2),'*g'); // original
plot(z(:,1),z(:,2),'or'); // filtred
set(gca(),"grid",[1 1])
 
x=log(z); //Считываются прологарифмированные координаты для x и y 
n=size(z,1)
mX2=sum(x(:,1).^2)/n
mX1=sum(x(:,1))/n
mY1=sum(x(:,2))/n
mX1Y1=sum(x(:,1).*x(:,2))/n
//disp([mX2, mX1, mY1, mX1Y1])
 
Kx = [1 mX1; mX1 mX2]; //disp(Kx)
Ky = [mY1; mX1Y1];     //disp(Ky)
a1=Kx\Ky
a1(1)=exp(a1(1))
k=0.0:0.5:5
PTC=a1(1)*k^a1(2);
plot2d(k,PTC);
// проверка на встроенной функции
c=[0;0];
function [ZR]=G(c,z)
    ZR = z(2)-c(1)*(z(1)^c(2))
endfunction
//[a2,err]=datafit(G,Fd',c); // оригинальные точки
[a2,err]=datafit(G,z',c);  //отфильтрованные точки
//disp(err);
PTC1=a2(1)*k^a2(2)
plot2d(k,PTC1,5);
disp([a1 a2])

@olbond · 11.12.2018, 03:50

Все можно сделать гораздо проще. Поскольку область определения апроксимирующей функции известна, можно линеаризовать функцию использую десятичный логарифм:

$https://www.cyberforum.ru/cgi-bin/latex.cgi?log10\left(y\right)=log10\left(a\right)+b*log10\left(x\right)$

и затем применить линейный МНК (также известный как нормальное уравнение):

$https://www.cyberforum.ru/cgi-bin/latex.cgi?\binom{a}{b}=\left(x^T*x\right)^{-1}*x^T*y$

Из набора экспериментальных точек надо использовать только положительные:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
Fd=read('C:\Users\Username\Desktop\probe.txt',-1,2);
 
ind = find(Fd(:,1)>0 & Fd(:,2)>0);
 
x = log10(Fd(ind,1));
 
y = log10(Fd(ind,2));
 
X = [ones(x),x];
 
theta = pinv(X'*X)*X'*y;
 
a = 10^theta(1);
b = theta(2);
 
plot2d(Fd(:,1),Fd(:,2),style=-9);xgrid(color('grey75'))
plot2d([0:0.1:5],a*[0:0.1:5].^b, style=2);

Новые блоги и статьи Все статьи Все блоги /
Access VikBal 11.12.2025 Помогите пожалуйста !! Как объединить 2 одинаковые БД Access с разными данными.	Новый ноутбук volvo 07.12.2025 Всем привет. По скидке в "черную пятницу" взял себе новый ноутбук Lenovo ThinkBook 16 G7 на Амазоне: Ryzen 5 7533HS 64 Gb DDR5 1Tb NVMe 16" Full HD Display Win11 Pro	Музыка, написанная Искусственным Интеллектом volvo 04.12.2025 Всем привет. Некоторое время назад меня заинтересовало, что уже умеет ИИ в плане написания музыки для песен, и, собственно, исполнения этих самых песен. Стихов у нас много, уже вышли 4 книги, еще 3. . .	От async/await к виртуальным потокам в Python IndentationError 23.11.2025 Армин Ронахер поставил под сомнение async/ await. Создатель Flask заявляет: цветные функции - провал, виртуальные потоки - решение. Не threading-динозавры, а новое поколение лёгких потоков. Откат?. . .	Поиск "дружественных имён" СОМ портов Argus19 22.11.2025 Поиск "дружественных имён" СОМ портов На странице: https:/ / norseev. ru/ 2018/ 01/ 04/ comportlist_windows/ нашёл схожую тему. Там приведён код на С++, который показывает только имена СОМ портов, типа,. . .
Сколько Государство потратило денег на меня, обеспечивая инсулином. Programma_Boinc 20.11.2025 Сколько Государство потратило денег на меня, обеспечивая инсулином. Вот решила сделать интересный приблизительный подсчет, сколько государство потратило на меня денег на покупку инсулинов. . . .	Ломающие изменения в C#.NStar Alpha Etyuhibosecyu 20.11.2025 Уже можно не только тестировать, но и пользоваться C#. NStar - писать оконные приложения, содержащие надписи, кнопки, текстовые поля и даже изображения, например, моя игра "Три в ряд" написана на этом. . .	Мысли в слух kumehtar 18.11.2025 Кстати, совсем недавно имел разговор на тему медитаций с людьми. И обнаружил, что они вообще не понимают что такое медитация и зачем она нужна. Самые базовые вещи. Для них это - когда просто люди. . .	Создание Single Page Application на фреймах krapotkin 16.11.2025 Статья исключительно для начинающих. Подходы оригинальностью не блещут. В век Веб все очень привыкли к дизайну Single-Page-Application . Быстренько разберем подход "на фреймах". Мы делаем одну. . .	Фото: Daniel Greenwood kumehtar 13.11.2025

Аппроксимировать экспериментальные данные степенной функцией a*x^b методом наименьших квадратов

Решение