Here's my question: How does one _implement_ constrained fitting parameters?
I know different aspects of this have been addressed before in a few
places, but I have not been able to put the pieces together
successfully. Can anyone help me out?
I have, of course, found posts such as
>> snip
$a = 12.4 + 5*sin(~0)
should give variable that is bounded by 12.4+-5.
>> snip
but I'm not sure where or when such a variable should be defined, or
how to use this constrained variable as a fitting parameter.
This is what I'm trying to do:
Fit a peak using several functions (e.g. two Lorentzians and one user-defined)
Constrain the positions of every peak to fall within a certain range
Constrain the peak heights to be >= zero
So far I have tried to write some scripts defining each peak, and
restricting the parameters (height & center) of each respective peak
in different ways to try to find what works. In many cases a script
loads with no errors, but the resulting parameters still end up
outside of the defined range, so I'm obviously doing something wrong.
Here is an example showing a few of the things I have tried.
#---------- myscript.fit----------
# Define variable for center of 1st Lorentzian
$c1 = 1592+3*sin(~0)
#
# Define functions, with restrictions on peak positions
define Lorentz1(height, c1, hwhm) = height/(1+((x-c1)/hwhm)^2) #try
using pre-defined variable for fitting
define Lorentz2(height, center=1554+3*sin(~0), hwhm) =
height/(1+((x-center)/hwhm)^2) #try restricting normal fitting
parameter 'center'
define Userdef(height, center=1575+5*sin(~0), fwhm, qinv=-0.2) =
height*((1+(x-center)*qinv/(fwhm))^2)/(1+((x-center)/fwhm)^2)
#
# Ensure all peak heights are non-negative (and try imposing
constraints after function definitions)
Lorentz1.height=(~{sqrt(%Lorentz1.height)})^2 #do I need the % ?
Lorentz2.height=(~{sqrt(Lorentz2.height)})^2 #or not?
Userdef.height=(~{sqrt(1000)})^2 # can I just start with an arbitrary
value that, when varied, will remain positive?
#----------end----------
Any help or advice is greatly appreciated.
Many thanks in advance,
Erik
On Tue, May 27, 2008 at 5:26 PM, Erik Einarsson
<er...@photon.t.u-tokyo.ac.jp> wrote:
>
> I have, of course, found posts such as
>>> snip
> $a = 12.4 + 5*sin(~0)
> should give variable that is bounded by 12.4+-5.
>>> snip
> but I'm not sure where or when such a variable should be defined, or
> how to use this constrained variable as a fitting parameter.
You can use the bounded variable as a parameter of peak, e.g.:
guess Lorentzian center=$a
You can also define a Loretzian-like function that has height >= 0:
define LorentzianPositive(sqrt_height=sqrt(height), center, hwhm) =
Lorentzian(sqrt_height^2, center, hwhm)
or even define a function like this:
define Lx(sqrt_height=sqrt(height), center, cx=0, hwhm) =
Lorentzian(sqrt_height^2, {center}+2*sin(cx), hwhm)
Marcin
--
Marcin Wojdyr | http://www.unipress.waw.pl/~wojdyr/
>
>You can use the bounded variable as a parameter of peak, e.g.:
>guess Lorentzian center=$a
>
I tried this, but often got an error because there was no maximum
within the specified range. I want to define the peak locations
because I'm trying to fit what appears in the spectrum as two broad
peaks with a total of five peaks, thus some must be small, and to the
side of the main peak.
>
>or even define a function like this:
>define Lx(sqrt_height=sqrt(height), center, cx=0, hwhm) =
>Lorentzian(sqrt_height^2, {center}+2*sin(cx), hwhm)
>
I think this is close to what I'm looking for, but doesn't this define
only the constraint box while still leaving the center free?
-Erik
I don't know how you add peaks, but you can always set $a to be a center.
>>or even define a function like this:
>>define Lx(sqrt_height=sqrt(height), center, cx=0, hwhm) =
>>Lorentzian(sqrt_height^2, {center}+2*sin(cx), hwhm)
>>
> I think this is close to what I'm looking for, but doesn't this define
> only the constraint box while still leaving the center free?
Yes, the center is free in this box.
>define Lx(sqrt_height=sqrt(height), center, cx=0, hwhm) =
>Lorentzian(sqrt_height^2, {center}+2*sin(cx), hwhm)
>
I just replace "{center}" with a number, around which the peak
position is constrained by (in this case) 2.
Many thanks!
I think this shouldn't be necessary. The center parameter won't change.
If for some reasons it doesn't work, you may remove {} from the
definition and explicitly use constant number as a parameter (e.g. 2,
not ~2).
Using different definition for every peak will be tedious.