The math
The correlation coefficient between an array of boolean vectors is a value between 0 and 1 defined by the following equations, covariance and standard deviation:
$\operatorname{cov} (X,Y)=\frac{1}{n}\sum_{i=1}^n (x_i-\overline{X})(y_i-\overline{Y})$
$\sigma_X=\sqrt{\frac{1}{n}\sum_{i=1}^n (x_i-\overline{X})^2}$
$\sigma_Y=\sqrt{\frac{1}{n}\sum_{i=1}^n (y_i-\overline{Y})^2}$
Where $\overline{X}$ is the average value of a set, and $X$ is the set of $x$ values, which would be the first array in your variables [0,1,0,1] , and $Y$ is the set of $y$ values, the second array [0,0,1,1]. Note the array indexing is 1-4 for the equations, but will not be so in computer code.
The correlation coefficient is the final equation:
$\rho_{X,Y}= \left|\frac{\operatorname{cov}(X,Y)}{\sigma_X \sigma_Y}\right|$
$\rho_{X,Y}= \left|\frac{\sum_{i=1}^n (x_i-\overline{X})(y_i-\overline{Y})}{\sqrt{\sum_{i=1}^n (x_i-\overline{X})^2} \ \cdot \ \sqrt{\sum_{i=1}^n (y_i-\overline{Y})^2}}\right|$
In the scenario where $\sigma_X$ and $\sigma_Y$ are both 0, the vectors are most likely completely correlated, and a different comparison must be made. If one of them is 0, the results are probably not correlated. You will probably not run into these situations if you are analyzing an s-box for bit independence.
Your example
Because your dataset is small, it is easy to use as an example, but is still a bad example, because the correlation coefficient is 0:
The standard dev of your X and Y sets is both 0.5, but the covariance is 0:
$0 = \left|\frac{0}{0.5 * 0.5}\right|$
Now, if you change your X,Y to [0,1,0,1],[1,0,1,0], you can see that they are completely correlated [y = not(x)], and the math agrees, as the covariance calculation is -0.25:
$1 = \left|\frac{-0.25}{0.5 * 0.5}\right|$
Optimization and code
You may notice that standard deviation looks like the square root of covariance between the same set, and for this calculation it is. You may also notice you cancel out the 1/N in the final equation. The intermediate calculations of average values and subtractions are also the same, so you can combine them all into a single calculation. Careful optimization will also help compensate for any floating point rounding errors.
Optimized pseudocode for the final equation on arrays of X and Y would look something like this:
mean_X = ArithMean(X)
mean_Y = ArithMean(Y)
For i = 1 to N
val_x = X(i) - mean_X
val_y = Y(i) - mean_Y
deltasum_c = deltasum_c + (val_x * val_y)
deltasum_x = deltasum_x + (val_x * val_x)
deltasum_y = deltasum_y + (val_y * val_y)
Next i
CoVariance = deltasum_c
StdDeviation_xy = Sqrt(deltasum_x * deltasum_y)
If StdDeviation_xy <> 0 Then COR = Abs(CoVariance / StdDeviation_xy)
Elseif deltasum_x = deltasum_y Then COR = 1
Else COR = 0