<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Mahalanobis Distance
This notebook shows how I think we should do Mahalanobis distance for the SECIM project. From JMP Website
Step1:  Generation of a simulated data set
As an example I will simulate a data set
Step2:  Calculate Mahalanobis distance
There is a Mahalanobis function in scipy, but from what I can tell it only works on 1d arrays and not matrices. I think it will be easiest to o the calculation by hand.
<ASSISTANT_TASK:>
Python Code:
import pandas as pd
import numpy as np
import scipy as sp
import scipy.stats as stats
import matplotlib.pyplot as plt
import cPickle as pickle
import os
%matplotlib inline
Explanation: Mahalanobis Distance
This notebook shows how I think we should do Mahalanobis distance for the SECIM project. From JMP Website:
The Mahalanobis distance takes into account the correlation structure of the data and the individual scales. For each value, the Mahalanobis distance is denoted $M_i$ and is computed as
$M_i = \sqrt{(Y_i - \bar Y)^T S^{-1}(Y_i - \bar Y)}$
where:
* $Y_i$ is the data for the ith row
* Y is the row of means
* S is the estimated covariance matrix for the data
End of explanation
covSim = np.array([[1.0, .8, .2, .2],
                   [.8, 1.0, .3, .3],
                   [.3, .3, 1.0, .8],
                   [.2, .2, .8, 1.0]])
np.random.seed(111)
datSim = np.random.multivariate_normal([2, 3, 8, 9], covSim, size=1000)
dfSim = pd.DataFrame(data=datSim, columns=['sample1', 'sample2', 'sample3', 'sample4'])
# Save for comparing in sas
dfSim.to_csv('/home/jfear/tmp/dfSim.csv', index=False)
 
dfSim.head()
Explanation: Generation of a simulated data set
As an example I will simulate a data set
End of explanation
# Calculate the covaranice matrix from the data
covHat = dfSim.cov()
covHat
# Get the inverse of the covarance matrix
covHatInv = np.linalg.inv(covHat)
covHatInv
# Calculate the column means
colMean = dfSim.mean(axis=0)
colMean
# Subtract the mean from each value
dfSimCenter = (dfSim - colMean).T
dfSimCenter.head()
# Calculate the mahalanobis distance
MD = np.sqrt(np.dot(np.dot(dfSimCenter.T, covHatInv), dfSimCenter))
MDval = np.diag(MD)
plt.scatter(x=range(len(MDval)), y=MDval)
plt.axhline(np.percentile(MDval, 95), ls='--', lw=2)
plt.show()
Explanation: Calculate Mahalanobis distance
There is a Mahalanobis function in scipy, but from what I can tell it only works on 1d arrays and not matrices. I think it will be easiest to o the calculation by hand.
End of explanation
<END_TASK>