Combined model behavior and VoigtModel

[alexserebryanskiy]

Hello,

I'm trying to use BayesicFitting in spectrum analysis. The model contains several (many) Voigt profiles and coded as VoigtModel() + VoigtModel() + .... The model also has one more model of Astropy model type. But it seems that combined model recognize only one VoitModel + Astropy Model. Moreover, I cannot specify parameters for VoigtModel (although in help there's a statement that I can do it via model.setParameters). What am I doing wrong?

[dokester]

Hi Alex,

I am not sure what is happening. I have no problem making a 3 line Voigt model plus an Astropy polynomial background.

m = VoigtModel( ) + VoigtModel( ) + VoigtModel( )
pm = modeling.models.Polynomial1D( 1 )
m += AstropyModel( pm )
params = numpy.arange( 14, dtype=float )
m.parameters = params
print( m )
print( m.parameters )

Which produces as output:

Voigt: +
Voigt: +
Voigt: +
AstropyModel( Polynomial1D )
[ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13.]

That seems all OK.

I see one option where this might go wrong. The astropy.modeling.model needs to be a FittableModel. Otherwise it misses some necessary items, that I could not work around.

Would it be possible to introduce your astropy-model as a UserModel(). That will give you more freedom. See eg. in the examples directory: Bessel-J0-UserModel.ipynb or in the test directory TestUserModel.py.

Please note that for fitting a spectrum BayesicFitting provides a RepeatingModel() which repeats the same model several times but also can keep some parameters at the same value for all lines. Another option is CombiModel() which is meant for line multiplets with even more options to fix similar parameters. Both RepeatingModel and CombiModel can be combined with any other Model, eg for providing a background. Even AstropyModel or UserModel should not be a problem.

I hope this helps you.

Do.

[alexserebryanskiy]

Dear Do,

Thank you very much for the prompt reply and very valuable information regarding your software. I have tried CombiModel() for several Voigt profiles and then added another CombiModel() for the set of Gaussian profiles. Seems working just fine because the number of parameters is equal to what is expected (in my case 189). But when I add yet another model (say Polynomial(3) or AstropyModel(SBPL1D)) the number of parameters of the resulting model is reduced by the amount of additional model parameters (become 185 or 186, depending on the model) !

best regards,

Alex

On Tue, Feb 6, 2024 at 5:43 PM Do Kester @.***> wrote:

Hi Alex,

I am not sure what is happening. I have no problem making a 3 line Voigt model plus an Astropy polynomial background.

m = VoigtModel( ) + VoigtModel( ) + VoigtModel( ) pm = modeling.models.Polynomial1D( 1 ) m += AstropyModel( pm ) params = numpy.arange( 14, dtype=float ) m.parameters = params print( m ) print( m.parameters )

Which produces as output:

Voigt: + Voigt: + Voigt: + AstropyModel( Polynomial1D ) [ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13.]

That seems all OK.

I see one option where this might go wrong. The astropy.modeling.model needs to be a FittableModel. Otherwise it misses some necessary items, that I could not work around.

Would it be possible to introduce your astropy-model as a UserModel(). That will give you more freedom. See eg. in the examples directory: Bessel-J0-UserModel.ipynb or in the test directory TestUserModel.py.

Please note that for fitting a spectrum BayesicFitting provides a RepeatingModel() which repeats the same model several times but also can keep some parameters at the same value for all lines. Another option is CombiModel() which is meant for line multiplets with even more options to fix similar parameters. Both RepeatingModel and CombiModel can be combined with any other Model, eg for providing a background. Even AstropyModel or UserModel should not be a problem.

I hope this helps you.

Do.

— Reply to this email directly, view it on GitHub https://github.com/dokester/BayesicFitting/issues/21#issuecomment-1929345001, or unsubscribe https://github.com/notifications/unsubscribe-auth/AZHTSEMTQZXGNPCX65GDOZLYSIJNBAVCNFSM6AAAAABC3QAXOKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSMRZGM2DKMBQGE . You are receiving this because you authored the thread.Message ID: @.***>

[dokester]

Hi Alex,

I tried some things but I cannot reproduce your problem. Here is my trial, which worked as expected.

nv = 50 
vm = VoigtModel()
cm = CombiModel( vm, nv, addCombi={2:[0]*nv, 3:[0]*nv} )
print( cm )
print( cm.npars )

ng = 29
gm = GaussModel()
cm += CombiModel( gm, ng )
print( cm )
print( cm.npars )

cm += PolynomialModel( 1 )
print( cm )
print( cm.npars )

Which resulted in:

Combi of 50 times Voigt
102
Combi of 50 times Voigt +
Combi of 29 times Gauss
189
Combi of 50 times Voigt +
Combi of 29 times Gauss +
Polynomial: f( x:p ) = p_189 + p_190 * x
191

In the CombiModel of Voigt I made the GaussWidth and the LorentzWidth the same for all lines. Hence the 102 parameters for it (50 for amplitude; 50 for position and 2 for all widths). I could do similar things for the Combi of the gausses.

So please send me some lines of code that demonstrate where it fails.

Do.

[alexserebryanskiy]

Dear Do,

Thanks a lot for your help. I think I found a flaw in my code. Here is what is going on: when a combined model is created as you explained, i.e. cm=VoigtModel(); cm+=GaussModel(); cm+=PolynomialModel() then it is working fine as expected, but if I do it this way: m1 = VoigtModel(); m2 = GaussModel(); cm = m1 + m2 + PolynomialModel() then I have a problem (number of parameters is wrong)! It looks like I have to be very careful with model manipulations.

with Best regards,

Alex

On Wed, Feb 7, 2024 at 6:13 PM Do Kester @.***> wrote:

Hi Alex,

I tried some things but I cannot reproduce your problem. Here is my trial, which worked as expected.

nv = 50 vm = VoigtModel() cm = CombiModel( vm, nv, addCombi={2:[0]nv, 3:[0]nv} ) print( cm ) print( cm.npars )

ng = 29 gm = GaussModel() cm += CombiModel( gm, ng ) print( cm ) print( cm.npars )

cm += PolynomialModel( 1 ) print( cm ) print( cm.npars )

Which resulted in:

Combi of 50 times Voigt 102 Combi of 50 times Voigt + Combi of 29 times Gauss 189 Combi of 50 times Voigt + Combi of 29 times Gauss + Polynomial: f( x:p ) = p_189 + p_190 * x 191

In the CombiModel of Voigt I made the GaussWidth and the LorentzWidth the same for all lines. Hence the 102 parameters for it (50 for amplitude; 50 for position and 2 for all widths). I could do similar things for the Combi of the gausses.

So please send me some lines of code that demonstrate where it fails.

Do.

— Reply to this email directly, view it on GitHub https://github.com/dokester/BayesicFitting/issues/21#issuecomment-1931920472, or unsubscribe https://github.com/notifications/unsubscribe-auth/AZHTSEMLVVRXNTUXVSUMUZDYSNVWBAVCNFSM6AAAAABC3QAXOKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSMZRHEZDANBXGI . You are receiving this because you authored the thread.Message ID: @.***>

[alexserebryanskiy]

Sorry, just another question. It seems that I cannot set limits on the parameters of the combined model. When I tried to do that, I got the message, "The (compound) model does not have 5 parameters".

On Wed, Feb 7, 2024 at 6:13 PM Do Kester @.***> wrote:

Hi Alex,

I tried some things but I cannot reproduce your problem. Here is my trial, which worked as expected.

nv = 50 vm = VoigtModel() cm = CombiModel( vm, nv, addCombi={2:[0]nv, 3:[0]nv} ) print( cm ) print( cm.npars )

ng = 29 gm = GaussModel() cm += CombiModel( gm, ng ) print( cm ) print( cm.npars )

cm += PolynomialModel( 1 ) print( cm ) print( cm.npars )

Which resulted in:

Combi of 50 times Voigt 102 Combi of 50 times Voigt + Combi of 29 times Gauss 189 Combi of 50 times Voigt + Combi of 29 times Gauss + Polynomial: f( x:p ) = p_189 + p_190 * x 191

In the CombiModel of Voigt I made the GaussWidth and the LorentzWidth the same for all lines. Hence the 102 parameters for it (50 for amplitude; 50 for position and 2 for all widths). I could do similar things for the Combi of the gausses.

So please send me some lines of code that demonstrate where it fails.

Do.

— Reply to this email directly, view it on GitHub https://github.com/dokester/BayesicFitting/issues/21#issuecomment-1931920472, or unsubscribe https://github.com/notifications/unsubscribe-auth/AZHTSEMLVVRXNTUXVSUMUZDYSNVWBAVCNFSM6AAAAABC3QAXOKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSMZRHEZDANBXGI . You are receiving this because you authored the thread.Message ID: @.***>

[dokester]

Hi Alex,

In both of your examples I cannot find fault. Lets look at the first one.

m1 = VoigtModel() ; m2 = GaussModel() ; cm = m1 + m2 + PolynomialModel( 1 )
print( cm )
print( cm.npars, len( cm.parameters ) )

Which results in:

Voigt: +
Gauss: f( x:p ) = p_4 * exp( -0.5 * ( ( x - p_5 ) / p_6 )^2 ) +
Polynomial: f( x:p ) = p_7 + p_8 * x
9 9

Indeed 9 is what we expect: 4 for Voigt + 3 for Gauss + 2 for Poly.

For the 2nd example i have

nv = 5 
vm = VoigtModel()
cm = CombiModel( vm, nv )

prior = UniformPrior( limits=[-10,20] )
for k in range( cm.npars ) :  
    cm.setPrior( k, prior=prior )

# printclass is a method in Tools.py which prints all attributes of a class
printclass( cm )
print( "Nr of priors is ", len( cm.priors ) )

This results in

+++++++++++++++++++++++++++++++++++++++++++++++++++++++
Combi of 5 times Voigt
+++++++++++++++++++++++++++++++++++++++++++++++++++++++
_head           Combi of 5 times Voigt
_next           None
_npchain        20
_operation      0
addindex        []
addvalue        []
deep            1
deltaP          [1.e-05]
expandindex     [ 0  1  2  3 ... 16 17 18 19] 20
fixed           None
mlist           []
model           Voigt + Voigt + Voigt + Voigt + Voigt
mulindex        []
mulvalue        []
ndim            1
nmp             4
nonZero         []
npbase          20
npcmax          20
npmax           20
nrepeat         5
parNames        []
parameters      [1. 0. 1. 1. ... 1. 0. 1. 1.] 20
parlist         None
posIndex        []
priors          [UniformPrior UniformPrior UniformPrior UniformPrior UniformPrior UniformPrior UniformPrior UniformPrior... ] 20
select          [ 0  1  2  3 ... 16 17 18 19] 20
stdevs          None
tiny            1e-20
xUnit           
yUnit           
Nr of priors is  20

As you can see the model has 20 priors, one for each of its parameters. In this case they as all the same.

Maybe you tried to put priors on the VoigtModel itself, which of course only has 4 parameters. If you put them on the CombiModel it is OK.

Success Do.

PS. As you can see I found the button to make some parts of the text into visual code. 😏

[dokester]

Hi Alex,

Spurred by your spectral fitting problem, I dusted off a piece of code that I had lying around which actually does that. I will comment it where needed. Lets assume that the variables freq and flux, contain arrays of frequencies and fluxes of the same length.

First we define the model.

vgt = VoigtModel( )
model = RepeatingModel( 1, vgt, same=[2,3] )
model.growPrior = ExponentialPrior( scale = 20 )
lolim = [ 0.0, 1105, 0.0001, 0.0001]
hilim = [10.0, 1106, 0.05, 0.05]
model.setLimits( lowLimits=lolim, highLimits=hilim )

# Fix p[1] at -1.0 and let p[0] vary a few percent to fit the continuum.
em = ExpModel( fixed={1:-1.0} )
em.setLimits( lowLimits=[0.99], highLimits=[1.01])
model |= em

RepeatingModel is just what it says: it repeats the same model (here vgt) an unknown number of times. Using NestedSampler we automatically decide how many lines we need. For the definition, we start with 1. The same=[2,3] ensures that all parameters 2 and 3 (lorentzwidth and gausswidth) are the same (if there are more than 1 line) The growprior is a prior that governs the growth of the model: more or less lines. The scale should be about the number you would expect. The next 3 lines define and set the limits (low and high) for the parameters in voigt, resp amplitude, frequency, g-width, and l-width. You will have to adapt them to your case, but always avoid 0 in the widths. The model does not like it. Note that all lines use the same priors. The model defines emission lines. If this is the case for you, skip the next lines. I had absorption lines. We have to pipe (like a UNIX pipe) the result through an ExpModel where we fix the second parameter (decay) to -1. The amplitude may vary in a limited range around 1, to accommodate for deviations of the continuum from 1. Here too you might have to adapt it to your case.

Now we set up NestedSampler

ns = NestedSampler( freq, model, flux, verbose=2, limits=[0.001, 0.1] )
evidence = ns.sample( plot=True )
sl = ns.samples
yfit = ns.yfit

The limits are for the JeffreysPrior on the scale of the Gaussian error distribution. The sample method sets NestedSampler churning. In my case with 2000 data points and eventually 27 lines took a few hours. So try something simple first. At the end, it will plot the results. The variable sl will contain a list of weighted samples drawn from the posterior and yfit is the average fit over all samples. Look at other examples of what you can do with them.

Hope this helps you to get started.

Do.

[alexserebryanskiy]

Dear Do,

Thank you so much for your time to assist me in mastering BayessicFitting! I really appreciate it. I'll learn what you suggested today and will try to use it. In my research I'm trying to analyse Seyferts spectra to decide which model (not physical but spectral one) is the best using Evidence parameter. I used to use emcee for Bayesian inference but unfortunately they do not provide a way to compute the evidence. It will take much longer for me to develop such a tool and looks like you BaesicFitting is that tool. Unfortunately, I did not find the appropriate example of how to use it for my particular task and that is why I decided to bother you (sorry for that) to try it as soon as possible. Seyfert spectra are quite complicated for different reasons (dynamics, active nuclear etc., wide and narrow line blending etc.) and your explanation is very valuable to get the idea of what can be done with BayesicFitting. Thanks again!

with Best Regards,

Alex

On Fri, Feb 9, 2024 at 2:22 AM Do Kester @.***> wrote:

Hi Alex,

Spurred by your spectral fitting problem, I dusted off a piece of code that I had lying around which actually does that. I will comment it where needed. Lets assume that the variables freq and flux, contain arrays of frequencies and fluxes of the same length.

First we define the model.

vgt = VoigtModel( ) model = RepeatingModel( 1, vgt, same=[2,3] ) model.growPrior = ExponentialPrior( scale = 20 ) lolim = [ 0.0, 1105, 0.0001, 0.0001] hilim = [10.0, 1106, 0.05, 0.05] model.setLimits( lowLimits=lolim, highLimits=hilim )

Fix p[1] at -1.0 and let p[0] vary a few percent to fit the continuum.

em = ExpModel( fixed={1:-1.0} ) em.setLimits( lowLimits=[0.99], highLimits=[1.01]) model |= em

RepeatingModel is just what it says: it repeats the same model (here vgt) an unknown number of times. Using NestedSampler we automatically decide how many lines we need. For the definition, we start with 1. The same=[2,3] ensures that all parameters 2 and 3 (lorentzwidth and gausswidth) are the same (if there are more than 1 line) The growprior is a prior that governs the growth of the model: more or less lines. The scale should be about the number you would expect. The next 3 lines define and set the limits (low and high) for the parameters in voigt, resp amplitude, frequency, g-width, and l-width. You will have to adapt them to your case, but always avoid 0 in the widths. The model does not like it. Note that all lines use the same priors. The model defines emission lines. If this is the case for you, skip the next lines. I had absorption lines. We have to pipe (like a UNIX pipe) the result through an ExpModel where we fix the second parameter (decay) to -1. The amplitude may vary in a limited range around 1, to accommodate for deviations of the continuum from 1. Here too you might have to adapt it to your case.

Now we set up NestedSampler

ns = NestedSampler( freq, model, flux, verbose=2, limits=[0.001, 0.1] ) evidence = ns.sample( plot=True ) sl = ns.samples yfit = ns.yfit

The limits are for the JeffreysPrior on the scale of the Gaussian error distribution. The sample method sets NestedSampler churning. In my case with 2000 data points and eventually 27 lines took a few hours. So try something simple first. At the end, it will plot the results. The variable sl will contain a list of weighted samples drawn from the posterior and yfit is the average fit over all samples. Look at other examples of what you can do with them.

Hope this helps you to get started.

Do.

— Reply to this email directly, view it on GitHub https://github.com/dokester/BayesicFitting/issues/21#issuecomment-1934875068, or unsubscribe https://github.com/notifications/unsubscribe-auth/AZHTSEPSAN7M2YGMDEUTIU3YSUXZTAVCNFSM6AAAAABC3QAXOKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSMZUHA3TKMBWHA . You are receiving this because you authored the thread.Message ID: @.***>

[alexserebryanskiy]

Dear Do,

Here is the code with slight modification of yours to replicate the "problem" (mostly my misunderstanding the package's features):

m1 = VoigtModel() ; m2 = GaussModel() ; cm = m1 + m2 + PolynomialModel( 1 ) print( cm ) print( cm.npars, len( cm.parameters ) )

nv = 5 vm = VoigtModel() cm = CombiModel( vm, nv ) cm += PolynomialModel(3) prior = UniformPrior( limits=[-10,20] ) for k in range( cm.npars ) : cm.setPrior( k, prior=prior )

printclass is a method in Tools.py which prints all attributes of a class

printclass( cm ) print( "Nr of priors is ", len( cm.priors ) )

which produce the following results:

Voigt: + Gauss: f( x:p ) = p_4 exp( -0.5 ( ( x - p_5 ) / p_6 )^2 ) + Polynomial: f( x:p ) = p_7 + p_8 x 9 9 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ Combi of 5 times Voigt + Polynomial: f( x:p ) = p_20 + p_21 x + p_22 x^2 + p_23 x^3 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ _head Combi of 5 times Voigt + Polynomial: f( x:p ) = p_20 + p_21 x + p_22 x^2 + p_23 x^3 _next Polynomial: f( x:p ) = p_0 + p_1 x + p_2 x^2 + p_3 x^3 _npchain 24 _operation 0 addindex [] addvalue [] deep 1 deltaP [1.e-05] expandindex [ 0 1 2 3 ... 16 17 18 19] 20 fixed None mlist [] model Voigt + Voigt + Voigt + Voigt + Voigt mulindex [] mulvalue [] ndim 1 nmp 4 nonZero [] npbase 20 npcmax 20 npmax 20 nrepeat 5 parNames [] parameters [1. 0. 1. 1. ... 0. 0. 0. 0.] 24 parlist None posIndex [] priors [UniformPrior UniformPrior UniformPrior UniformPrior UniformPrior UniformPrior UniformPrior UniformPrior... ] 20 select [ 0 1 2 3 ... 16 17 18 19] 20 stdevs None tiny 1e-20 xUnit yUnit Nr of priors is 20

As you can see the number of parameters is still 20, not 23 as I expected. And now if I try :

lolim = -1 * np.ones(cm.npars,) hilim = np.ones(cm.npars,)

cm.setLimits(lolim,hilim)

I will get that error message: The (compound) model does not have 5 parameters

Regards,

Alex

On Thu, Feb 8, 2024 at 5:28 PM Do Kester @.***> wrote:

Hi Alex,

In both of your examples I cannot find fault. Lets look at the first one.

m1 = VoigtModel() ; m2 = GaussModel() ; cm = m1 + m2 + PolynomialModel( 1 ) print( cm ) print( cm.npars, len( cm.parameters ) )

Which results in:

Voigt: + Gauss: f( x:p ) = p_4 exp( -0.5 ( ( x - p_5 ) / p_6 )^2 ) + Polynomial: f( x:p ) = p_7 + p_8 * x 9 9

Indeed 9 is what we expect: 4 for Voigt + 3 for Gauss + 2 for Poly.

For the 2nd example i have

nv = 5 vm = VoigtModel() cm = CombiModel( vm, nv )

prior = UniformPrior( limits=[-10,20] ) for k in range( cm.npars ) : cm.setPrior( k, prior=prior )

printclass is a method in Tools.py which prints all attributes of a class

printclass( cm ) print( "Nr of priors is ", len( cm.priors ) )

This results in

+++++++++++++++++++++++++++++++++++++++++++++++++++++++ Combi of 5 times Voigt +++++++++++++++++++++++++++++++++++++++++++++++++++++++ _head Combi of 5 times Voigt _next None _npchain 20 _operation 0 addindex [] addvalue [] deep 1 deltaP [1.e-05] expandindex [ 0 1 2 3 ... 16 17 18 19] 20 fixed None mlist [] model Voigt + Voigt + Voigt + Voigt + Voigt mulindex [] mulvalue [] ndim 1 nmp 4 nonZero [] npbase 20 npcmax 20 npmax 20 nrepeat 5 parNames [] parameters [1. 0. 1. 1. ... 1. 0. 1. 1.] 20 parlist None posIndex [] priors [UniformPrior UniformPrior UniformPrior UniformPrior UniformPrior UniformPrior UniformPrior UniformPrior... ] 20 select [ 0 1 2 3 ... 16 17 18 19] 20 stdevs None tiny 1e-20 xUnit yUnit Nr of priors is 20

As you can see the model has 20 priors, one for each of its parameters. In this case they as all the same.

Maybe you tried to put priors on the VoigtModel itself, which of course only has 4 parameters. If you put them on the CombiModel it is OK.

Success Do.

PS. As you can see I found the button to make some parts of the text into visual code. 😏

— Reply to this email directly, view it on GitHub https://github.com/dokester/BayesicFitting/issues/21#issuecomment-1933880285, or unsubscribe https://github.com/notifications/unsubscribe-auth/AZHTSEOO4CHE3JD3LS76PV3YSSZGRAVCNFSM6AAAAABC3QAXOKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSMZTHA4DAMRYGU . You are receiving this because you authored the thread.Message ID: @.***>

[dokester]

Hi Alex,

You were right. You found a corner of the package that I did not test well enough. And of course there were some errors.

It goes too far to explain where it exactly failed, but it has to do with how compound models are partitioned into simpler ones and the distribution of priors over them.

I have repaired it in my private sources and it will appear in the next release. That might take a while. There is a lot of testing to do and other hoops to jump throughs, before you can publish a new release. And I have to do it all on my own. So bear with me.

However not all is bad as you can program around it by setting the limits to the constituent models and then combine them. But only for RepeatingModel.

nv = 2
ng = 3
vm = VoigtModel()
hilim = 1 + numpy.arange( vm.npars )
lolim = - hilim
vm.setLimits( lolim, hilim )

cm = RepeatingModel( nv, vm, dynamic=False )

gm = GaussModel()
hilim = 1 + numpy.arange( gm.npars ) / 5
lolim = - hilim
gm.setLimits( lolim, hilim )

cm += RepeatingModel( ng, gm, dynamic=False )

pm = PolynomialModel( 2 )
pm.setLimits( lolim-5, hilim+5 )
cm += pm

for k in range( cm.npars ) :
    print( "prior   %2d" % k, cm.getPrior( k ) )

You need to turn the dynamic behaviour off in at least one Repeating model as a compound model can only have one dynamic constituent. However CombiModel is not dynamic at all, so you maybe dont need it.

Please note that the order of parameters in RepeatingModel( nrepeat, model ) is the reverse of that in CombiModel( model, ncombi ). Murphy's law is unavoidable, unless you really take care.

Thanks for spotting the failures; it really helps when people complain and I appreciate it.

Do.

[alexserebryanskiy]

Dear Do,

Thank you for your guidance in BayesicFitting. I understand what it means to develop and then try to maintain something. It is always difficult. I appreciate you spending your time replying to me. And I'm not complaining for sure ;) I use this opportunity to learn something new and valuable.

In one of my previous versions of code using BayessicFitting I did exactly as you suggested to place the limits. And, I think, CombiModel also works when using the same way to set the limits. I'm running the code with CombiModel right now. However, I think I do not fully understand the keyword "limits" in NestedSampler. I probably need to read your paper one more time.

with best regards,

Alex

On Sat, Feb 10, 2024 at 2:35 AM Do Kester @.***> wrote:

Hi Alex,

You were right. You found a corner of the package that I did not test well enough. And of course there were some errors.

It goes too far to explain where it exactly failed, but it has to do with how compound models are partitioned into simpler ones and the distribution of priors over them.

I have repaired it in my private sources and it will appear in the next release. That might take a while. There is a lot of testing to do and other hoops to jump throughs, before you can publish a new release. And I have to do it all on my own. So bear with me.

However not all is bad as you can program around it by setting the limits to the constituent models and then combine them. But only for RepeatingModel.

nv = 2 ng = 3 vm = VoigtModel() hilim = 1 + numpy.arange( vm.npars ) lolim = - hilim vm.setLimits( lolim, hilim )

cm = RepeatingModel( nv, vm, dynamic=False )

gm = GaussModel() hilim = 1 + numpy.arange( gm.npars ) / 5 lolim = - hilim gm.setLimits( lolim, hilim )

cm += RepeatingModel( ng, gm, dynamic=False )

pm = PolynomialModel( 2 ) pm.setLimits( lolim-5, hilim+5 ) cm += pm

for k in range( cm.npars ) : print( "prior %2d" % k, cm.getPrior( k ) )

You need to turn the dynamic behaviour off in at least one Repeating model as a compound model can only have one dynamic constituent. However CombiModel is not dynamic at all, so you maybe dont need it.

Please note that the order of parameters in RepeatingModel( nrepeat, model ) is the reverse of that in CombiModel( model, ncombi ). Murphy's law is unavoidable, unless you really take care.

Thanks for spotting the failures; it really helps when people complain and I appreciate it.

Do.

— Reply to this email directly, view it on GitHub https://github.com/dokester/BayesicFitting/issues/21#issuecomment-1936568904, or unsubscribe https://github.com/notifications/unsubscribe-auth/AZHTSEIXALZEXK4WNURANRLYS2CCVAVCNFSM6AAAAABC3QAXOKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSMZWGU3DQOJQGQ . You are receiving this because you authored the thread.Message ID: @.***>

[dokester]

Hi Alex,

The method model.setLimits( lolim, hilim ) sets a UniformPrior with limits to the parameters. It is nothing more than

for k, lo, hi in enumerate( zip( lolim, hilim ) ) :
     model.setPrior( k, UniformPrior( limits=[lo, hi] ) )

However, when I read my manual again, I saw that Priors are attributes of simple models. They should be assigned before the models are used in a compound model. I.e. chained together or used in a CombiModel or RepeatingModel. In the discussions we had here, I had more or less forgotten what I wrote in the manual, and why. Please read the manual, here on github: BayesicFitting/docs/manual.md, especially the section on Priors.

For CombiModel and RepeatingModel it means that the number of priors is determined by the encapsulated model (4 for Voigt, 3 for Gauss, etc.), independent of the number of repeats of the model. I don't see any use for different priors on the different components of a spectrum. Moreover, RepeatingModel is a dynamic model, that NestedSampler can grow or shrink until it finds the best fit. Repeating the priors is then a necessity.

So maybe it should be forbidden to assign priors to a compound model. On the other hand it is not my policy to forbid things just because I dont see the the benefit of it, as long as I can guarantee that it is in principle OK. And there it fails. Because RepeatingModel and CombiModel need less priors than there are parameters and on the outside of a model it is not directly clear which models it contains. So I decided to issue a serious warning when a user is setting priors cq. limits to a compound model. In the next release.

What do you think about this.

Do.

[alexserebryanskiy]

Dear Do,

Thank you for the clarification. I'll be back to this issue tomorrow. Meanwhile I built toy-model and tried to fitt it . The code and results are attached. Something is wrong, but I cannot figure out what I'm doing wrong.

best regards,

Alex

------ code -----

from BayesicFitting import (VoigtModel, CombiModel, GaussModel, PolynomialModel, LevenbergMarquardtFitter, AstropyModel, NestedSampler, RepeatingModel, ExponentialPrior, ExpModel) from BayesicFitting import formatter as fmt import numpy as np import matplotlib.pyplot as plt from astropy.modeling.models import Voigt1D, Gaussian1D from astropy.modeling.models import Linear1D as L1D

x = np.arange(1000, 5000, 0.1)

NV = 2 # number of Voigt profiles NG = 2 # number of Gauss profiles

vm = VoigtModel() lolim = np.array([0.8, 2000, 4, 4]) hilim = np.array([2.5, 4000., 60., 70.]) vm.setLimits( lowLimits=lolim, highLimits=hilim ) model = CombiModel(vm, NV)

model = RepeatingModel(NV, vm, dynamic = False)

gm = GaussModel() lolim = np.array([0.8, 2000, 4]) hilim = np.array([2.5, 4000., 70.]) gm.setLimits( lowLimits=lolim, highLimits=hilim ) model += CombiModel(gm, NG)

model += RepeatingModel(NG, gm, dynamic = False)

pm = AstropyModel(L1D(slope=0.0001, intercept=0.5)) #PolynomialModel(1) pm.setLimits(lowLimits=[-1, 0.1], highLimits=[1, 0.8])

model.priors model += pm

for k in range(model.npars): print(f" prior {k}: {model.getPrior(k)}")

Model_Astropy = Voigt1D(x_0=2050., amplitude_L=1., fwhm_L=35, fwhm_G=23, name='v1') Model_Astropy += Voigt1D(x_0=2550, amplitude_L=2, fwhm_L=40, fwhm_G=30, name='v2') Model_Astropy += Gaussian1D(amplitude=2, mean=3050, stddev=20, name='g1') Model_Astropy += Gaussian1D(amplitude=1.5, mean=3350, stddev=22, name='g2') Model_Astropy += L1D(slope=0.0001, intercept=0.5)

rng = np.random.default_rng(seed=42) Noise = rng.normal(0., 0.02, x.shape)

Obs = Model_Astropy(x) + Noise

fig,ax = plt.subplots(ncols=1, nrows=1, tight_layout=True, figsize=(10,6)) ax.plot(x, Model_Astropy(x)) ax.plot(x, Obs, '.')

model.parameters = np.array([2040., 0.8, 30., 30., 2600., 1., 30., 30., 1., 3100, 30., 2., 3400., 30., 0.0, 0.7])

fitter = LevenbergMarquardtFitter(x, model) param = fitter.fit(Obs, plot=True)

print( "Parameters :", fmt( param, max=None ) ) print( "StDevs :", fmt( fitter.stdevs, max=None ) ) print( "Chisq :", fmt( fitter.chisq ) ) print( "Scale :", fmt( fitter.scale ) ) print( "Evidence :", fmt( fitter.getEvidence( limits=[-100,100] ) ) )

ns = NestedSampler( x, model, Obs, verbose=2, limits=[0.01, 1.0] ) evidence = ns.sample( plot=True )

----- results ----- (figure is attached)

Parameters [ 2.016 3049.952 17.633 4.002 0.815 2050.079 10.133 16.872 1.497 3349.995 21.966 1.473 2550.006 27.441 0.000 0.511 0.028] Engines success reject failed best calls GalileanEngine 172263 55971 66412 0 21600 ChordEngine 147994 2188605 3311 0 21600 Calls to LogL 2690527 to dLogL 55971 Samples 21600 Evidence 37126.295 +- 0.609

On Fri, Feb 16, 2024 at 1:39 AM Do Kester @.***> wrote:

Hi Alex,

The method model.setLimits( lolim, hilim ) sets a UniformPrior with limits to the parameters. It is nothing more than

for k, lo, hi in enumerate( zip( lolim, hilim ) ) : model.setPrior( k, UniformPrior( limits=[lo, hi] ) )

However, when I read my manual again, I saw that Priors are attributes of simple models. They should be assigned before the models are used in a compound model. I.e. chained together or used in a CombiModel or RepeatingModel. In the discussions we had here, I had more or less forgotten what I wrote in the manual, and why. Please read the manual, here on github: BayesicFitting/docs/manual.md, especially the section on Priors.

For CombiModel and RepeatingModel it means that the number of priors is determined by the encapsulated model (4 for Voigt, 3 for Gauss, etc.), independent of the number of repeats of the model. I don't see any use for different priors on the different components of a spectrum. Moreover, RepeatingModel is a dynamic model, that NestedSampler can grow or shrink until it finds the best fit. Repeating the priors is then a necessity.

So maybe it should be forbidden to assign priors to a compound model. On the other hand it is not my policy to forbid things just because I dont see the the benefit of it, as long as I can guarantee that it is in principle OK. And there it fails. Because RepeatingModel and CombiModel need less priors than there are parameters and on the outside of a model it is not directly clear which models it contains. So I decided to issue a serious warning when a user is setting priors cq. limits to a compound model. In the next release.

What do you think about this.

Do.

— Reply to this email directly, view it on GitHub https://github.com/dokester/BayesicFitting/issues/21#issuecomment-1947113958, or unsubscribe https://github.com/notifications/unsubscribe-auth/AZHTSELLSPF4DA376XZTYRDYTZP7BAVCNFSM6AAAAABC3QAXOKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSNBXGEYTGOJVHA . You are receiving this because you authored the thread.Message ID: @.***>

[dokester]

Hi Alex,

Your code did not work for me either. So I made some adaptations, keeping the essentials. It took me a while to do so. Here is what I found.

First of all, I want to note that the Astropy (AP) VoigtModel has not the same definition as the one in BayesicFiiting (BF). The BF model is made such that when the Lorentz-width is going to 0, the VoigtModel approaches a GaussModel with same amplitude, center and width and when the Gauss-width goes to 0 it approaches a LorentzModel. In AP the widths are FWHM instead of sigma. Also the amplitude has another value. In BF the max of the model is 1.0 so the amplitude is indeed the top. It is a matter of choices and not very relevant to fitting.

When running your code, it got stuck in some local minimum, not the one we are looking for. We are searching in a 17 dimensional space for a tiny sweet spot where the global maximum of the likelihood is located. If we are only a small bit off, especially with the line centers, the model does not fit at all. The likelihood landscape is riddled with small ripples with no easy way to travel to the global maximum. Most of the time, such behaviour can be ameliorated by restricting the available space for some of the parameters and/or widening it for others. All to locate the sweet spot before it gets frozen in a local minimum.

The main culprit I think is the slope of the Polynomial for the background. It can vary between [-1,1]. Both limits are obviously very wrong: The background would go from +-1000 to +-5000, while all of the data lie between 0 and 3. If the outcome is so wrong that you would not believe it, you should adapt the Prior. I set a UniformPrior for the slope with limits=[-0.01,0.01]. And for the intercept i set a UniformPrior with limits=[-1,1].

Then we have the prior on the amplitude of the lines. The low limit is set at 0.8. If you try to place a line in a random spot on the frequency axis, you are already lucky to hit the wing, let alone the center. From the wing, the algorithm can find a path to the center. For that to happen the wing position has to be accepted first. An acceptable amplitude in the wing might be much less that 0.8. So with this low limit, valuable starting positions might be rejected. I reset the low limit for the amplitude to 0.0, giving it more space.

As I wrote before, the GaussModel is an approximation in the VoigtModel. So I used 4 VoigtModel instead of 2 Voigt and 2 Gauss. With this setup, when we hit a line, the algorithm can adapt the widths to form either a Gauss or a fat-winged Voigt. Every hit is OK. But when you hit a Voigtline with a GaussModel, it can never leave, because it is mostly OK. It freezes in a possibly almost global maximum. Not the one we are looking for. Using 4 VoigtModels increases to number of parameters to 19. It is a larger space but easier to search.

With these adaptations I rerun your example. I decreased to number of point by factor 100, to speed things up. And I still got stuck. Two lines were found, the other two were in limbo. At that moment I was still using the ComboModel. I thought that a dynamic RepeatingModel should have a better chance of locating the lines. The "dynamic" means that it can have a varying number of models. NestedSampling automatically finds the optimal number. Lines in limbo are simply rejected as the same model without the limbo is more probable. (Occam's Razor). And that simpler model can try to grow with a new line at a hopefully more fortunate location.

With this model I am now running the example again. It has not yet finished but it looks like it goes well. I expect it to finish today. Then I will attach the code and the plot here.

As I always say: Fitting an equal amount of science, craft and art. There are no simple recipes that always work.

Best Regards Do.

[alexserebryanskiy]

Dear Do,

Thank you so much for your explanation and the time you devoted to my problem! That's a lot of information. I'm reading your email and will be back to you with my reply as soon as possible. Meanwhile, I also tried to use ZEUS, NESTLE, and DYNESTY with some success.

Thank you!

Best regards,

Alex

On Fri, Feb 23, 2024 at 6:12 PM Do Kester @.***> wrote:

Hi Alex,

Your code did not work for me either. So I made some adaptations, keeping the essentials. It took me a while to do so. Here is what I found.

First of all, I want to note that the Astropy (AP) VoigtModel has not the same definition as the one in BayesicFiiting (BF). The BF model is made such that when the Lorentz-width is going to 0, the VoigtModel approaches a GaussModel with same amplitude, center and width and when the Gauss-width goes to 0 it approaches a LorentzModel. In AP the widths are FWHM instead of sigma. Also the amplitude has another value. In BF the max of the model is 1.0 so the amplitude is indeed the top. It is a matter of choices and not very relevant to fitting.

When running your code, it got stuck in some local minimum, not the one we are looking for. We are searching in a 17 dimensional space for a tiny sweet spot where the global maximum of the likelihood is located. If we are only a small bit off, especially with the line centers, the model does not fit at all. The likelihood landscape is riddled with small ripples with no easy way to travel to the global maximum. Most of the time, such behaviour can be ameliorated by restricting the available space for some of the parameters and/or widening it for others. All to locate the sweet spot before it gets frozen in a local minimum.

The main culprit I think is the slope of the Polynomial for the background. It can vary between [-1,1]. Both limits are obviously very wrong: The background would go from +-1000 to +-5000, while all of the data lie between 0 and 3. If the outcome is so wrong that you would not believe it, you should adapt the Prior. I set a UniformPrior for the slope with limits=[-0.01,0.01]. And for the intercept i set a UniformPrior with limits=[-1,1].

Then we have the prior on the amplitude of the lines. The low limit is set at 0.8. If you try to place a line in a random spot on the frequency axis, you are already lucky to hit the wing, let alone the center. From the wing, the algorithm can find a path to the center. For that to happen the wing position has to be accepted first. An acceptable amplitude in the wing might be much less that 0.8. So with this low limit, valuable starting positions might be rejected. I reset the low limit for the amplitude to 0.0, giving it more space.

As I wrote before, the GaussModel is an approximation in the VoigtModel. So I used 4 VoigtModel instead of 2 Voigt and 2 Gauss. With this setup, when we hit a line, the algorithm can adapt the widths to form either a Gauss or a fat-winged Voigt. Every hit is OK. But when you hit a Voigtline with a GaussModel, it can never leave, because it is mostly OK. It freezes in a possibly almost global maximum. Not the one we are looking for. Using 4 VoigtModels increases to number of parameters to 19. It is a larger space but easier to search.

With these adaptations I rerun your example. I decreased to number of point by factor 100, to speed things up. And I still got stuck. Two lines were found, the other two were in limbo. At that moment I was still using the ComboModel. I thought that a dynamic RepeatingModel should have a better chance of locating the lines. The "dynamic" means that it can have a varying number of models. NestedSampling automatically finds the optimal number. Lines in limbo are simply rejected as the same model without the limbo is more probable. (Occam's Razor). And that simpler model can try to grow with a new line at a hopefully more fortunate location.

With this model I am now running the example again. It has not yet finished but it looks like it goes well. I expect it to finish today. Then I will attach the code and the plot here.

As I always say: Fitting an equal amount of science, craft and art. There are no simple recipes that always work.

Best Regards Do.

— Reply to this email directly, view it on GitHub https://github.com/dokester/BayesicFitting/issues/21#issuecomment-1961220077, or unsubscribe https://github.com/notifications/unsubscribe-auth/AZHTSEMYQTQ5DS22B5ALJFDYVCBTVAVCNFSM6AAAAABC3QAXOKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSNRRGIZDAMBXG4 . You are receiving this because you authored the thread.Message ID: @.***>

[dokester]

I think we can close this discussion