|
|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
想通过matlab拟合得到二元函数y=f(x1,x2)的表达式,但是这里面18组初始的b值是参考别人的模板,因为是用另外的数据生成的,请问怎么用程序实现我这组数据里面的18组初始b值呢?比如用nlinfit,或者cftool,请大佬多多指教! 之前用cftool拟合后得到的poly23多项式函数的函数值域与已知数据相差几个数量级,函数表达式明显是错误的
8 u, Z% n0 L3 P3 Z代码如下:
& ~8 [. }. d# `7 D* Y- _clear,clc, y" }% j( c, U: M; M7 \
x1=[1.2:.1:1.7]'; V' a1 N0 G, M+ T) b) y
x2=[1000 5000 7500 10000 15000];
7 c3 a. _. Y% _% my=[14964 13166 14235 15550 19200
7 m$ p7 Y0 o4 h- f4 C7 p1 J 13479 13479 13090 14235 17014
9 `- I: Z7 [9 K( j& | 13747 13747 12750 13139 14991
$ X: C9 v2 Q8 V9 G6 g 13917 13917 13019 12553 13967
' m- }* Q, S6 T" ]& _" U( }2 X3 W 14065 14065 13386 12821 13188
* f" Y, \- C- `" X4 x 14306 14306 13578 13040 12728];& r! m6 [) e1 s
n1=length(x1);n2=length(x2);
3 y5 ?* x2 o& |5 X( H9 g. Fx1=x1*ones(1,n2);x1=x1( ;7 O1 _* e. a; P' F7 ~
x2=ones(n1,1)*x2;x2=x2( ;5 u: x" v8 V) P5 W# w4 o
y=y( ; X=[x1,x2]; n=length(y);) x$ J% c- X* p7 X( p
str=num2str([1:n]');" P6 \* R1 T4 k0 J
fx1=@(b,x1,x2)(b(1)+b(2)*x2+b(3)*x2.^2+b(4)*x1.*x2+b(5)*x1.^3+b(6)*x2.^3+b(7)*x1.^4+b(8)*x2.^4+b(9)*x1.*x2.^3+b(10)*x1.^2.*x2.^2+b(11)*x1.^5+b(12)*x2.^5)./(1+b(13)*exp(b(14)*x1+b(15)*x2+b(16)*x1.^2+b(17)*x2.^2+b(18)*x1.*x2));% M, q B* n3 X7 q
fx2=@(b,X)(b(1)+b(2)*X(:,2)+b(3)*X(:,2).^2+b(4)*X(:,1).*X(:,2)+b(5)*X(:,1).^3+b(6)*X(:,2).^3+b(7)*X(:,1).^4+b(8)*X(:,2).^4+b(9)*X(:,1).*X(:,2).^3+b(10)*X(:,1).^2.*X(:,2).^2+b(11)*X(:,1).^5+b(12)*X(:,2).^5)./(1+b(13)*exp(b(14)*X(:,1)+b(15)*X(:,2)+b(16)*X(:,1).^2+b(17)*X(:,2).^2+b(18)*X(:,1).*X(:,2)));8 ]' w0 _% R1 D1 m0 ]: G0 h: |
b=[345.45 -191.53 41.80 -2.6793 83.3156 -4.4545 -156.1311 0.2333 -0.0086582 0.251804 81.565 -0.0048047 989422 -146.607 4.003 95.35213571 -0.25589 2.752];
2 C9 M- C1 p, q4 E# V! bfor l=1:5
# F) K s( A1 x% b* Z: S3 s, E b=lsqcurvefit(fx2,b,X,y);
2 V0 T; a, K% f7 ]% h b=nlinfit(X,y,fx2,b);
4 m! X! X& Q4 G! Jend
/ P7 Z9 L7 {! T2 B6 h! u; B2 Pb
8 d( Y2 s( C7 f! K+ S. Rfigure(1),clf' e/ A4 P8 k6 R& r- P. T: d
plot3(x1,x2,y,'o')
( _, \2 f3 K! }: hstem3(x1,x2,y,'filled')* A' w5 Q/ |* E4 s3 o4 T2 @/ Z
text(x1,x2,y+.01,str)0 I' M: E+ o. d( V
hold on" u# |5 W8 {8 H/ j
[x11,x22]=meshgrid(min(x1):range(x1)/80:max(x1),min(x2):range(x2)/80:max(x2));1 C% z6 E! n! U) {) b% g- [: ]. d$ y
yhat=fx1(b,x11,x22);
$ w8 l: }/ O; u5 G# EsuRF(x11,x22,yhat)8 Z% ?' Q0 z. p! v
shading interp
2 U" G+ I1 M' @1 R* }alpha(.8)
" g' M: c2 ^6 C4 a: `axis tight, f" C: u8 s/ }- s1 G" B2 m) I8 Q+ A' d* `* S
xlabel('α');ylabel('C'),zlabel('V1')
) x* P" E- W5 k; ESSy=var(y)*(n-1)
. o$ }8 D, J( T) P% Dy1=fx1(b,x1,x2);: u4 S9 k8 G! h& p, W
RSS=(y-y1)'*(y-y1)8 l/ k- j2 W$ d: U" G7 ]
rsquare=(SSy-RSS)/SSy
7 W6 D) z* T$ @5 ~MSe=RSS/(n-length(b))8 K! w4 U8 T) n/ |3 R
; ^& y8 O0 W: l6 u' P1 ^7 M# ]) D |
|