Fit with constraints within a procedure

Share icon

Tue, 03/08/2016 - 03:38 am

Hi
I have plenty of data to fit with a power law.
I read on the forum that the built-in fonctions can't be constrained, so I write a user-defined function and I try to introduce constraints. But it doesn't work.

Without constraints, the macro is OK, and performs well.

But when i add the final command /C=const, the procedure is stopped. I certainly did something wrong, but ....
Thanks for any advice.
and Have a nice day !

Hélène


Here is the code (I also attached the procedure version)

macro traittspabsCbaselinePower(namebase,namemask,numdep,numfin)
string namebase,namemask
variable numdep,numfin
Prompt namebase,"Nom de base du spectre à analyser. Ex Abs_Cer_T"
Prompt namemask,"Nom de la wave mask"
Prompt numdep,"n° wave depart"
Prompt numfin,"n° wave fin"

silent 1
variable curA,curB
variable i=numdep
string namespectre,nomgraph
string namespectrecor,nameBLref
nomgraph=namebase+"_BLpow"
Const="contraintes"
Make/O/T/N=4 Const
Make/O/D/N=3 W_coef
W_coef={0,2e-5,11}

Const= {"K0 > -2e-3","K0 < 2e-3","K1 > 0","K2 > 0"}

// on va prendre le spectre 1 comme réferérence et définir la ligne de base de ref avec.
i=1
namespectre=namebase+num2str(i)
namespectrecor=namespectre+"_corBLpow"
nameBLref="BLref"
duplicate/o $namespectre $namespectrecor,$nameBLref
FuncFit/NTHR=0/TBOX=0 FctPower, W_coef , $namespectre /M=$namemask /D /C=Const
//CurveFit/NTHR=0/TBOX=0 poly 4, $namespectre /M=$namemask
$nameBLref=(W_coef[0]+W_coef[1]*x^W_coef[2])
WaveStats/Q/R=(xcsr(A),xcsr(B)) $namespectrecor
$namespectrecor-=V_avg
 FitPower.ipf
PASQUIER wrote:
I read on the forum that the built-in fonctions can't be constrained

Can you provide a link? It is not true, mostly. Linear fits (line, poly and poly2D) cannot be constrained, but a Power fit can.
Quote:
//CurveFit/NTHR=0/TBOX=0 poly 4, $namespectre /M=$namemask

I understand now- "Power" refers to a function like y = a*x^b. You wish to fit a polynomial with constraints. For that, indeed, you need a user-defined fit function.

Well, I just looked at the code in your procedure file:
Quote:
return w[0]+w[1]*x^(w[2])

That is a Power function. You can use the built-in Power function and still apply constraints.
Quote:
so I write a user-defined function and I try to introduce constraints. But it doesn't work.

Without constraints, the macro is OK, and performs well.

But when i add the final command /C=const, the procedure is stopped. I certainly did something wrong, but ....

When you say, "the procedure is stopped," what do you mean? Does it exit early? Is there an error message?
My guess is that the problem is that you have set the constraint on K1 and K2 to exactly zero. If the constraint becomes active, that will result in a Singular Matrix error. You need to choose some positive, non-zero number that is still "small" within the context of your problem. So you might use
Const= {"K0 > -2e-3","K0 < 2e-3","K1 > 1e-10","K2 > 0.01"}
I chose those numbers based on your initialization of the coefficient wave. And that points to another possible problem:
Quote:
W_coef={0,2e-5,11}

The initialization for w_coef[2] represents a *very* large exponent. That is bound to cause numerical instability.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

March 8, 2016 at 10:05 am - Permalink

Dear John

many thanks for your help.
You are right : the problem comes from the constraints values. Now it is working.

Concerning my sentence " I read on the forum that the built-in fonctions can't be constrained", in fact I have too briefly summarised your remark posted in a previous topic : http://www.igorexchange.com/node/3076
In fact in this topic you specified that only built-in Linear fits (line, poly and poly2D) can't be constrained. Mea Culpa !

Thank you again for your help

Hélène

March 9, 2016 at 12:57 am - Permalink