constraining parameters in multi-Gaussian 1D fitting

[Jess]

Hi,

I am attempting to do multi-Gaussian fitting of 1D profiles using Craig
Markwardt's mpfitfun, and the simple multi_gauss script kindly offered
in earlier advice to Ben Tupper on non-interactive multi Gaussian
fitting. I believe I should be able to make it work, as for each
profile I know the number of gaussians involved, 3-parameter Gaussians
are sufficient, and each of those parameters are fixed or heavily
constrained. I am using a starting parameter array P of size P =
fltarr(n_gauss), and parinfo to apply my constraints.

One constraint I am unable yet to do is: I = would like to be able to
tie the peak flux of the Gaussians such that the peak flux of last
Gaussian is always greater than that of the first Gaussian.
I tried using
parinfo((n_gauss-1)*3).tied = 'GT P[0]'

However the tied structure of parinfo doesn't seem to be meant to
accept operators like GT,LT, etc. The alternative is to use the limits
structure and saying that the lower limit of peakflux on last gaussian
must be greater than upper limit of peak flux on first gaussian with:
parinfo((n_gauss-1)*3).limited(0) = 1
parinfo((n_gauss-1)*3).limits(0) = parinfo(0).limits(1)

However this requires assigning rather tight bounds to P[0] which I
really dont know well. Is there a smarter way I can do this using
'tied' or another structure in parinfo?

Thanks,
Jess

[Craig Markwardt]

Post by Jess
One constraint I am unable yet to do is: I = would like to be able to
tie the peak flux of the Gaussians such that the peak flux of last
Gaussian is always greater than that of the first Gaussian.
I tried using
parinfo((n_gauss-1)*3).tied = 'GT P[0]'
However the tied structure of parinfo doesn't seem to be meant to
accept operators like GT,LT, etc. ...

True. MPFIT's TIED fields are limited to equality constraints only.

...

Post by Jess
However this requires assigning rather tight bounds to P[0] which I
really dont know well. Is there a smarter way I can do this using
'tied' or another structure in parinfo?

One approach is to fit one gaussian at a time starting with the
strongest one and working your way down to the fainter ones. This
assumes that the peaks are well enough separated that they can be fit
in succession. If the peaks are blended, well, that's a tough
situation.

Long ago I had somebody trying to do something similar, i.e. fit an
arbitrary number of gaussians with arbitrary centroids, widths and
amplitudes. I warned him that I thought it would turn into an
unconstrained mess, it did, and I think he eventually gave up. Sorry,
but least squares fitting is not mind reading. The more external
constraints you can apply, based on knowledge of the problem, the
better.

Craig

[Jess]

Hi Craig (and everyone else),

Thanks for clarifying the use of "tied". In that case I will have to
use limits. Perhaps I can re-fit a velocity profile altering my parinfo
limits, if the first multi-Gaussian fit of that profile doesnt yield
the right relationship between the peak fluxes of each Gaussian.

Unfortunately I can't sequentially fit single Gaussians, as the
velcoity profile is a sum of overlapping Gaussians. Each profile
comprises emission from a line-of-sight through a edge-on galaxy disk
in cylindrical rotation V(R), with gas distribution assumed to be
azimuthally symmetric, monotonicly decreasing radial flux profile,
isothermal and with isotropic velocity dispersion.

Thus I assume to know the following constraints:
- number of gaussians in any profile
- "good" parameter estimates of all but one gaussian in any profile.
- The centroid of the unknown gaussian must be higher than all the
others.
- The peak flux of the unknown gaussian should be equal or higher than
all the others.
- For the other Gaussians in a profile, each comprising flux from a
radial bin higher than R at the line of nodes, the parameters held
fixed or strongly limited from having already fitted flux at those
radii when analysing velocity profile of sightlines further out. I
start at the outermost sightline, where I am only fittin gas at Rmax,
then fit the next innermost using the kinematics at the higher R
sightline as parameter estimates.
- In practise, the centroids of higher R gaussians can be projected to
the new sightline and held fixed as they are found accurately in the
outer sightline, while the peak flux and velocity dispersion are
vulnerable to noise and biased by the telescope spatial beam at the
outermost sightlines, so must be bounded for the multi-gaussian fit of
the next few sightlines, until they stabilise.

I will try to generate reasonable bounding limits to supply for these
high R gaussians, such that the output obeys reasonable assumptions of
the input galaxy:
- that radial flux profile is flat or decreasing,
- the centroids of the gaussians at lower radii in line-of-sight
profile is always higher than the Gaussian centroids from gas at higher
radii projected onto the line of sight.

If this works I should find that the 3 parameter fits of higher radius
Gaussian gas components in a sightline profile should stabilise after a
few sightline profiles, as I fit that gas component in successive
sightlines. Thus the multi-gaussian fit to each profile should have few
unknowns, 3 due to the unknown gaussian at the line-of-nodes radius,
the dispersion and peak flux of the next couple of gaussians, say 3-4.
So I never have more than 8-11 unknowns, though less in the outermost
profiles with fewer gaussians. One advantage of an exponential
decreasing radial flux profile is that the poor fits to the highest r
components, while they affect the velocity profile fit of the next
innermost sightline, they become unimportant, providing its fitting the
line-of-node gaussian well.

What do you think Craig? Is the solution too degenerate to be solved
with curvefit algorithms like mpfit, despite the heavy constraints I am
trying to place on most of the parameters in the multi-gaussian curve
fit?

Thanks for your advice,
Jess

[Craig Markwardt]

Post by Jess
Hi Craig (and everyone else),
Thanks for clarifying the use of "tied". In that case I will have to
use limits. Perhaps I can re-fit a velocity profile altering my parinfo
limits, if the first multi-Gaussian fit of that profile doesnt yield
the right relationship between the peak fluxes of each Gaussian.
Unfortunately I can't sequentially fit single Gaussians, as the
velcoity profile is a sum of overlapping Gaussians. Each profile
comprises emission from a line-of-sight through a edge-on galaxy disk
in cylindrical rotation V(R), with gas distribution assumed to be
azimuthally symmetric, monotonicly decreasing radial flux profile,
isothermal and with isotropic velocity dispersion.

... extensive description deleted ...

Post by Jess
What do you think Craig? Is the solution too degenerate to be solved
with curvefit algorithms like mpfit, despite the heavy constraints I am
trying to place on most of the parameters in the multi-gaussian curve
fit?

Sorry, I got a bit lost in that description.

If possible, I would try to fit the "new" peak while holding the other
peaks absolutely constant, and then only after the "new" peak fit has
converged, allow more of the parameters to vary. You want to be as
close as possible to the optimal solution before you allow the most
parameters to vary.

Happy fitting,
Craig

[JD Smith]

Post by Craig Markwardt

Post by Jess
One constraint I am unable yet to do is: I = would like to be able to
tie the peak flux of the Gaussians such that the peak flux of last
Gaussian is always greater than that of the first Gaussian.
I tried using
parinfo((n_gauss-1)*3).tied = 'GT P[0]'
However the tied structure of parinfo doesn't seem to be meant to
accept operators like GT,LT, etc. ...

True. MPFIT's TIED fields are limited to equality constraints only.

What if you availed yourself of the ITERPROC procedure to enforce the
constraint, dragging the fitter (kicking and screaming if necessary)
back into line if it attempts to step out? Any reason this wouldn't
work?

JD

[Craig Markwardt]

Post by JD Smith

Post by Craig Markwardt

Post by Jess
One constraint I am unable yet to do is: I = would like to be able to
tie the peak flux of the Gaussians such that the peak flux of last
Gaussian is always greater than that of the first Gaussian.
I tried using
parinfo((n_gauss-1)*3).tied = 'GT P[0]'
However the tied structure of parinfo doesn't seem to be meant to
accept operators like GT,LT, etc. ...

True. MPFIT's TIED fields are limited to equality constraints only.

What if you availed yourself of the ITERPROC procedure to enforce the
constraint, dragging the fitter (kicking and screaming if necessary)
back into line if it attempts to step out? Any reason this wouldn't
work?

It might work, it might not. I suspect that in general the fitter
might get stuck. For example, if we *started* out by thinking the
tallest peak was on the left part of the curve - and made a model
function to match that - but the truth is that a peak on the right is
truly he tallest, then we can enforce all the constraints in the world
and won't come out with the best fit.

Now, with the additional info that the original poster provided, this
will probably not be the case, so s/he might be alright.

Craig