1
$\begingroup$

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  
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Neal Jul 21 '17 at 14:17
  • $\begingroup$ 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. $\endgroup$ – Kromanen Jul 21 '17 at 14:55
  • $\begingroup$ 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. $\endgroup$ – Kromanen Jul 21 '17 at 15:19
  • $\begingroup$ 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. $\endgroup$ – user121049 Jul 21 '17 at 15:22
  • $\begingroup$ @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? $\endgroup$ – Kromanen Jul 21 '17 at 16:38
1
$\begingroup$

The given data is graphically represented on the next figure :

enter image description here

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

enter image description here

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 :

enter image description here

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

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Kromanen Jul 23 '17 at 22:48
  • $\begingroup$ 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? $\endgroup$ – Kromanen Jul 23 '17 at 22:55
  • $\begingroup$ 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. $\endgroup$ – JJacquelin Jul 24 '17 at 6:14
  • $\begingroup$ 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. $\endgroup$ – JJacquelin Jul 24 '17 at 6:35
  • $\begingroup$ 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. $\endgroup$ – JJacquelin Jul 24 '17 at 6:43
1
$\begingroup$

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")
    curve(pnorm(x), col="red", lty="dashed", add=T)
  qqnorm(x, main="(b) Normal Q-Q Plot")
par(mfrow=c(1,1))

enter image description here

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)

enter image description here

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Kromanen Jul 22 '17 at 13:55
  • $\begingroup$ 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. $\endgroup$ – BruceET Jul 22 '17 at 16:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.