The form of the equation that I used to fit the characteristics curve just involved variables in the exponents. It is by no way unique as can by the fact that a multiplying current can be canceled out by a constant in an exponent. The fit proved challenging since my computer kept encountering singularities and it was difficult to find a starting point for a solution.
I was able to get a fairly good fit by fitting the upper and lower branches of the curve separately. The lower part is just the regular current for a reversed diode and an exponential curve is a good approximation. The upper branch includes in addition the breakdown current for the diode so one can fit the difference between the measured values and the extrapolated fit for the lower part of the curve. I also tried a linear estimate of the zero* of the variance but that proved rather slow. A quadratic estimate of the zero of the variance proved to be quicker after one got away from β = 0 where there was a singularity.
The fit shown started with a linear estimate of the zero and then switched to a quadratic estimate for faster convergence. It confirms the presence of two currents.
*edit: It was the zero of the least squares fit for which the variance is a minimum.
**edit: I corrected for the voltage drop across the milliammeter and got a steeper slope above 5.1 V and the rms deviation, δ'_rms, reduced to 0.320 mA which probably means that there is better agreement with the assumption of two currents. The diode may also have some internal resistance.