<jupyter_start><jupyter_text>normal PCA（linear projection）<jupyter_code>from sklearn.decomposition import PCA

# PCA assumes that the origin is centered, the data set is centralized automatically
pca1 = PCA(n_components=2)
# reduce the data set down to 2-dimensional space
X_reduced1 = pca1.fit_transform(X_processed)
X_reduced1.shape
# principal component's weight coefficient, 
# it can obtain the expression of principal component in feature space
pca1.components_
# The larger these absolute values are, the more a specific feature contributes to that principal component
abs(pca1.components_)
# the explained variance ratio of the data set along the principal component axis 
pca1.explained_variance_ratio_
# missing explained variance ratio
1 - pca1.explained_variance_ratio_.sum()
# the projection of a data set onto a hyperplane that is defined by the principal component，2272 x 2
X_reduced1
import matplotlib.pyplot as plt

plt.figure(1)

legend = []
# Draw a data point projection plot that is projected to a two-dimensional plane using normal PCA
for i, label in enumerate(tag):
    plt.scatter(X_reduced1[:, 0][labels==label], X_reduced1[:, 1][labels==label])
    legend.append("Label: {}".format(label))
plt.legend(legend)
plt.xlabel("$PC1$")
plt.ylabel("$PC2$")
plt.title("2-dimension PCA")
plt.grid()
save_fig("2D_PCA")
pca1.components_.shape
X_reduced1
# restore the original data set, that is reconstructed data
X_recovered = pca1.inverse_transform(X_reduced1)
X_recovered.shape
# Compare the original data set with the recovered data set after dimension reduction,
# False indicates that these are different, that is, recovered data set losed data
np.allclose(X_recovered, X_processed)
from sklearn.metrics import mean_squared_error

# mean square distance between reconstructed data and the original data
mean_squared_error(X_recovered, X_processed)
# don't set the number of dimensions to be reduced,
# that is to find all principal component axes which exist in the feature space
pca2 = PCA()
pca2.fit(X_processed)
# variance interpretation proportion on each principal component axis
pca2.explained_variance_ratio_
cumsum = np.cumsum(pca2.explained_variance_ratio_)
# calculate the minimum dimension required when the variance interpretation proportion is retained at 95 %
d = np.argmax(cumsum>=0.95) + 1
d
plt.figure(figsize=(6,4))
plt.plot(cumsum, linewidth=3)
plt.xlabel("Dimensions")
plt.ylabel("Explained Variance")
plt.grid(True)
save_fig("explained_variance_plot")
plt.show()
import matplotlib.pyplot as plt
import numpy as np

def biplot(reduced_data, labels, pc, variable):
    """plot componential biplot for two principle components
    
    :param reduced_data: data processed by PCA
    :param labels: labels of original dataset 
    :param pc: all the principle components
    :param variable: the name of the variables of the data set
    """
    plt.figure(1, figsize=(14, 10))
    
    legend = []  # 
    classes = np.unique(labels)  # label type
    n = pc.shape[1]
    # colors = ['g', 'r', 'y']
    # markers = ['o', '^', 'x']
    
    x = reduced_data[:, 0]  # variable contributions for PC1
    y = reduced_data[:, 1]  # variable contributions for PC2
    scalex = 1.0/(x.max() - x.min())
    scaley = 1.0/(y.max() - y.min())
    
    # Draw a data point projection plot that is projected to 
    # a two-dimensional plane using normal PCA
    for i, label in enumerate(classes):
        plt.scatter(x[labels==label] * scalex,
                    y[labels==label] * scaley,
                    linewidth=0.01) 
        # hyperparameter in plt.scatter(): c=colors[i], marker=markers[i]
        legend.append("Label: {}".format(label))
     
    plt.legend(legend)
    
    # plot arrows as the variable contribution,
    # each variable has a score for PC1 and for PC2 respectively
    for i in range(n):
        plt.arrow(0, 0, pc[0, i], pc[1, i], color='k', alpha=0.7,
                  linewidth=1, )
        plt.text(pc[0, i]*1.01, pc[1, i]*1.01, variable[i],
                 ha='center', va='center', color='k', fontsize=12)
        
    plt.xlabel("$PC1$")
    plt.ylabel("$PC2$")
    plt.title("Componential biplot")
    plt.grid()
    save_fig("Componential biplot")
biplot(X_reduced1, labels, pca1.components_, variable)
from sklearn.decomposition import PCA

# PCA assumes that the origin is centered, the data set is centralized automatically
pca3 = PCA(n_components=3)
# reduce the data set down to 3-dimensional space
X_reduced3 = pca3.fit_transform(X_processed)
X_reduced3.shape
import matplotlib.pyplot as plt
import numpy as np

def componential_plot_static_3d(reduced_data, labels, pc, variable):
    """draw a static componential plot in 3d for three principle components
    
    :param reduced_data: data processed by PCA
    :param labels: labels of original dataset 
    :param pc: all the principle components
    :param variable: the name of the variables of the data set
    """
    
    fig = plt.figure(1, figsize=(14, 12))
    ax = plt.axes(projection='3d')
    
    legend = []  # 
    classes = np.unique(labels)  # label type
    n = pc.shape[1]
    # colors = ['g', 'r', 'y']
    # markers = ['o', '^', 'x']
    
    x = reduced_data[:, 0]  # variable contributions for PC1
    y = reduced_data[:, 1]  # variable contributions for PC2 
    z = reduced_data[:, 2]  # variable contributions for PC3
    scalex = 1.0/(x.max() - x.min())
    scaley = 1.0/(y.max() - y.min())
    scalez = 1.0/(z.max() - z.min())
    
    # Draw a data point projection plot that is projected to 
    # a three-dimensional space using normal PCA
    for i, label in enumerate(classes):
        ax.scatter3D(x[labels==label] * scalex,
                     y[labels==label] * scaley,
                     z[labels==label] * scalez,
                     linewidth=0.01) 
        # hyperparameter in plt.scatter(): c=colors[i], marker=markers[i]
        legend.append("Label: {}".format(label))
    
    ax.legend(legend)
    
    # the initial angle to draw the 3d plot
    azim = -60  # azimuth
    elev = 30   # elevation
    ax.view_init(elev, azim)  # set the angles
    
    # plot arrows as the variable contribution,
    # each variable has a score for PC1, for PC2 and for PC3 respectively
    for i in range(n):
        ax.quiver(0, 0, 0, pc[0, i], pc[1, i], pc[2, i], color='k', alpha=0.7,
                  linewidth=1, arrow_length_ratio=0.05)
        ax.text(pc[0, i]*1.1, pc[1, i]*1.1, pc[2, i]*1.1, variable[i],
                ha='center', va='center', color='k', fontsize=12)
        
    ax.set_xlabel("$PC1$")
    ax.set_ylabel("$PC2$")
    ax.set_zlabel("$PC3$")
    plt.title("Componential Plot in 3 Dimension")
    plt.grid()
    save_fig("Componential_Plot_in_3_Dimension")
componential_plot_static_3d(X_reduced3, labels, pca3.components_, variable)
import matplotlib.pyplot as plt
from matplotlib import animation 
import numpy as np

def componential_plot_animated_3d(reduced_data, labels, pc, variable):
    """draw a animated componential plot in 3d for three principle components
    
    :param reduced_data: data processed by PCA
    :param labels: labels of original dataset 
    :param pc: all the principle components
    :param variable: the name of the variables of the data set
    """
    
    fig = plt.figure(1, figsize=(14, 12))
    ax = plt.axes(projection='3d')
    
    legend = []  # 
    classes = np.unique(labels)  # label type
    n = pc.shape[1]
    # colors = ['g', 'r', 'y']
    # markers = ['o', '^', 'x']
    
    x = reduced_data[:, 0]  # variable contributions for PC1
    y = reduced_data[:, 1]  # variable contributions for PC2 
    z = reduced_data[:, 2]  # variable contributions for PC3
    scalex = 1.0/(x.max() - x.min())
    scaley = 1.0/(y.max() - y.min())
    scalez = 1.0/(z.max() - z.min())
    
    def init():
        #  A function used to draw a clear and initial frame
        
        # draw a data point projection plot that is projected to 
        # a three-dimensional space using normal PCA
        for i, label in enumerate(classes):
            ax.scatter3D(x[labels==label] * scalex,
                         y[labels==label] * scaley,
                         z[labels==label] * scalez,
                         linewidth=0.01) 
            # hyperparameter in plt.scatter(): c=colors[i], marker=markers[i]
            legend.append("Label: {}".format(label))

        ax.legend(legend)

        # plot arrows as the variable contribution,
        # each variable has a score for PC1, for PC2 and for PC3 respectively
        for i in range(n):
            ax.quiver(0, 0, 0, pc[0, i], pc[1, i], pc[2, i], color='k', alpha=0.7,
                      linewidth=1, arrow_length_ratio=0.05)
            ax.text(pc[0, i]*1.1, pc[1, i]*1.1, pc[2, i]*1.1, variable[i],
                    ha='center', va='center', color='k', fontsize=12)

        ax.set_xlabel("$PC1$")
        ax.set_ylabel("$PC2$")
        ax.set_zlabel("$PC3$")
        plt.title("Compositional Plot in 3 Dimension")
        plt.grid()
    
    def rotate(angle): 
        ax.view_init(elev=10, azim=angle) 
    
    # animate, frames=values will be passed to rotate, interval means the the delay between frames in milliseconds
    rot_animation = animation.FuncAnimation(fig, rotate, init_func=init,
                                            frames=np.arange(0,362,2), interval=50) 
    # save the animat
    rot_animation.save('./images/animated_3d_composiontion_plot.gif', dpi=80, writer='imagemagick') 
componential_plot_animated_3d(X_reduced3, labels, pca3.components_, variable)<jupyter_output><empty_output><jupyter_text>Randomized PCA<jupyter_code># quickly finds an approximation of the first d principal components
# by default, svd_solver is actually set to "auto"
rnd_pca = PCA(n_components=2, svd_solver="randomized")
X_reduced2 = rnd_pca.fit_transform(X_processed)<jupyter_output><empty_output><jupyter_text>Incremental PCA<jupyter_code>from sklearn.decomposition import IncrementalPCA

# for large training sets, the data can be splited into mini-batches and fed to IPCA 
# one mini-batch at a time
n_batches = 100
inc_pca = IncrementalPCA(n_components=2)
for X_batch in np.array_split(X_processed, n_batches):
    inc_pca.partial_fit(X_batch)
    
X_reduced4 = inc_pca.transform(X_processed)
# compare the mean value of the training set used for PCA and INC_PCA in each dimension and see if they are equal
np.allclose(pca1.mean_, inc_pca.mean_)
# compare the dimensionless data of PCA and INC_PCA and see if they are same
np.allclose(X_reduced1, X_reduced4)
# the time required to test the dimension reduction of incremental PCA
%timeit IncrementalPCA(n_components=5, batch_size=100).fit(X_processed)
# the time required to test the dimension reduction of ordinary PCA
%timeit PCA(n_components=5).fit(X_processed)
# the time required to test the dimension reduction of random PCA
%timeit PCA(n_components=5, svd_solver="randomized", random_state=42)<jupyter_output>2.52 µs ± 111 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)<jupyter_text>PCA kernel（non-linear projection）<jupyter_code>from sklearn.decomposition import KernelPCA

# call the PCA with RBF kernel for nonlinear projection
rbf_pca = KernelPCA(n_components=2, kernel="rbf", gamma=0.04)
X_reduced5 = rbf_pca.fit_transform(X_processed)
# reduced to two dimensional data set
X_reduced5
import matplotlib.pyplot as plt

lin_pca = KernelPCA(n_components=2, kernel="linear", fit_inverse_transform=True)
rbf_pca = KernelPCA(n_components=2, kernel="rbf", gamma=0.0433, fit_inverse_transform=True)
poly_pca = KernelPCA(n_components=2, kernel="poly", gamma=0.001, coef0=1, fit_inverse_transform=True)
sig_pca = KernelPCA(n_components=2, kernel="sigmoid", gamma=0.001, coef0=1, fit_inverse_transform=True)

model = [(221, lin_pca, "Linear kernel"), 
         (222, poly_pca, "Polynomial kernel, $\gamma=10^{-3}, r=1$"),
         (223, sig_pca, "Sigmoid kernel, $\gamma=10^{-3}, r=1$"), 
         (224, rbf_pca, "RBF kernel, $\gamma=0.0433$")]

def plot_2D_reduced(model, X, labels, tag, fig_name=True):
    """
    2D mapping of data points using dimensionality reduction method
    
    :param model: dimension reduction algorithm
    :param X: data sets only contain eigenvalues
    :param labels: data sets only contain label values
    :param tag: true category tag type
    :param fig_name: figure name
    """
    
    plt.figure(figsize=(11, 8))
    legend = []
    for subplot, pca, title in model:
        X_reduced = pca.fit_transform(X)
        plt.subplot(subplot)
        plt.title(title, fontsize=14)
        for i, label in enumerate(tag):
            plt.scatter(X_reduced[:, 0][labels==label], X_reduced[:, 1][labels==label])
            legend.append("Label: {}".format(label))
        if subplot == 221 or subplot == 223:
            plt.ylabel("$x_2$", fontsize=18, rotation=0)
        if subplot == 223 or subplot == 224:
            plt.xlabel("$x_1$", fontsize=18)
        plt.grid(True)
        plt.legend(legend)
    save_fig(fig_name)
    plt.show()

plot_2D_reduced(model, X_processed, labels, tag, fig_name="kernel_pca_plot")
from sklearn.model_selection import GridSearchCV
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import Pipeline

# pipeline first reduct dimension then classificate（can be changed to regression model）
clf = Pipeline([
        ("kpca", KernelPCA(n_components=2)),
        ("log_reg", LogisticRegression(solver="liblinear", multi_class="auto")),
    ])

# parameters of PCA kernels
param_grid = {
        "kpca__gamma": np.linspace(0.03, 0.05, 10),
        "kpca__kernel": ["rbf", "sigmoid", "linear", "poly"],
        "kpca__coef0": np.linspace(0, 1, 10)
    }

# look for optimal parameters automatically, and the score is the accuracy of the classification
grid_search = GridSearchCV(clf, param_grid, cv=3, scoring="accuracy", n_jobs=-1, verbose=1)
grid_search.fit(X_processed, labels)
# the optimal model
grid_search.best_estimator_
# parameters of the optimal model
grid_search.best_params_
# optimal model score
grid_search.best_score_
# dimensionality reduction using the optimal model
lin_pca_new = KernelPCA(n_components=2, kernel="linear", gamma=0.03, coef0=0, fit_inverse_transform=True)
X_reduced_new = lin_pca_new.fit_transform(X_processed)
# reconstruct pre-image data，only if fit_inverse_transform=True，inverse_transform method can be used
X_preimage = lin_pca_new.inverse_transform(X_reduced_new)
X_preimage.shape
from sklearn.metrics import mean_squared_error

# the square distance between the image data before reconstruction and the original data
mean_squared_error(X_processed, X_preimage)<jupyter_output><empty_output><jupyter_text>manifold learning（nonlinear dimension reduction）<jupyter_code>from sklearn.manifold import LocallyLinearEmbedding, Isomap, TSNE, MDS

# four common manifold learning algorithms
lle = LocallyLinearEmbedding(n_components=2, n_neighbors=10, random_state=42)
mds = MDS(n_components=2, random_state=42)
isomap = Isomap(n_components=2)
tsne = TSNE(n_components=2, random_state=42)

model_manifold = [(221, lle, "LocallyLinearEmbedding, n_neighbors=10"), 
                  (222, mds, "Multidimensional Scaling"),
                  (223, isomap, "Isomap"), 
                  (224, tsne, "t-Distributed Stochastic Neighbor Embedding")]

plot_2D_reduced(model_manifold, X_processed, labels, tag, fig_name="manifold_plot")<jupyter_output>Saving figure manifold_plot
{'fraction_non_alphanumeric': 0.1101486748545572, 'fraction_numerical': 0.024886877828054297, 'mean_word_length': 3.654332129963899, 'pattern_counts': {'":': 3, '<': 19, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '403568', 'license': 'MIT', 'n_tokens_mistral': 5309, 'n_tokens_neox': 4946, 'n_words': 1441, 'path': 'english/dimensionality_reduction/dimensionality_reduction.ipynb', 'repo_name': 'ZJUEarthData/ML_explore', 'size': '15470'}