pearson.c File Reference

Pearson's correlation coefficient and regression line. More...

#include "libtpcmodel.h"

Go to the source code of this file.

Functions
int	pearson (double *x, double *y, int nr, double *k, double *kSD, double *b, double *bSD, double *r, double *ySD)
int	pearson2 (double *x, double *y, char *is, int nr, double *k, double *kSD, double *b, double *bSD, double *r, double *ySD)
int	pearson3 (double *x, double *y, int nr, double *k, double *kSD, double *b, double *bSD, double *r, double *ySD)
int	pearson4 (double *x, double *y, int nr, double start, double end, double *k, double *kSD, double *b, double *bSD, double *r, double *ySD)
	Calculate slope and intercept of a line and Pearson's correlation coefficient.
int	best_pearson (double *x, double *y, int nr, int min_nr, int *first, int *last, double *k, double *kSD, double *b, double *bSD, double *r, double *ySD)
int	mean (double *x, double *y, int nr, double *xmean, double *xsd, double *ymean, double *ysd)
int	regr_line (double *x, double *y, int n, double *m, double *c)
int	highest_slope (double *x, double *y, int n, int slope_n, double *m, double *c, double *xi, double *xh)
int	highest_slope_after (double *x, double *y, int n, int slope_n, double x_start, double *m, double *c, double *xi, double *xh)

Detailed Description

Pearson's correlation coefficient and regression line.

Author

Vesa Oikonen

Definition in file pearson.c.

Function Documentation

best_pearson()

int best_pearson	(	double *	x,
		double *	y,
		int	nr,
		int	min_nr,
		int *	first,
		int *	last,
		double *	k,
		double *	kSD,
		double *	b,
		double *	bSD,
		double *	r,
		double *	ySD )

Find the best linear fit to double data (x[], y[]) with nr points.

Data may contain NaN's, which are not used.

See also

highest_slope, regr_line, pearson

Returns

Returns the nr of points actually used, or 0, in case of error.

Parameters

x	Data x values
y	Data y values
nr	Number of data sample values
min_nr	Minimum nr of data points to use
first	Index [0..last-2] of the first point to use initially, and after fitting
last	Index [first+1..nr-1] of the last point to use initially, and after fitting
k	slope
kSD	S.D. of slope
b	y axis intercept
bSD	S.D. of y axis intercept
r	Pearson's correlation coefficient, r
ySD	Residual variance of y values

Definition at line 252 of file pearson.c.

  {
  int i, n, m, x1, x2, b1, b2, bm;
  double *nx, *ny;
  double tk, tkSD, tb, tbSD, tr, tySD;
 
 
  if(0) {
    fprintf(stdout, "best_pearson(x, y, %d, %d, %d, %d, k, kSD, b, bSD, r, ySD)\n",
      nr, min_nr, *first, *last);
    fflush(stdout);
  }
  /* Remove NaN's and those outside range first-last */
  /* Allocate memory for new data pointers          */
  nx=(double*)calloc(nr, sizeof(double)); if(nx==NULL) return 0;
  ny=(double*)calloc(nr, sizeof(double)); if(ny==NULL) {free((char*)nx); return 0;}
  /* Copy data to pointers */
  if(*last>nr-1) *last=nr-1;
  if(*first<0) *first=0;
  for(i=*first, n=0; i<=*last; i++)
    if(!isnan(x[i]) && !isnan(y[i])) {nx[n]=x[i]; ny[n]=y[i]; n++;}
 
  /* Check that we have enough points */
  if(n<2 || n<min_nr) {free((char*)nx); free((char*)ny); return 0;}
  if(n==min_nr) {
    i=pearson(nx, ny, n, k, kSD, b, bSD, r, ySD);
    free((char*)nx); free((char*)ny);
    if(i) return 0; else return n;
  }
 
  /* OK, let's start  */
  *k=*kSD=*b=*bSD=*r=*ySD=0.; b1=b2=bm=0;
  for(x1=0, m=n; n-x1>min_nr; x1++) {
    /*printf("I x1=%i x2=%i m=%i n=%i\n", x1, x2, m, n);*/
    for(x2=n-1, m=x2-x1+1; m>=min_nr; x2--, m--) {
      if(pearson(nx+x1, ny+x1, m, &tk, &tkSD, &tb, &tbSD, &tr, &tySD))
        continue;
      if((tr>*r) ||
          (tr==*r && m>bm) ||
          (tr==*r && m==bm && x1>b1) ||
          (tr==*r && m==bm && x1==b1 && tk>*k) ) {
        bm=m; b1=x1; b2=x2;
        *k=tk; *kSD=tkSD; *b=tb; *bSD=tbSD; *r=tr; *ySD=tySD;
      }
      /*printf("Fit range: %2i-%2i  ;  best fit: %2i-%2i\n", x1, x2, b1, b2);*/
    }
  }
 
  /* Set first&last  */
  for(i=*first; i<=*last; i++)
    if(x[i]==nx[b1] && y[i]==ny[b1]) {*first=i; break;}
  for(i=*last; i>=*first; i--)
    if(x[i]==nx[b2] && y[i]==ny[b2]) {*last=i; break;}
  /*printf("FIRST=%i  LAST=%i\n", *first, *last);*/
 
  free((char*)nx); free((char*)ny);
  return bm;
}

highest_slope()

int highest_slope	(	double *	x,
		double *	y,
		int	n,
		int	slope_n,
		double *	m,
		double *	c,
		double *	xi,
		double *	xh )

Finds the regression line with the highest slope for x,y data.

See also

highest_slope_after, regr_line, best_pearson

Returns

Return 0 if ok.

Parameters

x	An array of x axis values
y	An array of y axis values
n	The number of values in x and y arrays
slope_n	The number of samples used to fit the line
m	Pointer where calculated slope is written; NULL if not needed
c	Pointer where calculated y axis intercept is written; NULL if not needed
xi	Pointer where calculated x axis intercept is written; NULL if not needed
xh	Pointer where the place (x) of the highest slope is written; NULL if not needed

Definition at line 424 of file pearson.c.

  {
  int i, ret, i_at_max=0;
  double slope, ic, max_slope, ic_at_max=0.0;
 
  /* Check the data */
  if(x==NULL || y==NULL) return(1);
  if(n<2) return(2);
  if(slope_n>n) return(3);
 
  /* Compute */
  max_slope=-1.0E200;
  for(i=0; i<=n-slope_n; i++) {
    ret=regr_line(x+i, y+i, slope_n, &slope, &ic); if(ret) continue;
    if(slope>max_slope) {max_slope=slope; ic_at_max=ic; i_at_max=i;}
  } /* next slope */
  if(max_slope==-1.0E200) return(10);
  if(m!=NULL) *m=max_slope;
  if(c!=NULL) *c=ic_at_max;
  if(xi!=NULL) {if(max_slope!=0.0) *xi=-ic_at_max/max_slope; else *xi=0.0;}
  if(xh!=NULL) {
    *xh=0.0; for(i=i_at_max; i<i_at_max+slope_n; i++) *xh+=x[i];
    *xh/=(double)slope_n;
  }
 
  return(0);
}

Referenced by dftFixPeak().

highest_slope_after()

int highest_slope_after	(	double *	x,
		double *	y,
		int	n,
		int	slope_n,
		double	x_start,
		double *	m,
		double *	c,
		double *	xi,
		double *	xh )

Finds the regression line with the highest slope for x,y data after specified x.

See also

highest_slope, regr_line

Returns

Return 0 if ok.

Parameters

x	An array of x axis values.
y	An array of y axis values.
n	The number of values in x and y arrays.
slope_n	The number of samples used to fit the line.
x_start	Estimation start x value, samples with smaller x are ignored; can usually be set to zero.
m	Pointer where calculated slope is written; NULL if not needed.
c	Pointer where calculated y axis intercept is written; NULL if not needed.
xi	Pointer where calculated x axis intercept is written; NULL if not needed.
xh	Pointer where the place (x) of the highest slope is written; NULL if not needed.

Definition at line 475 of file pearson.c.

  {
  int i, ret, i_at_max=0;
  double slope, ic, max_slope, ic_at_max=0.0;
 
  /* Check the data */
  if(x==NULL || y==NULL) return(1);
  if(n<2) return(2);
  if(slope_n>n) return(3);
 
  /* Compute */
  max_slope=-1.0E200;
  for(i=0; i<=n-slope_n; i++) if(x[i]>=x_start) {
    ret=regr_line(x+i, y+i, slope_n, &slope, &ic); if(ret) continue;
    if(slope>max_slope) {max_slope=slope; ic_at_max=ic; i_at_max=i;}
  } /* next slope */
  if(max_slope==-1.0E200) return(10);
  if(m!=NULL) *m=max_slope;
  if(c!=NULL) *c=ic_at_max;
  if(xi!=NULL) {if(max_slope!=0.0) *xi=-ic_at_max/max_slope; else *xi=0.0;}
  if(xh!=NULL) {
    *xh=0.0; for(i=i_at_max; i<i_at_max+slope_n; i++) *xh+=x[i];
    *xh/=(double)slope_n;
  }
 
  return(0);
}

mean()

int mean	(	double *	x,
		double *	y,
		int	nr,
		double *	xmean,
		double *	xsd,
		double *	ymean,
		double *	ysd )

Calculates the mean and SD of data. Data (y data) may contain NaN's.

See also

dmean, dmedian, dmean_nan, regr_line, pearson

Returns

Returns !=0 in case of an error.

Parameters

x	Data x values
y	Data y values
nr	Number of data sample values
xmean	Calculated x mean
xsd	Calculated SD of x mean
ymean	Calculated y mean
ysd	Calculated SD of y mean

Definition at line 341 of file pearson.c.

  {
  int   i, n;
  double sumsqr, sqrsum;
 
  for(i=0, sumsqr=sqrsum=0.0, n=0; i<nr; i++) if(!isnan(x[i]) && !isnan(y[i])) {
    sumsqr+=x[i]*x[i]; sqrsum+=x[i]; n++;}
  if(n<=0) return -1;
  *xmean=sqrsum/(double)n;
  if(n==1) *xsd=0.0; else {
    sqrsum*=sqrsum;
    *xsd=sqrt( (sumsqr - sqrsum/(double)n) / (double)(n-1) );
  }
  for(i=0, sumsqr=sqrsum=0.0, n=0; i<nr; i++) if(!isnan(x[i]) && !isnan(y[i])) {
    sumsqr+=y[i]*y[i]; sqrsum+=y[i]; n++;}
  if(n<=0) return -1;
  *ymean=sqrsum/(double)n;
  if(n==1) *ysd=0.0; else {
    sqrsum*=sqrsum;
    *ysd=sqrt( (sumsqr - sqrsum/(double)n) / (double)(n-1) );
  }
  return 0;
}

Referenced by dftMeanTAC(), doubleMatchRel(), imgsegmClusterExpand(), least_trimmed_square(), mertwiRandomExponential(), noiseSD4SimulationFromDFT(), and resMean().

pearson()

int pearson	(	double *	x,
		double *	y,
		int	nr,
		double *	k,
		double *	kSD,
		double *	b,
		double *	bSD,
		double *	r,
		double *	ySD )

Calculate slope and intercept of a line and Pearson's correlation coefficient.

See also

regr_line, best_pearson, highest_slope, pearson2, pearson3, pearson4, mean, imgsegmPearson, nlopt1D

Returns

If successful, function returns value 0.

Parameters

x	data x values
y	data y values
nr	number of data sample values
k	slope
kSD	S.D. of slope
b	y axis intercept
bSD	S.D. of y axis intercept
r	Pearson's correlation coefficient, r
ySD	Residual variance of y values

Definition at line 14 of file pearson.c.

  {
  int i;
  double e, f, g, meanx, meany, kk, bb, ke, be, rr, sumsdcy=0.0;
  double sumx=0.0, sumy=0.0, sumsx=0.0, sumsdx=0.0, sumsdy=0.0, sumdxdy=0.0;
  double sumxy=0.0, sumsy=0.0;
 
 
  if(0) {
    fprintf(stdout, "pearson(x[], y[], %d, *k, *kSD, *b, *bSD, *r, *ySD)\n", nr);
    fflush(stdout);
  }
  /* Check that there is some data */
  if(x==NULL || y==NULL || nr<2) return(1);
 
  /* If only 2 datapoints */
  if(nr==2) {
    f=x[1]-x[0]; if(fabs(f)<1.0E-50) return(1);
    *k=(y[1]-y[0])/f; *b=y[0]-(*k)*x[0];
    *kSD=*bSD=*ySD=0.; *r=1.;
    return(0);
  }
 
  /* Calculate (x,y) sums and means */
  for(i=0; i<nr; i++) {
    sumx+=x[i]; sumy+=y[i];
    sumsx+=x[i]*x[i]; sumsy+=y[i]*y[i];
    sumxy+=x[i]*y[i];
  }
  meanx=sumx/(double)nr;
  meany=sumy/(double)nr;
  /* and then based on means */
  for(i=0; i<nr; i++) {
    f=x[i]-meanx; sumsdx+=f*f;
    g=y[i]-meany; sumsdy+=g*g;
    sumdxdy+=f*g;
  }
  if(sumsdx<1.0e-50 || sumsdy<1.0e-50) return(3);
  /* Regression coefficient */
  kk=sumdxdy/sumsdx; *k=kk;
  /* Intercept with y axis */
  bb=(sumsdx*sumy - sumx*sumdxdy)/((double)nr*sumsdx); *b=bb;
  /* Errors */
  for(i=0; i<nr; i++) {
    f=kk*x[i]+bb-y[i];
    sumsdcy+=f*f;
  }
  /* Deviation of y values */
  if(sumsdcy<=1.0e-12) e=0.0; else e=sqrt(sumsdcy/(double)(nr-2)); *ySD=e;
  /* SD of slope and intercept */
  ke=e/sqrt(sumsdx); be=e/sqrt((double)nr-sumx*sumx/sumsx);
  *kSD=ke; *bSD=be;
  be=sqrt(sumsdcy/(double)(nr-2))/sqrt((double)nr-(sumx*sumx)/sumsx);
  /* Pearson's correlation coefficient */
  rr=(sumxy-((sumx*sumy)/(double)nr)) /
     sqrt((sumsx-sumx*sumx/(double)nr)*(sumsy-sumy*sumy/(double)nr));
  /* Correct for small sample size */
  if(nr>4) rr*=1.0+(1.0-rr*rr)/(double)(2*(nr-4));
  *r=rr;
  if(0) {
    fprintf(stdout, "k=%14.7e +- %14.7e\n", kk, ke);
    fprintf(stdout, "b=%14.7e +- %14.7e\n", bb, be);
    fprintf(stdout, "r=%14.7e ySD=%14.7e\n", rr, e);
    fflush(stdout);
  }
 
  return(0);
}

Referenced by best_logan_regr(), best_pearson(), dft_end_line(), extrapolate_monoexp(), pearson2(), pearson3(), and pearson4().

pearson2()

int pearson2	(	double *	x,
		double *	y,
		char *	is,
		int	nr,
		double *	k,
		double *	kSD,
		double *	b,
		double *	bSD,
		double *	r,
		double *	ySD )

Calculate slope and intercept of a line and Pearson's correlation coefficient.

Array char is[] specifies whether single (x,y) points are used in the fit.

See also

regr_line, best_pearson, pearson, highest_slope

Returns

If successful, function returns value 0.

Parameters

x	data x values
y	data y values
is	Switch values: 0=do not use this point
nr	number of data sample values
k	slope
kSD	S.D. of slope
b	y axis intercept
bSD	S.D. of y axis intercept
r	Pearson's correlation coefficient, r
ySD	Residual variance of y values

Definition at line 109 of file pearson.c.

  {
  int i, j;
  double *nx, *ny;
 
  /* Allocate memory for new data pointers  */
  nx=(double*)calloc(nr, sizeof(double)); if(nx==NULL) return 1;
  ny=(double*)calloc(nr, sizeof(double)); if(ny==NULL) {free((char*)nx); return 1;}
 
  /* Copy data to pointers */
  for(i=0, j=0; i<nr; i++) if(is[i]) {nx[j]=x[i]; ny[j]=y[i]; j++;}
 
  /* Use pearson() */
  i=pearson(nx, ny, j, k, kSD, b, bSD, r, ySD); free((char*)nx); free((char*)ny);
  return (i);
}

pearson3()

int pearson3	(	double *	x,
		double *	y,
		int	nr,
		double *	k,
		double *	kSD,
		double *	b,
		double *	bSD,
		double *	r,
		double *	ySD )

Calculate slope and intercept of a line and Pearson's correlation coefficient.

Data points may contain NaN's.

See also

regr_line, pearson, best_pearson, highest_slope

Returns

If successful, function returns value 0.

Parameters

x	data x values
y	data y values
nr	number of data sample values
k	slope
kSD	S.D. of slope
b	y axis intercept
bSD	S.D. of y axis intercept
r	Pearson's correlation coefficient, r
ySD	Residual variance of y values

Definition at line 154 of file pearson.c.

  {
  int i, j;
  double *nx, *ny;
 
  /* Allocate memory for new data pointers  */
  nx=(double*)calloc(nr, sizeof(double)); if(nx==NULL) return 1;
  ny=(double*)calloc(nr, sizeof(double)); if(ny==NULL) {free((char*)nx); return 1;}
 
  /* Copy data to pointers  */
  for(i=0, j=0; i<nr; i++)
    if(!isnan(x[i]) && !isnan(y[i])) {nx[j]=x[i]; ny[j]=y[i]; j++;}
 
  /* Use pearson()  */
  i=pearson(nx, ny, j, k, kSD, b, bSD, r, ySD); free((char*)nx); free((char*)ny);
  return (i);
}

pearson4()

int pearson4	(	double *	x,
		double *	y,
		int	nr,
		double	start,
		double	end,
		double *	k,
		double *	kSD,
		double *	b,
		double *	bSD,
		double *	r,
		double *	ySD )

Calculate slope and intercept of a line and Pearson's correlation coefficient.

Data points may contain NaN's. Fit start and end times are specified.

See also

regr_line, best_pearson, highest_slope

Returns

If successful, function returns value 0.

Parameters

x	data x values
y	data y values
nr	number of data sample values
start	fit start time
end	fit end time
k	slope
kSD	S.D. of slope
b	y axis intercept
bSD	S.D. of y axis intercept
r	Pearson's correlation coefficient, r
ySD	Residual variance of y values

Definition at line 198 of file pearson.c.

  {
  int i, j;
  double *nx, *ny;
 
  if(0) {
    fprintf(stdout, "pearson4(x[], y[], %d, %g, %g, *k, *kSD, *b, *bSD, *r, *ySD)\n", nr,start,end);
    fflush(stdout);
  }
  /* Allocate memory for new data pointers */
  if((nx=(double*)calloc(nr, sizeof(double)))==NULL) return -3;
  if((ny=(double*)calloc(nr, sizeof(double)))==NULL) {free((char*)nx); return -3;}
 
  /* Copy data to pointers */
  for(i=0, j=0; i<nr; i++)
    if(x[i]>=start && x[i]<=end && !isnan(x[i]) && !isnan(y[i])) {
      nx[j]=x[i]; ny[j]=y[i]; j++;}
 
  /* Use pearson() */
  i=pearson(nx, ny, j, k, kSD, b, bSD, r, ySD);
  free((char*)nx); free((char*)ny);
  return i;
}

regr_line()

int regr_line	(	double *	x,
		double *	y,
		int	n,
		double *	m,
		double *	c )

Calculates regression line slope (m) and y axis intercept.

Data (x and y data) may contain NaN's.

See also

highest_slope, pearson, imgsegmPearson

Returns

Returns 0 if ok.

Parameters

x	An array of x axis values.
y	An array of y axis values.
n	The number of values in x and y arrays.
m	Pointer where calculated slope is written.
c	Pointer where calculated y axis intercept is written.

Definition at line 387 of file pearson.c.

  {
  double xsum=0.0, ysum=0.0, x2sum=0.0, xysum=0.0, delta;
  int i, nn=0;
 
  /* Check the data */
  if(x==NULL || y==NULL) return(1);
  if(n<2) return(2);
 
  /* Compute */
  for(i=0, nn=0; i<n; i++) {
    if(isnan(x[i]) || isnan(y[i])) continue;
    xsum+=x[i]; ysum+=y[i]; x2sum+=x[i]*x[i]; xysum+=x[i]*y[i]; nn++;
  }
  if(nn<2) return(2);
  delta=(double)nn*x2sum - xsum*xsum; if(delta==0.0) return(3);
  *m=((double)nn*xysum - xsum*ysum)/delta;
  *c=(x2sum*ysum - xsum*xysum)/delta;
  return(0);
}

Referenced by highest_slope(), and highest_slope_after().