Ok, I see I read the figure 3b wrong (only glanced at it before). I extracted the data from their extended figure 5a: dV sigma R_major -4.16 16.41

Ok, I see I read the figure 3b wrong (only glanced at it before). I extracted the data from their extended figure 5a: dV sigma R_major -4.16 16.41 -3.98 10.55 3.41 -1.97 15.03 7.25 -1.85 -2.24 5.12 -0.74 1.17 5.75 0.29 -14.39 5.97 1.16 -40.39 11.3 1.52 2.24 6.61 1.58 -3.52 13.43 4.49 -1.39 10.23 6.77 Then ran this simulation in R: # Used for biweight midvariance calculation (r.bw) require(asbio) # Utils: ## Simulation function: ### "For a given value of sigma_test we generated 1000 samples of 10 velocities, ### distributed according to a Gaussian of width sigma_test. The ten velocities ### in each sample were then perturbed with errors, drawn from Gaussians with ### widths equal to the empirically-determined uncertainties in the measured ### dispersions. run_sim = function(nSim, n, sigma_test, err){ replicate(nSim, rnorm(n, 0, sigma_test) + rnorm(n, 0, err)) } ## Clean up the biweight midvariance output r.bw2 = function(x){ as.numeric(sqrt(r.bw(x))) } # Scraped data roughly matches what is reported in paper ## MAD (~4.7) s_obs_nmad = mad(dat$dV) ## biweight (~8.4) s_obs_bi = r.bw2(dat$dV) ## RMS (~14.3) s_obs_rms = sqrt(mean(dat$dV^2)) # Run simulation to test sigma values in .1 km/s intervals ## "Using the biweight estimator, “measured” dispersions sigma_obs,test were calculated ## for all samples. If the value 8.4 is within the 16th – 84th percentile of the ## distribution of sigma_obs,test then sigma_test is within the + /-1 sigma uncertainty ## on sigma_intr. This method gives sigma_intr = 3.2 (+5.5, -3.2) kms^-1." sigmas = seq(.1, 30, by = .1) nSim = 1e4 out = as.data.frame(matrix(nrow = length(sigmas), ncol = 2)) colnames(out) = c("sigma_test", "s_obs_bi_percentile") for(i in 1:length(sigmas)){ # Calculate biweights for each sim result sim_res = run_sim(nSim, 10, sigmas[i], dat$sigma) sigma_obs_test = apply(sim_res, 2, function(x) r.bw2(x)) # Output sigma_test and the percentile of 8.4 out[i, ] = c(sigmas[i], ecdf(sigma_obs_test)(8.4)) # Print progress print(out[i, ]) } # Check which sigma_test values are within the accepted percentile range out$check = out$s_obs_bi_percentile > .16 & out$s_obs_bi_percentile < .84 # Plot the results png("results.png", width = 960, height = 960) plot(out[, 1], out[, 2], xlab = "sigma_test", ylab = "Percentile(8.4)", panel.first = grid()) dev.off() I get that their observed sigma_biweight = 8.4 has percentile of p = .0146 for sigma_test = 20 (https://image.ibb.co/ncRF27/results.png https://image.ibb.co/ncRF27/results.png). In their terminology that would mean it is just within the 99% CI, ie had odds of happening of ~ 1/70 if MOND were correct. I believe I replicated their results since it had percentile of p = .6336 for sigma_test = 0.1 and p = .1607 for sigma_test = 8.6, which matches to their sigma_int = 3.2 (+5.5, -3.2) using 16th - 84th percentiles result. Where did you get 1/100,000?