Linear Discriminant Analysis biplot I am not sure about plotting a biplot for LDA. Suppose I have a 3 feature data set with 3 classes. Then, I perform LDA to reduce the dimensionality to 2. How can I create a biplot for the LDA?
My understanding of biplot (I may be completely wrong) is simply a projection of the original features in the hyperplane (or line) of the dimensionality reduction. So, for example, using Python scikit-learn, can I simply perform the following?
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
lda = LinearDiscriminantAnalysis(n_components=2)
X_lda = lda.fit_transform(X_std,y) #X_std is input data matrix X standardized by Standardscaler, y is a vector of target values
org_features = np.identity(3)
proj_features = lda.transform(org_features)

Then, just connect 0,0 to the 3 points in proj_features to get the 3 arrows corresponding to the 3 original features?

Update:
I am still looking for confirmation of the code. For now, I am just interpreting the three points obtained from the product of 3D identity matrix and the LDA projection matrix as the change in location in the 2D LDA plane as I change the value of each feature by 1 Standard Deviation (all features were standardized by StandardScaler) while keeping the rest of the features constant.
 A: I am going to provide an example but keep in mind that it is not 100% clear that the lda.scalings_ stores the eigenvectors (loadings of the variables). This is something that I am currently trying to verify.
Example using iris data and sklearn:
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis


iris = datasets.load_iris()
X = iris.data
y = iris.target
#In general a good idea is to scale the data
scaler = StandardScaler()
scaler.fit(X)
X=scaler.transform(X)

lda = LinearDiscriminantAnalysis()
lda.fit(X,y)

x_new = lda.transform(X)   

def myplot(score,coeff,labels=None):
    xs = score[:,0]
    ys = score[:,1]
    n = coeff.shape[0]
    scalex = 1.0/(xs.max() - xs.min())
    scaley = 1.0/(ys.max() - ys.min())
    #plt.scatter(xs * scalex,ys * scaley, c = y) #with scaling
    plt.scatter(xs ,ys, c = y) #without scaling
    for i in range(n):
        plt.arrow(0, 0, coeff[i,0], coeff[i,1],color = 'r',alpha = 0.5)
        if labels is None:
            plt.text(coeff[i,0]* 1.15, coeff[i,1] * 1.15, "Var"+str(i+1), color = 'g', ha = 'center', va = 'center')
        else:
            plt.text(coeff[i,0]* 1.15, coeff[i,1] * 1.15, labels[i], color = 'g', ha = 'center', va = 'center')
#plt.xlim(-1,1) # use these then you use:plt.scatter(xs * scalex,ys * scaley, c = y)
#plt.ylim(-1,1)
plt.xlabel("PC{}".format(1))
plt.ylabel("PC{}".format(2))
plt.grid()

#Call the function. 

# Important: here I think that lda.scalings_ contains the 2 eigenvectors (loadings of the variables). The shape is [n_features,n_components] so [4,2] in our case. So in the myplot function I plot for each variable i, the values that are in [i,0] and [i,1]. All these assuming that the lda.scalings_ contain the eigenvectors.

myplot(x_new[:,0:2], lda.scalings_) 
plt.show()

Results

