mpfit: multivariate fit [message #53931] |
Mon, 07 May 2007 21:19  |
Dave[3]
Messages: 12 Registered: November 2006
|
Junior Member |
|
|
Hi all,
I'm trying to use MPFIT to numerically estimate a coordinate
transformation matrix that relates two sensors. One set is
uncalibrated, and the other has a known calibration. So, I have a set
of observed vectors (xyz_obs) and a set of known vectors (xyz_known)
and I'm trying to estimate (in a least squares sense) the
transformation matrix T that relates them. Judiciously choosing the
data in the fictional example below, I expect the transformation
matrix to be 2 * identity(3,3).
When I execute the code below, I get:
IDL> .r foo
% Compiled module: TRANS.
% Compiled module: $MAIN$.
Iter 1 CHI-SQUARE = 1278958.0 DOF = 291
P(0) = 1.00000
P(1) = 0.00000
P(2) = 0.00000
P(3) = 0.00000
P(4) = 1.00000
P(5) = 0.00000
P(6) = 0.00000
P(7) = 0.00000
P(8) = 1.00000
% MPFIT: Error detected while calling MPFIT_FDJAC2:
% MPFIT: Out of range subscript encountered: FJAC.
% MPFIT: Error condition detected. Returning to MAIN level.
Any ideas on what I'm doing wrong here?
Thanks!
Dave
%%%
% Contents of foo.pro
%%%
function trans, K, X=x, Y=y, err=err, forward=fw
model = K ## x
if keyword_set(fw) then return, model else return, (y-model)/err
end
; MAIN
; Attempt to estimate the transformation matrix given a set
; of observed Cartesian vectors and a set of known cartesian
; vectors.
n = 100 ; number of 'observations'
v = [1.0d, 0.15, 0.5] ; template vector
xyz_obs = dblarr(n, 3) ; observations
for i=0, n-1 do $
xyz_obs[i,*] = reform( v+0.01*randomn(seed,3), 1, 3)
xyz_known = dblarr(n, 3) ; known values (trivial scaling)
for i=0, n-1 do $
xyz_known[i,*] = reform( v*2.0, 1, 3)
; Estimate the transformation matrix, T
T0 = identity(3, /DOUBLE) ; initial guess transformation matrix
f = {x: xyz_obs, y: xyz_known, err: 0.01}
T = mpfit('trans', T0, functargs=f, COVAR=S2)
end
|
|
|
|
|
Re: mpfit: multivariate fit [message #54032 is a reply to message #53953] |
Sat, 12 May 2007 17:00   |
Craig Markwardt
Messages: 1869 Registered: November 1996
|
Senior Member |
|
|
hradilv <hradilv@yahoo.com> writes:
> On May 11, 3:05 am, Craig Markwardt
> <craigm...@REMOVEcow.physics.wisc.edu> wrote:
>>
>> David Fanning is right. You can't simply recode the function like you
>> did. According to the documentation:
>>
>
> Agreed - I was just posting to point out _where_ the problem was, not
> _what_ the problem was.
>
> I'm glad you cleared it up - and I'm very gratefule, Craig, for your
> excellent programs. I use them just about everyday!
OK. Thanks for your kind words!
Craig
--
------------------------------------------------------------ --------------
Craig B. Markwardt, Ph.D. EMAIL: craigmnet@REMOVEcow.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response
------------------------------------------------------------ --------------
|
|
|
Re: mpfit: multivariate fit [message #54105 is a reply to message #53931] |
Mon, 21 May 2007 16:52  |
Richard French
Messages: 173 Registered: December 2000
|
Senior Member |
|
|
Thanks, Craig -
I appreciate the clarification, and I'm glad the solution was so simple.
This probably deserves some clarification in the MPFIT documentation as
well, since it differs from the conventions of all other least squares
fitting routines I've used (such as the IMSL routine, or CURVEFI), and the
documentation in MPFIT.PRO for MYFUNCT seems to me to say the opposite of
what you did to fix the problem:
The confusion stems from just what the 'model' really means:
It states in the documentation that
Model=F(x,p)
dp(*,i) = FGRAD(x,p,i)
...
Return,(y-model)/err
Where FGRAD() is a user function which must compute the derivative of the
model with respect to the parameter P(i) at X..... DP(I,J) is the derivative
of the I'th point with respect to the J'th parameter."
Since MYFUNCT returns (y-model)/err, and dp is stated as being the
derivative array with respect to 'model' (not with respect to
(y-model)/err), it seems to me that the documentation is saying that the
'model' function with respect to which one should take derivatives is
model=F(x,P), not (y-F(x,P))/err.
Dick
On 5/21/07 2:13 AM, in article m28xbiembq.fsf@phloem.local, "Craig
Markwardt" <craigmnet@REMOVEcow.physics.wisc.edu> wrote:
>
> Richard, you sent me your code off-line.
>
> Yes, there was an obvious problem, which is that your derivatives were
> computed from d(MODEL)/d(P[i]).
>
> However, the "user" function for MPFIT is always (DATA-MODEL)/SIGMA [*].
> This means that the derivative is actually [ -d(MODEL)/d(P[i]) ]/SIGMA.
> When I inserted a negative sign in your derivatives, the fit then
> turned out fine.
>
> Craig
>
> [*] - This is true whether the user function is regular or "(EXTERNAL)".
>
> "Richard G. French" <rfrench@wellesley.edu> writes:
>> On the subject of MPFIT, I'm trying to implement the '(EXTERNAL)' function
>> option (the case where the user does not supply a function to be called by
>> MPFIT itself, but is expected to supply function values and a jacobian
>> matrix in each call to MPFIT). I'm trying to do this with a very simple
>> polynomial function for which the partial derivatives are very easy to
>> compute, but the fitted parameters are not the correct values, the number of
>> degrees of freedom returned is a negative number, and in general I am
>> clearly not getting things to work.
>>
>> In trying to debug this, I found what looks like an undocumented keyword:
>> EXTERNAL_INIT whose purpose is a bit obscure to me. I tried fiddling with
>> this but it did not fix the problem
>>
>> Does anyone (Craig or someone else) have a simple example of MPFIT that uses
>> the EXTERNAL option? If so, I'd be grateful if you would send it along. If
>> not, then I can post my failing example so that someone can tell me what I
>> am doing wrong.
>>
>> Thanks very much.
>>
>> Dick French
>> rfrench@REMOVEwellesley.edu
>>
|
|
|
Re: mpfit: multivariate fit [message #54124 is a reply to message #53931] |
Sun, 20 May 2007 23:13  |
Craig Markwardt
Messages: 1869 Registered: November 1996
|
Senior Member |
|
|
Richard, you sent me your code off-line.
Yes, there was an obvious problem, which is that your derivatives were
computed from d(MODEL)/d(P[i]).
However, the "user" function for MPFIT is always (DATA-MODEL)/SIGMA [*].
This means that the derivative is actually [ -d(MODEL)/d(P[i]) ]/SIGMA.
When I inserted a negative sign in your derivatives, the fit then
turned out fine.
Craig
[*] - This is true whether the user function is regular or "(EXTERNAL)".
"Richard G. French" <rfrench@wellesley.edu> writes:
> On the subject of MPFIT, I'm trying to implement the '(EXTERNAL)' function
> option (the case where the user does not supply a function to be called by
> MPFIT itself, but is expected to supply function values and a jacobian
> matrix in each call to MPFIT). I'm trying to do this with a very simple
> polynomial function for which the partial derivatives are very easy to
> compute, but the fitted parameters are not the correct values, the number of
> degrees of freedom returned is a negative number, and in general I am
> clearly not getting things to work.
>
> In trying to debug this, I found what looks like an undocumented keyword:
> EXTERNAL_INIT whose purpose is a bit obscure to me. I tried fiddling with
> this but it did not fix the problem
>
> Does anyone (Craig or someone else) have a simple example of MPFIT that uses
> the EXTERNAL option? If so, I'd be grateful if you would send it along. If
> not, then I can post my failing example so that someone can tell me what I
> am doing wrong.
>
> Thanks very much.
>
> Dick French
> rfrench@REMOVEwellesley.edu
>
--
------------------------------------------------------------ --------------
Craig B. Markwardt, Ph.D. EMAIL: craigmnet@REMOVEcow.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response
------------------------------------------------------------ --------------
|
|
|
Re: mpfit: multivariate fit [message #54126 is a reply to message #54032] |
Sun, 20 May 2007 20:25  |
Richard French
Messages: 173 Registered: December 2000
|
Senior Member |
|
|
On the subject of MPFIT, I'm trying to implement the '(EXTERNAL)' function
option (the case where the user does not supply a function to be called by
MPFIT itself, but is expected to supply function values and a jacobian
matrix in each call to MPFIT). I'm trying to do this with a very simple
polynomial function for which the partial derivatives are very easy to
compute, but the fitted parameters are not the correct values, the number of
degrees of freedom returned is a negative number, and in general I am
clearly not getting things to work.
In trying to debug this, I found what looks like an undocumented keyword:
EXTERNAL_INIT whose purpose is a bit obscure to me. I tried fiddling with
this but it did not fix the problem
Does anyone (Craig or someone else) have a simple example of MPFIT that uses
the EXTERNAL option? If so, I'd be grateful if you would send it along. If
not, then I can post my failing example so that someone can tell me what I
am doing wrong.
Thanks very much.
Dick French
rfrench@REMOVEwellesley.edu
|
|
|