More problems with Curvefit [message #35619] |
Mon, 30 June 2003 07:35  |
heather.williams
Messages: 3 Registered: July 2002
|
Junior Member |
|
|
Good afternoon, everyone. I'm having some problems using Curvefit (in
IDL 5.4) to fit my data. I've reviewed the messages which have already
been posted on this subject, and had a good look at the MPCURVEFIT
substitute, but am no wiser.
My code looks like this:
PRO data_fit
;Define the vectors of tabulated:
F18_x = FLTARR(12) & F18_x(0) = 9.055227833
F18_x(1) = 9.908886278 & F18_x(2) = 11.86860889
F18_x(3) = 13.281685 & F18_x(4) = 16.69834393
F18_x(5) = 19.52864256 & F18_x(6) = 23.17273836
F18_x(7) = 28.51793219 & F18_x(8) = 31.23624055
F18_x(9) = 33.53401408 & F18_x(10) = 38.12262897
F18_x(11) = 39.15701348
F18_y = FLTARR(12) & F18_y(0) = 0.108598707
F18_y(1) = 0.329883541 & F18_y(2) = 0.504690343
F18_y(3) = 0.685805013 & F18_y(4) = 0.780161321
F18_y(5) = 0.87284238 & F18_y(6) = 0.890067419
F18_y(7) = 0.907523914 & F18_y(8) = 0.98011631
F18_y(9) = 0.943832957 & F18_y(10) = 0.966238284
F18_y(11) = 1
X = FLTARR(12) & X(*) = F18_x(*) - F18_x(0)
Y = FLTARR(12) & Y(*) = F18_y(*) - F18_x(0)
;Define a vector of weights:
W = 1.0
;Provide an initial guess of the function's parameters:
A = [1.0, 1.0]
;Compute the parameters a0 and a1:
yfit = CURVEFIT(X, Y, W, A, SIGMA_A, FUNCTION_NAME = 'fit_funct')
;Print the parameters, which are returned in A:
PRINT, A
END
PRO fit_funct, X, A, F, PDER
F = (1.0 - EXP(-A[0] * X)) + (1.0 - EXP(-A[1]*X))
; PDER's column dimension is equal to the number of
; elements in xi and its row dimension is equal to
; the number of parameters in the function F:
pder = FLTARR(N_ELEMENTS(X), 2)
; Compute the partial derivatives with respect to
; a0 and place in the first row of PDER:
pder[*, 0] = A[0] * EXP(-A[0] * X)
pder[*, 1] = A[1] * EXP(-A[1] * X)
END
Which looks alright, if not particularly elegant, to me. However, when
I run it, I get this error message (which relates to the line
beginning y_fit =) :
% Operands of matrix multiply have incompatible dimensions: <FLOAT
Array[2, 12]>, <FLOAT Array[1, 2]>.
% Error occurred at: CURVEFIT 269
O:\Rsi\Idl54\lib\curvefit.pro
% DATA_FIT 21 H:\PhD
IDL\Progs\data_fit.pro
% $MAIN$
How do I avoid this error and get the fit to work?
Thanks for your help,
Heather Williams
PhD Student, Manchester PET Centre
Manchester, UK
|
|
|
|
Re: More problems with Curvefit [message #35706 is a reply to message #35619] |
Tue, 01 July 2003 04:54  |
heather.williams
Messages: 3 Registered: July 2002
|
Junior Member |
|
|
Thanks to everyone for your help - it's now working happily, but even
with the recommended changes, this function doesn't seem to match the
data! I shall go in search of a more appropriate model, but any
suggestions from the previous respondents would be welcome.
- Heather
|
|
|
Re: More problems with Curvefit [message #35711 is a reply to message #35619] |
Mon, 30 June 2003 11:43  |
Craig Markwardt
Messages: 1869 Registered: November 1996
|
Senior Member |
|
|
Paul van Delst <paul.vandelst@noaa.gov> writes:
> Try doing
> w = make_array(12, value = 1.0 )
> or
> w = replicate( 1.0, 12 )
> to actually give you a vector for w. The IDL documentation is a wee bit misleading here as
> the actual words say "For no weighting, set Weightsi = 1.0" where the "i" suffix is an
> indicator that Weights is an array but it's not entirely clear (to me at least.)
As a side note, MPCURVEFIT, and the other "MP" fitting programs can
accept either a scalar or a vector for the uncertainties or weights.
Either way it does the "right thing." And of course the MP fitting
programs do not necessarily require the user to compute function
derivatives.
To Heather Williams, the original poster, you should be aware that
your fitting function is degenerate, since it contains a linear
combination of the *same* basis function. This leads to a singular
normal matrix, and takes a canned routine like CURVEFIT on a trip to
la-la land. This effect is further exacerbated by the choice of
initial conditions (both coefficients equal). MPCURVEFIT is more
robust in this sense, since it will do a singular value decomposition
in the face of a singular matrix.
Happy fitting!
Craig
--
------------------------------------------------------------ --------------
Craig B. Markwardt, Ph.D. EMAIL: craigmnet@cow.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response
------------------------------------------------------------ --------------
|
|
|
Re: More problems with Curvefit [message #35715 is a reply to message #35619] |
Mon, 30 June 2003 10:17  |
the_cacc
Messages: 104 Registered: October 2001
|
Senior Member |
|
|
1) You should put 'fit_funct' at the top of the file so it compiles
first.
2) Change "W = 1.0" to "W = replicate(1.0,12)"
3) pder should be the partial deriv wrt the parameters A - ie.
pder[*,0] = X * EXP(-A[0]*X) and pder[*,1] = X * EXP(-A[1]*X).
4) Looking at your data (PLOT,x,y), it seems that the choice of
fitting function is unlikely to fit well... ever. You need to allow
for a constant offset - the value 9 looks "right". If you can be sure
it is 9, then simply change y to y+9 in the call to curvefit.
5) To see the fit add this line at the end of your code: PLOT,x,y &
fit_funct,x,a,f & OPLOT,x,f-9.
Hope this gives you somewhere to start. I recommend changing your fit
function to deal with the offset issue. More generally, I would also
recommend looking at AMOEBA for fitting a small number of parameters -
very stable, doesn't need derivatives and not fussy about
discontinuities etc.
Ciao.
|
|
|
Re: More problems with Curvefit [message #35716 is a reply to message #35619] |
Mon, 30 June 2003 09:19  |
Paul Van Delst[1]
Messages: 1157 Registered: April 2002
|
Senior Member |
|
|
Heather Williams wrote:
>
> Good afternoon, everyone. I'm having some problems using Curvefit (in
> IDL 5.4) to fit my data. I've reviewed the messages which have already
> been posted on this subject, and had a good look at the MPCURVEFIT
> substitute, but am no wiser.
>
> My code looks like this:
>
> PRO data_fit
>
> ;Define the vectors of tabulated:
>
> F18_x = FLTARR(12) & F18_x(0) = 9.055227833
> F18_x(1) = 9.908886278 & F18_x(2) = 11.86860889
> F18_x(3) = 13.281685 & F18_x(4) = 16.69834393
> F18_x(5) = 19.52864256 & F18_x(6) = 23.17273836
> F18_x(7) = 28.51793219 & F18_x(8) = 31.23624055
> F18_x(9) = 33.53401408 & F18_x(10) = 38.12262897
> F18_x(11) = 39.15701348
> F18_y = FLTARR(12) & F18_y(0) = 0.108598707
> F18_y(1) = 0.329883541 & F18_y(2) = 0.504690343
> F18_y(3) = 0.685805013 & F18_y(4) = 0.780161321
> F18_y(5) = 0.87284238 & F18_y(6) = 0.890067419
> F18_y(7) = 0.907523914 & F18_y(8) = 0.98011631
> F18_y(9) = 0.943832957 & F18_y(10) = 0.966238284
> F18_y(11) = 1
>
> X = FLTARR(12) & X(*) = F18_x(*) - F18_x(0)
> Y = FLTARR(12) & Y(*) = F18_y(*) - F18_x(0)
>
> ;Define a vector of weights:
> W = 1.0
Try doing
w = make_array(12, value = 1.0 )
or
w = replicate( 1.0, 12 )
to actually give you a vector for w. The IDL documentation is a wee bit misleading here as
the actual words say "For no weighting, set Weightsi = 1.0" where the "i" suffix is an
indicator that Weights is an array but it's not entirely clear (to me at least.)
This made your code work for me.....but I got the following:
% CURVEFIT: Failed to converge- CHISQ increasing without bound.
1.00000 1.00000
% Program caused arithmetic error: Floating overflow
Changing the line
yfit = CURVEFIT(X, Y, W, A, SIGMA_A, FUNCTION_NAME = 'fit_funct')
to use Craig's MPCURVEFIT
yfit = MPCURVEFIT(X, Y, W, A, SIGMA_A, FUNCTION_NAME = 'fit_funct')
gave me the following:
IDL> data_fit
% Compiled module: MPCURVEFIT.
% Compiled module: MPFIT.
Iter 1 CHI-SQUARE = 1215.4340 DOF = 10
P(0) = 1.00000
P(1) = 1.00000
Iter 2 CHI-SQUARE = 1212.4130 DOF = 10
P(0) = 0.887302
P(1) = 0.887302
Iter 3 CHI-SQUARE = 1208.7012 DOF = 10
P(0) = 0.776462
P(1) = 0.776462
Iter 4 CHI-SQUARE = 1203.9254 DOF = 10
P(0) = 0.665368
P(1) = 0.665368
Iter 5 CHI-SQUARE = 1189.1635 DOF = 10
P(0) = 0.448375
P(1) = 0.448375
Iter 6 CHI-SQUARE = 945.99725 DOF = 10
P(0) = 0.0271755
P(1) = 0.0271756
Iter 7 CHI-SQUARE = 440.82562 DOF = 10
P(0) = -0.0570307
P(1) = -0.0570309
-0.0570307 -0.0570309
% Program caused arithmetic error: Floating underflow
% Program caused arithmetic error: Floating overflow
% Program caused arithmetic error: Floating illegal operand
This output just *has* to be more meaningful wrt diagnostics than the curvefit output.
paulv
--
Paul van Delst
CIMSS @ NOAA/NCEP/EMC
Ph: (301)763-8000 x7748
Fax:(301)763-8545
|
|
|