# Non-linear regression for cumulative distribution function

I have twenty probability distributions based on a simulation. The corresponding cumulative distribution plot for one distribution looks like this:

Simulated result

I believe that most of the results will look something like this, not necessarily symmetric though.

I would like to have an approximate simple analytic formula for this curve which has small errors around the percentiles 10% and 90%, and for the median. It would be nice if the function could be defined with as few points as possible, maybe the three points.

I have constructed an approximation with a 6th order polynomial. I am not very satisfied with this as I need too many points (101) to have small errors in the extremes, as shown in the picture below.

Polynomial regression

I also tried with a sigmoid function

$$\frac{1}{1+e^{k(x-x_{\text{med}})}}$$

where $k$ is a constant and $x_{\text{med}}$ is the median of the curve. It looks smooth, but I had to adjust the constant $k$ manually. And as far as I can see, this function will be symmetric. The simulated distributions may have skewness.

Hoping for some suggestions for possible functions. Thanks in advance.

EDIT:
As Neal points out, I can determine the constant $k$ by calculating the slope at the median. The problem with the sigmoid function is that it doesn't handle skewness. User121049 proposes that the generalized logistic function

$$Y(t) = A + \frac{K-A}{(C+Qe^{-B(t-M)})^{1/\nu}}$$

might be an option, but I have problems determining the constants.

I added the data for two of the distributions for the percentiles 0% to 100% with 1% interval. I multiplied the numbers with 100 to make them better looking for sharing. In the final model I do not want to extract these many percentiles per distribution. Hopefully it would be enough with less than 10, or to use some other parameters like the variance. The reason: I have many distributions and this is an Excel model, so extracting 101 points will make the model slow. I only use the whole data set beneath to test the approximated formula. I need a compact formula for the CDF to use further on in the model AND to share in digestible written format as a report.

D1
-48,223
-40,862
-38,091
-35,840
-34,064
-32,759
-30,909
-29,986
-28,598
-27,683
-26,887
-26,058
-25,345
-24,582
-23,994
-23,660
-23,121
-22,636
-21,747
-21,274
-20,811
-20,507
-19,718
-19,488
-18,993
-18,397
-17,898
-17,199
-16,847
-16,526
-16,046
-15,593
-15,127
-14,655
-14,173
-13,687
-13,262
-12,844
-12,439
-12,073
-11,706
-11,197
-10,648
-10,277
-9,810
-9,531
-9,091
-8,885
-8,555
-8,209
-7,819
-7,470
-7,027
-6,726
-6,317
-6,032
-5,556
-4,919
-4,615
-4,008
-3,806
-3,414
-3,118
-2,732
-2,030
-1,313
-0,803
-0,538
-0,205
0,387
0,660
1,032
1,450
1,938
2,291
2,746
3,576
4,084
4,513
5,117
6,077
6,750
7,555
8,138
9,032
9,693
10,085
10,591
11,264
11,901
12,716
13,125
13,820
14,536
15,390
16,456
17,485
19,564
21,245
25,268
39,824

D2
-31,925
-25,450
-21,897
-19,410
-18,221
-17,111
-16,028
-15,207
-14,121
-12,840
-11,846
-11,422
-10,586
-9,718
-8,817
-8,351
-7,719
-6,897
-6,573
-6,211
-5,832
-5,208
-4,715
-4,403
-3,975
-3,400
-2,779
-2,239
-1,715
-1,299
-0,811
-0,385
-0,016
0,341
0,789
1,153
1,681
2,177
2,432
2,698
3,049
3,484
3,824
4,151
4,634
5,021
5,347
5,756
6,166
6,704
7,048
7,394
7,856
8,271
8,564
9,166
9,594
10,066
10,422
10,761
11,438
11,750
12,152
12,543
12,834
13,448
13,780
14,226
14,562
15,033
15,378
15,818
16,517
16,828
17,565
18,025
18,345
19,148
19,482
20,175
20,481
21,044
21,592
22,309
23,259
23,677
24,891
25,624
26,056
26,751
27,312
27,963
28,931
29,793
31,111
32,043
33,324
34,444
37,477
40,485
52,978

• Sigmoid looks like the best initial choice. You can compute the constant $k$ analytically from the estimated slope at the median or the estimated initial exponential growth rate at the left tail. Might be able to say more about better estimating functions if you described your simulation? – Neal Jul 21 '17 at 14:17
• Thanks @Neal. I worked out the derivative from the median, which returned $k=-4f'$. With the slope from p50 and 051 I got the number for $k$. The error is greatest at p92 (3%), which is unfortunate. – Kromanen Jul 21 '17 at 14:55
• In my model I have three uniform distributions on the intervals [-1, -1/3), [-1/3,1/3] and (1/3,1]. I also have $n$ choice distributions with three probabilities each. They return a value from the uniform distributions per iteration. All of these $n$ values between -1 and 1 are multiplied with weights 1, 2 or 3, summed together and divided by the sum of the weights. This ensures that we have a number between -1 and 1. The distribution for this calculation is the one I have presented above. I also have different combination of weights which return other distributions for the same simulation. – Kromanen Jul 21 '17 at 15:19
• There is the generalised logistic curvehttps://en.wikipedia.org/wiki/Generalised_logistic_function , you may be able to get away with not generalising it this much. Also Gompertz curve. – user121049 Jul 21 '17 at 15:22
• @user121049 I looked into it, but I have trouble determining C, Q and $\nu$. I have all the percentiles of my simulation, so I should be able to use the formula. Do you have any suggestions? – Kromanen Jul 21 '17 at 16:38

The given data is graphically represented on the next figure :

First, we will look how the logistic function fit to the given data : $$(x_1,y_1),(x_2,y_2),...,(x_i,y_i),...,(x_n,y_n)$$ $$y_i\simeq\frac{1}{1-e^{k(x_i-x_{\text{m}})}} \tag 1$$ $$x_i\simeq x_{\text{m}}+\frac{1}{k}\ln\left(\frac{1}{y_i}-1\right)$$ We compute $$\quad (z_1),(z_2),...,(z_i),...,(z_n)\quad$$ with : $$\quad z_i=\ln\left(\frac{1}{y_i}-1\right)$$ Then, the points $$(x_i,z_i)$$ are plotted on the next figure (respectively in BLUE and RED for the two given examples).

If the function $$(1)$$ was perfectly convenient, the points would have been on a straight line. We observe a non-negligible deviation for small and large values of $$x$$.

This draw to add a corrective term which has to be a symmetrical odd function. The simplest one has the form $$\quad \alpha (x-c)^3\quad$$ where $$\alpha$$ is a small coefficient to be determined. $$c$$ is directly found in the given data for $$y(c)=0.5$$ Don't confuse $$c$$ with $$x_m$$ above, even if they are on same order of magnitude.

For the first data : $$c=-7.819$$ and for the second : $$c=7.048$$ (no adjustment is necessary).

We see on the figure that with this kind of corrective term, the points (plotted in BLACK) can become nicely aligned.

In fact, the proposed function is : $$y(x)\simeq\frac{1}{1-e^{k(x-x_{\text{m}})+\alpha(x-c)^3}} \tag 2$$ where there are three parameters to adjust : $$k$$ , $$x_{\text{m}}$$ and $$\alpha$$.

What is more, the computation of those three parameters is very easy, in fact a simple linear regression ( no need for recursive calculus, no initial guess).

Consider the data : $$(x_1,z_1),(x_2,z_2),...,(x_i,z_i),...,(x_n,z_n)$$ and the linear relationship (with the above known value of $$c$$ ) : $$z=kx+\beta+\alpha (x-c)^3$$ where $$\beta=-kx_{\text{m}} \quad\to\quad x_{\text{m}}=-\frac{\beta}{k}$$

An usual linear regression for $$k$$ , $$\beta$$ , $$\alpha$$ leads to the wanted parameters of equation (2).

The result is shown on the next figure :

Of course, no need to take all the digits given by the computer. Only three or four significant digits are largely sufficient.

Note : The value of $$c$$ is not critical. It comes from the $$50^{th}$$ point in the given data. But one can take any other point around. For example for the first data, instead of $$-7.819$$ one can take $$7$$ or $$8$$ without a signifiant change of the final fitting.

Note: In the regression calculus, the points $$(x_0,y_0=0)$$ and $$(x_{100},y_{100}=1)$$ are excluded since they are obviously deviant for finite value of $$x$$.

• Excellent! An elegant solution. A few questions: Why do you choose c to be the median? Could it be any value in the interval of the dataset? I see that I had a typo in the question (minus-sign before the exponential function in the logistic function). You also have the minus-sign, but it seems that you change the sign in the derivation of the linear function to the correct one. Is this correct? – Kromanen Jul 23 '17 at 22:48
• Do you know how robust this is for different distributions? I am using 11 points (from p5 to p95) to determine the parameters. Is there any guidelines for number of point and which ones to use? – Kromanen Jul 23 '17 at 22:55
• It is not necessary that $c$ be exactly the median. On figure 2 one see that the deviation from the straight line occurs almost on a symmetrical manner for large negative and positive $x$. So the corrective function should be almost symmetric. That is why $\alpha(x-c)^3$ is convenient. But if $c$ is not exactly the median (and in fact $c$ is not the median since we take a value $c$ before computing the median $x_m$), then later, the value computed for $\alpha$ will be slightly different, which counterbalances and leads to quite the same fitting. – JJacquelin Jul 24 '17 at 6:14
• I don't understand why you use 11 points (from p5 to p95). Why not using all ? The more the number of points is large, the best is the accuracy. So my guideline is to take all experimental points, of course except the points obviously out of range or the points which are known to have a too low accuracy in measurements. Note that the above method is especially valuable because it well takes account of the points in the ranges of large negative and large positive $x$. Instead of the fictional examples with 100 points, you should have given the examples with 11 points that you really use. – JJacquelin Jul 24 '17 at 6:35
• About the sign of $k$ : It doesn't matter. If you put $+$ in the initial formula (as you did), the resulting value is negative (see it on figure 3). If you put $-$ in the formula, with exactly the same procedure, the resulting value will be positive. So, don't worry about that, we get the same result in both cases. – JJacquelin Jul 24 '17 at 6:43

A session in R statistical software illustrates my Comment on the Kolmogorov-Smirnov test and related topics. Similar functions and procedures are available in other software packages.

(1) Simulate $n = 1000$ values known to be from $\mathsf{Norm}(0,1).$ Then look at: (a) An empirical CDF (ECDF) to visually assess fit of that ECDF to the standard normal CDF (dotted red). (b) A normal probability plot (also called Q-Q plot) in which the vertical scale of the ECDF is transformed so that normal data should fall very nearly along a straight line. [Retain the set.seed statement to repeat exactly the same simulated sample; omit it for a fresh simulation. It is best to ignore quirky behavior of a few points at the extremes.]

set.seed(1234);  x = rnorm(1000)
par(mfrow=c(1,2))
plot(ecdf(x), main="(a) ECDF")
qqnorm(x, main="(b) Normal Q-Q Plot")
par(mfrow=c(1,1))


Roughly speaking, the Kolmogorov-Smirnov test statistic finds the maximum distance $D$ between the ECDF and the CDF. (A suitable transformation makes the distance independent of the mean and variance of the normal distribution.)

ks.test(x, pnorm)

One-sample Kolmogorov-Smirnov test

data:  x
D = 0.026371, p-value = 0.4901
alternative hypothesis: two-sided


A P-value smaller than 0.05 indicates evidence of a poor fit of the ECDF to the hypothetical CDF.

(2) Simulate n = 1000 values known not to be normal. Here $T_i \sim \mathsf{T}(df=3)$ and $X_i = T_i/\sqrt{3}$ so that $E(X_i) = 0$ and $SD(X_i) = 1.$ Observe the results of the same procedures as above for non-normal data.

The following code is used for simulation. The rest of the code is as in (1).

set.seed(1235);  x = rt(1000, 3)/sqrt(3)


Notice the poor fit of the ECDF to the standard normal density (dotted red). Also, the distinctly nonlinear Q-Q plot.

ks.test(x, pnorm)

One-sample Kolmogorov-Smirnov test

data:  x
D = 0.098214, p-value = 8.367e-09
alternative hypothesis: two-sided


The very small P-value indicates that there is almost no chance that actual standard normal data could have produced these 1000 observations.

Notes: In (2) I chose data that would make an ECDF that looks 'sigmoidal' in shape, but is clearly not the standard normal CDF. Less extreme departures from normality would be more difficult to detect graphically, but, with reasonably large sample sizes, the Kolmogorov-Smirnov test has good power to detect nonnormal data.

The K-S test in R requires one to know the mean and SD of the hypothetical normal distribution and is limited to a few thousand observations.

Another good test of normality is the Shapiro-Wilk test (shapiro.test in R), which is based on different theoretical methods.

• Thanks, very interesting. From now on I will definitely use Q-Q plots to see if my theoretical formula matches the data. I am not quite sure how I can use this to create a formula, though. And I am using Excel with RiskAmp as add-in to build my model. No way around that, unfortunately. – Kromanen Jul 22 '17 at 13:55
• Part of the message of my Answer is that it is hard to imagine you can find a truly useful rational approx to the CDF of a simulated distribution using only a few hundred points. You might get a more useful answer if you describe the underlying purpose of what you are doing. Perhaps on our statistics site because this project seems pretty far from traditional mathematical statistics. – BruceET Jul 22 '17 at 16:09