时间戳

首页 > 解决方案 > 如何找到方程拟合参数的统计量

matlab - 如何找到方程拟合参数的统计量

问题描述

例如：a、b和cn 是三个常数，需要在特定方程中使用数据拟合方法计算。

如何计算自定义方程（例如二次高原方程）的参数的统计数据（平均值、标准差、方差、偏度值和学生 t 检验值）？

例子：

x=[0,40,80,100,120,150,170,200], 
y=[1865,2855,3608,4057,4343,4389,4415,4478]

y=a*(x+n)^2+b*(x+n)+c, x < xc(Ymax) ....(1)   y=yp, x >= xc(Ymax) ....(2)

我已经通过给定的代码拟合了这个方程：

yf = @(b,x) b(1).*(x+n).^2+b(2)*(x+n)+b(3); B0 = [0.006; 21; 1878];
[Bm,normresm] = fminsearch(@(b) norm(y - yf(b,x)), B0); a=Bm(1);
b=Bm(2); c=Bm(3); xc=(-b/(2*a))-n; p=p=a*(xc+n)^2+b*(xc+n)+c;  
if (x < xc) 
    yfit = a.*(x+n).^2+ b*(x+n)+c;  
else 
    yfit = p;  
end 
plot(x,yfit,'*') 
hold on; plot(x,y); hold off

注意：我已经使用了该polyfit命令，它很有帮助并为我提供了结果。但是，我真的觉得它不合适，因为没有自定义方程的选项。我可以通过任何代码找到这些统计信息吗？

标签： matlabcurve-fittingequation

解决方案

问题 1、2 和 4) 如果您对方程系统有先前的了解，好的做法是设置接近最终结果的初始值：

你所拥有的是一个超定的线性方程组。

y(1) = a*x(1)^2 + b*x(1) + c
y(2) = a*x(2)^2 + b*x(2) + c
y(3) = a*x(3)^2 + b*x(3) + c
…
y(n) = a*x(n)^2 + b*x(n) + c

or in general:
y = A*X, where

A = [a; b; c]
X = [x(1)^2 x(1) 1;
     x(2)^2 x(2) 1;
     x(3)^2 x(2) 1;
     ...
     x(n)^2 x(n) 1]

拟合超定系统（因为它没有解）的一种常见做法是“最小二乘拟合”（mldivide，\（链接））

x=[0; 40; 80; 100; 120; 150; 170; 200];
y=[1865; 2855; 3608; 4057; 4343; 4389; 4415; 4478];

X = [x.^2 x ones(numel(x),1)];

A = y\X;

a0=A(1); %- initial value for a
b0=A(2); %- initial value for b
c0=A(3); %- initial value for c

当您自定义 X 和 A 时，您可以自定义方程式

但您也可以将初始值设置为 1，它对结果的影响应该可以忽略不计。更多与问题 4相关

a0=1;
b0=1; 
c0=1;

或随机值

rng(10);
A = rand(3,1);

a0=A(1);
b0=A(2); 
c0=A(3);

问题 3 - 统计 如果您需要对优化过程的监控进行更多控制，请使用更通用的匿名函数形式（在下面的代码中> myfun）来保存参数的所有中间值（a_iter、b_iter、c_iter）

function Fiting_ex()
global a_iter b_iter c_iter
a_iter = 0;
b_iter = 0;
c_iter = 0;
x=[0; 40; 80; 100; 120; 150; 170; 200];
y=[1865; 2855; 3608; 4057; 4343; 4389; 4415; 4478];

X = [x.^2 x ones(numel(x),1)];

A = y\X;

a0=A(1);
b0=A(2); 
c0=A(3);

B0 = [a0; b0; c0];

[Bm,normresm] = fminsearch(@(b) myfun(b,x,y),B0);

a=Bm(1);
b=Bm(2); 
c=Bm(3); 

xc=-b/(2*a); 
p=c-(b^2/(4*a));
yfit = zeros(numel(x),1);
for i=1:numel(x)
    if (x(i) < xc)
        yfit(i) = a.*x(i).^2+ b*x(i)+c;
    else
        yfit(i) = p;
    end
end

plot(x,yfit,'*') 
hold on; 
plot(x,y); 
hold off

% Statistic on optimization process
a_mean = mean(a_iter(2:end)); % mean value
a_var = var(a_iter(2:end)); % variance
a_std = std(a_iter(2:end)); % standard deviation

function f = myfun(Bm, x, y)
global a_iter b_iter c_iter
a_iter = [a_iter Bm(1)];
b_iter = [b_iter Bm(2)]; 
c_iter = [c_iter Bm(3)];

yf = Bm(1)*(x).^2+Bm(2)*(x)+Bm(3); 
a=Bm(1);
b=Bm(2); 
c=Bm(3); 
xc=-b/(2*a); 
p=c-(b^2/(4*a));
yfit = zeros(numel(x),1);
for i=1:numel(x)
    if (x(i) < xc)
        yfit(i) = a.*x(i).^2+ b*x(i)+c;
    else
        yfit(i) = p;
    end
end
f = norm(y - yfit);

推荐阅读