Prompted by change from scipy.math.factorial (v1.8.1) to scipy.special.factorial (v1.9.2) (see #15)
Other scipy functions may have changed or moved around as well

For large shape parameters, the current implementation for the beta pdf utilizes gamma functions which tend to blow up to inf when a, b are ~ 150 and larger. For example, GammaFunction(150) ~ 1e260, which pretty quickly crosses machine precision while multiplying/dividing such terms.

.... causes error for large datasets where one would like to check for large degree_max.

@JenniferBurt

Because the normal PDF is not being sampled finely enough, if we have measurement errors that are too small then the PDFs are 0 everywhere (giving us nans).

I have set a threshold of 50 sigma at the moment, but need to carefully check if this suffices for period, insolation, etc.. Maybe add some descriptive flags?

For > 2 dimensions, the current weight estimation is erroneous if unequal number of degrees are used for different dimensions.

While the other predicting functions for M-dwarfs work, this one is giving some weird issues. Double-check.

https://github.com/shbhuk/mrexo/blob/92ecf369325dd22683c37d85f8224d5ed220bd05/mrexo/datasets/MdwarfRuns/MdwarfPrediction.py#L92C1-L92C43

https://docs.python.org/3/library/decimal.html

To avoid the floating point accuracy limited non zero weights from biasing the AIC.

The _find_indv_pdf() function tends to zero when the measurement uncertainty is too small. This is especially a problem for the period dimension where the sigmas are typically < 1e-3, 1e-4. This is because at such small sigma values we need to assume only a beta distribution, and not a product of beta and normal.

Therefore we need to figure out how to handle this.. Perhaps an if statement?

Code up, and Test that this works.

There is a bug in the code when the code is run with no measurement uncertainties, i.e. if Mass_sigma or Radius_sigma = None, then the results should not make sense.

Write function to save all the outputs such as joint dist, etc..

In the AIC function we need to estimate the number of weights to penalize the log-likelihood and calculate the AIC.

While we have tried a few different thresholds to calculate weights (nonzero, w > 1e-8), we need to find something more robust

Can try this 'Design Effect' Estimated Sample Size method?
https://en.wikipedia.org/wiki/Design_effect#Effective_sample_size

@hematthi

Seeing that the corner weights and joint distribution are too high, and also asymmetric (i.e. the upper values are larger than lower)

Conditional dist plots
Joint dist plots