|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
想通过matlab拟合得到二元函数y=f(x1,x2)的表达式,但是这里面18组初始的b值是参考别人的模板,因为是用另外的数据生成的,请问怎么用程序实现我这组数据里面的18组初始b值呢?比如用nlinfit,或者cftool,请大佬多多指教! 之前用cftool拟合后得到的poly23多项式函数的函数值域与已知数据相差几个数量级,函数表达式明显是错误的
$ v. S5 S' h" H; W1 x; ?9 L0 L代码如下:6 W" F8 R |( X! u( N& D
clear,clc
4 A* N$ H5 }; b b4 H- ex1=[1.2:.1:1.7]';
6 A4 \6 Z3 X6 C( ?; q# bx2=[1000 5000 7500 10000 15000];
$ l, M# g5 ^# q1 d. Xy=[14964 13166 14235 15550 19200, d2 t/ _) v& w1 p$ i1 E
13479 13479 13090 14235 17014$ l- ?% U$ E; n+ y1 V
13747 13747 12750 13139 14991 [9 x. _. v2 c) _& n
13917 13917 13019 12553 13967
! P i$ R# d# J9 P* U 14065 14065 13386 12821 13188
9 @9 a* F4 w" X5 w s 14306 14306 13578 13040 12728];
! ^2 ?' V9 g( X) h, q2 {n1=length(x1);n2=length(x2);8 y4 J" d1 v) |! e. c" U t; G
x1=x1*ones(1,n2);x1=x1( ;) J& H6 E; ]1 D! R: V; x8 T
x2=ones(n1,1)*x2;x2=x2( ;% [/ T4 d/ K! n1 \
y=y( ; X=[x1,x2]; n=length(y);& w% w6 \5 D0 p' j* f* c0 s- y* J
str=num2str([1:n]');( `8 s! R8 m/ U
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));$ t& U* f8 C) Z+ f) t- j1 j
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)));; }3 q& n& s4 e1 j& [- E& w9 g
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];
* n) \9 `7 t8 x2 Efor l=1:53 s- B+ A9 f/ I# N) ^# T) i
b=lsqcurvefit(fx2,b,X,y);
" q4 V& F# e( B% U" U# o0 l6 U% k- S b=nlinfit(X,y,fx2,b);9 w$ z4 E: t* `) y2 @& S4 I
end* f/ A7 w" W+ ?, \
b
% N8 f4 u' Y' R/ E9 ?6 Dfigure(1),clf
) i" Z5 @* p# V2 _3 B0 Tplot3(x1,x2,y,'o')8 L/ R7 [7 L# f0 d' m" e" k
stem3(x1,x2,y,'filled')
' s2 O1 c7 y3 w: d8 j: N4 r, k# xtext(x1,x2,y+.01,str)1 f5 f5 Y X6 o, O$ I0 ^+ ?+ E# X
hold on/ ?: U6 O' I+ u9 t" g
[x11,x22]=meshgrid(min(x1):range(x1)/80:max(x1),min(x2):range(x2)/80:max(x2));
* p8 I$ }8 |) b) A" W$ H iyhat=fx1(b,x11,x22);
' L! l. j/ S: A7 |! @$ `suRF(x11,x22,yhat)+ p' L1 |( V" ?/ [! o; P4 Q. {+ j
shading interp
% ?( i) Z; y, Z% E9 g" O, i7 Qalpha(.8) E i8 _4 m2 x8 i2 m
axis tight
4 [. j: N3 e* U& m9 ixlabel('α');ylabel('C'),zlabel('V1')0 e4 f; x, k3 G; A. ]! v/ X; O% }+ M
SSy=var(y)*(n-1)
2 H9 r; [6 o, H# O; R8 i) Ey1=fx1(b,x1,x2);3 L# U3 ?; i0 H+ \, J. u
RSS=(y-y1)'*(y-y1)
+ e" J: O8 s, z9 b7 trsquare=(SSy-RSS)/SSy0 L+ Y3 j5 G# t0 z
MSe=RSS/(n-length(b))
9 T6 s" p5 k$ c" D
4 O; r: R, j% S. v5 k3 i/ [2 a |
|