|
|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
想通过matlab拟合得到二元函数y=f(x1,x2)的表达式,但是这里面18组初始的b值是参考别人的模板,因为是用另外的数据生成的,请问怎么用程序实现我这组数据里面的18组初始b值呢?比如用nlinfit,或者cftool,请大佬多多指教! 之前用cftool拟合后得到的poly23多项式函数的函数值域与已知数据相差几个数量级,函数表达式明显是错误的
6 Z/ ^2 i3 B* W5 u8 U( g代码如下:
# J5 x6 a1 S8 e& p8 n% @$ qclear,clc; A4 f/ i& _: _4 z% `8 Z/ W) q
x1=[1.2:.1:1.7]';% M6 ]& O6 S: M# ?6 f
x2=[1000 5000 7500 10000 15000];) V F& e7 ^0 j
y=[14964 13166 14235 15550 192004 v# y5 a' P' F' M0 Z5 E7 C
13479 13479 13090 14235 17014# [, D% v% r9 w0 c% W% n
13747 13747 12750 13139 149919 F, ^: e7 l3 v+ S
13917 13917 13019 12553 13967
# M( ^: e9 |( g, `* ~. { 14065 14065 13386 12821 13188
' @9 n1 C* _5 v6 D3 x! N! x0 d 14306 14306 13578 13040 12728];/ m5 @( D0 j6 j n3 I! B' u
n1=length(x1);n2=length(x2);
9 ` h: x2 L; W6 c4 j% g# r7 ix1=x1*ones(1,n2);x1=x1( ;
, ^4 W$ p( V Q" g3 z* a- Mx2=ones(n1,1)*x2;x2=x2( ;% F. S7 Y: k# C! x+ v3 e3 T% E
y=y( ; X=[x1,x2]; n=length(y);
8 k, `# ]/ ~" xstr=num2str([1:n]');
4 l8 }: G9 Y. Q6 o V- C) afx1=@(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));3 X8 V3 A8 H, @& l# ^! a
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)));. j/ l; Y$ m& q1 ~
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];
& j. @8 g3 Y# [9 mfor l=1:5. y9 l; \7 z# ^8 W7 M h2 S
b=lsqcurvefit(fx2,b,X,y);" a8 v$ r$ L4 {, {. A/ u9 R
b=nlinfit(X,y,fx2,b);7 m7 E& R- N8 ~3 X$ _
end5 D/ i1 J1 Q, w, y# j0 W
b
6 C# M6 N' S( f+ P" e9 p+ ^figure(1),clf7 E! l/ `; R- u
plot3(x1,x2,y,'o')
, }; V* c/ U" X/ C; X+ \stem3(x1,x2,y,'filled')
: ]& P7 o: B1 rtext(x1,x2,y+.01,str)
" J1 i3 U! n: W! Y; dhold on v% k" i0 M' I" S8 x: I
[x11,x22]=meshgrid(min(x1):range(x1)/80:max(x1),min(x2):range(x2)/80:max(x2));
! u& ^( P! I3 s3 r8 V# nyhat=fx1(b,x11,x22);
2 R3 X+ V) r9 W1 g" N7 MsuRF(x11,x22,yhat)
6 o# j4 ~' \; i0 u7 Mshading interp
; W Z- @7 e6 H* talpha(.8)# F- q6 r- S/ _2 M& [6 b9 a3 u
axis tight
# ]& L" p h3 b8 M/ S9 hxlabel('α');ylabel('C'),zlabel('V1')! E- g, {3 g: s% a6 S5 R, \8 n
SSy=var(y)*(n-1)$ V2 a) |& K$ }2 t: F
y1=fx1(b,x1,x2);
2 p% _: w, q1 |4 v2 G0 T# D- ORSS=(y-y1)'*(y-y1)
$ p& \/ j9 j! B/ R: Crsquare=(SSy-RSS)/SSy
: o/ g+ ]6 e- UMSe=RSS/(n-length(b))8 p. p7 s1 z/ F5 W
7 v% x: [1 ~ O* {5 R- I% v
|
|