Correlation is the process of quantifying the relationship between two sets of values, and in this post I will be writing code in C to calculate possibly the best-known type of correlation - the Pearson Correlation Coefficient.

## An Overview of Correlation

Consider the following three data sets and their graphs, or, more accurately, scatter plots. The x-axis represents an independent variable, and the y-axis represents a dependent variable.

## Data Set 1

x | y |
---|---|

10 | 32 |

20 | 44 |

40 | 68 |

45 | 74 |

60 | 92 |

65 | 98 |

75 | 110 |

80 | 116 |

In this data set all the points are on a straight line - there is clearly a relationship between them.

## Data Set 2

x | y |
---|---|

10 | 40 |

20 | 40 |

40 | 60 |

45 | 80 |

60 | 90 |

65 | 110 |

75 | 100 |

80 | 130 |

Here the data is a bit more scattered but there is still a clear relationship, subject to some minor fluctuations.

## Data Set 3

x | y |
---|---|

10 | 100 |

20 | 10 |

40 | 130 |

45 | 90 |

60 | 40 |

65 | 80 |

75 | 180 |

80 | 50 |

In this last set of data the points appear completely random, with little or no relationship between x and corresponding y values.

## Quantifying the Relationships

Relationships within bivariate data such as the examples above can be precisely quantified by calculating the Pearson Correlation Coefficient, usually referred to as ρ, the Greek letter rho. This gives us a number between -1 and 1, -1 being a perfect negative correlation and 1 being a perfect positive correlation. Values in between represent various degrees of imperfect correlation, with 0 being none at all.

As rough estimates, in **Data Set 1** we would expect ρ to be 1, in **Data Set 2** we would expect ρ to be around 0.9, and in **Data Set 3** ρ would be close to 0.

Full details of the Pearson Correlation Coefficient are on Wikipedia. In short the formula is

Formula for Pearson Correlation Coefficient

ρX,Y = cov(X,Y) / (σX * σY)

where cov is the covariance of the data paired data, and σ is the standard deviation of the data. These will be explored as we get to them.

## Caveat

There is a famous phrase in statistics: "correlation does not imply causation". A strong correlation between two variables might just be by chance, even where common sense might make you believe the association could be genuine. A strong correlation should only ever be used as a basis for further investigation, particularly in this age of Big Data number crunching where seemingly unrelated data sets can be compared quickly and cheaply.

## Coding

For this project we will write a function to create sample arrays of data corresponding to that shown above, and then write the code necessary to calculate ρ for each of them.

Create a new folder and then create the following empty files:

- main.c
- data.h
- data.c
- pearsoncorrelation.h
- pearsoncorrelation.c

Open main.c and type or paste the following. You can also download the source code from the Downloads page if you prefer.

main.c

#include<stdio.h>

#include<stdlib.h>

#include<stdbool.h>

#include"data.h"

#include"pearsoncorrelation.h"

//--------------------------------------------------------

// FUNCTION PROTOTYPES

//--------------------------------------------------------

void printdata(double* independent, double* dependent, int size);

//--------------------------------------------------------

// FUNCTION main

//--------------------------------------------------------

int main(void)

{

puts("---------------------------------------\n| code-in-c.com - Pearson Correlation |\n---------------------------------------\n");

double independent[8];

double dependent[8];

double rho;

for(int dataset = 1; dataset <=3; dataset++)

{

if(populatedata(independent, dependent, dataset))

{

printf("Data Set %d\n-----------\n", dataset);

printdata(independent, dependent, 8);

rho = pearson_correlation(independent, dependent, 8);

printf("Pearson Correlation Coefficient rho = %lf\n\n", rho);

}

}

return EXIT_SUCCESS;

}

//--------------------------------------------------------

// FUNCTION printdata

//--------------------------------------------------------

void printdata(double* independent, double* dependent, int size)

{

for(int i = 0; i < size; i++)

{

printf("%3.0lf\t%3.0lf\n", independent[i], dependent[i]);

}

}

In main we declare a pair of double arrays to hold our bivariate data. The size of these is hard-coded in as this is only testing and demonstration code. The actual function which calculates the correlation will be general-purpose, and will take data of any size. We also declare the variable rho to take the correlation coefficient.

The populatedata function can populate arrays with any of the three sample data sets above, so we for-loop from 1 to 3 to calculate the correlation for each in turn. The data set number and the data itself is printed by calling the very trivial printdata function, after which we call the not-so-trivial pearson_correlation function and print the result.

We will now write the code to populate the arrays. Open data.h and paste in the following.

data.h

//--------------------------------------------------------

// FUNCTION PROTOTYPES

//--------------------------------------------------------

bool populatedata(double independent[8], double dependent[8], int dataset);

Now open data.c and enter this.

data.c

#include<string.h>

#include<stdio.h>

#include<stdbool.h>

//--------------------------------------------------------

// FUNCTION populatedata

//--------------------------------------------------------

bool populatedata(double independent[8], double dependent[8], int dataset)

{

bool success = true;

switch(dataset)

{

case 1:

memcpy(independent, (double[8]){10,20,40,45,60,65,75,80}, sizeof(double[8]));

memcpy(dependent, (double[8]){32,44,68,74,92,98,110,116}, sizeof(double[8]));

break;

case 2:

memcpy(independent, (double[8]){10,20,40,45,60,65,75,80}, sizeof(double[8]));

memcpy(dependent, (double[8]){40,40,60,80,90,110,100,130}, sizeof(double[8]));

break;

case 3:

memcpy(independent, (double[8]){10,20,40,45,60,65,75,80}, sizeof(double[8]));

memcpy(dependent, (double[8]){100,10,130,90,40,80,180,50}, sizeof(double[8]));

break;

default:

puts("Invalid dataset requested");

success = false;

break;

}

return success;

}

Basically all we are doing here is copying hard-coded data into the arrays supplied as arguments and returning true. If the dataset index is not 1, 2 or 3 we return false. Note that we checked this return value in main.

Now the moment you've all been waiting for - we will now actually code up the pearson_correlation function. First paste the prototype into pearsoncorrelation.h.

pearsoncorrelation.h

//--------------------------------------------------------

// FUNCTION PROTOTYPES

//--------------------------------------------------------

double pearson_correlation(double* independent, double* dependent, int size);

Then open pearsoncorrelation.c and paste all the following code into it.

pearsoncorrelation.c

#include<stdlib.h>

#include<stdbool.h>

#include<stdio.h>

#include<string.h>

#include<math.h>

#include"pearsoncorrelation.h"

//--------------------------------------------------------

// FUNCTION PROTOTYPES

//--------------------------------------------------------

static double arithmetic_mean(double* data, int size);

static double mean_of_products(double* data1, double* data2, int size);

static double standard_deviation(double* data, int size);

//--------------------------------------------------------

// FUNCTION pearson_correlation

//--------------------------------------------------------

double pearson_correlation(double* independent, double* dependent, int size)

{

double rho;

// covariance

double independent_mean = arithmetic_mean(independent, size);

double dependent_mean = arithmetic_mean(dependent, size);

double products_mean = mean_of_products(independent, dependent, size);

double covariance = products_mean - (independent_mean * dependent_mean);

// standard deviations of independent values

double independent_standard_deviation = standard_deviation(independent, size);

// standard deviations of dependent values

double dependent_standard_deviation = standard_deviation(dependent, size);

// Pearson Correlation Coefficient

rho = covariance / (independent_standard_deviation * dependent_standard_deviation);

return rho;

}

//--------------------------------------------------------

// FUNCTION arithmetic_mean

//--------------------------------------------------------

static double arithmetic_mean(double* data, int size)

{

double total = 0;

// note that incrementing total is done within the for loop

for(int i = 0; i < size; total += data[i], i++);

return total / size;

}

//--------------------------------------------------------

// FUNCTION mean_of_products

//--------------------------------------------------------

static double mean_of_products(double* data1, double* data2, int size)

{

double total = 0;

// note that incrementing total is done within the for loop

for(int i = 0; i < size; total += (data1[i] * data2[i]), i++);

return total / size;

}

//--------------------------------------------------------

// FUNCTION standard_deviation

//--------------------------------------------------------

static double standard_deviation(double* data, int size)

{

double squares[size];

for(int i = 0; i < size; i++)

{

squares[i] = pow(data[i], 2);

}

double mean_of_squares = arithmetic_mean(squares, size);

double mean = arithmetic_mean(data, size);

double square_of_mean = pow(mean, 2);

double variance = mean_of_squares - square_of_mean;

double std_dev = sqrt(variance);

return std_dev;

}

As you can see much of the arithmetic is farmed out to separate functions:

- arithmetic_mean
- mean_of_products
- standard_deviation

Just to recap, the formula for the Pearson Correlation Coefficient is

Formula for Pearson Correlation Coefficient

ρX,Y = cov(X,Y) / (σX * σY)

This is implemented near the end of the pearson_correlation function, most of which is concerned with calculating the three components of the formula:

- The covariance of the bivariate data
- The standard deviation of the independent data
- The standard deviation of the dependent data

The code is hopefully self-documenting enough not to need an explanation of the actual calculation but I will add a couple of brief notes about covariance and standard deviation.

Covariance is a rather complex concept which you might like to look into on Wikipedia or elsewhere, but an easy way to remember how to calculate it is **mean of products minus product of means**. The **mean of products** is the mean of each pair of values multiplied together, and the **product of means** is simply the means of the two sets of data multiplied together.

Standard deviation can be thought of roughly as the average amount by which the individual values vary from the mean. (It is not that exactly, which is a different figure known as mean of absolute deviation or MAD, but the concepts are similar.) The standard deviation is the square root of the variance which uses squared values to eliminate negative differences. An easy way to remember how to calculate a variance is **mean of squares minus square of mean**. The **mean of squares** is simple enough, just square all the values and calculate their mean; and the **square of mean** is self-explanatory. Once you have calculated the variance don't forget to calculate its square root to get the standard deviation.

## Compile and Run at Last

We have now finished writing code so can actually compile and run it. In the terminal type

Compiling

gcc main.c data.c pearsoncorrelation.c -std=c11 -lm -o main

Note all three .c files are included, and we also need -lm as we're using math.h.

Now run the program with

Running

./main

The output is:

Program Output

---------------------------------------

| code-in-c.com - Pearson Correlation |

---------------------------------------

Data Set 1

-----------

10 32

20 44

40 68

45 74

60 92

65 98

75 110

80 116

Pearson Correlation Coefficient rho = 1.000000

Data Set 2

-----------

10 40

20 40

40 60

45 80

60 90

65 110

75 100

80 130

Pearson Correlation Coefficient rho = 0.960215

Data Set 3

-----------

10 100

20 10

40 130

45 90

60 40

65 80

75 180

80 50

Pearson Correlation Coefficient rho = 0.207800

Earlier I gave rough estimates of ρ for the three data sets as 1, 0.9 and 0. As you can see they are actually 1, 0.960215 and 0.2078.

## Next Step: Regression

If you have discovered a strong correlation in your data your next step will probably be to put numbers to that relationship by working out the formula of the straight line which best fits the data. That process is called linear regression and will be the subject of a future post.

Please let me know any ideas you might have for improvements to this code in Comments below, and follow Code in C to keep up with future posts.

Follow @code_in_c