Expression

[ZIHABAKE Daniel]

Hello everyone,

i'm modeling emission lines by setting expression to relate one amplitude to the other.

since the flux of one line is equal the  third of the other, i am expressing amplitude in terms of flux. But when i use quad in integrating it brings an error of not being defined yet it is imported earlier at the start.

How can i handle this error?

the attached is the data.

thank you

Daniel.

from numpy import loadtxt
import numpy as np
import matplotlib.pyplot as plt
from lmfit import Parameters
from lmfit.models import GaussianModel, LinearModel
from scipy.integrate import quad
import math as mt

params={'legend.fontsize':'18','axes.labelsize':'18',
'axes.titlesize':'18','xtick.labelsize':'18',
'ytick.labelsize':'18','lines.linewidth':2,'axes.linewidth':2,'animation.html': 'html5'}
plt.rcParams.update(params)
plt.rcParams.update({'figure.max_open_warning': 0})
plt.rcParams['figure.figsize'] = (12,8)


data = loadtxt('Data.txt')
x = data[:, 0]
y = data[:,1]

#selected region for fitting
mask = (x>= 6680) & (x <= 6770)
mask = (x>= 6400) & (x <= 6800)
xarray, yarray = x[mask], y[mask]
x=xarray

y=yarray
#  Note: use builtin models
gmodel = ( GaussianModel(prefix='g2_')+GaussianModel(prefix='g1_')+GaussianModel(prefix='g4_')+GaussianModel(prefix='g3_')+GaussianModel(prefix='g7_')+LinearModel(prefix='off_')+GaussianModel(prefix='g8_'))

#  Note: many fewer bounds:
params = Parameters()

params.add('off_intercept', value=1)
params.add('off_slope', value=0)
params.add('g2_amplitude', value=3, vary= True)
params.add('g2_center', value=6731, min=6725, max=6740)
params.add('g2_sigma',  value= 3.6021,min=0)
params.add('g1_amplitude', value=1.2, vary=True)
params.add('g1_center', value=6719, min=6710, max=6724)
params.add('g1_sigma', expr='g2_sigma')

params.add('g3_amplitude', value=1.6)
params.add('g3_center', value=6585, min=6570, max=6595)
params.add('g3_sigma', expr='g2_sigma')

params.add('g7_amplitude', value=5)
params.add('g7_center', value=6567, min=6555, max=6575)
params.add('g7_sigma', expr='g2_sigma')

params.add('g4_center', value=6549, min=6540, max=6555)
params.add('g4_sigma', expr='g2_sigma')
params.add('g4_amplitude', expr='quad(1/(g3_sigma*np.sqrt(2*np.pi))*mt.exp(-(x-g3_center)**2/2*g3_sigma**2), 6500, 6650)/quad(3*mt.exp(-(x-g4_center)**2/2*g4_sigma**2), 6500, 6650)')
params.add('g8_amplitude', value=7, min=1)
params.add('g8_center', value=6565, min=6554, max=6575)
params.add('g8_sigma', value=7, min=0)


result = gmodel.fit(y, params, x=x)
print(result.fit_report())

comps = result.eval_components(x=x)

fig, axes = plt.subplots(2, 1, gridspec_kw={'hspace': 0, 'height_ratios': [3,1]}, figsize=(12, 8))

axes[0].set_title('MCG_00_29_023')
axes[0].plot(x, y, 'k', label='observed')
axes[0].set_ylim(ymin=-0.1, ymax=2)
axes[0].plot(x, result.best_fit, 'r--', label='best fit')
axes[0].plot(x, comps['g1_']+comps['off_'], 'g-', label='g1')
axes[0].plot(x, comps['g2_']+comps['off_'], 'g-', label='g2')

axes[0].plot(x, comps['g3_']+comps['off_'], 'g-', label='g3')
axes[0].plot(x, comps['g4_']+comps['off_'], 'g-', label='g4')
axes[0].plot(x, comps['g7_']+comps['off_'], 'g-', label='g7')
axes[0].plot(x, comps['g8_']+comps['off_'], 'b-', label='g8')
axes[0].plot(x, comps['off_'], 'k-', label='offset')

residual=result.best_fit-y
axes[1].set_ylim(ymin=-0.15, ymax=0.15)
axes[1].plot(x,residual)
axes[1].set_xlabel(r'wavelength [$\rm \AA$]')
axes[0].set_ylabel(r'flux [$\rm erg / (cm^{2} \, s \, \AA)$]')
plt.savefig('Ha_Data.pdf')
plt.show()

[Matt Newville]

Hi Daniel, 

On Thu, Aug 4, 2022 at 4:23 AM ZIHABAKE Daniel <dziha...@gmail.com> wrote:

Hello everyone,

i'm modeling emission lines by setting expression to relate one amplitude to the other.

since the flux of one line is equal the  third of the other, i am expressing amplitude in terms of flux. But when i use quad in integrating it brings an error of not being defined yet it is imported earlier at the start.

How can i handle this error?

Wouldn't you just say that 

   params.add('g4_amplitude', expr='g3_amplitude/3')

?

Amplitude *is* the strength of that unit-normalized Gaussian: the area of the curve.  

Otherwise, you would need to add the `quad` function (FWIW, `exp` is already available) into the `asteval` interpreter used for the Parameter expressions.  But I think you do not need to do that. 

--Matt 

[ZIHABAKE Daniel]

Thank you Matt for the explanation,

About Area we tried to compare amplitude but we found it not giving correct measurements.

Let me try with integration.

Thank you once again

Daniel.

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CA%2B7ESbotPn0dcbcTqWdc5ZX%2BK0W4WGs93ichzV_-OH%2BHvdnkUw%40mail.gmail.com.

[Matt Newville]

On Thu, Aug 4, 2022 at 11:55 AM ZIHABAKE Daniel <dziha...@gmail.com> wrote:

Thank you Matt for the explanation,

About Area we tried to compare amplitude but we found it not giving correct measurements.

Hmm, that seems odd. To be clear, "amplitude" is the factor multiplying a unit-normalized Gaussian.  It will be the area under the curve.

Let me try with integration.

That would definitely take some effort.  In the end, it should give the same result, to numerical precision.   Well, maybe except for the fact that your formula would not give the ratio of the area of two Gaussian. Maybe these are "typos", but you left out both amplitudes and are not normalizing consistently.  You would also be integrating over a limited range, but if sigma is close to your initial values, that would be +/-5*sigma or more, so it might not matter.  

But, really, "amplitude" will be the area under each Gaussian.  That's built into the definition. 

--Matt 

[ZIHABAKE Daniel]

Dear Matt,

On Fri, 5 Aug 2022 at 04:36, Matt Newville <newv...@cars.uchicago.edu> wrote:

On Thu, Aug 4, 2022 at 11:55 AM ZIHABAKE Daniel <dziha...@gmail.com> wrote:

Thank you Matt for the explanation,

About Area we tried to compare amplitude but we found it not giving correct measurements.

Hmm, that seems odd. To be clear, "amplitude" is the factor multiplying a unit-normalized Gaussian.  It will be the area under the curve.

The reason why we prefer integration instead of using amplitude is that "two lines may have the same amplitude but with different widths. Therefore , they will also have different area".

Let me try with integration.

That would definitely take some effort.  In the end, it should give the same result, to numerical precision.   Well, maybe except for the fact that your formula would not give the ratio of the area of two Gaussian.

The ratio of the area of two Gaussian in that formula is three, and from there i have to express one amplitude in terms of the rest parameters in that formula.

Maybe these are "typos", but you left out both amplitudes and are not normalizing consistently.  You would also be integrating over a limited range, but if sigma is close to your initial values, that would be +/-5*sigma or more, so it might not matter.  

the limit of integration is from 6500 to 6650.

But, really, "amplitude" will be the area under each Gaussian.  That's built into the definition. 

--Matt 

Than you .

Daniel

You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CA%2B7ESbp9zA%3DUBat2b0OETcc6PpE8ke8urB%2B5WF3%2B%3DDGJt-BF2A%40mail.gmail.com.

[Matt Newville]

On Fri, Aug 5, 2022 at 1:30 AM ZIHABAKE Daniel <dziha...@gmail.com> wrote:

Dear Matt,

On Fri, 5 Aug 2022 at 04:36, Matt Newville <newv...@cars.uchicago.edu> wrote:

On Thu, Aug 4, 2022 at 11:55 AM ZIHABAKE Daniel <dziha...@gmail.com> wrote:

Thank you Matt for the explanation,

About Area we tried to compare amplitude but we found it not giving correct measurements.

Hmm, that seems odd. To be clear, "amplitude" is the factor multiplying a unit-normalized Gaussian.  It will be the area under the curve.

The reason why we prefer integration instead of using amplitude is that "two lines may have the same amplitude but with different widths. Therefore , they will also have different area".

That is not correct.  Amplitude is the area of the curve, independent of the value of sigma. 

 

Let me try with integration.

That would definitely take some effort.  In the end, it should give the same result, to numerical precision.   Well, maybe except for the fact that your formula would not give the ratio of the area of two Gaussian.

The ratio of the area of two Gaussian in that formula is three, and from there i have to express one amplitude in terms of the rest parameters in that formula.

Um, no.  The integral over x of  

      [1/(sqrt(2*pi)*sigma)] * exp[-(x-center)**2 / (2*sigma**2)]

is 1. 

The formula you posted leaves out the "1/(sqrt(2*pi)*sigma)" when "integrating" peak 4.  The numerator (integral over peak 3) will be 1, and the denominator will be "sqrt(2*pi)*sigma/3" for peak 4.   That is not a sensible value for an amplitude. 

--Matt 

[ZIHABAKE Daniel]

Ooh, i understand you very well. The mathematics in my fitting is that "the ratio of the flux of the two Gaussian lines is three. Therefore, by equating 3 with that ratio quad(1/(g3_sigma*np.sqrt(2*np.pi))*mt.exp(-(x-g3_center)**2/2*g3_sigma**2), 6500, 6650)/quad((1/g4_sigma)*np.sqrt(2*np.pi)*mt.exp(-(x-g4_center)**2/2*g4_sigma**2), 6500, 6650) 

from this equation, i had to get amplitude 4 which is equal to 1/g4_sigma)*np. sqrt(2*np.pi), and the remaining expression should replace amplitude 4 in setting the initial parameters. The reason why three is appearing in expression.

the only variable is x others are initial guesses.

Simply the message is that we equated the two areas of Gaussian function using the ratio between which is 3. Then from there, we solved amplitude 4 which we wanted to set the initial parameter for but in terms of area, and that amplitude is that missing term in the denominator.

i hope my message is clear.

ratio and solved  " 

will be "sqrt(2*pi)*sigma/3" for peak 4.   That is not a sensible value for an amplitude. 

--Matt 

I appreciate the help.

Daniel

You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CA%2B7ESbqjrzPpa9uyLO4oN%3DExbM0t-OFw9W2rL%2BCLcqP99TybRQ%40mail.gmail.com.

[Matt Newville]

On Fri, Aug 5, 2022 at 8:30 AM ZIHABAKE Daniel <dziha...@gmail.com> wrote:

On Fri, 5 Aug 2022 at 15:39, Matt Newville <newv...@cars.uchicago.edu> wrote:

On Fri, Aug 5, 2022 at 1:30 AM ZIHABAKE Daniel <dziha...@gmail.com> wrote:

Dear Matt,

On Fri, 5 Aug 2022 at 04:36, Matt Newville <newv...@cars.uchicago.edu> wrote:

On Thu, Aug 4, 2022 at 11:55 AM ZIHABAKE Daniel <dziha...@gmail.com> wrote:

Thank you Matt for the explanation,

About Area we tried to compare amplitude but we found it not giving correct measurements.

Hmm, that seems odd. To be clear, "amplitude" is the factor multiplying a unit-normalized Gaussian.  It will be the area under the curve.

The reason why we prefer integration instead of using amplitude is that "two lines may have the same amplitude but with different widths. Therefore , they will also have different area".

That is not correct.  Amplitude is the area of the curve, independent of the value of sigma. 

 

Let me try with integration.

That would definitely take some effort.  In the end, it should give the same result, to numerical precision.   Well, maybe except for the fact that your formula would not give the ratio of the area of two Gaussian.

The ratio of the area of two Gaussian in that formula is three, and from there i have to express one amplitude in terms of the rest parameters in that formula.

Um, no.  The integral over x of  

      [1/(sqrt(2*pi)*sigma)] * exp[-(x-center)**2 / (2*sigma**2)]

is 1. 

The formula you posted leaves out the "1/(sqrt(2*pi)*sigma)" when "integrating" peak 4.  The numerator (integral over peak 3) will be 1, and the denominator

Ooh, i understand you very well. The mathematics in my fitting is that "the ratio of the flux of the two Gaussian lines is three. Therefore, by equating 3 with that ratio quad(1/(g3_sigma*np.sqrt(2*np.pi))*mt.exp(-(x-g3_center)**2/2*g3_sigma**2), 6500, 6650)/quad((1/g4_sigma)*np.sqrt(2*np.pi)*mt.exp(-(x-g4_center)**2/2*g4_sigma**2), 6500, 6650) 

from this equation, i had to get amplitude 4 which is equal to 1/g4_sigma)*np. sqrt(2*np.pi), and the remaining expression should replace amplitude 4 in setting the initial parameters. The reason why three is appearing in expression.

the only variable is x others are initial guesses.

Simply the message is that we equated the two areas of Gaussian function using the ratio between which is 3. Then from there, we solved amplitude 4 which we wanted to set the initial parameter for but in terms of area, and that amplitude is that missing term in the denominator.

i hope my message is clear.

Last time, I promise.  If you want the area of Gaussian #4 to be 1/3 the area of Gaussian #3, then

    params.add('g4_amplitude', expr='g3_amplitude/3')

is exactly what you want. If your expression uses sigma it is either reducible to an expression that does not use sigma or is wrong.

--Matt