comp.lang.idl-pvwave archive
Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado

Home » Public Forums » archive » mpfit: multivariate fit
Show: Today's Messages :: Show Polls :: Message Navigator
E-mail to friend 
Switch to threaded view of this topic Create a new topic Submit Reply
mpfit: multivariate fit [message #53931] Mon, 07 May 2007 21:19 Go to next message
Dave[3] is currently offline  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 #53952 is a reply to message #53931] Fri, 11 May 2007 06:45 Go to previous messageGo to next message
Vince Hradil is currently offline  Vince Hradil
Messages: 574
Registered: December 1999
Senior Member
On May 11, 8:44 am, hradilv <hrad...@yahoo.com> wrote:
> 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!

Umm... that's grateful, not gratefule...
Re: mpfit: multivariate fit [message #53953 is a reply to message #53931] Fri, 11 May 2007 06:44 Go to previous messageGo to next message
Vince Hradil is currently offline  Vince Hradil
Messages: 574
Registered: December 1999
Senior Member
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!
Re: mpfit: multivariate fit [message #54032 is a reply to message #53953] Sat, 12 May 2007 17:00 Go to previous messageGo to next message
Craig Markwardt is currently offline  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 Go to previous message
Richard French is currently offline  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 Go to previous message
Craig Markwardt is currently offline  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 Go to previous message
Richard French is currently offline  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
  Switch to threaded view of this topic Create a new topic Submit Reply
Previous Topic: Re: Idlgrsymbol size?
Next Topic: Re: write a root group to HDF5 file

-=] Back to Top [=-
[ Syndicate this forum (XML) ] [ RSS ] [ PDF ]

Current Time: Wed Oct 08 15:12:02 PDT 2025

Total time taken to generate the page: 0.00577 seconds