1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
|
function [p,err]=fit_dat(G,p0,Z,W,pmin,pmax,DG)
//
// [p,err]=fit_dat(G,p0,Z [,W] [,pmin,pmax] [,DG])
//
// Function used for fitting data to a model.
// For a given function G(p,z), this function finds the best vector
// of parameters p for approximating G(p,z_i)=0 for a set of measurement
// vectors z_i. Vector p is found by minimizing
// G(p,z_1)'WG(p,z_1)+G(p,z_2)'WG(p,z_2)+...+G(p,z_n)'WG(p,z_n)
//
// G: Scilab function (e=G(p,z), e: nex1, p: npx1, z: nzx1)
// p0: initial guess (size npx1)
// Z: matrix [z_1,z_2,...z_n] where z_i (nzx1) is the ith measurement
// W: weighting matrix of size nexne (optional; default 1)
// pmin: lower bound on p (optional; size npx1)
// pmax: upper bound on p (optional; size npx1)
// DG: partial of G wrt p (optional; S=DG(p,z), S: nexnp)
//
// Examples
//
//deff('y=FF(x)','y=a*(x-b)+c*x.*x')
//X=[];Y=[];
//a=34;b=12;c=14;for x=0:.1:3, Y=[Y,FF(x)+100*(rand-.5)];X=[X,x];end
//Z=[Y;X];
//deff('e=G(p,z)','a=p(1),b=p(2),c=p(3),y=z(1),x=z(2),e=y-FF(x)')
//[p,err]=fit_dat(G,[3;5;10],Z)
//xset('window',0)
//xbasc();
//plot2d(X',Y',-1)
//plot2d(X',FF(X)',5,'002')
//a=p(1),b=p(2),c=p(3);plot2d(X',FF(X)',12,'002')
//
//a=34;b=12;c=14;
//deff('s=DG(p,z)','y=z(1),x=z(2),s=-[x-p(2),-p(1),x*x]')
//[p,err]=fit_dat(G,[3;5;10],Z,DG)
//xset('window',1)
//xbasc();
//plot2d(X',Y',-1)
//plot2d(X',FF(X)',5,'002')
//a=p(1),b=p(2),c=p(3);plot2d(X',FF(X)',12,'002')
//
//
// Copyright INRIA
[lhs,rhs]=argn(0)
boun=%f
if rhs==3 then
W=1;GR=%f //Jacobian not provided
elseif rhs==4 then
if type(W)==11|type(W)==13 then
DG=W;W=1;GR=%t //Jacobian provided
else
GR=%f //Jacobian not provided
end
elseif rhs==5 then
if type(pmin)==1 then
pmax=pmin,pmin=W,W=1,GR=%f,boun=%t
else
DG=pmin;GR=%t
end
elseif rhs==6 then
boun=%t
if type(pmax)==1 then
GR=%f
else
DG=pmax,pmax=pmin,pmin=W,W=1,GR=%t
end
elseif rhs==7 then
GR=%t;
else
error('wrong number of arguments')
end
if size(W,1)~=size(W,2) then error('Weighting matrix must be square');end
nz=size(Z,2);mz=size(p0,'*');
deff('e=GG(p,Z)','e=0;for i=1:nz,e=e+G(p,Z(:,i))''*W*G(p,Z(:,i)),end')
if ~GR then
deff('g=DGG(ps,Z)',['g=[]';
'nrm=norm(ps);if nrm<1 then nrm=1,end';
'[d,nu]=colcomp(rand(mz,mz));d=10*%eps*nrm*d';
'for v=d,e=GG(ps+v),g=[g,e-f],end;g=g/d;'])
else
deff('g=DGG(ps,Z)','g=0;for i=1:nz,g=g+2*DG(ps,Z(:,i))''*W*G(ps,Z(:,i)),end')
end
deff('[f,g,ind]=costf(p,ind)',['f=GG(p,Z);';
'g=DGG(p,Z);';
'disp(''error= ''+string(f));'])
if ~boun then
[err,p]=optim(costf,p0)
else
if or(p0<pmin)|or(p0>pmax) then
error('initial guess not feasible')
end
[err,p]=optim(costf,'b',pmin,pmax,p0)
end
|
