I'd like to highlight a problem I encountered while trying to fit a multi-parameter null hypothesis to a dataset. I've resolved this issue in a somewhat messy manner in my fork of Bioverse, but there is probably a much cleaner way how to do it that would be of interest to the rest of the users of Bioverse.
When I was trying to fit a null hypothesis with 3 parameters I encountered an error that occurred in the calculation of the optimum parameters used to calculate the AIC and BIC in hypothesis.py. See the line below:
Here the optimum parameters for the null hypothesis are being computed as the average of the array of measurements. At first I didn't understand why the average of Y was being used for theta_opt_null, but I realized that for the very specific case of a null hypothesis defined by the function f_null (i.e. a flat function of theta), the mean of Y would be a reasonable estimate of the optimum value of theta. For any other null hypothesis, this value for theta_opt_null is incorrect. When theta for the null hypothesis is more than one parameter, the current calculation only give a scalar value which causes a crash in the computation of the AIC.
I resolved this problem in my fork by using the output of the dynesty fit for the null hypothesis to calculate the maximum likelihood parameters for the model. I modified the sample_posterior_dynesty function to return the dynesty chain for the null hypothesis, when previously it was only returning chain for the alternative hypothesis. I also modified the dynesty fitting for the null hypothesis which previously seemed fine tuned to a 1 parameter model and had more live points and stricter convergence criteria than the alternative hypothesis.
I use the output of the null hypothesis dynesty chain to calculate the maximum likelihood parameters and obtain the AIC and BIC. I also modified the output results dictionary to include the results from fitting the null hypothesis as it was helpful to diagnose whether the dynesty chains were converging.
Currently it only works for dynesty fitting. I'm realizing that the function for emcee fitting in sample_posterior_emcee isn't fitting for the null hypothesis, only the alternative hypothesis. I think the code approximating theta_opt_null as the average of Y may be a relic of trying to get dAIC with emcee without actually using the results of fitting the null hypothesis model.
I'd be interested to hear what you all think would be the best way to resolve this problem in the code. My fix works for dynesty, but maybe there is a simple way to get it to work for the other methods.
I'd like to highlight a problem I encountered while trying to fit a multi-parameter null hypothesis to a dataset. I've resolved this issue in a somewhat messy manner in my fork of Bioverse, but there is probably a much cleaner way how to do it that would be of interest to the rest of the users of Bioverse.
When I was trying to fit a null hypothesis with 3 parameters I encountered an error that occurred in the calculation of the optimum parameters used to calculate the AIC and BIC in
hypothesis.py. See the line below:https://github.com/danielapai/bioverse/blob/c5b1bafe7ca70b54e2040bb7764098c8ead504e6/bioverse/hypothesis.py#L294
Here the optimum parameters for the null hypothesis are being computed as the average of the array of measurements. At first I didn't understand why the average of Y was being used for
theta_opt_null, but I realized that for the very specific case of a null hypothesis defined by the functionf_null(i.e. a flat function of theta), the mean of Y would be a reasonable estimate of the optimum value of theta. For any other null hypothesis, this value fortheta_opt_nullis incorrect. When theta for the null hypothesis is more than one parameter, the current calculation only give a scalar value which causes a crash in the computation of the AIC.I resolved this problem in my fork by using the output of the dynesty fit for the null hypothesis to calculate the maximum likelihood parameters for the model. I modified the
sample_posterior_dynestyfunction to return the dynesty chain for the null hypothesis, when previously it was only returning chain for the alternative hypothesis. I also modified the dynesty fitting for the null hypothesis which previously seemed fine tuned to a 1 parameter model and had more live points and stricter convergence criteria than the alternative hypothesis.I use the output of the null hypothesis dynesty chain to calculate the maximum likelihood parameters and obtain the AIC and BIC. I also modified the output results dictionary to include the results from fitting the null hypothesis as it was helpful to diagnose whether the dynesty chains were converging.
Here is a link to my commit that resolved these problems in my fork: https://github.com/danielapai/bioverse/commit/6972c24e353fc9edfcb511035b40d166680fc473
Currently it only works for dynesty fitting. I'm realizing that the function for emcee fitting in
sample_posterior_emceeisn't fitting for the null hypothesis, only the alternative hypothesis. I think the code approximating theta_opt_null as the average of Y may be a relic of trying to get dAIC with emcee without actually using the results of fitting the null hypothesis model.I'd be interested to hear what you all think would be the best way to resolve this problem in the code. My fix works for dynesty, but maybe there is a simple way to get it to work for the other methods.