x ~ scale * N(theta, 1) or x ~ N(theta, scale)? 🤔

[niels-leif-bracher]

Hello 👋 ,
I was trying to understand the benchmark datasets used to evaluate BayesFlow's performance. One thing that I noted for a couple of datasets is the scaling factor in front of the distribution for sampling the observations. One example is linked below, which is the T.2 task from the sbib paper.
Generating points for the original task would look like this:
theta ~ scale * N(0,1) = N(0, scale * 1)
x ~ N(theta, scale * 1)
with scale being the same for both distributions i.e., scale=0.1

Now your implementation is:
theta ~ scale * N(0,1) (same as above)
x ~ scale * N(theta, 1) = N(scale * theta, scale * 1)
which is equivalent to having a different prior compared to the T.2 task for the same data distribution, i.e.:
theta ~ scale^2 * N(0, 1) = N(0, scale^2)
x ~ N(theta, scale * 1)

maybe im also missing something but i would be glad for clarification. I also wanted to ask what code you used for the multivariate dataset from your paper part 3.3?

https://github.com/stefanradev93/BayesFlow/blob/5ff1c1ea8ee2bda8cf333ae2890b95d05185af23/bayesflow/benchmarks/gaussian_linear_uniform.py#L85

[stefanradev93]

Hi Niels,

the task T2 in the BayesFlow repo is the same as the one in Lueckmann et al., only written slightly differently in code. :)
numpy parameterizes the normal with its scale, not its variance.
Also, I think you are referring to task T1, since T2 has a uniform prior over the locations.

Can you also please specify which paper you mean?

[stefanradev93]

Follow up: You are right in assuming that the default scaling factor in the BayesFlow differs from Lueckmann et al., since the 0.1 in the original paper is a variance and not a scale! Nice spot! :)

You can still adjust it with a keyword argument, but I will modify the code for out of the box compatibility.

[niels-leif-bracher]

Hi Stefan,
Thanks for your fast reply! Sorry for the confusion caused by mixing up the Tasks.

Indeed the task I described is T1 and not T2. Regarding the values from the paper, I think the implementation from SBIB and what they wrote in the paper is also a little different since, as you pointed out, they have 0.1 in the paper as a variance but implemented a scale (standard deviation, instead of variance) of 0.1 in their code.

Anyway, my point is that having a factor in front of a distribution is only equivalent to having that scale factor as std for the normal distribution when the mean is 0. I.e., in the case where you sample the observations:
x = scale * np.random.normal(theta)
then this is equivalent to:
x=np.random.normal(scale * theta, scale)
which is not equivalent to the original task:
x=np.random.normal(theta, scale)

Equivalently for T2:
your implementation:
theta ~ U(-1,1), x = scale * np.random.normal(theta)
such that the multiplication of the scale with the loc (theta) indirectly modifies the prior since:
theta ~ U(-1, 1), x = np.random.normal(scale * theta, scale) <-> theta ~ U(-scale, scale), x=np.random.normal(theta, scale)

Sorry for the confusion with variance and scale and T1 and T2, but I hope my point/question is now a little bit clearer.

Regarding the paper that mentions the multivariate dataset, I mean the BayesFlow paper section 3.3. Do you also provide the BayesFlow model settings somewhere for that evaluation task?

[stefanradev93]

You are exactly right. I wanted to avoid using the np.random.multivariate_normal and leverage the fact that the distros in tasks T1 and T2 factorize, but overlooked the scale <-> variance conversion. :-D

How do we "fix" this if the SBIB implementation / description differs as well? Should we call the optional argument diagonal_var and set it to 0.1?

Regarding the multivariate dataset, you can find the code in the super legacy repo;

https://github.com/stefanradev93/cINN/blob/master/MVN_Gaussian_cINN_Example.ipynb

However, it would be much easier if you use the D and n_obs arguments from the benchmark model T1 to get a similar model. We could also think of a T1 and T2 versions with a full covariance. I would be happy to hear your thoughts on that.

[niels-leif-bracher]

Ok, so to fix the indirect scaling of the prior when sampling the observations should be straightforward instead of:
x = scale * rng.normal(theta) do x = rng.normal(theta, scale) or x = theta + scale * rng.normal(size=theta.shape) (reparametrization trick)
scale is a scalar but rng.normal broadcasts theta and scale so there shouldn't be a problem and the resulting distribution for rvs x stays isotropic/spherical Gaussian.

Btw I found that the changes from above need to be changed for the following scripts that I looked at (might need also to be changed for other scripts; those are just the ones I looked at):
https://github.com/stefanradev93/BayesFlow/blob/5ff1c1ea8ee2bda8cf333ae2890b95d05185af23/bayesflow/benchmarks/gaussian_linear_uniform.py#L84-L85

https://github.com/stefanradev93/BayesFlow/blob/5ff1c1ea8ee2bda8cf333ae2890b95d05185af23/bayesflow/benchmarks/gaussian_linear.py#L82-L83

https://github.com/stefanradev93/BayesFlow/blob/5ff1c1ea8ee2bda8cf333ae2890b95d05185af23/bayesflow/benchmarks/gaussian_mixture.py#L93-L94

To fix the inconsistency sbib has in their implementation/paper, I guess I would stick to their implementation since the benchmark results in their paper are result of using a scale of 0.1, although 0.1 was given as variance (=> scale~~0.31...) in their appendix. If they want to fix it, I guess the easiest way would be to change the scale values (or description: scale = 0.1, x|theta ~ N(theta, scale^2 * 1) in the paper and releasing a new version instead of changing the scale in the code to match the description in the paper since that would mean they also need to rerun the experiments to yield the benchmark results for the new datasets.
regarding the links in your implementation, e.g.:
https://github.com/stefanradev93/BayesFlow/blob/5ff1c1ea8ee2bda8cf333ae2890b95d05185af23/bayesflow/benchmarks/gaussian_linear_uniform.py#L21
it might need to be updated once there eventually is a new version of the paper, or you can add an additional link to their implementation. I don't know of an elegant solution, unfortunately.

Let me know what you think!

Thank you for providing the link to the legacy repo. I am currently trying to reproduce your results from BaysFlow but am not using this repo here but an implementation for the flow using freia and to make sure everything works as expected I wanted to reproduce your results on a couple of benchmarks where numbers were provided :)

[stefanradev93]

Thanks for spotting the bugs! These are now fixed in the dev branch (soon to become main) and pointers to potential inconsistencies between paper and the two implementations have been added. I will resolve the issue upon a "thumbs up".