File: ecdf.c

package info (click to toggle)
fitspng 1.4-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 1,804 kB
  • sloc: sh: 1,232; ansic: 860; makefile: 14
file content (112 lines) | stat: -rw-r--r-- 2,593 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
/*

  FITSPNG    -- empirical distribution function
  Copyright (C) 2019  Filip Hroch, Masaryk University, Brno, CZ

  Fitspng is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  Fitspng is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with Fitspng.  If not, see <http://www.gnu.org/licenses/>.

*/

#include "fitspng.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <string.h>

#include <assert.h>

#define EPSILON 2*FLT_EPSILON
#define CMP(x,y) ((x) < (y) ? -1 : ((x) > (y) ? 1 : 0))

static int fcmp(const void *u, const void *v)
{
  return CMP(*(const float *) u, *(const float *) v);

  /*
  const float x = *(const float *) u;
  const float y = *(const float *) v;
  if( x < y )
    return -1;
  else if( x > y )
    return 1;
  else
    return 0;
  */
}

/* alternative with modified data */
int ecdf(long n, const float *data, float *xcdf, float *ycdf)
{
  assert(n > 0 && data);

  long m = n * sizeof(float);
  float *d;
  if( (d = malloc(m)) == NULL )
    perror("There is no room for data sort in CDF.");

  memcpy(d,data,m);
  qsort(d,n,sizeof(float),fcmp);

  /* remove duplicities */
  m = 1;
  for(long i = 1; i < n; i++) {
    if( fabsf(d[i] - d[m]) > EPSILON )
      d[m++] = d[i];
  }

  /* setup CDF */
  memcpy(xcdf,d,m*sizeof(float));

  float h = 1.0 / (float)(m + 1);
  for(long i = 0; i < m; i++) ycdf[i] = i*h;

  free(d);
  return m;
}

float quantile(long ncdf, const float *xcdf, const float *ycdf, float q)
{
  long n = ncdf;

  if( n == 0 )
    return 0.0;
  else if( n == 1 )
    return xcdf[0];

  if( q < ycdf[0] )
    return xcdf[0];
  else if( q > ycdf[n-1] )
    return(xcdf[n-1]);
  else {
    float h = 1.0 / (float)(n + 1);
    float r = q / h;
    int m = roundf(r);
    if( fabsf(m - r) < EPSILON )
      return(xcdf[m]);
    else {
      int low = r;
      int high = low + 1;

     float dy = ycdf[high] - ycdf[low];
      if( fabsf(dy) > EPSILON )
        // inverse by linear interpolation
        return((xcdf[high] - xcdf[low])/dy*(q - ycdf[low]) + xcdf[low]);
      else {
        // nearly singular
        return((xcdf[high] + xcdf[low]) / 2);
      }
    }
  }
}