Tentative fix for continuity SFH artefact

[joriswitstok]

Hi Adam,

Hope you’re well! I’m back with a bagpipes question – I have been getting stuck into non-parametric SFHs and was getting some nice fits, although I noticed the two most recent bins always seemed to be fixed at the same value, which seemed suspicious to me.

I had a look under the hood and I think I may have found a bug causing this behaviour (see below, it may well be I misunderstood your implementation). Let me know if you agree!

Cheers, Joris

---

As an example, below I was trying to manually reproduce a 6-bin SFH myself using some best-fit values of the ΔSFR parameters (printing several values for each step to debug). The inferred SFR is exactly the same in the last two bins, while I'm fairly sure this shouldn't be the case?

In [17]: log_sampling = 0.0025
    ...: log_age_max = np.log10(cosmo.age(0).value)+9. + 2*log_sampling
    ...: ages = 10**np.arange(6., log_age_max, log_sampling)
    ...: sfr = np.zeros_like(ages)
    ...: 
    ...: # As in continuity function
    ...: bin_edges = np.array([0.0, 3.0, 10.0, 25.60021148954903, 65.53708283096381, 167.77631808807664, 429.51092259926116])[::-1]*10**6
    ...: n_bins = len(bin_edges) - 1
    ...: dsfrs = [-0.01, 0.23, 0.43, 0.58, 3.6]
    ...: 
    ...: for i in range(1, n_bins+1):
    ...:     mask = (ages < bin_edges[i-1]) & (ages > bin_edges[i])
    ...:     print("ages", np.min(ages[mask]*1e-6), "to", np.max(ages[mask]*1e-6), "Myr")
    ...:     print(i, dsfrs[:i], 10**np.sum(dsfrs[:i]))
    ...:     sfr[mask] += 10**np.sum(dsfrs[:i])
ages 167.88040181239046 to 429.04218923929284 Myr
1 [-0.01] 0.9772372209558107
ages 65.6901116476696 to 166.91678072123963 Myr
2 [-0.01, 0.23] 1.6595869074375607
ages 25.703957827701682 to 65.3130552647899 Myr
3 [-0.01, 0.23, 0.43] 4.466835921509632
ages 10.000000000003599 to 25.55641901065433 Myr
4 [-0.01, 0.23, 0.43, 0.58] 16.982436524617444
ages 3.0026174208617835 to 9.942600739533136 Myr
5 [-0.01, 0.23, 0.43, 0.58, 3.6] 67608.29753919819
ages 1.0 to 2.98538261891847 Myr
6 [-0.01, 0.23, 0.43, 0.58, 3.6] 67608.29753919819

If I simply prepend dsfrs with a zero, it seems to solve the issue and (I think) each of the 5 ΔSFR parameters nicely describe the logarithmic ratio of SFRs in adjacent bins – I assume the code later takes care of normalisation using the total stellar mass:

In [18]: log_sampling = 0.0025
    ...: log_age_max = np.log10(cosmo.age(0).value)+9. + 2*log_sampling
    ...: ages = 10**np.arange(6., log_age_max, log_sampling)
    ...: sfr = np.zeros_like(ages)
    ...: 
    ...: # As in continuity function, with adaptation
    ...: bin_edges = np.array([0.0, 3.0, 10.0, 25.60021148954903, 65.53708283096381, 167.77631808807664, 429.51092259926116])[::-1]*10**6
    ...: n_bins = len(bin_edges) - 1
    ...: dsfrs = [0, -0.01, 0.23, 0.43, 0.58, 3.6]
    ...: 
    ...: for i in range(1, n_bins+1):
    ...:     mask = (ages < bin_edges[i-1]) & (ages > bin_edges[i])
    ...:     print("ages", np.min(ages[mask]*1e-6), "to", np.max(ages[mask]*1e-6), "Myr")
    ...:     print(i, dsfrs[:i], 10**np.sum(dsfrs[:i]))
    ...:     sfr[mask] += 10**np.sum(dsfrs[:i])
ages 167.88040181239046 to 429.04218923929284 Myr
1 [0] 1
ages 65.6901116476696 to 166.91678072123963 Myr
2 [0, -0.01] 0.9772372209558107
ages 25.703957827701682 to 65.3130552647899 Myr
3 [0, -0.01, 0.23] 1.6595869074375607
ages 10.000000000003599 to 25.55641901065433 Myr
4 [0, -0.01, 0.23, 0.43] 4.466835921509632
ages 3.0026174208617835 to 9.942600739533136 Myr
5 [0, -0.01, 0.23, 0.43, 0.58] 16.982436524617444
ages 1.0 to 2.98538261891847 Myr
6 [0, -0.01, 0.23, 0.43, 0.58, 3.6] 67608.29753919819

[ACCarnall]

Hi Joris,

Thanks for this, definitely possible that I've made some sort of mistake. Please could you upload a minimum working example code file demonstrating the issue? Just so I can understand exactly what the issue is and run some tests.

Cheers, Adam

[joriswitstok]

Sure! I've uploaded one into my forked repository (https://github.com/joriswitstok/bagpipes/blob/master/mwe.py).

Result before the fix:

After:

[ACCarnall]

Hi Joris,

Apologies for the long radio-silence on this, but now term is over I should have some time to investigate. I'll try to take a look later on this week. Just checking if anything has changed in the mean time?

Cheers, Adam

[joriswitstok]

Hi Adam,

No worries! There's no updates since last time, hopefully the MWE should clarify things.

Cheers, Joris