{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn.feature_selection import VarianceThreshold\n",
    "from sklearn.datasets import load_iris\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 属性降维"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCA降维"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(150, 4)\n",
      "[[5.1 3.5 1.4 0.2]\n",
      " [4.9 3.  1.4 0.2]\n",
      " [4.7 3.2 1.3 0.2]\n",
      " [4.6 3.1 1.5 0.2]\n",
      " [5.  3.6 1.4 0.2]]\n",
      "(150, 2)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "from sklearn.datasets import load_iris\n",
    "data = load_iris()\n",
    "X = data.data\n",
    "y = data.target\n",
    "print(X.shape)\n",
    "print(X[:5])\n",
    "# 初始化PCA，设置主成分数目\n",
    "pca = PCA(n_components=2)\n",
    "\n",
    "# 对数据进行PCA降维\n",
    "X_pca = pca.fit_transform(X)\n",
    "\n",
    "# 输出降维后的数据形状\n",
    "print(X_pca.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ICA独立成分分析"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(150, 4)\n",
      "(150, 2)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.decomposition import FastICA\n",
    "print(X.shape)\n",
    "# 初始化ICA，设置主成分数目\n",
    "ica = FastICA(n_components=2)\n",
    "ica.fit(X)\n",
    "X_ica = ica.transform(X)\n",
    "print(X_ica.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "LDA线性判别分析"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(150, 4)\n",
      "(150, 2)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n",
    "print(X.shape)\n",
    "lda = LinearDiscriminantAnalysis(n_components=2)\n",
    "lda.fit(X, y)\n",
    "X_kpca = lda.transform(X)\n",
    "print(X_kpca.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 特征选择\n",
    "### 过滤法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "方差过滤"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "选取的数据已经被清洗过，所以方差过滤的效果不会太好"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(150, 4)\n",
      "(150, 4)\n"
     ]
    }
   ],
   "source": [
    "data = load_iris()\n",
    "X = data.data\n",
    "y = data.target\n",
    "print(X.shape)\n",
    "selector = VarianceThreshold(1*10**-40)\n",
    "x_feature_selection = selector.fit_transform(X)\n",
    "print(x_feature_selection.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "相关性过滤 \\\n",
    "卡方过滤"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(150, 4)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=300 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier as RFC\n",
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn.feature_selection import SelectKBest\n",
    "from sklearn.feature_selection import chi2\n",
    "\n",
    "# 留下300个特征\n",
    "X_fschi = SelectKBest(chi2, k=300).fit_transform(X, y)\n",
    "print(X_fschi.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用学习曲线来获得一个最优的超参数K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=390 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=380 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=370 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=360 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=350 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=340 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=330 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=320 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=310 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=300 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=290 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=280 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=270 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=260 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=250 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=240 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=230 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=220 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n",
      "d:\\Python\\CondaEnvs\\DataAnalysis\\lib\\site-packages\\sklearn\\feature_selection\\_univariate_selection.py:779: UserWarning: k=210 is greater than n_features=4. All the features will be returned.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhVklEQVR4nO3df1BVdeL/8dcFL3C3gDAUf4DgtpPRkq6ig0Bmti7EpuL+Sttdy91y1lZLsmaSyq10kjZX98ck5I9odbYNZ9Jad6W+0q6WhrMEZYkW6obBKsTgKFiucIP39w+H+/HGz4s/eHN9PmbOTB7e5/J+z+nEs3MvR4cxxggAAMBiAX09AQAAgO4QLAAAwHoECwAAsB7BAgAArEewAAAA6xEsAADAegQLAACwHsECAACsN6CvJ3CxtLa26vjx4woNDZXD4ejr6QAAgB4wxuj06dMaNmyYAgI6v4/iN8Fy/PhxxcTE9PU0AABAL1RXVys6OrrTr/tNsISGhko6t+CwsLA+ng0AAOiJxsZGxcTEeH6Od8ZvgqXtbaCwsDCCBQCAfqa7j3PwoVsAAGA9ggUAAFiPYAEAANYjWAAAgPUIFgAAYD2CBQAAWI9gAQAA1iNYAACA9QgWAABgPYIFAABYj2ABAADWI1gAAID1CBYAAGA9n4PlnXfe0fTp0zVs2DA5HA69/vrr3R7z9ttvKzExUSEhIfrmN7+pF154od2YLVu26MYbb1RwcLBuvPFGvfbaa75ODQAA+Cmfg+XLL7/UmDFj9Pzzz/dofGVlpb7//e9r0qRJ+uCDD/TYY4/pwQcf1JYtWzxj9u7dq1mzZmnOnDn68MMPNWfOHN15553697//7ev0AACAH3IYY0yvD3Y49Nprr2nmzJmdjnn00Ue1bds2ffzxx5598+fP14cffqi9e/dKkmbNmqXGxka98cYbnjG33367IiIi9Morr/RoLo2NjQoPD1dDQ4PCwsJ6tyAAAHBZ9fTn9yX/DMvevXuVlpbmtS89PV2lpaVyu91djikuLu70dZuamtTY2Oi1AQAA/3TJg6W2tlZRUVFe+6KiovTVV1+pvr6+yzG1tbWdvm5OTo7Cw8M9W0xMzMWfPAAAsMJl+S0hh8Ph9ee2d6HO39/RmK/vO192drYaGho8W3V19UWcMQAAsMmAS/0NhgwZ0u5OSV1dnQYMGKBrr722yzFfv+tyvuDgYAUHB1/8CQMAAOtc8jssycnJKioq8tq3Y8cOjR8/Xk6ns8sxKSkpl3p6AACgH/D5DssXX3yhI0eOeP5cWVmpffv2aeDAgRoxYoSys7N17Ngxbdq0SdK53wh6/vnntXjxYs2bN0979+7Viy++6PXbP4sWLdItt9yi3/72t8rMzNTf/vY3vfXWW9qzZ89FWCIAAOjvfL7DUlpaqrFjx2rs2LGSpMWLF2vs2LH6zW9+I0mqqalRVVWVZ/zIkSNVWFioXbt26Tvf+Y6WL1+uP/3pT/rRj37kGZOSkqKCggK99NJLGj16tP785z9r8+bNSkpKutD1AQAAP3BBz2GxCc9hAQCg/7HmOSwAAAAXimABAADWI1gAAID1CBYAAGA9ggUAAFiPYAEAANYjWAAAgPUIFgAAYD2CBQAAWI9gAQAA1iNYAACA9QgWAABgPYIFAABYj2ABAADWI1gAAID1CBYAAGA9ggUAAFiPYAEAANYjWAAAgPUIFgAAYD2CBQAAWI9gAQAA1iNYAACA9QgWAABgPYIFAABYj2ABAADWI1gAAID1CBYAAGA9ggUAAFiPYAEAANYjWAAAgPUIFgAAYD2CBQAAWI9gAQAA1iNYAACA9QgWAABgPYIFAABYj2ABAADWI1gAAID1CBYAAGA9ggUAAFiPYAEAANYjWAAAgPUIFgAAYD2CBQAAWI9gAQAA1iNYAACA9QgWAABgPYIFAABYj2ABAADWI1gAAID1CBYAAGA9ggUAAFiPYAEAANYjWAAAgPUIFgAAYD2CBQAAWI9gAQAA1iNYAACA9XoVLLm5uRo5cqRCQkKUmJio3bt3dzl+zZo1io+Pl8vl0qhRo7Rp06Z2Y/7whz9o1KhRcrlciomJ0UMPPaSzZ8/2ZnoAAMDPDPD1gM2bNysrK0u5ublKTU3V2rVrlZGRoYMHD2rEiBHtxufl5Sk7O1vr16/XhAkTVFJSonnz5ikiIkLTp0+XJL388stasmSJ8vPzlZKSokOHDmnu3LmSpN///vcXtkIAANDvOYwxxpcDkpKSNG7cOOXl5Xn2xcfHa+bMmcrJyWk3PiUlRampqVq5cqVnX1ZWlkpLS7Vnzx5J0sKFC/Xxxx/rn//8p2fMww8/rJKSkm7v3rRpbGxUeHi4GhoaFBYW5suSAABAH+npz2+f3hJqbm5WWVmZ0tLSvPanpaWpuLi4w2OampoUEhLitc/lcqmkpERut1uSdPPNN6usrEwlJSWSpE8//VSFhYW64447Op1LU1OTGhsbvTYAAOCffAqW+vp6tbS0KCoqymt/VFSUamtrOzwmPT1dGzZsUFlZmYwxKi0tVX5+vtxut+rr6yVJs2fP1vLly3XzzTfL6XTquuuu05QpU7RkyZJO55KTk6Pw8HDPFhMT48tSAABAP9KrD906HA6vPxtj2u1rs3TpUmVkZGjixIlyOp3KzMz0fD4lMDBQkrRr1y4988wzys3N1fvvv6+tW7fqH//4h5YvX97pHLKzs9XQ0ODZqqure7MUAADQD/gULJGRkQoMDGx3N6Wurq7dXZc2LpdL+fn5OnPmjI4ePaqqqirFxcUpNDRUkZGRks5FzZw5c3Tffffppptu0g9+8AOtWLFCOTk5am1t7fB1g4ODFRYW5rUBAAD/5FOwBAUFKTExUUVFRV77i4qKlJKS0uWxTqdT0dHRCgwMVEFBgaZNm6aAgHPf/syZM55/bhMYGChjjHz8TDAAAPBDPv9a8+LFizVnzhyNHz9eycnJWrdunaqqqjR//nxJ596qOXbsmOdZK4cOHVJJSYmSkpJ08uRJrV69WuXl5dq4caPnNadPn67Vq1dr7NixSkpK0pEjR7R06VLNmDHD87YRAAC4cvkcLLNmzdKJEye0bNky1dTUKCEhQYWFhYqNjZUk1dTUqKqqyjO+paVFq1atUkVFhZxOp6ZMmaLi4mLFxcV5xjzxxBNyOBx64okndOzYMQ0aNEjTp0/XM888c+ErBAAA/Z7Pz2GxFc9hAQCg/7kkz2EBAADoCwQLAACwHsECAACsR7AAAADrESwAAMB6BAsAALAewQIAAKxHsAAAAOsRLAAAwHoECwAAsB7BAgAArEewAAAA6xEsAADAegQLAACwHsECAACsR7AAAADrESwAAMB6BAsAALAewQIAAKxHsAAAAOsRLAAAwHoECwAAsB7BAgAArEewAAAA6xEsAADAegQLAACwHsECAACsR7AAAADrESwAAMB6BAsAALAewQIAAKxHsAAAAOsRLAAAwHoECwAAsB7BAgAArEewAAAA6xEsAADAegQLAACwHsECAACsR7AAAADrESwAAMB6BAsAALAewQIAAKxHsAAAAOsRLAAAwHoECwAAsB7BAgAArEewAAAA6xEsAADAegQLAACwHsECAACsR7AAAADrESwAAMB6BAsAALAewQIAAKxHsAAAAOsRLAAAwHoECwAAsB7BAgAArEewAAAA6/UqWHJzczVy5EiFhIQoMTFRu3fv7nL8mjVrFB8fL5fLpVGjRmnTpk3txpw6dUoLFizQ0KFDFRISovj4eBUWFvZmegAAwM8M8PWAzZs3KysrS7m5uUpNTdXatWuVkZGhgwcPasSIEe3G5+XlKTs7W+vXr9eECRNUUlKiefPmKSIiQtOnT5ckNTc363vf+54GDx6sV199VdHR0aqurlZoaOiFrxAAAPR7DmOM8eWApKQkjRs3Tnl5eZ598fHxmjlzpnJyctqNT0lJUWpqqlauXOnZl5WVpdLSUu3Zs0eS9MILL2jlypX65JNP5HQ6e7WQxsZGhYeHq6GhQWFhYb16DQAAcHn19Oe3T28JNTc3q6ysTGlpaV7709LSVFxc3OExTU1NCgkJ8drncrlUUlIit9stSdq2bZuSk5O1YMECRUVFKSEhQStWrFBLS0unc2lqalJjY6PXBgAA/JNPwVJfX6+WlhZFRUV57Y+KilJtbW2Hx6Snp2vDhg0qKyuTMUalpaXKz8+X2+1WfX29JOnTTz/Vq6++qpaWFhUWFuqJJ57QqlWr9Mwzz3Q6l5ycHIWHh3u2mJgYX5YCAAD6kV596NbhcHj92RjTbl+bpUuXKiMjQxMnTpTT6VRmZqbmzp0rSQoMDJQktba2avDgwVq3bp0SExM1e/ZsPf74415vO31ddna2GhoaPFt1dXVvlgIAAPoBn4IlMjJSgYGB7e6m1NXVtbvr0sblcik/P19nzpzR0aNHVVVVpbi4OIWGhioyMlKSNHToUF1//fWegJHOfS6mtrZWzc3NHb5ucHCwwsLCvDYAAOCffAqWoKAgJSYmqqioyGt/UVGRUlJSujzW6XQqOjpagYGBKigo0LRp0xQQcO7bp6am6siRI2ptbfWMP3TokIYOHaqgoCBfpggAAPyQz28JLV68WBs2bFB+fr4+/vhjPfTQQ6qqqtL8+fMlnXur5u677/aMP3TokP7yl7/o8OHDKikp0ezZs1VeXq4VK1Z4xtx///06ceKEFi1apEOHDmn79u1asWKFFixYcBGWCAAA+jufn8Mya9YsnThxQsuWLVNNTY0SEhJUWFio2NhYSVJNTY2qqqo841taWrRq1SpVVFTI6XRqypQpKi4uVlxcnGdMTEyMduzYoYceekijR4/W8OHDtWjRIj366KMXvkIAANDv+fwcFlvxHBYAAPqfS/IcFgAAgL5AsAAAAOsRLAAAwHoECwAAsB7BAgAArEewAAAA6xEsAADAegQLAACwHsECAACsR7AAAADrESwAAMB6BAsAALAewQIAAKxHsAAAAOsN6OsJ2MwYo/+5W/p6GgAAWMHlDJTD4eiT702wdOF/7hbd+Jv/19fTAADACgeXpesbQX2TDrwlBAAArMcdli64nIE6uCy9r6cBAIAVXM7APvveBEsXHA5Hn936AgAA/4e3hAAAgPUIFgAAYD2CBQAAWI9gAQAA1iNYAACA9QgWAABgPYIFAABYj2ABAADWI1gAAID1CBYAAGA9ggUAAFiPYAEAANYjWAAAgPUIFgAAYD2CBQAAWI9gAQAA1iNYAACA9QgWAABgPYIFAABYj2ABAADWI1gAAID1CBYAAGA9ggUAAFiPYAEAANYjWAAAgPUIFgAAYD2CBQAAWI9gAQAA1iNYAACA9QgWAABgPYIFAABYj2ABAADWI1gAAID1CBYAAGA9ggUAAFiPYAEAANYjWAAAgPUIFgAAYD2CBQAAWI9gAQAA1utVsOTm5mrkyJEKCQlRYmKidu/e3eX4NWvWKD4+Xi6XS6NGjdKmTZs6HVtQUCCHw6GZM2f2ZmoAAMAPDfD1gM2bNysrK0u5ublKTU3V2rVrlZGRoYMHD2rEiBHtxufl5Sk7O1vr16/XhAkTVFJSonnz5ikiIkLTp0/3GvvZZ5/pkUce0aRJk3q/IgAA4HccxhjjywFJSUkaN26c8vLyPPvi4+M1c+ZM5eTktBufkpKi1NRUrVy50rMvKytLpaWl2rNnj2dfS0uLJk+erF/84hfavXu3Tp06pddff73H82psbFR4eLgaGhoUFhbmy5IAAEAf6enPb5/eEmpublZZWZnS0tK89qelpam4uLjDY5qamhQSEuK1z+VyqaSkRG6327Nv2bJlGjRokO69994ezaWpqUmNjY1eGwAA8E8+BUt9fb1aWloUFRXltT8qKkq1tbUdHpOenq4NGzaorKxMxhiVlpYqPz9fbrdb9fX1kqR3331XL774otavX9/jueTk5Cg8PNyzxcTE+LIUAADQj/TqQ7cOh8Prz8aYdvvaLF26VBkZGZo4caKcTqcyMzM1d+5cSVJgYKBOnz6tn//851q/fr0iIyN7PIfs7Gw1NDR4turq6t4sBQAA9AM+feg2MjJSgYGB7e6m1NXVtbvr0sblcik/P19r167V559/rqFDh2rdunUKDQ1VZGSkPvroIx09etTrA7itra3nJjdggCoqKnTddde1e93g4GAFBwf7Mn0AANBP+XSHJSgoSImJiSoqKvLaX1RUpJSUlC6PdTqdio6OVmBgoAoKCjRt2jQFBATohhtu0P79+7Vv3z7PNmPGDE2ZMkX79u3jrR4AAOD7rzUvXrxYc+bM0fjx45WcnKx169apqqpK8+fPl3TurZpjx455nrVy6NAhlZSUKCkpSSdPntTq1atVXl6ujRs3SpJCQkKUkJDg9T2uueYaSWq3HwAAXJl8DpZZs2bpxIkTWrZsmWpqapSQkKDCwkLFxsZKkmpqalRVVeUZ39LSolWrVqmiokJOp1NTpkxRcXGx4uLiLtoiAACAf/P5OSy24jksAAD0P5fkOSwAAAB9gWABAADWI1gAAID1CBYAAGA9ggUAAFiPYAEAANYjWAAAgPUIFgAAYD2CBQAAWI9gAQAA1iNYAACA9QgWAABgPYIFAABYj2ABAADWI1gAAID1CBYAAGA9ggUAAFiPYAEAANYjWAAAgPUIFgAAYD2CBQAAWI9gAQAA1iNYAACA9QgWAABgPYIFAABYj2ABAADWI1gAAID1CBYAAGA9ggUAAFiPYAEAANYjWAAAgPUIFgAAYD2CBQAAWI9gAQAA1iNYAACA9QgWAABgPYIFAABYj2ABAADWI1gAAID1CBYAAGA9ggUAAFiPYAEAANYjWAAAgPUIFgAAYD2CBQAAWI9gAQAA1iNYAACA9QgWAABgPYIFAABYj2ABAADWI1gAAID1CBYAAGA9ggUAAFiPYAEAANYjWAAAgPUIFgAAYD2CBQAAWI9gAQAA1iNYAACA9QgWAABgvV4FS25urkaOHKmQkBAlJiZq9+7dXY5fs2aN4uPj5XK5NGrUKG3atMnr6+vXr9ekSZMUERGhiIgITZ06VSUlJb2ZGgAA8EM+B8vmzZuVlZWlxx9/XB988IEmTZqkjIwMVVVVdTg+Ly9P2dnZeuqpp3TgwAE9/fTTWrBggf7+9797xuzatUt33XWXdu7cqb1792rEiBFKS0vTsWPHer8yAADgNxzGGOPLAUlJSRo3bpzy8vI8++Lj4zVz5kzl5OS0G5+SkqLU1FStXLnSsy8rK0ulpaXas2dPh9+jpaVFERERev7553X33Xf3aF6NjY0KDw9XQ0ODwsLCfFkSAADoIz39+e3THZbm5maVlZUpLS3Na39aWpqKi4s7PKapqUkhISFe+1wul0pKSuR2uzs85syZM3K73Ro4cGCnc2lqalJjY6PXBgAA/JNPwVJfX6+WlhZFRUV57Y+KilJtbW2Hx6Snp2vDhg0qKyuTMUalpaXKz8+X2+1WfX19h8csWbJEw4cP19SpUzudS05OjsLDwz1bTEyML0sBAAD9SK8+dOtwOLz+bIxpt6/N0qVLlZGRoYkTJ8rpdCozM1Nz586VJAUGBrYb/9xzz+mVV17R1q1b292ZOV92drYaGho8W3V1dW+WAgAA+gGfgiUyMlKBgYHt7qbU1dW1u+vSxuVyKT8/X2fOnNHRo0dVVVWluLg4hYaGKjIy0mvs7373O61YsUI7duzQ6NGju5xLcHCwwsLCvDYAAOCffAqWoKAgJSYmqqioyGt/UVGRUlJSujzW6XQqOjpagYGBKigo0LRp0xQQ8H/ffuXKlVq+fLnefPNNjR8/3pdpAQAAPzfA1wMWL16sOXPmaPz48UpOTta6detUVVWl+fPnSzr3Vs2xY8c8z1o5dOiQSkpKlJSUpJMnT2r16tUqLy/Xxo0bPa/53HPPaenSpfrrX/+quLg4zx2cq6++WldfffXFWCcAAOjHfA6WWbNm6cSJE1q2bJlqamqUkJCgwsJCxcbGSpJqamq8nsnS0tKiVatWqaKiQk6nU1OmTFFxcbHi4uI8Y3Jzc9Xc3Kwf//jHXt/rySef1FNPPdW7lQEAAL/h83NYbMVzWAAA6H8uyXNYAAAA+gLBAgAArEewAAAA6xEsAADAegQLAACwHsECAACsR7AAAADrESwAAMB6BAsAALAewQIAAKxHsAAAAOsRLAAAwHoECwAAsB7BAgAArEewAAAA6xEsAADAegQLAACwHsECAACsR7AAAADrESwAAMB6BAsAALAewQIAAKxHsAAAAOsRLAAAwHoECwAAsB7BAgAArEewAAAA6xEsAADAegQLAACwHsECAACsR7AAAADrESwAAMB6BAsAALAewQIAAKxHsAAAAOsRLAAAwHoECwAAsB7BAgAArEewAAAA6xEsAADAegQLAACwHsECAACsR7AAAADrESwAAMB6BAsAALAewQIAAKxHsAAAAOsRLAAAwHoECwAAsB7BAgAArEewAAAA6xEsAADAegQLAACwHsECAACsR7AAAADrESwAAMB6BAsAALAewQIAAKxHsAAAAOsN6OsJXCzGGElSY2NjH88EAAD0VNvP7baf453xm2A5ffq0JCkmJqaPZwIAAHx1+vRphYeHd/p1h+kuafqJ1tZWHT9+XKGhoXI4HH09nUuusbFRMTExqq6uVlhYWF9P57K5UtctsfYrce1X6rol1n4lrd0Yo9OnT2vYsGEKCOj8kyp+c4clICBA0dHRfT2Nyy4sLOyK+Bf6667UdUus/Upc+5W6bom1Xylr7+rOShs+dAsAAKxHsAAAAOsRLP1UcHCwnnzySQUHB/f1VC6rK3XdEmu/Etd+pa5bYu1X6tq74jcfugUAAP6LOywAAMB6BAsAALAewQIAAKxHsAAAAOsRLBbJycnRhAkTFBoaqsGDB2vmzJmqqKjwfN3tduvRRx/VTTfdpKuuukrDhg3T3XffrePHj3u9zq233iqHw+G1zZ49+3Ivp8e6W7ckzZ07t92aJk6c6DWmqalJDzzwgCIjI3XVVVdpxowZ+u9//3s5l+Kznqz96+tu21auXOkZ09/OeV5enkaPHu15MFZycrLeeOMNz9eNMXrqqac0bNgwuVwu3XrrrTpw4IDXa/TH8y11vXZ/vcbbdHfe/fU6727d/niNXxIG1khPTzcvvfSSKS8vN/v27TN33HGHGTFihPniiy+MMcacOnXKTJ061WzevNl88sknZu/evSYpKckkJiZ6vc7kyZPNvHnzTE1NjWc7depUXyypR7pbtzHG3HPPPeb222/3WtOJEye8Xmf+/Plm+PDhpqioyLz//vtmypQpZsyYMearr7663EvqsZ6s/fw119TUmPz8fONwOMx//vMfz5j+ds63bdtmtm/fbioqKkxFRYV57LHHjNPpNOXl5cYYY5599lkTGhpqtmzZYvbv329mzZplhg4dahobGz2v0R/PtzFdr91fr/E23Z13f73Ou1u3P17jlwLBYrG6ujojybz99tudjikpKTGSzGeffebZN3nyZLNo0aLLMMNLo6N133PPPSYzM7PTY06dOmWcTqcpKCjw7Dt27JgJCAgwb7755qWc7kXVk3OemZlpbrvtNq99/f2cG2NMRESE2bBhg2ltbTVDhgwxzz77rOdrZ8+eNeHh4eaFF14wxvjP+W7TtvaO+OM1fr7z136lXOfGdH3O/fUav1C8JWSxhoYGSdLAgQO7HONwOHTNNdd47X/55ZcVGRmpb3/723rkkUc8f5t1f9DZunft2qXBgwfr+uuv17x581RXV+f5WllZmdxut9LS0jz7hg0bpoSEBBUXF1+eiV8E3Z3zzz//XNu3b9e9997b7mv99Zy3tLSooKBAX375pZKTk1VZWana2lqvcxkcHKzJkyd7zqW/nO+vr70j/niNS52v3d+v8+7OuT9e4xeL3/zlh/7GGKPFixfr5ptvVkJCQodjzp49qyVLluinP/2p11+Q9bOf/UwjR47UkCFDVF5eruzsbH344YcqKiq6XNPvtc7WnZGRoZ/85CeKjY1VZWWlli5dqttuu01lZWUKDg5WbW2tgoKCFBER4fV6UVFRqq2tvdzL6JWenPONGzcqNDRUP/zhD73298dzvn//fiUnJ+vs2bO6+uqr9dprr+nGG2/0/OCJioryGh8VFaXPPvtMkvr9+e5s7V/nj9d4V2v35+u8p+fcn67xi66P7/CgE7/+9a9NbGysqa6u7vDrzc3NJjMz04wdO9Y0NDR0+VqlpaVGkikrK7sUU72oult3m+PHjxun02m2bNlijDHm5ZdfNkFBQe3GTZ061fzqV7+6JHO92Hqy9lGjRpmFCxd2+1r94Zw3NTWZw4cPm/fee88sWbLEREZGmgMHDph3333XSDLHjx/3Gn/fffeZ9PR0Y0z/P9+drf18/nqN92TtbfzpOu/puv3pGr/YeEvIQg888IC2bdumnTt3Kjo6ut3X3W637rzzTlVWVqqoqKjbv3583LhxcjqdOnz48KWa8kXR3brPN3ToUMXGxnrWNGTIEDU3N+vkyZNe4+rq6tr9n7qNerL23bt3q6KiQvfdd1+3r9cfznlQUJC+9a1vafz48crJydGYMWP0xz/+UUOGDJGkdv/HfP657O/nu7O1t/HXa1zqfu3n86frvCfr9rdr/GIjWCxijNHChQu1detW/etf/9LIkSPbjWn7D9nhw4f11ltv6dprr+32dQ8cOCC3262hQ4deimlfsJ6s++tOnDih6upqz5oSExPldDq9bo/W1NSovLxcKSkpl2zuF8qXtb/44otKTEzUmDFjun1d2895R4wxampq8tz2Pv9cNjc36+233/acy/56vjvTtnbJP6/xrpy/9q/zl+u8Ix2t29+v8QvWh3d38DX333+/CQ8PN7t27fL61bUzZ84YY4xxu91mxowZJjo62uzbt89rTFNTkzHGmCNHjpinn37avPfee6aystJs377d3HDDDWbs2LHW/tpfd+s+ffq0efjhh01xcbGprKw0O3fuNMnJyWb48OHtfs01OjravPXWW+b99983t912m/W/7tjd2ts0NDSYb3zjGyYvL6/da/THc56dnW3eeecdU1lZaT766CPz2GOPmYCAALNjxw5jzLlfaw4PDzdbt241+/fvN3fddVeHv9bc3863MV2v3V+v8TZdrd2fr/Pu/n03xv+u8UuBYLGIpA63l156yRhjTGVlZadjdu7caYwxpqqqytxyyy1m4MCBJigoyFx33XXmwQcfbPcsA5t0t+4zZ86YtLQ0M2jQION0Os2IESPMPffcY6qqqrxe53//+59ZuHChGThwoHG5XGbatGntxtimu7W3Wbt2rXG5XB0+d6E/nvNf/vKXJjY21gQFBZlBgwaZ7373u17/8W5tbTVPPvmkGTJkiAkODja33HKL2b9/v9dr9MfzbUzXa/fXa7xNV2v35+u8u3/fjfG/a/xScBhjzOW4kwMAANBbfIYFAABYj2ABAADWI1gAAID1CBYAAGA9ggUAAFiPYAEAANYjWAAAgPUIFgAAYD2CBQAAWI9gAQAA1iNYAACA9QgWAABgvf8PFAFQL6qNGWMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "score = []\n",
    "for i in range(390, 200, -10):\n",
    "    x_fschi = SelectKBest(chi2, k=i).fit_transform(X, y)\n",
    "    once = cross_val_score(\n",
    "        RFC(n_estimators=10, random_state=0), x_fschi, y, cv=5).mean()\n",
    "    score.append(once)\n",
    "plt.plot(range(390, 200, -10), score)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 10.81782088   3.7107283  116.31261309  67.0483602 ]\n",
      "[4.47651499e-03 1.56395980e-01 5.53397228e-26 2.75824965e-15]\n"
     ]
    }
   ],
   "source": [
    "chivalue, pvalues_chi = chi2(X, y)\n",
    "print(chivalue)\n",
    "print(pvalues_chi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3\n"
     ]
    }
   ],
   "source": [
    "k = chivalue.shape[0] - (pvalues_chi > 0.05).sum()\n",
    "print(k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "F检验"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.int64(4)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.feature_selection import f_classif\n",
    "F, pvalues_f = f_classif(X, y)\n",
    "k = F.shape[0] - (pvalues_f > 0.05).sum()\n",
    "k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "互信息法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.int64(4)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.feature_selection import mutual_info_regression as MIC\n",
    "result = MIC(X, y)\n",
    "k = X.shape[1] - sum(result <= 0)\n",
    "k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Embedding嵌入法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.feature_selection import SelectFromModel\n",
    "from sklearn.ensemble import RandomForestClassifier as RFC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3.68063403e-04, 1.55507755e-04, 3.40179337e-06, ...,\n",
       "       4.08596449e-06, 0.00000000e+00, 3.36483008e-06])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rfc = RFC()\n",
    "rfc.fit(X, y)\n",
    "fi = rfc.feature_importances_\n",
    "fi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.        , 0.0002761 , 0.00055221, 0.00082831, 0.00110442,\n",
       "       0.00138052, 0.00165663, 0.00193273, 0.00220884, 0.00248494,\n",
       "       0.00276104, 0.00303715, 0.00331325, 0.00358936, 0.00386546,\n",
       "       0.00414157, 0.00441767, 0.00469378, 0.00496988, 0.00524598])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "threshold = np.linspace(0, fi.max(), 20)\n",
    "threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "score = []\n",
    "for i in threshold:\n",
    "    x_embedded = SelectFromModel(rfc, threshold=i).fit_transform(X, y)\n",
    "    once = cross_val_score(RFC(n_estimators=10, random_state=0),\n",
    "                           x_embedded, y, cv=5).mean()\n",
    "    score.append(once)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAgO0lEQVR4nO3de3DU1f3/8deybC7VJI1GI5hAom0hFLGQ0JBEvg6jTUwFoZcROi1KL0xxsCXETiUqxcJIrBR6GUiEQCxMW8NUxNIaW9JWbgabJoKK0KCDmAwmzYTRBOXbJC7n94e/7NdtArKBsO9sn4+Zz4z5cD675xwZ85zPXvQ455wAAAAMGxbuCQAAAHwcggUAAJhHsAAAAPMIFgAAYB7BAgAAzCNYAACAeQQLAAAwj2ABAADmDQ/3BC6WM2fO6O2331ZcXJw8Hk+4pwMAAM6Dc06nTp3SyJEjNWzY2e+jREywvP3220pNTQ33NAAAwAA0NzcrJSXlrH8eMcESFxcn6cMFx8fHh3k2AADgfHR2dio1NTXwe/xsIiZYel8Gio+PJ1gAABhiPu7tHLzpFgAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5BAsAADCPYAEAAOYRLAAAwDyCBQAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5IQfLnj17NGPGDI0cOVIej0fPPPPMx16ze/duZWZmKiYmRtddd50ef/zxPmO2bdumcePGKTo6WuPGjdP27dtDnRoAAIhQIQfL+++/rxtvvFFr1649r/FvvvmmvvjFL2rq1Kk6cOCAHnjgAX3/+9/Xtm3bAmP279+v2bNna+7cuXr55Zc1d+5c3Xnnnfr73/8e6vQAAEAE8jjn3IAv9ni0fft2zZo166xj7r//fu3YsUNHjhwJnFuwYIFefvll7d+/X5I0e/ZsdXZ26rnnnguMue2225SYmKgnn3zyvObS2dmphIQEdXR0KD4+fmALAgAAl9T5/v4e9Pew7N+/X/n5+UHnCgoKVF9fr56ennOOqa2tPevjdnV1qbOzM+gAAACRadCDpbW1VcnJyUHnkpOT9cEHH6i9vf2cY1pbW8/6uKWlpUpISAgcqampF3/yAADAhEvyKSGPxxP0c++rUB8939+Y/zz3USUlJero6Agczc3NF3HGAADAkuGD/QTXXHNNnzslbW1tGj58uK688spzjvnPuy4fFR0drejo6Is/YQAAYM6g32HJyclRTU1N0LmdO3cqKytLPp/vnGNyc3MHe3oAAGAICPkOy3vvvac33ngj8PObb76pgwcP6oorrtCoUaNUUlKiEydOaMuWLZI+/ETQ2rVrVVxcrPnz52v//v3atGlT0Kd/Fi1apP/5n//RT37yE82cOVO///3v9Ze//EX79u27CEsEAABDXch3WOrr6zVx4kRNnDhRklRcXKyJEyfqRz/6kSSppaVFTU1NgfHp6emqrq7Wrl279LnPfU4rVqzQL3/5S33lK18JjMnNzVVVVZWeeOIJTZgwQb/61a+0detWZWdnX+j6AABABLig72GxhO9hAQBg6DHzPSwAAAAXimABAADmESwAAMA8ggUAAJhHsAAAAPMIFgAAYB7BAgAAzCNYAACAeQQLAAAwj2ABAADmESwAAMA8ggUAAJhHsAAAAPMIFgAAYB7BAgAAzCNYAACAeQQLAAAwj2ABAADmESwAAMA8ggUAAJhHsAAAAPMIFgAAYB7BAgAAzCNYAACAeQQLAAAwj2ABAADmESwAAMA8ggUAAJhHsAAAAPMIFgAAYB7BAgAAzCNYAACAeQQLAAAwj2ABAADmESwAAMA8ggUAAJhHsAAAAPMIFgAAYB7BAgAAzCNYAACAeQQLAAAwj2ABAADmESwAAMA8ggUAAJhHsAAAAPMIFgAAYB7BAgAAzCNYAACAeQQLAAAwj2ABAADmESwAAMA8ggUAAJhHsAAAAPMIFgAAYB7BAgAAzCNYAACAeQQLAAAwj2ABAADmESwAAMC8AQVLWVmZ0tPTFRMTo8zMTO3du/ec49etW6eMjAzFxsZqzJgx2rJlS58xP//5zzVmzBjFxsYqNTVVixcv1r///e+BTA8AAESY4aFesHXrVhUVFamsrEx5eXlav369CgsLdfjwYY0aNarP+PLycpWUlKiiokKTJ09WXV2d5s+fr8TERM2YMUOS9Jvf/EZLlixRZWWlcnNzdfToUc2bN0+S9LOf/ezCVggAAIY8j3POhXJBdna2Jk2apPLy8sC5jIwMzZo1S6WlpX3G5+bmKi8vT6tWrQqcKyoqUn19vfbt2ydJuvfee3XkyBH99a9/DYy57777VFdX97F3b3p1dnYqISFBHR0dio+PD2VJAAAgTM7393dILwl1d3eroaFB+fn5Qefz8/NVW1vb7zVdXV2KiYkJOhcbG6u6ujr19PRIkm666SY1NDSorq5OknTs2DFVV1fr9ttvP+tcurq61NnZGXQAAIDIFFKwtLe3y+/3Kzk5Oeh8cnKyWltb+72moKBAGzduVENDg5xzqq+vV2VlpXp6etTe3i5JmjNnjlasWKGbbrpJPp9P119/vaZNm6YlS5acdS6lpaVKSEgIHKmpqaEsBQAADCEDetOtx+MJ+tk51+dcr6VLl6qwsFBTpkyRz+fTzJkzA+9P8Xq9kqRdu3bpkUceUVlZmV566SU9/fTT+uMf/6gVK1acdQ4lJSXq6OgIHM3NzQNZCgAAGAJCCpakpCR5vd4+d1Pa2tr63HXpFRsbq8rKSp0+fVrHjx9XU1OT0tLSFBcXp6SkJEkfRs3cuXP1ne98RzfccIO+9KUvaeXKlSotLdWZM2f6fdzo6GjFx8cHHQAAIDKFFCxRUVHKzMxUTU1N0Pmamhrl5uae81qfz6eUlBR5vV5VVVVp+vTpGjbsw6c/ffp04J97eb1eOecU4nuCAQBABAr5Y83FxcWaO3eusrKylJOTow0bNqipqUkLFiyQ9OFLNSdOnAh818rRo0dVV1en7OxsvfPOO1qzZo0OHTqkzZs3Bx5zxowZWrNmjSZOnKjs7Gy98cYbWrp0qe64447Ay0YAAOC/V8jBMnv2bJ08eVLLly9XS0uLxo8fr+rqao0ePVqS1NLSoqampsB4v9+v1atXq7GxUT6fT9OmTVNtba3S0tICYx566CF5PB499NBDOnHihK666irNmDFDjzzyyIWvEAAADHkhfw+LVXwPCwAAQ8+gfA8LAABAOBAsAADAPIIFAACYR7AAAADzCBYAAGAewQIAAMwjWAAAgHkECwAAMI9gAQAA5hEsAADAPIIFAACYR7AAAADzCBYAAGAewQIAAMwjWAAAgHkECwAAMI9gAQAA5hEsAADAPIIFAACYR7AAAADzCBYAAGAewQIAAMwjWAAAgHkECwAAMI9gAQAA5hEsAADAPIIFAACYR7AAAADzCBYAAGAewQIAAMwjWAAAgHkECwAAMI9gAQAA5hEsAADAPIIFAACYR7AAAADzCBYAAGAewQIAAMwjWAAAgHkECwAAMI9gAQAA5hEsAADAPIIFAACYR7AAAADzCBYAAGAewQIAAMwjWAAAgHkECwAAMI9gAQAA5hEsAADAPIIFAACYR7AAAADzCBYAAGAewQIAAMwjWAAAgHkECwAAMI9gAQAA5hEsAADAPIIFAACYR7AAAADzBhQsZWVlSk9PV0xMjDIzM7V3795zjl+3bp0yMjIUGxurMWPGaMuWLX3GvPvuu1q4cKFGjBihmJgYZWRkqLq6eiDTAwAAEWZ4qBds3bpVRUVFKisrU15entavX6/CwkIdPnxYo0aN6jO+vLxcJSUlqqio0OTJk1VXV6f58+crMTFRM2bMkCR1d3frC1/4gq6++mo99dRTSklJUXNzs+Li4i58hQAAYMjzOOdcKBdkZ2dr0qRJKi8vD5zLyMjQrFmzVFpa2md8bm6u8vLytGrVqsC5oqIi1dfXa9++fZKkxx9/XKtWrdI///lP+Xy+AS2ks7NTCQkJ6ujoUHx8/IAeAwAAXFrn+/s7pJeEuru71dDQoPz8/KDz+fn5qq2t7fearq4uxcTEBJ2LjY1VXV2denp6JEk7duxQTk6OFi5cqOTkZI0fP14rV66U3+8/61y6urrU2dkZdAAAgMgUUrC0t7fL7/crOTk56HxycrJaW1v7vaagoEAbN25UQ0ODnHOqr69XZWWlenp61N7eLkk6duyYnnrqKfn9flVXV+uhhx7S6tWr9cgjj5x1LqWlpUpISAgcqampoSwFAAAMIQN6063H4wn62TnX51yvpUuXqrCwUFOmTJHP59PMmTM1b948SZLX65UknTlzRldffbU2bNigzMxMzZkzRw8++GDQy07/qaSkRB0dHYGjubl5IEsBAABDQEjBkpSUJK/X2+duSltbW5+7Lr1iY2NVWVmp06dP6/jx42pqalJaWpri4uKUlJQkSRoxYoQ+85nPBAJG+vB9Ma2treru7u73caOjoxUfHx90AACAyBRSsERFRSkzM1M1NTVB52tqapSbm3vOa30+n1JSUuT1elVVVaXp06dr2LAPnz4vL09vvPGGzpw5Exh/9OhRjRgxQlFRUaFMEQAARKCQXxIqLi7Wxo0bVVlZqSNHjmjx4sVqamrSggULJH34Us1dd90VGH/06FH9+te/1uuvv666ujrNmTNHhw4d0sqVKwNj7rnnHp08eVKLFi3S0aNH9eyzz2rlypVauHDhRVgiAAAY6kL+HpbZs2fr5MmTWr58uVpaWjR+/HhVV1dr9OjRkqSWlhY1NTUFxvv9fq1evVqNjY3y+XyaNm2aamtrlZaWFhiTmpqqnTt3avHixZowYYKuvfZaLVq0SPfff/+FrxAAAAx5IX8Pi1V8DwsAAEPPoHwPCwAAQDgQLAAAwDyCBQAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5BAsAADCPYAEAAOYRLAAAwDyCBQAAmEewAAAA8wgWAABgHsECAADMGx7uCVjmnNP/9vjDPQ0AAEyI9Xnl8XjC8twEyzn8b49f437053BPAwAAEw4vL9AnosKTDrwkBAAAzOMOyznE+rw6vLwg3NMAAMCEWJ83bM9NsJyDx+MJ260vAADwf3hJCAAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5BAsAADCPYAEAAOYRLAAAwDyCBQAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5BAsAADCPYAEAAOYRLAAAwDyCBQAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5BAsAADCPYAEAAOYRLAAAwDyCBQAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5BAsAADCPYAEAAOYRLAAAwDyCBQAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5BAsAADCPYAEAAOYRLAAAwLwBBUtZWZnS09MVExOjzMxM7d2795zj161bp4yMDMXGxmrMmDHasmXLWcdWVVXJ4/Fo1qxZA5kaAACIQMNDvWDr1q0qKipSWVmZ8vLytH79ehUWFurw4cMaNWpUn/Hl5eUqKSlRRUWFJk+erLq6Os2fP1+JiYmaMWNG0Ni33npLP/jBDzR16tSBrwgAAEQcj3POhXJBdna2Jk2apPLy8sC5jIwMzZo1S6WlpX3G5+bmKi8vT6tWrQqcKyoqUn19vfbt2xc45/f7dfPNN+ub3/ym9u7dq3fffVfPPPPMec+rs7NTCQkJ6ujoUHx8fChLAgAAYXK+v79Dekmou7tbDQ0Nys/PDzqfn5+v2trafq/p6upSTExM0LnY2FjV1dWpp6cncG758uW66qqr9O1vf/u85tLV1aXOzs6gAwAARKaQgqW9vV1+v1/JyclB55OTk9Xa2trvNQUFBdq4caMaGhrknFN9fb0qKyvV09Oj9vZ2SdILL7ygTZs2qaKi4rznUlpaqoSEhMCRmpoaylIAAMAQMqA33Xo8nqCfnXN9zvVaunSpCgsLNWXKFPl8Ps2cOVPz5s2TJHm9Xp06dUrf+MY3VFFRoaSkpPOeQ0lJiTo6OgJHc3PzQJYCAACGgJDedJuUlCSv19vnbkpbW1ufuy69YmNjVVlZqfXr1+tf//qXRowYoQ0bNiguLk5JSUl65ZVXdPz48aA34J45c+bDyQ0frsbGRl1//fV9Hjc6OlrR0dGhTB8AAAxRId1hiYqKUmZmpmpqaoLO19TUKDc395zX+nw+paSkyOv1qqqqStOnT9ewYcM0duxYvfrqqzp48GDguOOOOzRt2jQdPHiQl3oAAEDoH2suLi7W3LlzlZWVpZycHG3YsEFNTU1asGCBpA9fqjlx4kTgu1aOHj2quro6ZWdn65133tGaNWt06NAhbd68WZIUExOj8ePHBz3HJz/5SUnqcx4AAPx3CjlYZs+erZMnT2r58uVqaWnR+PHjVV1drdGjR0uSWlpa1NTUFBjv9/u1evVqNTY2yufzadq0aaqtrVVaWtpFWwQAAIhsIX8Pi1V8DwsAAEPPoHwPCwAAQDgQLAAAwDyCBQAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5BAsAADCPYAEAAOYRLAAAwDyCBQAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5BAsAADCPYAEAAOYRLAAAwDyCBQAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5BAsAADCPYAEAAOYRLAAAwDyCBQAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5BAsAADCPYAEAAOYRLAAAwDyCBQAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5BAsAADCPYAEAAOYRLAAAwDyCBQAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5BAsAADCPYAEAAOYRLAAAwDyCBQAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5BAsAADCPYAEAAOYRLAAAwDyCBQAAmEewAAAA8wYULGVlZUpPT1dMTIwyMzO1d+/ec45ft26dMjIyFBsbqzFjxmjLli1Bf15RUaGpU6cqMTFRiYmJuvXWW1VXVzeQqQEAgAgUcrBs3bpVRUVFevDBB3XgwAFNnTpVhYWFampq6nd8eXm5SkpK9PDDD+u1117Tj3/8Yy1cuFB/+MMfAmN27dqlr33ta3r++ee1f/9+jRo1Svn5+Tpx4sTAVwYAACKGxznnQrkgOztbkyZNUnl5eeBcRkaGZs2apdLS0j7jc3NzlZeXp1WrVgXOFRUVqb6+Xvv27ev3Ofx+vxITE7V27Vrddddd5zWvzs5OJSQkqKOjQ/Hx8aEsCQAAhMn5/v4O6Q5Ld3e3GhoalJ+fH3Q+Pz9ftbW1/V7T1dWlmJiYoHOxsbGqq6tTT09Pv9ecPn1aPT09uuKKK846l66uLnV2dgYdAAAgMoUULO3t7fL7/UpOTg46n5ycrNbW1n6vKSgo0MaNG9XQ0CDnnOrr61VZWamenh61t7f3e82SJUt07bXX6tZbbz3rXEpLS5WQkBA4UlNTQ1kKAAAYQgb0pluPxxP0s3Ouz7leS5cuVWFhoaZMmSKfz6eZM2dq3rx5kiSv19tn/GOPPaYnn3xSTz/9dJ87Mx9VUlKijo6OwNHc3DyQpQAAgCEgpGBJSkqS1+vtczelra2tz12XXrGxsaqsrNTp06d1/PhxNTU1KS0tTXFxcUpKSgoa+9Of/lQrV67Uzp07NWHChHPOJTo6WvHx8UEHAACITCEFS1RUlDIzM1VTUxN0vqamRrm5uee81ufzKSUlRV6vV1VVVZo+fbqGDfu/p1+1apVWrFihP/3pT8rKygplWgAAIMIND/WC4uJizZ07V1lZWcrJydGGDRvU1NSkBQsWSPrwpZoTJ04Evmvl6NGjqqurU3Z2tt555x2tWbNGhw4d0ubNmwOP+dhjj2np0qX67W9/q7S0tMAdnMsvv1yXX375xVgnAAAYwkIOltmzZ+vkyZNavny5WlpaNH78eFVXV2v06NGSpJaWlqDvZPH7/Vq9erUaGxvl8/k0bdo01dbWKi0tLTCmrKxM3d3d+upXvxr0XMuWLdPDDz88sJUBAICIEfL3sFjF97AAADD0DMr3sAAAAIQDwQIAAMwjWAAAgHkECwAAMI9gAQAA5hEsAADAPIIFAACYR7AAAADzCBYAAGAewQIAAMwjWAAAgHkECwAAMI9gAQAA5hEsAADAPIIFAACYR7AAAADzCBYAAGAewQIAAMwjWAAAgHkECwAAMI9gAQAA5hEsAADAPIIFAACYR7AAAADzCBYAAGAewQIAAMwjWAAAgHkECwAAMI9gAQAA5hEsAADAPIIFAACYR7AAAADzCBYAAGAewQIAAMwjWAAAgHkECwAAMI9gAQAA5hEsAADAPIIFAACYR7AAAADzCBYAAGAewQIAAMwjWAAAgHkECwAAMI9gAQAA5hEsAADAPIIFAACYR7AAAADzCBYAAGAewQIAAMwjWAAAgHkECwAAMI9gAQAA5hEsAADAPIIFAACYR7AAAADzCBYAAGAewQIAAMwjWAAAgHnDwz2Bi8U5J0nq7OwM80wAAMD56v293ft7/GwiJlhOnTolSUpNTQ3zTAAAQKhOnTqlhISEs/65x31c0gwRZ86c0dtvv624uDh5PJ6L9ridnZ1KTU1Vc3Oz4uPjL9rj/rdjXwcH+zo42NfBwb4OjqG2r845nTp1SiNHjtSwYWd/p0rE3GEZNmyYUlJSBu3x4+Pjh8S/+KGGfR0c7OvgYF8HB/s6OIbSvp7rzkov3nQLAADMI1gAAIB5BMvHiI6O1rJlyxQdHR3uqUQU9nVwsK+Dg30dHOzr4IjUfY2YN90CAIDIxR0WAABgHsECAADMI1gAAIB5BAsAADAv4oOlrKxM6enpiomJUWZmpvbu3XvO8bt371ZmZqZiYmJ03XXX6fHHH+8zZtu2bRo3bpyio6M1btw4bd++/YKfd6gJx77u2bNHM2bM0MiRI+XxePTMM89czCWZEI59LS0t1eTJkxUXF6err75as2bNUmNj40VdV7iFY1/Ly8s1YcKEwJd35eTk6Lnnnruo6wq3cP33tVdpaak8Ho+KiooudCmmhGNfH374YXk8nqDjmmuuuajrumAuglVVVTmfz+cqKirc4cOH3aJFi9xll13m3nrrrX7HHzt2zH3iE59wixYtcocPH3YVFRXO5/O5p556KjCmtrbWeb1et3LlSnfkyBG3cuVKN3z4cPfiiy8O+HmHmnDta3V1tXvwwQfdtm3bnCS3ffv2wV7qJRWufS0oKHBPPPGEO3TokDt48KC7/fbb3ahRo9x777036Gu+FMK1rzt27HDPPvusa2xsdI2Nje6BBx5wPp/PHTp0aNDXfCmEa1971dXVubS0NDdhwgS3aNGiwVrmJReufV22bJn77Gc/61paWgJHW1vboK83FBEdLJ///OfdggULgs6NHTvWLVmypN/xP/zhD93YsWODzn33u991U6ZMCfx85513uttuuy1oTEFBgZszZ86An3eoCde+flQkBouFfXXOuba2NifJ7d69O9QlmGRlX51zLjEx0W3cuDGU6ZsVzn09deqU+/SnP+1qamrczTffHFHBEq59XbZsmbvxxhsvcPaDK2JfEuru7lZDQ4Py8/ODzufn56u2trbfa/bv399nfEFBgerr69XT03POMb2POZDnHUrCta+RztK+dnR0SJKuuOKKkNdhjZV99fv9qqqq0vvvv6+cnJyBLseMcO/rwoULdfvtt+vWW2+90KWYEu59ff311zVy5Eilp6drzpw5Onbs2IUu6aKK2GBpb2+X3+9XcnJy0Pnk5GS1trb2e01ra2u/4z/44AO1t7efc0zvYw7keYeScO1rpLOyr845FRcX66abbtL48eMHuhwzwr2vr776qi6//HJFR0drwYIF2r59u8aNG3ehywq7cO5rVVWVXnrpJZWWll6MpZgSzn3Nzs7Wli1b9Oc//1kVFRVqbW1Vbm6uTp48eTGWdlFEzP+t+Ww8Hk/Qz865Puc+bvx/nj+fxwz1eYeacO1rpAv3vt5777165ZVXtG/fvpDmbV249nXMmDE6ePCg3n33XW3btk133323du/eHRHRIl36fW1ubtaiRYu0c+dOxcTEXNDcLQvH39fCwsLAP99www3KycnR9ddfr82bN6u4uDj0RQyCiA2WpKQkeb3ePlXa1tbWpzR7XXPNNf2OHz58uK688spzjul9zIE871ASrn2NdBb29Xvf+5527NihPXv2KCUl5UKWY0a49zUqKkqf+tSnJElZWVn6xz/+oV/84hdav379Ba0r3MK1rw0NDWpra1NmZmbgz/1+v/bs2aO1a9eqq6tLXq/3gtcXLuH++/pRl112mW644Qa9/vrrA1nKoIjYl4SioqKUmZmpmpqaoPM1NTXKzc3t95qcnJw+43fu3KmsrCz5fL5zjul9zIE871ASrn2NdOHcV+ec7r33Xj399NP629/+pvT09IuxJBOs/X11zqmrqyvUZZgTrn295ZZb9Oqrr+rgwYOBIysrS1//+td18ODBIR0rkq2/r11dXTpy5IhGjBgxkKUMjkv5Dt9LrffjYZs2bXKHDx92RUVF7rLLLnPHjx93zjm3ZMkSN3fu3MD43o+HLV682B0+fNht2rSpz8fDXnjhBef1et2jjz7qjhw54h599NGzfqz5bM871IVrX0+dOuUOHDjgDhw44CS5NWvWuAMHDkTcx8Uv9b7ec889LiEhwe3atSvoI42nT5++dIsfROHa15KSErdnzx735ptvuldeecU98MADbtiwYW7nzp2XbvGDKFz7+p8i7VNC4drX++67z+3atcsdO3bMvfjii2769OkuLi7O1O+tiA4W55xbt26dGz16tIuKinKTJk0K+qjm3Xff7W6++eag8bt27XITJ050UVFRLi0tzZWXl/d5zN/97nduzJgxzufzubFjx7pt27aF9LyRIBz7+vzzzztJfY677757MJYYFuHY1/72VJJ74oknBmOJYRGOff3Wt74VeM6rrrrK3XLLLRETK73C9d/Xj4q0YHEuPPs6e/ZsN2LECOfz+dzIkSPdl7/8Zffaa68NyvoGyuPc/393DgAAgFER+x4WAAAQOQgWAABgHsECAADMI1gAAIB5BAsAADCPYAEAAOYRLAAAwDyCBQAAmEewAAAA8wgWAABgHsECAADMI1gAAIB5/w936uDbfLJeZAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.96\n"
     ]
    }
   ],
   "source": [
    "plt.plot(threshold, score)\n",
    "plt.show()\n",
    "X_embedded = SelectFromModel(rfc, threshold=0.00067).fit_transform(X, y)\n",
    "X_embedded.shape\n",
    "print(cross_val_score(rfc, X_embedded, y, cv=5).mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 因子分析"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 使用sklearn实现因子分析"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn.decomposition import FactorAnalysis\n",
    "from sklearn.datasets import load_iris\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.metrics import accuracy_score\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   sepal length (cm)  sepal width (cm)  petal length (cm)  petal width (cm)\n",
      "0                5.1               3.5                1.4               0.2\n",
      "1                4.9               3.0                1.4               0.2\n",
      "2                4.7               3.2                1.3               0.2\n",
      "3                4.6               3.1                1.5               0.2\n",
      "4                5.0               3.6                1.4               0.2\n"
     ]
    }
   ],
   "source": [
    "# 加载示例数据集\n",
    "data = load_iris()\n",
    "X = data.data\n",
    "y = data.target\n",
    "feature_names = data.feature_names\n",
    "\n",
    "# 将数据转换为DataFrame\n",
    "df = pd.DataFrame(X, columns=feature_names)\n",
    "print(df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHFCAYAAAAOmtghAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAABR5UlEQVR4nO3deVxU9eL/8dcMO7K4gqCouKGIqGkqllumpqUtdu/ttnerm9lm5je1Tc26qO2r5r22WNeuv8JdS63ElUpLRUFNcwEVRFwAUdY5vz9MikQEHDjDzPv5eJxHcvicmfecx8l5e1aLYRgGIiIiIk7CanYAEREREXtSuRERERGnonIjIiIiTkXlRkRERJyKyo2IiIg4FZUbERERcSoqNyIiIuJUVG5ERETEqajciIiIiFNRuRGRSvvhhx+4+eabadasGV5eXgQHBxMTE8NTTz1ldrRK+/jjj7FYLCWTu7s7TZs25b777uPw4cMl4+Lj47FYLMTHx1f6PTZu3MikSZM4deqU/YKLyEWp3IhIpSxbtoxevXqRnZ3N9OnTWblyJW+99RZXXXUV8+bNMztelX300UckJCSwatUqHnzwQT7//HN69+5Nbm7uZb/2xo0bmTx5ssqNSA1xNzuAiNQu06dPJzw8nBUrVuDu/vtfIbfddhvTp0+3y3ucOXMGX19fu7xWRUVFRdGtWzcA+vfvT3FxMVOmTGHhwoXccccdNZpFRC6P9tyISKUcP36chg0blio251mtF/6VMnfuXGJiYvDz88PPz4/OnTsze/bskt/369ePqKgo1q5dS69evfD19eUf//gHANnZ2YwdO5bw8HA8PT1p0qQJo0ePvmBvimEYvP/++3Tu3BkfHx/q1avHrbfeyr59+6r8OXv27AnAwYMHyx23ePFiYmJi8PX1xd/fn4EDB5KQkFDy+0mTJvF///d/AISHh5cc/qrK4S0RqRiVGxGplJiYGH744Qcef/xxfvjhBwoLCy869oUXXuCOO+4gNDSUjz/+mAULFnDPPfdcUBjS0tK48847uf3221m+fDmjRo3izJkz9O3bl08++YTHH3+cr776inHjxvHxxx8zfPhwDMMoWf6hhx5i9OjRXHvttSxcuJD333+fpKQkevXqxdGjR6v0Offu3QtAo0aNLjpm7ty53HjjjQQEBPD5558ze/ZsTp48Sb9+/Vi/fj0ADzzwAI899hgA8+fPJyEhgYSEBK644ooq5RKRCjBERCohMzPTuPrqqw3AAAwPDw+jV69eRmxsrJGTk1Mybt++fYabm5txxx13lPt6ffv2NQDj22+/LTU/NjbWsFqtxqZNm0rN//LLLw3AWL58uWEYhpGQkGAAxmuvvVZqXGpqquHj42M8/fTT5b7/Rx99ZADG999/bxQWFho5OTnG0qVLjUaNGhn+/v5Genq6YRiGsXr1agMwVq9ebRiGYRQXFxuhoaFGx44djeLi4pLXy8nJMYKCgoxevXqVzHvllVcMwNi/f3+5WUTEPrTnRkQqpUGDBqxbt45NmzYxdepUbrzxRn755RcmTJhAx44dyczMBGDVqlUUFxfzyCOPXPI169WrxzXXXFNq3tKlS4mKiqJz584UFRWVTIMHDy51WGfp0qVYLBbuvPPOUuMaN25Mp06dKnz4p2fPnnh4eODv788NN9xA48aN+eqrrwgODi5z/O7duzly5Ah33XVXqcNxfn5+jBgxgu+//54zZ85U6L1FxL50QrGIVEm3bt1KTsAtLCxk3LhxvPHGG0yfPp3p06dz7NgxAJo2bXrJ1woJCblg3tGjR9m7dy8eHh5lLnO+RB09ehTDMC5aQlq2bFmhzzNnzhzat2+Pu7s7wcHBZWb6o+PHj180e2hoKDabjZMnT9b4idEionIjInbg4eHBxIkTeeONN9ixYwfw+7kqhw4dIiwsrNzlLRbLBfMaNmyIj48PH374YZnLNGzYsOS/FouFdevW4eXldcG4suaVpX379iVlrSIaNGgAnDtf6M+OHDmC1WqlXr16FX49EbEflRsRqZS0tLQy91bs3LkTOLfXAmDQoEG4ubkxY8YMYmJiKv0+N9xwA//6179o0KAB4eHh5Y6bOnUqhw8f5q9//Wul36eqIiIiaNKkCXPnzmXs2LElBS03N5e4uLiSK6jg94J19uzZGssn4spUbkSkUgYPHkzTpk0ZNmwY7dq1w2azsXXrVl577TX8/Px44oknAGjRogXPPPMMU6ZM4ezZs/z9738nMDCQ5ORkMjMzmTx5crnvM3r0aOLi4ujTpw9PPvkk0dHR2Gw2UlJSWLlyJU899RQ9evTgqquu4p///Cf33Xcfmzdvpk+fPtSpU4e0tDTWr19Px44defjhh+2+HqxWK9OnT+eOO+7ghhtu4KGHHiI/P59XXnmFU6dOMXXq1JKxHTt2BOCtt97innvuwcPDg4iICPz9/e2eS0TQ1VIiUjnz5s0zbr/9dqNNmzaGn5+f4eHhYTRr1sy46667jOTk5AvGz5kzx7jyyisNb29vw8/Pz+jSpYvx0Ucflfy+b9++RocOHcp8r9OnTxvPPfecERERYXh6ehqBgYFGx44djSeffLLkKqbzPvzwQ6NHjx5GnTp1DB8fH6NVq1bG3XffbWzevLncz3P+aqk/X5X1Z3++Wuq8hQsXGj169DC8vb2NOnXqGAMGDDA2bNhwwfITJkwwQkNDDavVWubriIj9WAzjDzeLEBEREanldCm4iIiIOBWVGxEREXEqKjciIiLiVFRuRERExKmo3IiIiIhTUbkRERERp+JyN/Gz2WwcOXIEf3//Mm/5LiIiIo7HMAxycnIIDQ0t9bDasrhcuTly5Mgln3MjIiIijik1NfWSD+R1uXJz/nbnqampBAQEmJxGREREKiI7O5uwsLAKPbbE5crN+UNRAQEBKjciIiK1TEVOKdEJxSIiIuJUVG5ERETEqajciIiIiFMxtdzMmDGD6OjokvNfYmJi+Oqrr8pdZs2aNXTt2hVvb29atmzJzJkzayitiIiI1AamlpumTZsydepUNm/ezObNm7nmmmu48cYbSUpKKnP8/v37GTp0KL1792bLli0888wzPP7448TFxdVwchEREXFUFsMwDLND/FH9+vV55ZVXuP/++y/43bhx41i8eDE7d+4smTdy5Ei2bdtGQkJChV4/OzubwMBAsrKydLWUiIhILVGZ72+HOeemuLiY//3vf+Tm5hITE1PmmISEBAYNGlRq3uDBg9m8eTOFhYVlLpOfn092dnapSURERJyX6eVm+/bt+Pn54eXlxciRI1mwYAGRkZFljk1PTyc4OLjUvODgYIqKisjMzCxzmdjYWAIDA0sm3Z1YRETEuZlebiIiIti6dSvff/89Dz/8MPfccw/JyckXHf/nm/ecP6p2sZv6TJgwgaysrJIpNTXVfuFFRETE4Zh+h2JPT09at24NQLdu3di0aRNvvfUWH3zwwQVjGzduTHp6eql5GRkZuLu706BBgzJf38vLCy8vL/sH/5Nim8GP+0+QkZNHkL833cPr42bVgzlFRERqmunl5s8MwyA/P7/M38XExLBkyZJS81auXEm3bt3w8PCoiXhl+npHGpOXJJOWlVcyLyTQm4nDIrkuKsS0XCIiIq7I1MNSzzzzDOvWrePAgQNs376dZ599lvj4eO644w7g3CGlu+++u2T8yJEjOXjwIGPGjGHnzp18+OGHzJ49m7Fjx5r1Efh6RxoPf/ZzqWIDkJ6Vx8Of/czXO9JMSiYiIuKaTN1zc/ToUe666y7S0tIIDAwkOjqar7/+moEDBwKQlpZGSkpKyfjw8HCWL1/Ok08+yXvvvUdoaChvv/02I0aMMCV/sc1g8pJkyrqW3gAswOQlyQyMbKxDVCIiIjXE4e5zU93seZ+bhF+P8/d/f3/JcZ8/2JOYVmWfEyQiIiKXVivvc1MbZeTkXXpQJcaJiIjI5VO5uQxB/t52HSciIiKXT+XmMnQPr09IoDcXO5vGwrmrprqH16/JWCIiIi5N5eYyuFktTBx27m7KZRUcA5g4LFInE4uIiNQglZvLdF1UCDPuvILGgRceenKzQouGdUxIJSIi4rp0tZSdlL5DsRez1+/nm50ZdAqry/yHe2nvjYiIyGWozPe3w92huLZys1pKXe4d3tCPH/atYVvqKT7asJ8Herc0MZ2IiIjr0GGpatI40Jtnrm8PwKsrd5Ny/IzJiURERFyDyk01uu3KMHq2rE9eoY0JCxJxsSOAIiIiplC5qUYWi4Wpt0Tj5W5lw97jfLH5kNmRREREnJ7KTTVr0bAOTw1qC8CUZclkZOtuxSIiItVJ5aYG/OOqcDo2CSQnr4gXFiWZHUdERMSpqdzUAHc3K9NGRONutfB1UjpfbU8zO5KIiIjTUrmpIZGhAYzs2wqA5xclkXWm0OREIiIizknlpgY9ek1rWjWqQ+bpfF5almx2HBEREaekclODvD3cmDYiGosFvvjpEOv3ZJodSURExOmo3NSwbi3qc3fP5gCMn5/ImYIikxOJiIg4F5UbE/zfde0IDfTm0MmzvLbyF7PjiIiIOBWVGxP4ebnz8i0dAfhww35+TjlpciIRERHnoXJjkv4RQdzSpQmGAePjEikospkdSURExCmo3Jjo+RsiaVDHk1+Onub9+L1mxxEREXEKKjcmqlfHk0nDOwDw3uq9/HI0x+REIiIitZ/KjcluiA7h2vZBFBYbPP1lIsU2PTlcRETkcqjcmMxisTDlpij8vdzZmnqKjzceMDuSiIhIraZy4wBCAn0YP7QdAK+u2E3qiTMmJxIREam9VG4cxN+vbEaP8PqcLSxmwvztGIYOT4mIiFSFyo2DsFotTB0RjZe7lfV7M/nyp0NmRxIREamVVG4cSHjDOjw5sC0AU5Ymk5GTZ3IiERGR2kflxsE8cHU4UU0CyM4rYtLiJLPjiIiI1DoqNw7G3c3KtBHRuFktLN+eztc70syOJCIiUquo3DigDqGBjOzbEoDnFyWRdabQ5EQiIiK1h8qNg3rsmja0bFSHYzn5/Gv5TrPjiIiI1BoqNw7K28ONaSOiAZi3OZUNezNNTiQiIlI7qNw4sCtb1Oeuns0BmDB/O2cKikxOJCIi4vhUbhzc09dFEBroTcqJM7y+8hez44iIiDg8lRsH5+/twcs3dwTgww372Zp6ytxAIiIiDk7lphbo3y6ImzqHYjNg3JeJFBTZzI4kIiLisFRuaokXhnWgfh1Pdh/NYeaaX82OIyIi4rBUbmqJ+nU8mTgsEoB3vtvDnqM5JicSERFxTCo3tcjwTqEMaBdEYbHB03GJFNv05HAREZE/U7mpRSwWCy/dHIWflztbUk4xJ+GA2ZFEREQcjspNLRMS6MP4Ie0AmP71blJPnDE5kYiIiGNRuamFbu/ejO7h9TlbWMwzC7ZjGDo8JSIicp7KTS1ktVqYektHPN2trNuTSdzPh82OJCIi4jBUbmqplo38GH1tGwCmLE3mWE6+yYlEREQcg8pNLfZg75Z0CA0g62whkxYnmR1HRETEIZhabmJjY7nyyivx9/cnKCiIm266id27d5e7THx8PBaL5YJp165dNZTacXi4WZk2Iho3q4Vl29NYkZRudiQRERHTmVpu1qxZwyOPPML333/PqlWrKCoqYtCgQeTm5l5y2d27d5OWllYytWnTpgYSO56oJoH8s09LAJ5fuIOss4UmJxIRETGXu5lv/vXXX5f6+aOPPiIoKIiffvqJPn36lLtsUFAQdevWrcZ0tccTA9qwYkc6+zJzmfrVTmJviTY7koiIiGkc6pybrKwsAOrXr3/JsV26dCEkJIQBAwawevXqi47Lz88nOzu71ORsvD3ciL3l3JPDP/8xlY17M01OJCIiYh6HKTeGYTBmzBiuvvpqoqKiLjouJCSEWbNmERcXx/z584mIiGDAgAGsXbu2zPGxsbEEBgaWTGFhYdX1EUzVo2UD7uzZDIDx87dztqDY5EQiIiLmsBgOcge4Rx55hGXLlrF+/XqaNm1aqWWHDRuGxWJh8eLFF/wuPz+f/PzfL5POzs4mLCyMrKwsAgICLju3I8nJK2TQG2tJy8rjn31a8szQ9mZHEhERsYvs7GwCAwMr9P3tEHtuHnvsMRYvXszq1asrXWwAevbsyZ49e8r8nZeXFwEBAaUmZ+Xv7cHLN5/b6/WfdfvYlnrK3EAiIiImMLXcGIbBo48+yvz58/nuu+8IDw+v0uts2bKFkJAQO6erna5pF8zwTqHYDBgXl0hBkc3sSCIiIjXK1KulHnnkEebOncuiRYvw9/cnPf3cfVoCAwPx8fEBYMKECRw+fJg5c+YA8Oabb9KiRQs6dOhAQUEBn332GXFxccTFxZn2ORzNxGGRrNtzjF3pOXyw5lceG+Cal8mLiIhrMnXPzYwZM8jKyqJfv36EhISUTPPmzSsZk5aWRkpKSsnPBQUFjB07lujoaHr37s369etZtmwZt9xyixkfwSE18PNi4rAOALzz3V72ZuSYnEhERKTmOMwJxTWlMick1WaGYfCPjzexevcxujavxxcPxWC1WsyOJSIiUiW17oRisT+LxcLLN3ekjqcbPx08yaffHzQ7koiISI1QuXFioXV9GD+kHQDTvt7FoZNnTE4kIiJS/VRunNwdPZpzZYt6nCko5pkFO3Cxo5AiIuKCVG6cnNVqYeqIaDzdraz95RgLthw2O5KIiEi1UrlxAa0a+fHEb5eDv7g0mczT+ZdYQkREpPZSuXER/+zTksiQAE6dKWTS4iSz44iIiFQblRsX4eFmZfqt0bhZLSxNTGNV8lGzI4mIiFQLlRsXEtUkkAd6n3vExXMLt5OdV2hyIhEREftTuXExT17blhYNfDmanU/s8l1mxxEREbE7lRsX4+3hxtQR0QB8/mMKCb8eNzmRiIiIfancuKCeLRtwe49mAEyYn0heYbHJiUREROxH5cZFjR/SjsYB3hw4foY3vvnF7DgiIiJ2o3LjogK8PXjppigA/r12H4mHTpkbSERExE5UblzYtZHBDOsUis2Ap79MpLDYZnYkERGRy6Zy4+ImDoukrq8Hu9JzmLV2n9lxRERELpvKjYtr6OfFxGGRALz1zR72Zpw2OZGIiMjlUbkRburchH4RjSgotjE+LhGbTU8OFxGR2kvlRrBYLLx0UxR1PN3YfPAkn/1w0OxIIiIiVaZyIwA0refL09e1A2DaV7s4fOqsyYlERESqRuVGStzVszndmtcjt6CYZxdsxzB0eEpERGoflRspYbVamDoiGk83K/G7j7Fo6xGzI4mIiFSayo2U0jrIj8cHtAZg8pIkMk/nm5xIRESkclRu5AIP9W1Fu8b+nDxTyOQlyWbHERERqRSVG7mAh5uV6bdGY7XAkm1H+Cb5qNmRREREKkzlRsoU3bQuD/ZuCcBzC3eQnVdociIREZGKUbmRixp9bVuaN/AlPTuPaV/tMjuOiIhIhajcyEX5eLoRe0tHAP77Qwrf7ztuciIREZFLU7mRcvVq1ZC/dw8DYML87eQVFpucSEREpHwqN3JJ44e0JzjAi/2Zubz5zR6z44iIiJRL5UYuKdDHg5duOnd46t/r9rHjcJbJiURERC5O5UYqZGBkMNdHh1BsM3j6y0QKi21mRxIRESmTyo1U2KRhHajr60FyWjb/XrfP7DgiIiJlUrmRCmvk78Xz10cC8OY3e/j12GmTE4mIiFxI5UYq5ZYrmtCnbSMKimyMj0vEZtOTw0VExLGo3EilWCwW/nVzFL6ebmw6cJL//phidiQREZFSVG6k0prW8+XpwREATF2+kyOnzpqcSERE5HcqN1Ild8W04IpmdcktKOa5hTswDB2eEhERx6ByI1XiZrUwbUQ0nm5WvtuVweJtR8yOJCIiAqjcyGVoE+zPo9e0BmDykmSOn843OZGIiIjKjVymkX1b0a6xPydyC3hxabLZcURERFRu5PJ4uluZNiIaqwUWbT3Cd7uOmh1JRERcnMqNXLZOYXW5/+pwAJ5dsIOcvEKTE4mIiCtTuRG7GDMwgmb1fUnLymP617vNjiMiIi5M5UbswsfTjam3nHty+KffH+TH/SdMTiQiIq5K5Ubsplfrhtx2ZRgA4+MSySssNjmRiIi4IpUbsasJQ9sT5O/Fvsxc3v52j9lxRETEBanciF0F+ngw5aYoAD5Yu48dh7NMTiQiIq7G1HITGxvLlVdeib+/P0FBQdx0003s3n3pk1HXrFlD165d8fb2pmXLlsycObMG0kpFDe7QmKEdG1NsMxgXl0hRsc3sSCIi4kJMLTdr1qzhkUce4fvvv2fVqlUUFRUxaNAgcnNzL7rM/v37GTp0KL1792bLli0888wzPP7448TFxdVgcrmUScM7EOjjQdKRbP69br/ZcURExIVYDAd64uGxY8cICgpizZo19OnTp8wx48aNY/HixezcubNk3siRI9m2bRsJCQmXfI/s7GwCAwPJysoiICDAbtnlQl/+dIixX2zDy93KV0/0pmUjP7MjiYhILVWZ72+HOucmK+vc+Rn169e/6JiEhAQGDRpUat7gwYPZvHkzhYUX3jwuPz+f7OzsUpPUjBFXNKF3m4bkF9kYP387NpvD9GgREXFiDlNuDMNgzJgxXH311URFRV10XHp6OsHBwaXmBQcHU1RURGZm5gXjY2NjCQwMLJnCwsLsnl3KZrFY+NfNHfH1dOPH/Sf4fFOK2ZFERMQFOEy5efTRR0lMTOTzzz+/5FiLxVLq5/NH1v48H2DChAlkZWWVTKmpqfYJLBUSVt+XsYMiAIhdvou0rLMmJxIREWfnEOXmscceY/HixaxevZqmTZuWO7Zx48akp6eXmpeRkYG7uzsNGjS4YLyXlxcBAQGlJqlZ9/RqQZdmdTmdX8RzC3bgQKd5iYiIEzK13BiGwaOPPsr8+fP57rvvCA8Pv+QyMTExrFq1qtS8lStX0q1bNzw8PKorqlwGN6uFaSOi8XCz8O2uDJYkppkdSUREnJip5eaRRx7hs88+Y+7cufj7+5Oenk56ejpnz/5+6GLChAncfffdJT+PHDmSgwcPMmbMGHbu3MmHH37I7NmzGTt2rBkfQSqobbA/j/ZvA8CkxUmcyC0wOZGIiDgrU8vNjBkzyMrKol+/foSEhJRM8+bNKxmTlpZGSsrvJ6KGh4ezfPly4uPj6dy5M1OmTOHtt99mxIgRZnwEqYSH+7UiItifE7kFTFmabHYcERFxUg51n5uaoPvcmGtr6ilueX8DNgM+uu9K+kcEmR1JRERqgVp7nxtxfp3D6nLfVefOrXp2/nZO5xeZnEhERJyNyo3UuKcGtSWsvg9HsvKY/vUus+OIiIiTUbmRGufr6c7UW6IB+PT7g2w6cMLkRCIi4kxUbsQUV7VuyF+7NcUwYFxcInmFxWZHEhERJ6FyI6Z5dmgkjfy92Hcsl3e/22t2HBERcRIqN2KaQF8PptzYAYCZa34l6UiWyYlERMQZqNyIqa6LCmFIVGOKbAbj4hIpKraZHUlERGo5lRsx3eQbOxDg7c6Ow9nMXr/f7DgiIlLLqdyI6YL8vXnuhkgAXl/1C/szc01OJCIitZnKjTiEv3RtytWtG5JfZGN8XCI2m0vdOFtEROxI5UYcgsVi4V83d8THw40f9p/gf5tSzY4kIiK1lMqNOIxmDXwZOzgCgNjlO0nPyjM5kYiI1EYqN+JQ7u3Vgs5hdcnJL+K5hTtwsee6ioiIHajciENxs1qYNiIaDzcL3+w8yrLtaWZHEhGRWsYu5aa4uJitW7dy8uRJe7ycuLiIxv6M6tcagImLkjiZW2ByIhERqU2qVG5Gjx7N7NmzgXPFpm/fvlxxxRWEhYURHx9vz3ziokb1b0XbYD+O5xYwZWmy2XFERKQWqVK5+fLLL+nUqRMAS5YsYf/+/ezatYvRo0fz7LPP2jWguCYvdzemjojGYoH5Ww4TvzvD7EgiIlJLVKncZGZm0rhxYwCWL1/OX/7yF9q2bcv999/P9u3b7RpQXNcVzepxX69wAJ5dsIPT+UUmJxIRkdqgSuUmODiY5ORkiouL+frrr7n22msBOHPmDG5ubnYNKK5t7OC2NK3nw+FTZ3l1xW6z44iISC1QpXJz33338de//pWoqCgsFgsDBw4E4IcffqBdu3Z2DSiuzdfTndhbOgLwScIBfjp4wuREIiLi6KpUbiZNmsR//vMf/vnPf7Jhwwa8vLwAcHNzY/z48XYNKNK7TSP+0rUphgFPf5lIXmGx2ZFERMSBWYzLvEtaXl4e3t7e9spT7bKzswkMDCQrK4uAgACz40gFZZ0pZMDra8g8nc9j17TmqUERZkcSEZEaVJnv7yrtuSkuLmbKlCk0adIEPz8/9u3bB8Dzzz9fcom4iD0F+now5cYOAMyI/5WdadkmJxIREUdVpXLz8ssv8/HHHzN9+nQ8PT1L5nfs2JH//Oc/dgsn8kdDOoYwuEMwRTaDcXGJFBXbzI4kIiIOqErlZs6cOcyaNYs77rij1NVR0dHR7Nq1y27hRP5syo1R+Hu7k3goi482HDA7joiIOKAqlZvDhw/TunXrC+bbbDYKCwsvO5TIxQQFePPc9e0BeG3Vbg5k5pqcSEREHE2Vyk2HDh1Yt27dBfO/+OILunTpctmhRMrz125h9GrVgLxCGxPmb9eTw0VEpBT3qiw0ceJE7rrrLg4fPozNZmP+/Pns3r2bOXPmsHTpUntnFCnFYrEw9ZZoBr25hoR9x5m3KZXbujczO5aIiDiIKu25GTZsGPPmzWP58uVYLBZeeOEFdu7cyZIlS0pu6CdSnZo18GXsb5eDv7x8J0ez80xOJCIijuKy73NT2+g+N86j2GZwy/sb2HYoi4GRwcy6qysWi8XsWCIiUg2q/T43Io7AzWph2q3RuFstrEo+yvLt6WZHEhERB1ClcmO1WnFzc7voJFJT2jUOYFT/c1fuTVy8g1NnCkxOJCIiZqvSCcULFiwo9XNhYSFbtmzhk08+YfLkyXYJJlJRj/RvxfLtaezNOM2UpTt57a+dzI4kIiImsus5N3PnzmXevHksWrTIXi9pdzrnxjn9dPAkt87ciGHAnH90p0/bRmZHEhEROzLtnJsePXrwzTff2PMlRSqka/N63BPTAoAJ87eTm19kbiARETGN3crN2bNneeedd2jatKm9XlKkUv5vcARN6vpw+NRZXlmx2+w4IiJikiqdc1OvXr1Sl9wahkFOTg6+vr589tlndgsnUhl1vNyJvaUjd3/4I58kHGBYp1C6Nq9ndiwREalhVSo3b7zxRqlyY7VaadSoET169KBePX2ZiHn6tG3EiCuaEvfzIcbFJbLs8avxctcVfCIirkQ38ROnc+pMAde+vobM0wU8PqANYwa2NTuSiIhcpsp8f1d4z01iYmKFA0RHR1d4rIi91fX1ZPLwKB6Z+zPvr97L0I6NaddYRVZExFVUuNx07twZi8VyyScwWywWiouLLzuYyOUY2rExgyKDWZl8lHFfJjJ/1FW4WfVoBhERV1DhcrN///7qzCFiVxaLhSk3RZGw7zjbDmXx0Yb9PNC7pdmxRESkBlS43DRv3rw6c4jYXXCAN88Obc/4+dt5deVuBkU2plkDX7NjiYhINavS1VLnJScnk5KSQkFB6ef5DB8+/LJCidjL364MY9HWIyTsO874+Yn894EeenK4iIiTq1K52bdvHzfffDPbt28vdR7O+S8NnXMjjsJisRB7S0eue2stG389zhebD/HXK8PMjiUiItWoSncofuKJJwgPD+fo0aP4+vqSlJTE2rVr6datG/Hx8XaOKHJ5WjSsU3I5+JRlyRzNzjM5kYiIVKcqlZuEhARefPFFGjVqhNVqxWq1cvXVVxMbG8vjjz9e4ddZu3Ytw4YNIzQ0FIvFwsKFC8sdHx8fj8ViuWDatWtXVT6GuJB/XBVOdNNAcvKKeGHRDrPjiIhINapSuSkuLsbPzw+Ahg0bcuTIEeDcSce7d1f8mT65ubl06tSJd999t1Lvv3v3btLS0kqmNm3aVGp5cT3ublamjYjG3WphRdJRvtqeZnYkERGpJlU65yYqKorExERatmxJjx49mD59Op6ensyaNYuWLSt+ue2QIUMYMmRIpd8/KCiIunXrVno5cW3tQwJ4uF8r3vluL88vSiKmVQPq+nqaHUtEROysSntunnvuOWw2GwAvvfQSBw8epHfv3ixfvpy3337brgHL0qVLF0JCQhgwYACrV6+u9vcT5/HoNa1p1agOmafzeXnZTrPjiIhINajSnpvBgweX/Llly5YkJydz4sSJC54Wbm8hISHMmjWLrl27kp+fz6effsqAAQOIj4+nT58+ZS6Tn59Pfn5+yc/Z2dnVlk8cn5e7G9NvjebWmQl88dMhhncOpXebRmbHEhERO6rSnptPPvmE3NzcUvPq169f7fcPiYiI4MEHH+SKK64gJiaG999/n+uvv55XX331osvExsYSGBhYMoWF6TJgV9e1eX3uiWkBwIT52zlTUGRuIBERsasqlZuxY8cSFBTEbbfdxtKlSykqMu/LoWfPnuzZs+eiv58wYQJZWVklU2pqag2mE0f1f4MjaFLXh0Mnz/Lqil/MjiMiInZUpXKTlpbGvHnzcHNz47bbbiMkJIRRo0axceNGe+e7pC1bthASEnLR33t5eREQEFBqEqnj5c7LN0cB8NHG/fycctLkRCIiYi9VKjfu7u7ccMMN/Pe//yUjI4M333yTgwcP0r9/f1q1alXh1zl9+jRbt25l69atwLmHc27dupWUlBTg3F6Xu+++u2T8m2++ycKFC9mzZw9JSUlMmDCBuLg4Hn300ap8DHFx/SKCuKVLEwwDxn2ZSH6R7qwtIuIMLuvZUgC+vr4MHjyYkydPcvDgQXburPgVKJs3b6Z///4lP48ZMwaAe+65h48//pi0tLSSogNQUFDA2LFjOXz4MD4+PnTo0IFly5YxdOjQy/0Y4qKevyGSNb8cY0/Gad5f/StP/nYnYxERqb0sxvkHQ1XSmTNnWLBgAf/973/55ptvCAsL4+9//zt33HEH7du3t3dOu8nOziYwMJCsrCwdohIAlmw7wmOfb8HDzcLSx3oT0djf7EgiIvInlfn+rtKem7///e8sWbIEX19f/vKXvxAfH0+vXr2qFFbEbDdEh7Bo6xG+2XmUcXGJxD3cCzernhwuIlJbVancWCwW5s2bx+DBg3F3v+wjWyKmslgsvHRTFD/sO87W1FN8vPEA918dbnYsERGpoiqdUDx37lyuv/56FRtxGo0DvZkw9Nzh1FdX7Cb1xBmTE4mISFVVuZ18++23fPvtt2RkZJQ8iuG8Dz/88LKDidS0264MY/G2w3y/7wQT5m/n0/u7V/uNKUVExP6qtOdm8uTJDBo0iG+//ZbMzExOnjxZahKpjaxWC1NvicbL3cr6vZl88dMhsyOJiEgVVGnPzcyZM/n444+566677J1HxFQtGtZhzMC2xH61i5eWJtMvohFB/t5mxxIRkUqo0p6bgoICXR0lTuv+q8Pp2CSQ7LwiJi5KMjuOiIhUUpXKzQMPPMDcuXPtnUXEIbi7WZk2Ihp3q4WvdqTz9Y40syOJiEglVOmwVF5eHrNmzeKbb74hOjoaDw+PUr9//fXX7RJOxCyRoQE81Lcl763+lecXJRHTsiGBvh6XXlBERExXpXKTmJhI586dAdixY0ep3+nqEnEWj13Thq92pLPvWC4vL09m+q2dzI4kIiIVUOXHL9RWevyCVMamAyf4y8wEAP77QA+uat3Q5EQiIq6pMt/fVTrn5ry9e/eyYsUKzp49C4CL9SRxAVe2qM/dMc0BGD8/kTMFRSYnEhGRS6lSuTl+/DgDBgygbdu2DB06lLS0cydcPvDAAzz11FN2DShitqeva0dooDepJ87y+spfzI4jIiKXUKVy8+STT+Lh4UFKSgq+vr4l8//2t7/x9ddf2y2ciCPw83Ln5Vs6AvDhhv1sTT1lbiARESlXlcrNypUrmTZtGk2bNi01v02bNhw8eNAuwUQcSf+IIG7u0gSbAeO+TKSgyHbphURExBRVKje5ubml9ticl5mZiZeX12WHEnFEz98QSf06nuw+msOM+F/NjiMiIhdRpXLTp08f5syZU/KzxWLBZrPxyiuv0L9/f7uFE3Ek9et4Mml4BwDeXb2HPUdzTE4kIiJlqdJ9bl555RX69evH5s2bKSgo4OmnnyYpKYkTJ06wYcMGe2cUcRjDokNYtOUw3+7K4Om4RL4c2Qs3q+7tJCLiSKq05yYyMpLExES6d+/OwIEDyc3N5ZZbbmHLli20atXK3hlFHIbFYuGlm6Pw83JnS8opPtl4wOxIIiLyJ7qJn0gVfPb9QZ5buAMfDzdWPtmHsPoXnoMmIiL2U5nv7yo/fqEsFosFb29vmjVrphOLxand3r0Zi7cd4cf9J3hmwXbm/KO7Hj0iIuIgqlRuOnfuXPIX+fkdP3/8i93Dw4O//e1vfPDBB3h7e9shpohjsVotTL2lI0PeWse6PZnE/XyYW7s2vfSCIiJS7ap0zs2CBQto06YNs2bNYtu2bWzdupVZs2YRERHB3LlzmT17Nt999x3PPfecvfOKOIyWjfwYfW1bAKYsTeZYTr7JiUREBKp4zk337t2ZMmUKgwcPLjV/xYoVPP/88/z4448sXLiQp556il9/daz7geicG7GnomIbN72/gR2Hs7m+Ywjv3XGF2ZFERJxStT84c/v27TRv3vyC+c2bN2f79u3AuUNX5585JeKs3N2sTBsRjZvVwrLtaaxISjc7koiIy6tSuWnXrh1Tp06loKCgZF5hYSFTp06lXbt2ABw+fJjg4GD7pBRxYB1CA3moT0sAnl+4g6yzhSYnEhFxbVU6ofi9995j+PDhNG3alOjoaCwWC4mJiRQXF7N06VIA9u3bx6hRo+waVsRRPT6gDV/vSGdfZi6xy3cydUS02ZFERFxWle9zc/r0aT777DN++eUXDMOgXbt23H777fj7+9s7o13pnBupLj/uP8FfP0gAYO4DPejVuqHJiUREnEdlvr91Ez8RO3pu4XY++z6FZvV9WTG6Dz6ebmZHEhFxCtVyE7/FixczZMgQPDw8WLx4cbljhw8fXtGXFXEq465rx7c7M0g5cYbXV+3m2esjzY4kIuJyKrznxmq1kp6eTlBQEFbrxc9DtlgsFBcX2y2gvWnPjVS373Yd5R8fb8ZqgQWjrqJTWF2zI4mI1HrVcim4zWYjKCio5M8Xmxy52IjUhGvaBXNj51BsBoyLS6SgyGZ2JBERl1KpS8GHDh1KVlZWyc8vv/wyp06dKvn5+PHjREZqN7zICzdEUr+OJ7vSc/hgjWPdyFJExNlVqtysWLGC/PzfbzE/bdo0Tpw4UfJzUVERu3fvtl86kVqqgZ8XE4edK/rvfLeXvRk5JicSEXEdlSo3fz49x8UutBKplOGdQrmmXRAFxTae/jKRYpv+fxERqQlVukOxiFyaxWLhpZui8PNy5+eUU3yacMDsSCIiLqFS5cZisWCxWC6YJyJlC63rw7gh5x5JMn3Fbg6dPGNyIhER51epxy8YhsG9996Ll5cXAHl5eYwcOZI6deoAlDofR0TOuaN7M5ZsPcKPB07wzIIdfHLflfpHgYhINarUHYrvu+++Co376KOPqhyouuk+N2KGX4+dZshb6ygosvHaXzoxomtTsyOJiNQqevxCOVRuxCzvrd7LKyt2U9fXg1VP9qWRv5fZkUREao1quYmfiFyef/ZpSWRIAKfOFDJpSZLZcUREnJbKjUgN8XCzMv3WaNysFpYlprEyKd3sSCIiTknlRqQGRTUJ5MHeLQF4ftEOsvMKTU4kIuJ8VG5Eatjoa9sQ3rAOR7PziV2+y+w4IiJOR+VGpIZ5e7gx9ZaOAHz+YwoJvx43OZGIiHNRuRExQY+WDbijRzMAJsxP5GxBscmJRESch8qNiEnGD2lH4wBvDhw/w5vf/GJ2HBERp2FquVm7di3Dhg0jNDQUi8XCwoULL7nMmjVr6Nq1K97e3rRs2ZKZM2dWf1CRauDv7cFLN0UB8O91+0g8dMrcQCIiTsLUcpObm0unTp149913KzR+//79DB06lN69e7NlyxaeeeYZHn/8ceLi4qo5qUj1uDYymGGdQrEZ8PSXiRQW28yOJCJS61Xq2VL2NmTIEIYMGVLh8TNnzqRZs2a8+eabALRv357Nmzfz6quvMmLEiGpKKVK9Jg6LZP2eY+xKz+GDNb/y6DVtzI4kIlKr1apzbhISEhg0aFCpeYMHD2bz5s0UFpZ9v5D8/Hyys7NLTSKOpKGfFxOHdQDg7W/3sjfjtMmJRERqt1pVbtLT0wkODi41Lzg4mKKiIjIzM8tcJjY2lsDAwJIpLCysJqKKVMqNnUPpF9GIgmIb4+MSsdlc6pFvIiJ2VavKDYDFYin18/nnfv55/nkTJkwgKyurZEpNTa32jCKVZbFYePnmjtTxdGPzwZN89sNBsyOJiNRatarcNG7cmPT00s/jycjIwN3dnQYNGpS5jJeXFwEBAaUmEUfUpK4P44a0A2DaV7s4fOqsyYlERGqnWlVuYmJiWLVqVal5K1eupFu3bnh4eJiUSsR+7uzRnG7N65FbUMyzC7aX7JkUEZGKM7XcnD59mq1bt7J161bg3KXeW7duJSUlBTh3SOnuu+8uGT9y5EgOHjzImDFj2LlzJx9++CGzZ89m7NixZsQXsTur1cLUEdF4ulmJ332M+T8fIuHX4yzaepiEX49TrHNxREQuyWKY+E/D+Ph4+vfvf8H8e+65h48//ph7772XAwcOEB8fX/K7NWvW8OSTT5KUlERoaCjjxo1j5MiRFX7P7OxsAgMDycrK0iEqcVjvrd7LKyt2Y7HAH/8PDQn0ZuKwSK6LCjEvnIiICSrz/W1quTGDyo3UBksTj/Do3C0XzD9/2vyMO69QwRERl1KZ7+9adc6NiCsothm8vGxnmb87/y+RyUuSdYhKROQiVG5EHMyP+0+QlpV30d8bQFpWHj/uP1FzoUREahGVGxEHk5Fz8WJTlXEiIq5G5UbEwQT5e1dwnFc1JxERqZ1UbkQcTPfw+oQEelP2Pbd/9/rKX0g+omeliYj8mcqNiINxs1qYOCwS4IKCc/5nTzcrmw6e5IZ31vHCoh2cOlNQoxlFRByZyo2IA7ouKoQZd15B48DSh6gaB3oz884riP+/flwfHYLNgDkJB+n/ajyf/5iiK6hERNB9bsyOI1KuYpvBj/tPkJGTR5C/N93D6+Nm/X1/zsa9mUxaksQvR08D0LFJIJNv7MAVzeqZFVlEpFroJn7lULkRZ1NYbGNOwkHeXPULOflFANzatSnjrmtHI510LCJOQjfxE3EhHm5W7r86nO/G9uPWrk0B+PKnQ1zzajwfrt9PYbHN5IQiIjVLe25EnMzPKSeZuCiJ7YezAGgb7Mek4R3o1aqhyclERKpOh6XKoXIjrqDYZvD/Nqcy/etdnDxTCMD10SE8O7Q9oXV9TE4nIlJ5Oiwl4uLcrBb+3r0Zq8f24+6Y5lgtsCwxjQGvreG91XvJLyo2O6KISLXRnhsRF5B8JJuJi3ew6cBJAFo08OWFYZFc0y7Y5GQiIhWjw1LlULkRV2UYBou2HuFfy3eSkZMPwIB2QTx/QyQtGtYxOZ2ISPlUbsqhciOu7nR+Ee98u4fZ6/dTZDPwdLPyzz4tGdW/Fb6e7mbHExEpk8pNOVRuRM7Zm3GayUuSWLcnE4DQQG+evT6SoR0bY7Fc6slWIiI1S+WmHCo3Ir8zDIOVyUd5cUkyh0+dBSCmZQMm39iBtsH+JqcTEfmdyk05VG5ELpRXWMyM+F+ZueZX8otsuFkt3BPTgtED2xDg7WF2PBERlZvyqNyIXFzqiTO8tCyZFUlHAWjo58n4Ie25pUsTrFYdqhIR86jclEPlRuTS1vxyjMmLk9iXmQtAl2Z1eXF4FB2bBpqcTERclcpNOVRuRCqmoMjGRxv28/a3e8gtKMZigduubMb/DY6gfh1Ps+OJiIvRHYpF5LJ5ult5qG8rvhvbj5s6h2IY8PmPKfR/NZ5PEw5QbHOpfxeJSC2iPTciUiE/7j/BC4t2sCs9B4DIkAAm39iBK1vUNzmZiLgCHZYqh8qNSNUVFduY+2MKr67YTXZeEQA3d2nChCHtCArwNjmdiDgzlZtyqNyIXL7jp/N5deVu/rcpFcOAOp5uPHFtG+7tFY6nu452i4j9qdyUQ+VGxH4SD53ihUVJbE09BUDLRnWYNKwDfdo2MjeYiDgdlZtyqNyI2JfNZhD38yGmfb2LzNMFAAzuEMxz10cSVt/X5HQi4ixUbsqhciNSPbLOFvLWN3v45LcrqbzcrTzcrxUj+7bC28PN7HgiUsup3JRD5Uakeu1Oz2HS4iQS9h0HoGk9H164IZKBkcF6IKeIVJnKTTlUbkSqn2EYLNuexsvLdpKWlQdAn7aNmDgsklaN/ExOJyK1kcpNOVRuRGrOmYIi3lu9l3+v3U9BsQ0PNwv/uDqcx65pg5+Xu9nxRKQWUbkph8qNSM07kJnLi0uT+W5XBgDBAV48M7Q9wzuF6lCViFSIyk05VG5EzPPtzqO8uDSZg8fPANC9RX0mDe9AZKj+XxSR8qnclEPlRsRceYXF/GfdPt5dvZe8QhtWC9zZszlPDYwg0NfD7Hgi4qBUbsqhciPiGA6fOsu/lu1k2fY0AOrX8eT/Bkfw125huFl1qEpESlO5KYfKjYhj2bg3k4mLk9iTcRqA6KaBTB7egS7N6pmcTEQcicpNOVRuRBxPYbGNOQkHeXPVL+Tkn3sg51+6NuXp69rRyN/L5HQi4ggq8/2tJ9yJiOk83Kzcf3U4347ty61dmwLwxU+HuObVeD5cv5+iYpvJCUWkNtGeGxFxOD8dPMnExTvYcTgbgIhgfyYN70BMqwYmJxMRs+iwVDlUbkRqh2KbwbxNqbyyYhcnzxQCcEN0CM9e356QQB+T04lITVO5KYfKjUjtcupMAa+t/IX//nAQmwE+Hm48ek1rHugdjpe7Hsgp4ipUbsqhciNSOyUdyWLioiQ2HzwJQIsGvkwc1oH+7YJMTiYiNUHlphwqNyK1l2EYLNx6mH8t38WxnHwABrQL4oVhkTRvUMfkdCJSnVRuyqFyI1L75eQV8s53e89dSWUz8HSz8s8+LRnVvxW+nnogp4gzUrkph8qNiPPYm3GayUuSWLcnE4DQQG+evT6SoR0b64GcIk6mVt3n5v333yc8PBxvb2+6du3KunXrLjo2Pj4ei8VywbRr164aTCwijqJ1kB9z/tGdmXd2pUldH45k5fHI3J+54z8/sOdojtnxRMQkppabefPmMXr0aJ599lm2bNlC7969GTJkCCkpKeUut3v3btLS0kqmNm3a1FBiEXE0FouF66Ia882YvjwxoA2e7lY2/nqcIW+tY8rSZLLzCs2OKCI1zNTDUj169OCKK65gxowZJfPat2/PTTfdRGxs7AXj4+Pj6d+/PydPnqRu3bpVek8dlhJxbqknzjBlaTIrk48C0NDPi/FD2nFLlyZY9UBOkVqrVhyWKigo4KeffmLQoEGl5g8aNIiNGzeWu2yXLl0ICQlhwIABrF69ujpjikgtE1bfl1l3d+OTf3SnZcM6ZJ7OZ+wX27h15kZ2HM4yO56I1ADTyk1mZibFxcUEBweXmh8cHEx6enqZy4SEhDBr1izi4uKYP38+ERERDBgwgLVr1170ffLz88nOzi41iYjz69u2EV+P7sP4Ie3w9XTj55RTDHt3Pc8s2M7J3AKz44lINTL9msk/X9FgGMZFr3KIiIggIiKi5OeYmBhSU1N59dVX6dOnT5nLxMbGMnnyZPsFFpFaw9Pdysi+rbipcxNiv9rJoq1HmPtDCssS0xg7OILbuzfDTYeqRJyOaXtuGjZsiJub2wV7aTIyMi7Ym1Oenj17smfPnov+fsKECWRlZZVMqampVc4sIrVT40Bv3rqtC/P+2ZN2jf3JOlvI8wt3MOyd9Ww+cMLseCJiZ6aVG09PT7p27cqqVatKzV+1ahW9evWq8Ots2bKFkJCQi/7ey8uLgICAUpOIuKYeLRuw9LGrmTy8AwHe7iSnZXPrzASenLeVjOw8s+OJiJ2YelhqzJgx3HXXXXTr1o2YmBhmzZpFSkoKI0eOBM7tdTl8+DBz5swB4M0336RFixZ06NCBgoICPvvsM+Li4oiLizPzY4hILeLuZuWeXi24ITqEV1bsZt7mVBZsOczKpHSeuLYN9/YKx9Pd9FuAichlMLXc/O1vf+P48eO8+OKLpKWlERUVxfLly2nevDkAaWlppe55U1BQwNixYzl8+DA+Pj506NCBZcuWMXToULM+gojUUg38vJg6Ipq/d2/GC4uT2JZ6in8t38W8TalMGt6B3m0amR1RRKpIj18QEZdnsxl8+fMhpn21i+O/XUl1XYfGPHdDe5rW8zU5nYiAni1VLpUbEbmYrLOFvPnNL8xJOEixzcDL3cqofq15qG9LvD3czI4n4tJUbsqhciMil7IrPZtJi5P4ft+5K6ma1vPhhRsiGRgZrAdyiphE5aYcKjciUhGGYbA0MY2Xl+0k/bcrqfq0bcTEYZG0auRncjoR16NyUw6VGxGpjNz8It5bvZf/rNtPQbENDzcL91/dkseuaU0dL9PvgyriMlRuyqFyIyJVsT8zl8lLkojffQyA4AAvnhnanuGdQnWoSqQGqNyUQ+VGRKrKMAy+3ZnBi0uTSTlxBoDu4fWZPLwD7UP094lIdVK5KYfKjYhcrrzCYv69dh/vxe8lr9CG1QJ39WzOmIERBPp6mB1PxCmp3JRD5UZE7OXwqbO8vCyZ5dvPPSOvfh1Pnh4cwV+7hWHVAzlF7ErlphwqNyJibxv2ZjJxcRJ7M04D0KlpIJNvjKJzWF1zg4k4EZWbcqjciEh1KCy28cnGA7z5zR5O5xcB8NduTXn6unY09PMyOZ1I7VeZ7289HU5ExA483Kw80Lsl343ty4grmgLw/zYfov+r8Xy0YT9FxTaTE4q4Du25ERGpBj8dPMHExUnsOJwNQESwP5OGdyCmVQOTk4nUTjosVQ6VGxGpKcU2g/9tSuGVFbs5daYQgBuiQ3j2+vaEBPqYnE6kdtFhKRERB+BmtXBHj+asfqofd/ZshtUCSxPTuObVNbwfv5f8omKzI4o4Je25ERGpITsOZzFxcRI/HTwJQHjDOrxwQyT92wWZnEzE8emwVDlUbkTETIZhsGDLYWK/2sWxnHwArm0fxPM3RNK8QR2T04k4LpWbcqjciIgjyMkr5J3v9vLh+v0U2Qw83a081Kclo/q1xsfTzex4Ig5H5aYcKjci4kj2ZuQwaXEy6/dmAtCkrg/PXt+eIVGN9UBOkT9QuSmHyo2IOBrDMFiRlM6UpTs5fOosAFe1bsCkYR1oE+xvcjoRx6ByUw6VGxFxVGcLipmx5ldmrvmVgiIb7lYL9/ZqwRPXtsHfWw/kFNemclMOlRsRcXQpx88wZVkyq5KPAtDQz4sJQ9pxc5cmeiCnuCyVm3Ko3IhIbRG/O4PJS5LZn5kLwBXN6vLijVFENQk0OZlIzVO5KYfKjYjUJvlFxXy4/gDvfLeHMwXFWCzw9+7N+L9BEdSr42l2PJEaozsUi4g4CS93Nx7u14rvnurH8E6hGAbM/SGF/q/F89n3Bym2udS/T0UqRHtuRERqke/3HWfS4iR2pecA0CE0gMnDO9CtRX2Tk4lULx2WKofKjYjUdkXFNv77QwqvrdxNdl4RALd0acL4Ie0ICvA2OZ1I9dBhKRERJ+buZuWeXi1YPbYft10ZhsUC87cc5prX1vDvtfsoLLaZHVHEVNpzIyJSy21LPcULi5PYlnoKgNZBfkwa1oGr2zQ0N5iIHemwVDlUbkTEGdlsBl/+dIhpX+/ieG4BAEOiGvPs9e1pWs/X5HQil0/lphwqNyLizLLOFvLGql/49Lcrqbw9rIzq15p/9mmJt8fvD+Qsthn8uP8EGTl5BPl70z28Pm66QaA4MJWbcqjciIgr2JWezcRFSfyw/wQAYfV9eOGGDlzbPogVSelMXpJMWlZeyfiQQG8mDovkuqgQsyKLlEvlphwqNyLiKgzDYGliGi8v20l69rkiExkSQHJa9gVjz++zmXHnFSo44pB0tZSIiGCxWBjWKZRvn+rLw/1a4W6lzGIDcP5fuZOXJOvGgFLrqdyIiDi5Ol7ujLuuHdNv7VTuOANIy8pjwc+HyDpTiIvt2Bcn4m52ABERqRkVPWF47JeJQCKe7laC/L1o5O9FkL8XQf7e5/4bcO7PjX77c4M6XjoZWRyKyo2IiIsI8q/Y3Yt9PaycKbRRUGTj0MmzHDp5ttzxVgs08DtfgH4rQQFevxWjP/7ZCy93t3JfS8QeVG5ERFxE9/D6hAR6k56VR1kHnCxA40Bv1o+7hsJiG8dy8snIyedYTh4ZOflkZOeT8Yc/Hzudz/HT+dgMOJaTz7GcfJIukSHQx6PU3p+SPUMB3qX2Evl5uWOxaG+QVI3KjYiIi3CzWpg4LJKHP/sZC5QqOOdrxMRhkbhZLbhZ3Qir70tY/fJvAFhUbONEbsG5wpOT91sB+v3Px07/VoRy8ikotpF1tpCss4XsyThd7uv6eLiV7PH54yGwRn6/F6Egfy/q+Xpi1SEx+RNdCi4i4mK+3pFW4/e5MQyDrLOFF90DlJGdV7Kn6HR+UYVf191qKdnbc/4Q2LkCVPocoYZ+Xni46Rqa2kz3uSmHyo2IiGPfofhMQdFF9wBl5Pxegk789piJirBYoL6vZ8khsN8LUOlzhIL8vfHx1HlBjkjlphwqNyIizqGgyEbm6d9KUHbeHwpQ6fOEMk/nU1SJe/f4e7nTqIxDYKX2Bvl7E+Cj84JqUmW+v3XOjYiI1Eqe7lZC6/oQWten3HE2m8GJMwUX7Pk5VsZ5QnmFNnLyi8g5VsS+Y7mXfH9dKu+YVG5ERMSpWa0WGvqdO+8mkov/i98wDHLyi0pOgP5jEcrI/m1P0G9/zs4r0qXyDkzlRkREhHOPqwjw9iDA24PWQX7ljs0rLL7kpfIZOfkcz9Wl8mZQuREREakkbw/7XCp//vCYvS6VL/mzSZfKO8qJ6io3IiIi1cTdzXpur0uANxB40XHlXSp//vDYHy+VP1tYzMHjZzh4/Ez571/GpfK/7wGy76XyZtxi4GJ0tZSIiEgtUtal8hk5F5agy7lU/oITpS9xqfzXO9J4+LOfL7jz9fl9NjPuvOKyC06tuhT8/fff55VXXiEtLY0OHTrw5ptv0rt374uOX7NmDWPGjCEpKYnQ0FCefvppRo4cWeH3U7kRERFX8OdL5TNy/ngYrPRNFIsv41L5hn6efPnTIXLyyr754h8f63E5h6hqzaXg8+bNY/To0bz//vtcddVVfPDBBwwZMoTk5GSaNWt2wfj9+/czdOhQHnzwQT777DM2bNjAqFGjaNSoESNGjDDhE4iIiDimql4q/8fzgC7nUvnzDCAtK48f958gplUDO3yySzN1z02PHj244oormDFjRsm89u3bc9NNNxEbG3vB+HHjxrF48WJ27txZMm/kyJFs27aNhISECr2n9tyIiIhU3h8vlf/jIbCNvx7nu10Zl1z+rds6c2PnJlV+/1qx56agoICffvqJ8ePHl5o/aNAgNm7cWOYyCQkJDBo0qNS8wYMHM3v2bAoLC/Hw8Lhgmfz8fPLz80t+zs7OtkN6ERER13KxS+U7hAZWqNwE+XtXZ7xSTHuKWGZmJsXFxQQHB5eaHxwcTHp6epnLpKenlzm+qKiIzMzMMpeJjY0lMDCwZAoLC7PPBxARERG6h9cnJNCbi51NY+HcVVPdw+vXWCbTH5H655sQGYZR7o2Jyhpf1vzzJkyYQFZWVsmUmpp6mYlFRETkPDerhYnDIgEuKDjnf544LLJG73djWrlp2LAhbm5uF+ylycjIuGDvzHmNGzcuc7y7uzsNGpR9kpKXlxcBAQGlJhEREbGf66JCmHHnFTQOLH3oqXGgt10uA68s08658fT0pGvXrqxatYqbb765ZP6qVau48cYby1wmJiaGJUuWlJq3cuVKunXrVub5NiIiIlIzrosKYWBkY92heMyYMdx1111069aNmJgYZs2aRUpKSsl9ayZMmMDhw4eZM2cOcO7KqHfffZcxY8bw4IMPkpCQwOzZs/n888/N/BgiIiLCuUNUNXW5d3lMLTd/+9vfOH78OC+++CJpaWlERUWxfPlymjdvDkBaWhopKSkl48PDw1m+fDlPPvkk7733HqGhobz99tu6x42IiIiUMP0OxTVN97kRERGpfSrz/W361VIiIiIi9qRyIyIiIk5F5UZEREScisqNiIiIOBWVGxEREXEqKjciIiLiVFRuRERExKmYehM/M5y/rU92drbJSURERKSizn9vV+T2fC5XbnJycgAICwszOYmIiIhUVk5ODoGBgeWOcbk7FNtsNo4cOYK/vz8Wi30f5pWdnU1YWBipqam6+/ElaF1VnNZVxWldVY7WV8VpXVVcda0rwzDIyckhNDQUq7X8s2pcbs+N1WqladOm1foeAQEB2vgrSOuq4rSuKk7rqnK0vipO66riqmNdXWqPzXk6oVhEREScisqNiIiIOBWVGzvy8vJi4sSJeHl5mR3F4WldVZzWVcVpXVWO1lfFaV1VnCOsK5c7oVhEREScm/bciIiIiFNRuRERERGnonIjIiIiTkXlRkRERJyKyk0FrV27lmHDhhEaGorFYmHhwoWXXGbNmjV07doVb29vWrZsycyZM6s/qIOo7PqKj4/HYrFcMO3atatmApskNjaWK6+8En9/f4KCgrjpppvYvXv3JZdzxW2rKuvKVbcrgBkzZhAdHV1yI7WYmBi++uqrcpdxxe0KKr+uXHm7+qPY2FgsFgujR48ud5wZ25XKTQXl5ubSqVMn3n333QqN379/P0OHDqV3795s2bKFZ555hscff5y4uLhqTuoYKru+ztu9ezdpaWklU5s2baopoWNYs2YNjzzyCN9//z2rVq2iqKiIQYMGkZube9FlXHXbqsq6Os/VtiuApk2bMnXqVDZv3szmzZu55ppruPHGG0lKSipzvKtuV1D5dXWeK25X523atIlZs2YRHR1d7jjTtitDKg0wFixYUO6Yp59+2mjXrl2peQ899JDRs2fPakzmmCqyvlavXm0AxsmTJ2skk6PKyMgwAGPNmjUXHaNt65yKrCttV6XVq1fP+M9//lPm77RdlVbeunL17SonJ8do06aNsWrVKqNv377GE088cdGxZm1X2nNTTRISEhg0aFCpeYMHD2bz5s0UFhaalMrxdenShZCQEAYMGMDq1avNjlPjsrKyAKhfv/5Fx2jbOqci6+o8V9+uiouL+d///kdubi4xMTFljtF2dU5F1tV5rrpdPfLII1x//fVce+21lxxr1nblcg/OrCnp6ekEBweXmhccHExRURGZmZmEhISYlMwxhYSEMGvWLLp27Up+fj6ffvopAwYMID4+nj59+pgdr0YYhsGYMWO4+uqriYqKuug4bVsVX1euvl1t376dmJgY8vLy8PPzY8GCBURGRpY51tW3q8qsK1ferv73v//x888/s2nTpgqNN2u7UrmpRhaLpdTPxm83g/7zfIGIiAgiIiJKfo6JiSE1NZVXX33V6f+yOO/RRx8lMTGR9evXX3Ksq29bFV1Xrr5dRUREsHXrVk6dOkVcXBz33HMPa9asueiXtitvV5VZV666XaWmpvLEE0+wcuVKvL29K7ycGduVDktVk8aNG5Oenl5qXkZGBu7u7jRo0MCkVLVLz5492bNnj9kxasRjjz3G4sWLWb16NU2bNi13rKtvW5VZV2Vxpe3K09OT1q1b061bN2JjY+nUqRNvvfVWmWNdfbuqzLoqiytsVz/99BMZGRl07doVd3d33N3dWbNmDW+//Tbu7u4UFxdfsIxZ25X23FSTmJgYlixZUmreypUr6datGx4eHialql22bNni9LvCDcPgscceY8GCBcTHxxMeHn7JZVx126rKuiqLK2xXF2MYBvn5+WX+zlW3q4spb12VxRW2qwEDBrB9+/ZS8+677z7atWvHuHHjcHNzu2AZ07araj1d2Ynk5OQYW7ZsMbZs2WIAxuuvv25s2bLFOHjwoGEYhjF+/HjjrrvuKhm/b98+w9fX13jyySeN5ORkY/bs2YaHh4fx5ZdfmvURalRl19cbb7xhLFiwwPjll1+MHTt2GOPHjzcAIy4uzqyPUCMefvhhIzAw0IiPjzfS0tJKpjNnzpSM0bZ1TlXWlatuV4ZhGBMmTDDWrl1r7N+/30hMTDSeeeYZw2q1GitXrjQMQ9vVH1V2XbnydvVnf75aylG2K5WbCjp/6d+fp3vuuccwDMO45557jL59+5ZaJj4+3ujSpYvh6elptGjRwpgxY0bNBzdJZdfXtGnTjFatWhne3t5GvXr1jKuvvtpYtmyZOeFrUFnrCDA++uijkjHats6pyrpy1e3KMAzjH//4h9G8eXPD09PTaNSokTFgwICSL2vD0Hb1R5VdV668Xf3Zn8uNo2xXFsP47cweERERESegE4pFRETEqajciIiIiFNRuRERERGnonIjIiIiTkXlRkRERJyKyo2IiIg4FZUbERERcSoqNyIiIuJUVG5EpMbde++9WCyWC6a9e/de1uv269eP0aNH2yekiNRaenCmiJjiuuuu46OPPio1r1GjRialKa2goABPT0+zY4hIFWnPjYiYwsvLi8aNG5ea3nrrLTp27EidOnUICwtj1KhRnD59utRyGzZsoG/fvvj6+lKvXj0GDx7MyZMnuffee1mzZg1vvfVWyZ6gAwcOALBmzRq6d++Ol5cXISEhjB8/nqKiopLX7NevH48++ihjxoyhYcOGDBw4EIBJkybRrFkzvLy8CA0N5fHHH6+x9SMiVadyIyIOw2q18vbbb7Njxw4++eQTvvvuO55++umS32/dupUBAwbQoUMHEhISWL9+PcOGDaO4uJi33nqLmJgYHnzwQdLS0khLSyMsLIzDhw8zdOhQrrzySrZt28aMGTOYPXs2L730Uqn3/uSTT3B3d2fDhg188MEHfPnll7zxxht88MEH7Nmzh4ULF9KxY8eaXiUiUgV6cKaI1Lh7772Xzz77DG9v75J5Q4YM4Ysvvig17osvvuDhhx8mMzMTgNtvv52UlBTWr19f5uv269ePzp078+abb5bMe/bZZ4mLi2Pnzp1YLBYA3n//fcaNG0dWVhZWq5V+/fqRlZXFli1bSpZ7/fXX+eCDD9ixYwceHh72+ugiUgO050ZETNG/f3+2bt1aMr399tusXr2agQMH0qRJE/z9/bn77rs5fvw4ubm5wO97bipj586dxMTElBQbgKuuuorTp09z6NChknndunUrtdxf/vIXzp49S8uWLXnwwQdZsGBBqUNZIuK4VG5ExBR16tShdevWJVNBQQFDhw4lKiqKuLg4fvrpJ9577z0ACgsLAfDx8an0+xiGUarYnJ8HlJpfp06dUmPCwsLYvXs37733Hj4+PowaNYo+ffqUZBERx6VyIyIOYfPmzRQVFfHaa6/Rs2dP2rZty5EjR0qNiY6O5ttvv73oa3h6elJcXFxqXmRkJBs3buSPR+A3btyIv78/TZo0KTeTj48Pw4cP5+233yY+Pp6EhAS2b99ehU8nIjVJ5UZEHEKrVq0oKirinXfeYd++fXz66afMnDmz1JgJEyawadMmRo0aRWJiIrt27WLGjBkl5+S0aNGCH374gQMHDpCZmYnNZmPUqFGkpqby2GOPsWvXLhYtWsTEiRMZM2YMVuvF/wr8+OOPmT17Njt27CjJ4+PjQ/Pmzat1PYjI5VO5ERGH0LlzZ15//XWmTZtGVFQU//3vf4mNjS01pm3btqxcuZJt27bRvXt3YmJiWLRoEe7u527ZNXbsWNzc3IiMjKRRo0akpKTQpEkTli9fzo8//kinTp0YOXIk999/P88991y5eerWrcu///1vrrrqqpI9RkuWLKFBgwbVtg5ExD50tZSIiIg4Fe25EREREaeiciMiIiJOReVGREREnIrKjYiIiDgVlRsRERFxKio3IiIi4lRUbkRERMSpqNyIiIiIU1G5EREREaeiciMiIiJOReVGREREnIrKjYiIiDiV/w/hJ+xD25Sr3wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "# 标准化数据\n",
    "scaler = StandardScaler()\n",
    "X_scaled = scaler.fit_transform(X)\n",
    "# 计算协方差矩阵的特征值\n",
    "cov_matrix = np.cov(X_scaled, rowvar=False)\n",
    "eigenvalues, _ = np.linalg.eig(cov_matrix)\n",
    "# 绘制碎石图\n",
    "plt.plot(range(1, len(eigenvalues) + 1), eigenvalues, marker='o')\n",
    "plt.title('Scree Plot')\n",
    "plt.xlabel('Factors')\n",
    "plt.ylabel('Eigenvalues')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Factor 1  Factor 2\n",
      "0 -1.328255 -0.559732\n",
      "1 -1.337759 -0.000107\n",
      "2 -1.402517  0.307654\n",
      "3 -1.300187  0.719597\n",
      "4 -1.333781 -0.363828\n",
      "                   Factor 1  Factor 2\n",
      "sepal length (cm)  0.880960 -0.447287\n",
      "sepal width (cm)  -0.416916 -0.553900\n",
      "petal length (cm)  0.999189  0.019153\n",
      "petal width (cm)   0.962289  0.058402\n"
     ]
    }
   ],
   "source": [
    "# 初始化因子分析模型，假设我们想要2个因子\n",
    "fa = FactorAnalysis(n_components=2, random_state=42)\n",
    "\n",
    "# 拟合模型并转换数据\n",
    "X_fa = fa.fit_transform(X_scaled)\n",
    "\n",
    "# 将结果转换为DataFrame\n",
    "df_fa = pd.DataFrame(X_fa, columns=['Factor 1', 'Factor 2'])\n",
    "print(df_fa.head())\n",
    "\n",
    "# 获取因子载荷矩阵\n",
    "loadings = pd.DataFrame(fa.components_.T, columns=[\n",
    "                        'Factor 1', 'Factor 2'], index=feature_names)\n",
    "print(loadings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 手动实现因子分析"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_factors = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 标准化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.10252609, -0.92662384, -2.22574466,  0.19894046, -1.05119277],\n",
       "       [-1.41317631, -0.09082591,  1.33010341, -0.83892082,  0.91055478],\n",
       "       [ 1.17516225,  0.235568  , -0.17837558,  0.46713812, -1.23929943],\n",
       "       [ 1.63974719,  0.75432769,  1.12727895, -1.82278903,  1.5553271 ],\n",
       "       [ 1.28722509, -0.140075  ,  0.59776196, -1.49282457,  0.41577788]])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = (X - np.mean(X, axis=0)) / np.std(X, axis=0)  # 标准化\n",
    "X[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 计算相关系数矩阵"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.01010101, -0.03578495, -0.00420521,  0.16590033, -0.03816365],\n",
       "       [-0.03578495,  1.01010101, -0.03823689, -0.09396191,  0.10202764],\n",
       "       [-0.00420521, -0.03823689,  1.01010101, -0.04038821, -0.14398062],\n",
       "       [ 0.16590033, -0.09396191, -0.04038821,  1.01010101,  0.02946688],\n",
       "       [-0.03816365,  0.10202764, -0.14398062,  0.02946688,  1.01010101]])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = np.cov(X, rowvar=False)  # 计算协方差矩阵\n",
    "R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 计算特征值和特征向量"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.23575513, 1.18666926, 0.80051743, 0.95759901, 0.86996422])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigenvalues, eigenvectors = np.linalg.eig(R)  # 计算特征值和特征向量\n",
    "eigenvalues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "选择最大的K个特征值对应的特征向量"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = eigenvalues.argsort()[::-1]  # 对特征值排序\n",
    "eigenvalues = eigenvalues[idx]\n",
    "eigenvectors = eigenvectors[:, idx]\n",
    "eigenvals = eigenvalues[:n_factors]\n",
    "eigenvecs = eigenvectors[:, :n_factors]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 计算因子载荷矩阵"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.56019332,  0.38473066],\n",
       "       [ 0.57375641,  0.10026824],\n",
       "       [-0.28962755, -0.62764749],\n",
       "       [-0.5377147 ,  0.54816824],\n",
       "       [ 0.468744  ,  0.57807272]])"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loadings = eigenvecs * np.sqrt(eigenvals)\n",
    "loadings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 计算因子得分"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-8.93672723e-01,  1.03620674e+00],\n",
       "       [ 9.97141351e-01, -1.11331623e+00],\n",
       "       [-1.05490003e+00,  1.07328147e-01],\n",
       "       [ 7.25810935e-01, -8.52306419e-02],\n",
       "       [ 1.86241992e-02, -3.97720830e-01],\n",
       "       [-1.63226641e+00, -5.25392421e-02],\n",
       "       [-1.38984979e-02,  2.35260990e+00],\n",
       "       [-9.51363586e-01,  1.24067084e+00],\n",
       "       [ 8.18615714e-01, -9.57352641e-02],\n",
       "       [ 7.52754346e-01,  1.62719775e+00],\n",
       "       [ 5.74344185e-01, -5.92751160e-01],\n",
       "       [-3.02722167e-01, -1.70621370e+00],\n",
       "       [ 4.68333006e-04, -2.71870263e-01],\n",
       "       [ 5.31842869e-01,  1.47632728e-01],\n",
       "       [-1.95138815e-01,  7.58449711e-01],\n",
       "       [ 5.52345186e-01, -4.67294798e-01],\n",
       "       [-1.22297039e+00,  3.15066081e-01],\n",
       "       [ 7.80695127e-01, -6.74109211e-01],\n",
       "       [-1.02751564e+00,  5.97553845e-01],\n",
       "       [ 1.31269371e+00, -4.97933850e-01],\n",
       "       [-2.72819389e-01,  1.11725993e+00],\n",
       "       [-2.93710916e-01, -1.19680312e+00],\n",
       "       [ 8.50033433e-01,  4.65341156e-01],\n",
       "       [ 3.24066720e-01, -1.65753521e-01],\n",
       "       [-3.06656709e-01, -1.32057761e+00],\n",
       "       [ 1.17653516e+00, -7.96624194e-01],\n",
       "       [-1.28555765e+00,  1.28712099e+00],\n",
       "       [-6.03963911e-01,  1.54517906e+00],\n",
       "       [ 2.83020117e-01, -7.89237756e-01],\n",
       "       [-3.23709733e-01, -7.30214797e-01],\n",
       "       [ 1.14063477e+00, -2.22028628e+00],\n",
       "       [ 9.68038954e-01, -1.05782929e+00],\n",
       "       [ 1.50019802e+00, -2.85673418e-01],\n",
       "       [-9.21343224e-01,  2.23288438e-01],\n",
       "       [ 1.24907535e+00,  9.50658122e-01],\n",
       "       [ 1.35939112e+00,  1.67495009e+00],\n",
       "       [ 1.66748321e-01,  9.15424211e-02],\n",
       "       [-2.85140804e-01, -8.11082327e-01],\n",
       "       [ 3.99217070e-01, -1.01381388e+00],\n",
       "       [-1.41742315e+00,  8.49382762e-01],\n",
       "       [ 8.05223101e-01,  2.10656220e+00],\n",
       "       [-2.52009816e+00, -6.25333192e-01],\n",
       "       [ 8.94683435e-01,  4.93239114e-01],\n",
       "       [ 5.12374470e-01,  3.09221863e-01],\n",
       "       [ 1.05602138e+00, -2.53139289e-01],\n",
       "       [ 2.75758732e-01,  1.59844718e+00],\n",
       "       [-1.47075199e+00, -1.54445603e+00],\n",
       "       [ 1.11459121e+00,  1.81063371e-01],\n",
       "       [ 9.69387403e-01,  4.99449607e-01],\n",
       "       [ 4.97758741e-01,  8.82457459e-01],\n",
       "       [ 2.00622662e+00, -4.63414284e-02],\n",
       "       [-1.36646548e+00,  2.51810110e-01],\n",
       "       [ 1.30425986e+00, -1.49409865e-01],\n",
       "       [-1.43666872e+00, -6.86471674e-01],\n",
       "       [-6.71157496e-01, -1.34827901e+00],\n",
       "       [-3.20384822e-01,  1.41801592e-01],\n",
       "       [-1.29282990e+00, -1.45192299e+00],\n",
       "       [-1.52617551e-01, -2.90664552e-01],\n",
       "       [-9.03950447e-02,  1.34660946e+00],\n",
       "       [-2.50055379e-01, -1.87806842e+00],\n",
       "       [-3.96924614e-01, -7.64384538e-01],\n",
       "       [ 6.17021169e-01,  1.21275406e+00],\n",
       "       [ 6.82038478e-01, -1.42894594e+00],\n",
       "       [ 7.98657949e-01,  1.49030922e+00],\n",
       "       [ 6.29827960e-02,  9.93270200e-01],\n",
       "       [ 1.79283264e-01,  9.65802957e-01],\n",
       "       [ 7.69451717e-01,  3.77052864e-01],\n",
       "       [-1.09712854e+00, -1.11093779e+00],\n",
       "       [-1.14222320e+00,  5.19056931e-02],\n",
       "       [ 2.91141448e-01,  2.58356237e-01],\n",
       "       [-1.46956821e+00, -6.97928256e-01],\n",
       "       [-1.76927687e+00,  8.65717500e-01],\n",
       "       [ 2.37036254e-01, -5.93131597e-01],\n",
       "       [ 1.29585399e+00, -4.83274429e-01],\n",
       "       [-3.53961529e-01,  1.11778855e+00],\n",
       "       [ 7.98554521e-03, -9.07973959e-02],\n",
       "       [-5.47522972e-01, -3.81394418e-01],\n",
       "       [-1.23246193e+00,  8.75478651e-01],\n",
       "       [-1.20580061e+00, -2.25957045e-01],\n",
       "       [-1.72918111e+00,  5.70980705e-01],\n",
       "       [-7.65670877e-02,  1.02450393e+00],\n",
       "       [ 3.17874827e-01, -7.72420483e-01],\n",
       "       [-3.20325625e-01, -3.78949840e-01],\n",
       "       [-2.92331315e-01,  5.23742287e-01],\n",
       "       [-9.04401483e-02, -3.73625851e-01],\n",
       "       [ 1.47471891e+00, -1.49899833e+00],\n",
       "       [ 5.84408624e-01, -2.59492950e+00],\n",
       "       [ 9.22530858e-01,  1.97120848e-01],\n",
       "       [ 4.55303500e-01,  2.07026581e+00],\n",
       "       [-2.42524297e+00, -2.81019055e-01],\n",
       "       [ 9.08144167e-01,  2.07879276e-01],\n",
       "       [ 1.28335616e+00, -7.70829240e-02],\n",
       "       [-2.43240104e+00,  1.20950063e-01],\n",
       "       [ 1.19014990e+00, -1.24287921e+00],\n",
       "       [-6.83817262e-01,  6.00775009e-01],\n",
       "       [-7.09405555e-01, -3.79417657e-01],\n",
       "       [ 1.28116603e-01,  1.10592246e+00],\n",
       "       [ 5.60044052e-01, -8.79909558e-01],\n",
       "       [ 2.03580911e+00,  1.38991865e+00],\n",
       "       [-1.69710363e-01, -1.34591399e+00]])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F = X @ loadings @ np.linalg.inv(loadings.T @ loadings)\n",
    "F"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "DataAnalysis",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.20"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}