|
|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
想通过matlab拟合得到二元函数y=f(x1,x2)的表达式,但是这里面18组初始的b值是参考别人的模板,因为是用另外的数据生成的,请问怎么用程序实现我这组数据里面的18组初始b值呢?比如用nlinfit,或者cftool,请大佬多多指教! 之前用cftool拟合后得到的poly23多项式函数的函数值域与已知数据相差几个数量级,函数表达式明显是错误的
1 \! O1 C' c! s. s. m" R) X3 T代码如下:) T: U0 @% X5 n- \0 B! n4 T
clear,clc
* m' D. u% t# px1=[1.2:.1:1.7]';* Y& `6 I; a1 D! e; {
x2=[1000 5000 7500 10000 15000];
o& M: m0 k+ g! _* z* P5 ~y=[14964 13166 14235 15550 192001 m. H" A* V* o+ b. w* ~
13479 13479 13090 14235 17014
0 K# \9 P, k2 e8 L$ [4 {" L 13747 13747 12750 13139 14991
6 A- a+ h8 a, U8 b) C% i& Y 13917 13917 13019 12553 13967$ G+ e2 S9 h! m7 n$ O8 N
14065 14065 13386 12821 13188
" q) ]1 l5 ~! E3 i- ]; I% z5 K8 A, K 14306 14306 13578 13040 12728];
% |1 x. v+ Q% p8 W* h( K& X/ m- wn1=length(x1);n2=length(x2);6 `& L# o3 _6 O
x1=x1*ones(1,n2);x1=x1( ;( u- H3 F1 M D* R3 X; e1 w
x2=ones(n1,1)*x2;x2=x2( ;
# w* ^$ [- \0 T e# m/ iy=y( ; X=[x1,x2]; n=length(y);
. z$ l- T# D' g* e$ F6 Q; Qstr=num2str([1:n]');$ b# ?$ P+ N/ F, h2 k
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));, W" W8 P( `( U# b# R
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)));" u! w2 W2 R$ h0 _
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];( [/ E- E; D/ L) A2 c B* v1 ^! G
for l=1:5. H/ _# A$ U+ Z7 v: O
b=lsqcurvefit(fx2,b,X,y);
* E8 Z- Y6 n( H) |3 r6 P, ` b=nlinfit(X,y,fx2,b);
0 P+ M/ x7 ^4 p( Y O# pend$ r! _ O$ C" p# d" o/ _! d, W
b
) I/ t; i9 _% A, w1 B! Ffigure(1),clf
. F( {0 [; U0 A- {+ }- M. Z! Cplot3(x1,x2,y,'o')
+ [! U4 v2 R' ~2 Y4 j" F7 Kstem3(x1,x2,y,'filled'): A; W* L1 ?+ f1 b! e' ~' p. F- _, ^
text(x1,x2,y+.01,str)0 A' v( D; J! ~
hold on
& ?, Y: i) ~. S3 c[x11,x22]=meshgrid(min(x1):range(x1)/80:max(x1),min(x2):range(x2)/80:max(x2));
3 T- Z* s! j4 _yhat=fx1(b,x11,x22);
/ U E* R3 e9 BsuRF(x11,x22,yhat)1 I# {- A' v# D/ N& w4 T* ?
shading interp. G% J) |* K$ e6 U% m
alpha(.8)0 M8 R$ N8 f9 x7 D# b* T$ O/ R
axis tight
% T2 o3 h* ]* z6 ^. rxlabel('α');ylabel('C'),zlabel('V1')9 j5 `/ U* K5 ?+ W/ F
SSy=var(y)*(n-1)$ d' W5 o& Y& O
y1=fx1(b,x1,x2);
! |! f& i& O5 q0 h( Y' JRSS=(y-y1)'*(y-y1). [6 B* {3 S& k; ~. d5 D1 J
rsquare=(SSy-RSS)/SSy1 ]5 v1 v2 W, `# Z1 }
MSe=RSS/(n-length(b))% l9 B" D, \/ p" z# U, H: b) N+ R! O5 [
# S9 f- u9 a/ n2 n: W) S |
|