Propagating asymmetric and large errors in a ratio of likelihood-fit parameters

Question

I'm currently trying to propagate the errors of two values:

$A=143^{+28.6}_{-26.6}$ , $B=19^{+10.9}_{-8.7}$ $\rightarrow$ $B/A = 0.133^{+?}_{-?}$

Both values A and B (#counts) are derived from a log-likelihood model fit on poisson data. The asymmetrical errors are due to the low number of events, which means that the Poisson-$\chi^2$ distribution (and thus the errors) is slightly asymmetric. I also know the numerical chi2 (likelihood) distribution for A and B.

Now, there are multiple major problems that I'm facing in the error propagation:

1) The errors are asymmetric. Of course you can propagate the upper/lower error for both separately, but there's no theoretical justification for that. However, sometimes even Barlow recommends it, despite that fact (http://www.slac.stanford.edu/BFROOT/www/Statistics/Report/report.pdf).

2) The errors are large (50%), which breaks assumptions for standard error propagation.

3) Even if you consider the errors as symmetric and small, calculating the ratio of two gaussian distributions can still result in an asymmetric ratio distribution.

The major struggle that I have right now is to incorporate all of those three issues into one method. If this is very time demanding, I could also assume the errors to be symmetric (because the difference isn't sooo big). However, even then, the problem of 2) and 3) still exist.

I've thought about combining the likelihood of A and B to C (C=A/B), but unfortunately that led to nothing.

Note: This is a crosspost from CrossValidated, because I feel that this community should be a better place for the problem (i.e. different error propagation is needed for a mean of some sample data in e.g. biology compared to the values based on a likelihood fit).

Answer 1

The issues you mention in the question are quite technical. So I disagree with your assessment that Physics SE is a better place for this question than Cross Validated SE.

The usual physics treatment of such issues is probably to make a reasonable approximation. The asymetry is small, and the errors are quoted to 3 sig figs, which does not seem justified. It is rarely appropriate to use more than 1 sig fig. From a physics point of view there is little (if anything) to be gained from using a statistically correct method for propagating the error. The final error, rounded to 1 sig fig, is unlikely to be any different whichever method is used.

I think your time would be better spent on trying to improve the accuracy of the result, by reviewing the experimental method, rather than trying to improve the accuracy of the error bounds.

Answer 2

From my perspective this is a perfect example when bootstrapping should be used. It relies only one few assumptions and yields great results if we are interested in the main body of the distribution. The main idea is that we treat the dataset as population and resample from it with replacement. By generating e.g. 2000 bootstrap samples we generate a bootstrap distribution of the quantity of interest -- in your case $A/B$.

Here R code:

generateFakeData = TRUE
if ( generateFakeData ){
    # Generate fake data, which is an asymmetric distribution around the average 
    # value. You don't need this step, because you posses the dataset. 
    nSim = 1e4
    
    # 1. Approx. the asymmetric standard deviation using a chi^2-distribution with 
    #    the correct degrees of freedoms.
    #       chi2.sd = sqrt(2*df) 
    #       => df = ...
    sd.A   = 27
    df.A   = round(sd.A^2/2)
    data.A = rchisq(n=nSim, df=df.A)
    # 2. Generate the fake dataset around the average value:
    muA    = 143
    data.A = data.A - mean(data.A) + muA
    data.A = abs(data.A) # ensure no neg. values
    mean(data.A)
    sd(data.A)
    
    # Same for B
    sd.B   = 10
    df.B   = round(sd.B^2/2)
    data.B = rchisq(n=nSim, df=df.B)
    # 2. Generate the fake dataset around the average value:
    muB    = 19
    data.B = data.B - mean(data.B) + muB
    data.B = abs(data.B) # ensure no neg. values
    mean(data.B)
    sd(data.B)
    
    ## Combine fake data into a dataset:
    DATA = data.frame(A = data.A, 
                      B = data.B)
}

## Bootstrapping:
# 1. Define a helper function:
bootFct = function(dataFrame, idx){
    # For each bootstrap sample, we only return the statistic of interest:
    # The code is "step-by-step", because not everybody is used to read code.
    A     = dataFrame[idx,1] 
    B     = dataFrame[idx,2]
    ratio = A/B
    out   = mean(ratio)
    return(out)
}
# 2. Use the helper function to obtain 2000 bootstrap samples:
bootSample = boot::boot(data=DATA, statistic=bootFct, R=2000)

## Plot a histogram of the bootstrap distribution:
library(ggplot2)
DF = data.frame(ratio = bootSample$t)
gg = ggplot(DF, aes(x=ratio)) +
     geom_histogram( color="black", alpha=0.3, bins=30)
print(gg)

## Calc CI assuming normal bootstrap distribution:
bootSample.SD     = sd(bootSample$t)
data.var.CI.upper = bootSample$t0 + 1.96*bootSample.SD
data.var.CI.lower = bootSample$t0 - 1.96*bootSample.SD
print(paste0('[', round.smartly(data.var.CI.lower), ', ', 
             round.smartly(data.var.CI.upper), ']'))

## Calc CI (non-parametric):
quantile(bootSample$t, probs=c(0.05, 0.95))