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

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).

• Quoting and combining error bars is not appropriate in this situation. You need to calculate the full probability distribution of C using a Monte Carlo simulation. Commented Jan 14, 2017 at 10:17
• Thank you for the comment! I guess that MC-generating A and B based on their probability (=exp(DeltaL)) and then calculating C should be sufficient?
– 0vbb
Commented Jan 15, 2017 at 2:58
• It would be helpful to have more details of your experiment and the meaning of A and B. Commented Jan 15, 2017 at 3:12
• I am counting events (1D energy histo) in a certain energy range. I know that the physical distribution of these events is the sum of two underlying physical processes. Thus I fit (Binned Likelihood) the counts in the data with two model MC-spectra in a joint fit. As a result, I get A and B based on the fit. However, the A and B that I want to compare (C=A/B) are from different, independent data sets.
– 0vbb
Commented Jan 15, 2017 at 3:20

## 2 Answers

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.

• Referring to the crosspost, I just realized that the appropriate error propagation depends on how you derive the initial values (A/B by a fit). I.e., I can't use Fieller's theorem to propagate large errors. Thank you also for the second point, although improving the accuracy will most likely result in waiting for 1 year of more data in my case.
– 0vbb
Commented Jan 15, 2017 at 3:07

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))