How to get the 'pblum' distribution data after MCMC？

[d11122]

Hello, I have a confusing question：
I want get the 'pblum' distribution like this:
lp=b.get_parameter('pblum@primary', ).get_distribution('emcee_sol'), but it doesn't work.
I can get other parameter（eg: 'incl') by incl=b.get_parameter('incl', component='binary', context='component').get_distribution('emcee_sol')
By the way, I want get pblum@primary and @secondary, but only 'pblum@primary' is involved in the MCMC calculation.

[kecnry]

Accessing the distribution on the pblum parameter that you did sample should work. If it doesn't, then that is likely a bug and may be related to #562.

Accessing the distribution of the pblum of the other star (assuming you're using something like pblum_mode='component-coupled') is unfortunately not currently supported. The relation between the pblums is not algebraic, but rather requires the mesh to be created and populated, so it's not possible in the current framework to propagate the distribution on the fly as the pblum values aren't stored anywhere. But it should be possible to do this manually by sampling from pblum@primary and calling b.compute_pblums. There is a proposed future extension which could support doing this within phoebe discussed here: #407 (comment).

I'll try to put together an example script which does the manual approach (and in doing so I might be able to reproduce what you're seeing with plbum@primary, in which case we can file an issue and hopefully get it fixed)... but I can't get to that right away. If you (or anyone) wants to try and gets a working example, please post it here and I'd be happy to include it in the example scripts on the website.

[d11122]

Oh, thanks! but I have got the distribution by lp1=b.get_parameter('pblum@primary@lc01@lc@dataset', ).get_distribution('emcee_sol'), so can I get the distribution of 'pblum@secondary stars '? If so, can you please provide some help (e.g. some pseudo-code)?

[kecnry]

Ok, great! Sure, I'll provide some steps that I think should work and try to write up an actual example script later... of if you get it working and are willing to share, feel free to post it back here and I'll adopt it for the website.

The easiest way would be to make use of the distl.function distribution. For that you need to write a function that takes inputs and returns the pblum@primary, so something like (note: this is pseudo-code and won't run as-is, but I tried to stick to the right syntax):

def compute_pblum_secondary(b, pblum_primary, other_distribution):
    b.set_value(qualifier='pblum', component='primary', value=pblum_primary)
    # set any other parameters that you sampled over that affect the secondary pblum (teff, etc)
    return b.compute_pblums()['pblum@secondary@lc01']

You can then define your propagated distribution as something like:

lp2 = distl.function(compute_pblum_secondary, args=(b, lp1, other_distribution), vectorized=False)

Note that whenever you sample from lp2 it will create the mesh under-the-hood (which is relatively expensive)... so if you want to plot, be careful to pass a sample size that is reasonable and won't take forever (1 million samples will take a while to propagate).

If your original MCMC run didn't contain a large number of samples, you could loop over those samples directly instead, compute the secondary luminosity for each of them, and build up a new distribution object from those or plot them manually. This will be a little more work to setup, but removes the random drawing and "downsampling" that you'll get from having distl do the looping and sampling for you.

[d11122]

Ok, thank you！I will try it.

[kecnry]

There is now a full example showing how to make use of this in the online documentation. Note that the pseudocode example above will not work as the bundle cannot be passed as an argument to the callable function sent to distl.