Skip to content

Instantly share code, notes, and snippets.

@lan2720
Created July 15, 2021 12:14
Show Gist options
  • Save lan2720/8b7fb498c61d08d175389440ff657257 to your computer and use it in GitHub Desktop.
Save lan2720/8b7fb498c61d08d175389440ff657257 to your computer and use it in GitHub Desktop.
Whitening.ipynb
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "d78d12f3-6463-4d6e-aaa8-6f2cdc5d0dcf",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "5128dd7c-dc4f-4b37-b299-3d57a0891f5e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"相关系数:\n",
" [[1. 0.97821977]\n",
" [0.97821977 1. ]]\n",
"特征值:\n",
" [450.18259359 3.1377546 ]\n",
"特征向量:\n",
" [[-0.44159531 -0.89721434]\n",
" [-0.89721434 0.44159531]]\n",
"(100, 2)\n"
]
},
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7f4dcb529e80>"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAIWCAYAAACLNPPAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABFO0lEQVR4nO3de3iU1b33/8/KJJDhIBFFgUGEthoPBAhEiyK0Sre0xUOM20MfW2utZbtrLW33QwXdVXp6oDs+Kj67tpdtPfSq/bVUMWqxRVvcFbXaBgOiIh5RGMCiEAphApNk/f6YmTBJ5jz3zNwz835dl5fknsnMymQy85m1vvd3GWutAAAA4JyKQg8AAACg1BCwAAAAHEbAAgAAcBgBCwAAwGEELAAAAIcRsAAAABxWWegBRDv66KPthAkTCj0MAACApNatW/eBtXZUrMtcFbAmTJig1tbWQg8DAAAgKWPMu/EuY4kQAADAYQQsAAAAhxGwAAAAHOaqGiwAAHBYMBjUtm3b1NnZWeihlLXq6mqNGzdOVVVVKX8PAQsAAJfatm2bhg8frgkTJsgYU+jhlCVrrT788ENt27ZNEydOTPn7WCIEAMClOjs7ddRRRxGuCsgYo6OOOirtWUQCFgAALka4KrxMfgcELAAAkJXPfvazam9vT3idm2++WX/605/yMyAXoAYLAABkxFora60ef/zxpNf93ve+l4cRuQczWAAAlIiWNr9mLlujiYtWaeayNWpp82d9m7fddpsmTZqkSZMm6Y477tCWLVtUW1urK6+8UpMmTdLWrVs1YcIEffDBB5Kk73//+6qtrdVZZ52lz33uc7r11lslSVdddZUefPBBSaGdW2655RZNmzZNdXV1eu2117Iep9swgwUAQAloafNr8cqNCgS7JUn+9oAWr9woSWqs92V0m+vWrdO9996rF154QdZaffzjH9cnPvEJvfHGG7r//vs1Y8aMPtf/+9//roceekgbNmxQMBjUtGnTNH369Ji3ffTRR+vFF1/UXXfdpVtvvVU///nPMxqjWzGDBQBACWhevbk3XEUEgt1qXr0549t85plndNFFF2no0KEaNmyYmpqatHbtWh1//PEDwpUkPfvss7rwwgtVXV2t4cOH6/zzz497201NTZKk6dOna8uWLRmP0a0cCVjGmBpjzIPGmNeMMZuMMWcYY0YaY540xrwR/v+RTtwXAAAYaHt7IK3j2Rg6dGjWtzF48GBJksfjUVdXV9a35zZOzWAtl/RHa+1JkqZI2iRpkaQ/W2tPkPTn8NcAACAHxtZ40zqeilmzZqmlpUUHDhxQR0eHHn74Yc2aNSvu9WfOnKnHHntMnZ2d2r9/v37/+99nfN/FLusaLGPMCEmzJV0lSdbaQ5IOGWMulPTJ8NXul/Q/km7I9v4AAMBAC+fW9qnBkiRvlUcL59ZmfJvTpk3TVVddpdNPP12SdM011+jII+MvSJ122mm64IILNHnyZB177LGqq6vTiBEjMr7/YmastdndgDFTJd0t6VWFZq/WSVogyW+trQlfx0jaE/k6noaGBtva2prVeAAAKBWbNm3SySefnPL1W9r8al69WdvbAxpb49XCubUZF7hnav/+/Ro2bJgOHDig2bNn6+6779a0adPyOoZciPW7MMass9Y2xLq+E2cRVkqaJul6a+0Lxpjl6rccaK21xpiYSc4YM1/SfEkaP368A8MBAKA8Ndb78h6o+ps/f75effVVdXZ26otf/GJJhKtMOBGwtknaZq19Ifz1gwoFrPeNMWOstTuMMWMk/SPWN1tr71ZoBkwNDQ3ZTacBAICC+vWvf13oIbhC1gHLWrvTGLPVGFNrrd0saY5Cy4WvSvqipGXh/z+S7X0BAOAmbliSgzs51Wj0ekkPGGMGSXpb0pcUOkNxhTHmy5LelXSpQ/cFAEDB5aKxJ0qHIwHLWrteUqwirzlO3D4AAG6TqLEnAQt0cgcAIAP5bOyJ4kPAAgAgA7lo7Ol2S5Ys6d28uZCGDRuW8PL29nbdddddad+ukz8fAQsAgAwsnFsrb5Wnz7FsG3uWK2utenp64n6drkwDlpMIWAAAZKCx3qelTXXy1XhlJPlqvFraVFfY+quXVki3T5KW1IT+/9KKrG/yhz/8oU488USdddZZ2rw5tHH0W2+9pU9/+tOaPn26Zs2apddee02S9P777+uiiy7SlClTNGXKFD333HOSpNtuu02TJk3SpEmTdMcdd0iStmzZotraWl155ZWaNGmS1q5d2+frrVu3qrm5WaeddpomT56sW265ZcDY9u/frzlz5mjatGmqq6vTI4+EGhYsWrRIb731lqZOnaqFCxdKUtzbivXzOcJa65r/pk+fbgEAQMirr76a+pU3/NbaHxxr7S1HHP7vB8eGjmeotbXVTpo0yXZ0dNi9e/faj370o7a5udmec8459vXXX7fWWvv888/bs88+21pr7aWXXmpvv/12a621XV1dtr29vfc29u/fb/ft22dPOeUU++KLL9p33nnHGmPsX//6V2utHfD16tWr7Ve+8hXb09Nju7u77bx58+xf/vIXa621Q4cOtdZaGwwG7d69e6211u7atct+9KMftT09Pfadd96xp556au/PEe+24v18scT6XUhqtXEyjVNtGgAAQCH9+XtSsF+BfTAQOj45s05Ja9eu1UUXXaQhQ4ZIki644AJ1dnbqueee0yWXXNJ7vYMHD0qS1qxZo1/+8peSJI/HoxEjRuiZZ57RRRddpKFDh0qSmpqatHbtWl1wwQU6/vjjNWPGjN7bif76iSee0BNPPKH6+npJodmqN954Q7Nnz+69vrVWN954o55++mlVVFTI7/fr/fffH/BzxLutffv2Dfj5nELAAgCgFOzdlt7xDPX09Kimpkbr16/P+rYioSvW19ZaLV68WP/2b/8W9/sfeOAB7dq1S+vWrVNVVZUmTJigzs7OAdeLd1uR5cpcoAYLAIBSMGJcesdTMHv2bLW0tCgQCGjfvn167LHHNGTIEE2cOFG/+93vJIXCy4YNGyRJc+bM0U9+8hNJUnd3t/bu3atZs2appaVFBw4cUEdHhx5++GHNmjUr6X3PnTtX99xzj/bv3y9J8vv9+sc/+u66t3fvXh1zzDGqqqrSU089pXfffVeSNHz4cO3bty/pbcX6+ZzCDBYAAKVgzs3SY1/vu0xY5Q0dz9C0adN02WWXacqUKTrmmGN02mmnSQrNHP37v/+7fvCDHygYDOryyy/XlClTtHz5cs2fP1+/+MUv5PF49JOf/ERnnHGGrrrqKp1++umSpGuuuUb19fXasmVLwvs+99xztWnTJp1xxhmSQq0ZfvWrX+mYY47pvc4VV1yh888/X3V1dWpoaNBJJ50kSTrqqKM0c+ZMTZo0SZ/5zGfU3Nwc87bi/XxOMKEaLXdoaGiwra2thR4GAACusGnTJp188smpf8NLK0I1V3u3hWau5tyccf0V+or1uzDGrLPWxtrJhhksAABKxuRLCVQuQQ0WAACAwwhYAAAADiNgAQAAOIyABQAA4DACFgAAgMMIWAAAICPDhg2LefynP/1p75Y59913n7Zv357PYbkCbRoAAICjrr322t5/33fffZo0aZLGjh1bwBHlHzNYAACUiFVvr9K5D56ryfdP1rkPnqtVb6/K6vaam5t15513SpK++c1v6pxzzpEU2tT5iiuukCTddNNNmjJlimbMmNG70fKSJUt066236sEHH1Rra6uuuOIKTZ06VYFAQOvWrdMnPvEJTZ8+XXPnztWOHTskSZ/85Cd1ww036PTTT9eJJ56otWvXZjX2QiNgAQBQAla9vUpLnluiHR07ZGW1o2OHljy3JKuQNWvWrN6g09raqv379ysYDGrt2rWaPXu2Ojo6NGPGDG3YsEGzZ8/Wz372sz7f/6//+q9qaGjQAw88oPXr16uyslLXX3+9HnzwQa1bt05XX321brrppt7rd3V16W9/+5vuuOMOffe738143G7AEiEAACVg+YvL1dnd2edYZ3enlr+4XPM+Mi+j25w+fbrWrVunf/7znxo8eLCmTZum1tZWrV27VnfeeacGDRqk8847r/e6Tz75ZMLb27x5s15++WX9y7/8i6TQhtBjxozpvbypqan3tpLtVeh2BCwAAErAzo6daR1PRVVVlSZOnKj77rtPZ555piZPnqynnnpKb775pk4++WRVVVXJGCNJ8ng86urqSnh71lqdeuqp+utf/xrz8sGDB6d8W27HEiEAACVg9NDRaR1P1axZs3Trrbdq9uzZmjVrln7605+qvr6+N1glM3z4cO3bt0+SVFtbq127dvUGrGAwqFdeeSWr8bkVAQsAgBKwYNoCVXuq+xyr9lRrwbQFWd3urFmztGPHDp1xxhk69thjVV1drVmzZqX8/VdddZWuvfZaTZ06Vd3d3XrwwQd1ww03aMqUKZo6daqee+65rMbnVsZaW+gx9GpoaLCtra2FHgYAAK6wadMmnXzyySlff9Xbq7T8xeXa2bFTo4eO1oJpCzKuv0JfsX4Xxph11tqGWNenBgsAgBIx7yPzCFQuwRIhAACAwwhYAAAADiNgAQDgYm6qlS5XmfwOCFgAALhUdXW1PvzwQ0JWAVlr9eGHH6q6ujr5laNQ5A4AgEuNGzdO27Zt065duwo9lLJWXV2tcePGpfU9BCwAAFwq0kkdxYeABQAoWS1tfjWv3qzt7QGNrfFq4dxaNdb7Cj0slAECFgCgJLW0+bV45UYFgt2SJH97QItXbpQkQhZyjiJ3AEBJal69uTdcRQSC3WpevblAI0I5YQYLAFCStrcH0joeC0uMyBQzWACAkjS2xpvW8f4iS4z+9oCsDi8xtrT5HRwlShUBCwBQkhbOrZW3ytPnmLfKo4Vza1P6fpYYkQ2WCAEAJSmylJfpEp8TS4woXwQsAEDJaqz3ZVwzNbbGK3+MMJXqEiPKG0uEAADEkO0SI8obM1gAAMSQ7RIjyhsBCwCAOLJZYkR5Y4kQAADAYcxgAQBchwafKHYELACAq7CHIEoBS4QAAFehwSdKAQELAOAqNPhEKSBgAQBcJds9BAE3IGABAFyFBp8oBRS5AwBchQafKAUELACA69DgE8WOJUIAAACHEbAAAAAcRsACAABwGAELAADAYQQsAAAAhxGwAAAAHEbAAgAAcBgBCwAAwGEELAAAAIfRyR0AUPJa2vxsvYO8ImABAEpaS5tfi1duVCDYLUnytwe0eOVGSSJkIWdYIgQAlLTm1Zt7w1VEINit5tWbCzQilAMCFgCgpG1vD6R1HHACAQsAUNLG1njTOg44gYAFAChpC+fWylvl6XPMW+XRwrm1BRoRygFF7gCAkhYpZOcsQuQTAQsAUPIa630EKuQVAQsAUHLoe4VCI2ABAEoKfa/gBhS5AwBKCn2v4AYELABASaHvFdyAgAUAKCn0vYIbUIMFAChasYrZF86t7VODJdH3CvnHDBYAoChFitn97QFZ9S1mX9pUJ1+NV0aSr8arpU11FLgjr5jBAgAUpUTF7M8uOodAhYJiBgsAUJQoZoebEbAAAEWJYna4GQELAFCU2MQZbkYNFgCgKLGJM9yMgAUAKAr9WzKcfdIoPfXaLsIVXMmxJUJjjMcY02aM+X3464nGmBeMMW8aY35rjBnk1H0BAMpLrJYMv3r+vQEtGlra/IUeKiDJ2RqsBZI2RX39I0m3W2s/JmmPpC87eF8AgDISqyVDf+w3CDdxJGAZY8ZJmifp5+GvjaRzJD0Yvsr9khqduC8AQPlJtfUCLRrgFk7NYN0h6duSesJfHyWp3VrbFf56myQWxgEAGUm19QItGuAWWQcsY8x5kv5hrV2X4ffPN8a0GmNad+3ale1wAAAlKFZLhv5o0QA3cWIGa6akC4wxWyT9RqGlweWSaowxkbMUx0mKWXlorb3bWttgrW0YNWqUA8MBAJSaxnrfgP0FPz9jPPsNwrWMtda5GzPmk5L+t7X2PGPM7yQ9ZK39jTHmp5Jestbelej7GxoabGtrq2PjAQAAyBVjzDprbUOsy3LZyf0GSd8yxrypUE3WL3J4XwAAAK7haKNRa+3/SPqf8L/flnS6k7cPAABQDNiLEAAAwGEELAAAAIcRsAAAABxGwAIAAHAYAQsAAMBhBCwAAACHEbAAAAAcRsACAABwGAELAADAYY52cgcAlJeWNr+aV2/W9vaAxtZ4tXBuLRsuAyJgAQAy1NLm1+KVGxUIdkuS/O0BLV65UZIIWSh7LBECADLSvHpzb7iKCAS71bx6c4FGBLgHAQsAkJHt7YG0jgPlhCVCAEBGxtZ45Y8RpsbWeFO+DWq4UKqYwQIAZGTh3Fp5qzx9jnmrPFo4tzal74/UcPnbA7I6XMPV0ubPwWiB/CJgAQAy0ljv09KmOvlqvDKSfDVeLW2qS3kGihoulDKWCAEAGWus92W8pEcNF0oZM1gAgIKIV6uVTg0X4FYELABAQWRbwwW4GUuEAIC8ij5zcIS3StVVFWo/EOQsQpQUAhYAIG/6d39vDwTlrfLo9sumEqxQUlgiBADkDWcOolwQsAAAecOZgygXBCwAQN5w5iDKBQELAJA3nDmIckGROwAgbyKF7Ow/iFJHwAIA5FWq3d/ZCBrFjIAFAHCd/u0cIhtBSyJkoShQgwUAcB3aOaDYEbAAAK5DOwcUOwIWAMB1aOeAYkfAAgAUTEubXzOXrdHERas0c9katbT5JdHOAcWPIncAQEGkUsjOWYQoVgQsAEDKnGydkKiQPdLKgUCFYkXAAgCkxOnWCRSyo5RRgwUASInTrRMoZEcpI2ABAFLi9IwThewoZQQsAEBKnJ5xaqz3aWlTnXw1XhlJvhqvljbVUXeFkkANFgAgJQvn1vapwZKyn3GikB2lioAFAEgJrROA1BGwAAApy8WMk5OtHwC3IGABAArG6dYPgFtQ5A4AKBinWz8AbkHAAgAUDM1GUaoIWACAgqHZKEoVAQsAUDA0G0WposgdAFAwtH5AqSJgAQAK2iqBZqMoRQQsAChz+W6VQN8rlANqsACgzOWzVUIkzPnbA7I6HOZa2vyO3xdQSAQsAChz+WyVQN8rlAsCFgCUuXy2SqDvFcoFAQsAylw+WyXQ9wrlgoAFAGWusd6npU118tV4ZST5arxa2lSXk8Jz+l6hXHAWIQAgrVYJ2ZwFSN8rlAsCFgAgZU60dKDvFcoBS4QAgJRxFiCQGgIWACBlnAUIpIaABQBIGWcBAqkhYAEAUsZZgEBqKHIHAKSMswCB1BCwAABp4SxAIDmWCAEAABxGwAIAAHAYAQsAAMBhBCwAAACHEbAAAAAcRsACAABwGAELAADAYQQsAAAAhxGwAAAAHEbAAgAAcBhb5QAA1NLmV/PqzfK3B+QxRt3Wysc+g0DGCFgAUOZa2vxavHKjAsFuSVK3tZIkf3tAi1dulCRCFpAmAhYAlLnm1Zt7w1V/gWC3mldv7r3e9vaAxjKzBSRFwAKAEhVZ9ksWira3BxLeTmQmKxLCmNkCkqPIHQBKUGTZz98ekNXhUNTS5h9w3bE13oS35TFmwAxX9MwWgIEIWABQgmIt+8ULRQvn1spb5Yl5O94qT29NVn/JZr6AckbAAoASFC/8xDreWO/T0qY6+cIzWR5jJEm+Gm+f4/0lm/kCyhk1WABQgsbWeOWPEabihaLGel9vPVV07Vbz6s06+6RRemidv8+MmLfKo4Vza3MzeKAEMIMFACXo7JNGyfQ7liwUtbT5NfW7T+gbv13fp3broXV+XTzdJ1+NV0aHZ7YocAfiYwYLAEpMS5tfD63zK7pyyki6eLpvQCiKbjBqJMWqtgoEu/XUa7v07KJzcjhqoLRkPYNljDnOGPOUMeZVY8wrxpgF4eMjjTFPGmPeCP//yOyHCwBIJlaBu5X0q+ff08xla3rPJIw+0zBynXgoaAfS48QSYZek/7DWniJphqTrjDGnSFok6c/W2hMk/Tn8NQAgxxKFoeh2DYkajPZHQTuQnqwDlrV2h7X2xfC/90naJMkn6UJJ94evdr+kxmzvCwCQXLIwFGnXkOqslFEomEXPfgFIzNEid2PMBEn1kl6QdKy1dkf4op2SjnXyvgAAsSXqaxUR6e6eisjSYaJmpQD6cixgGWOGSXpI0jestf+MvsxaaxVned8YM98Y02qMad21a5dTwwGAstW/r1Uska1z+gexyJmHvhqvjhxSNeD76OAOpMaRgGWMqVIoXD1grV0ZPvy+MWZM+PIxkv4R63uttXdbaxustQ2jRo1yYjgAUPYa6316dtE5uuOyqQNCVKRdQ3QQi7RfuP2yqdqybJ6eXXSO2g8EY942Be9Aclm3aTDGGEm/kLTJWntb1EWPSvqipGXh/z+S7X0BANITacsQb9Pn6Aaj/aXbrBTAYU70wZop6QuSNhpj1oeP3ahQsFphjPmypHclXerAfQEAYojuvp5OiEpk4dxaLV65kQ7uQAayDljW2mekAQ2DI+Zke/sAgMQi/awiQShSjC4po2AVHdZqhlRpcGWF9gaCA4IbgPjo5A4ARS5WP6tIMXq6Yah/WNtzIChvlUe3XzaVYAWkgYAFAEUuXtF5OsXo0Vvm9JdpWAPKGZs9A0CRi1d0nmoxekubXwsf3BAzXEUkugzAQAQsAChysfpZpVOM/t3HXlGwO9FOhKFCWxqMAqkjYAFACRhcefjl/MghVVraVJfykt6eOP2uollJ3/jterbLAVJEDRYAFLH+RemS1Bnsydn9ZXuGIlAumMECgAJqafNr5rI1mrhoVUazQ4nOIExVjXfgljiJsF0OkBwBCwAKJDL75G8PyCqzzZSdOINwyQWnqqoiXjvD9O4XQAgBCwAKxInZp1TOIEw2S9ZY71PzJVP67En4+Rnjk24WDSA+arAAoECcmH2KtZ2NJB041NUbpFLp8h5vO51YNV5slwMkR8ACgAJxYjPlSCha8ugrag8cPhtwz4GgFq/cqMGVFTFnyZY8+kpKRerJNosGEBsBCwAKxKnNlBvrfWpevblPwJJCQap/uIpoDwTV0uZPOWQRqID0UIMFAAXSWO/T0qa6PrVP6fSvipZJ0Tl9rYDcYQYLAArIqdmheMuNydDXCsgNZrAAoARkU3ROXyvAeQQsACgBjfU+HTkkdsPQI4dUDdirsD/6WgHOImABQIm45fxTY276PG/yGFVXJX65p68V4CxqsACgyLW0+XvbKIzwVqm6qkLtB4IaW+PV2SeN0kPr/HHPJpToawXkAgELAIpY/0ag7YGgvFUe3X7ZVDXW+zRz2ZqY4cpjjHqspa8VkCMELAAoYom222ms98WtreqxVu8sm5ePIQJliRosAChiybbbSWWvQgDOI2ABQBGLF5SspJnL1ujsk0bFLHyn5grILQIWALhUS5tfM5et0cRFq+J2XF84tzZuCwZ/e0APrfPr4uk+R7rFA0gdNVgA4EL9i9fjdVyP3ow5Vif3QLBbT722S88uOicPowYQwQwWALhQouL1/hrrfXp20TkycW6LJqJA/hGwAMCFkhWvx0JBO+AeBCwAcKFMwlKseiwK2oHCIGABgAtlEpYa631a2lRHQTvgAhS5A4ALRRevb28PpNxxvbHeR6ACXICABQAulW5Yit6TkC1wgMIiYAFACUi1rQOA/KAGCwBKQLy2Dv+xYkPMBqUAcouABQAlIF77hm5r9c3frteEBN3gATjPWGsLPYZeDQ0NtrW1tdDDAICiEF1zZYzUk8LLuVFon0IfNVpA1owx66y1DbEuowYLAIpQ/5qrVD8rR65GjRaQWywRAkARilVzla54W+8AyB4BCwCKkFP7C7JPIZAbBCwAKEKp7i8YbwPodG8HQHoIWABQhGJtpRPNW+XRHZdN1e2XTZUvHKL6hy32KQRyhyJ3AChC/bfSqRlSJWulvYHggC7ukf/T6R3IHwIWABSZ/kHp9sumphSU2KcQyB8CFgAUEbbEAYoDAQsAikiiLXGir8MyIFBYBCwAcIlUaqQSbYmz8HcbJCMFu0PtRJndAgqHswgBwAUiS3/+9oCsDoej/nsHJmqrEOyxveEqgmaiQGEQsADABeIt/fUPR8naM8RCM1Eg/whYAOAC8UJQ/+ON9T4tbaqTxyRrIXoYzUSB/CNgAYALxAtBsY431vv0fy+dMmAmq6rCqMrTN3jRTBQoDIrcAaAA+he0n33SKD20zt9nmTBROOrfaDRSFB/rGAXuQP4Za23ya+VJQ0ODbW1tLfQwACCn+veykkJh6uLpPj312i7CEVAkjDHrrLUNsS5jBgsA8ixeQftTr+3Ss4vOKdCoADiJGqxS9tIK6fZJ0pKa0P9fWlHoEQFQ6gXtAIoXAatUvbRCeuzr0t6tkmzo/499nZCF/CHgx5VOQTuA4kTAKlV//p4U7PdpOBgIHQdyjYCfUHQvqwsqntEzg76utwdfoSfNV3mMgBJBDVap2rstveOAkxIF/MmXFmZMLhIpXF+/6m59O/hzDTGHJElDAjtCQVTicQKKHDNYpWrEuPSOA04i4CfVWO/TkqEP9YarXsGAdq68URMXrdLMZWsGbJUDoDgQsErVnJulqn71HFXe0PFSRc2Pe5RhwG9p82vmsjXpBaM4gfMY+0HC/QgBuB8Bq1RNvlQ6/05pxHGSTOj/599ZussO1Py4S5kF/FQ3au4fwvZUHRPz9rbbo3r/zWbNQHEiYJWyyZdK33xZWtIe+n+phiuJon63KbOAn8pGzbFC2JKOi3XADurzfQfsIP1XV9/HifYNQPGhyB2lgZof95l8ackGqv5S6WsVK4Q90nOWbFD6duUKjTUfars9Sv/Vdake7Tmrz/Vo3wAUHwIWSsOIceHlwRjH3eSlFaFZtb3bQmObc3PZhJBSNrbGK3+MkBUdjOKFsEd7ztKjh86KeZnEZs38zaBYsUSI0lAMNT/UiRWvJCdQRPe1iugfjNKdhTKSfDVeLW2qK9/9CPmbQREjYLkRZ8OlrxhqfgpZJ8ZzKnMpvMk31vu0tKlOvhpv3GAUK4SZOHfpq/HqnWXz9Oyic4o7XGX7vKO2EkWMJUK3ibyYR15UIi/mkrvCghvluuYn26WKQtWJlfhz6pVXXtHNN98sn8+nj33sYzr++OM1ZswYjRkzRqNHj1ZVVVV2d5Bi09TGel/CMBS5rHn1Zm1vD2hsjVdnnzRKD63z96nNKpklQSeed9RWoogRsNyGDtju5MSbRaHqxOI9px6+Vlo5v+jrWsaPH68hQ4boD3/4g2pqaiRJlZWVOnTokC644ALdcsst2d2Bg2/ysUJYw/Ej+4SuhXNr05u1cmuNkhOvZbn+m3HrY4eSwBKh2/CJzZ2cWKooVJ1YvOeO7VZKdS35Wl7M8H6GDx+uu+66S2eeeaYk6dhjj9XRRx8tY4ystero6MhuXGk0TV319iqd9es5qruvTg0/P10PNNfKJvl5Gut9enbROZktCbq5Rinua9nW1H+/ufybcfNjh5JgrLWFHkOvhoYG29raWuhhFNbtk+J8Yjsu1MsKhbGkRlKsvxUT6jOWqkJ8Yo73nOov1nOs/8ydFHqDm/K/pDeecO7niHc/adTR7d+/X9dff71efvll1dTUaNeuXTLGaPjw4br66qt18cUXa/jw4Tkb26q3V+k7z9yioD3Ye6y6p0dLPtiteR0HtOqIGi0ffZx2Bv+p0UNHa/a42Xp629Pa2bFTo4eO1oJpCzTvI/MkhXpmpTSr5ebXi2TPu1R/v7n6m3HzY4eiYYxZZ61tiHkZActlHHijQQ4U84txrOdUTDHCYtw3SaM+gTPb56hDj29HR4euv/56/fWvf9WsWbO0fv16dXZ2yuPxaNiwYbrqqqt0ySWX6IgjjkhvfMne5F9aoXNbv6cdnoFl62OCXVqwp11Ljh6pzor4iwbVnmotOXOJgnunavHKjQPqsmKeTehU8O8v21Dz0grpDzdIgd2Jr1fIv59cPXYoKwSsiGJZby+WcZaTYg++0c8pUxFeHuzHO1K64Z2+x+K+CcWQzZulg292HR0deuKJJ3Teeedp586duueee7Rq1Sp1dnaqoqJCQ4cO1Re+8AVdPmmwav7+f7P/Ows/Nyb7jpY1Mc4LtFYjunu0t9Iz8LJ+xgwdo443F8XsqeWr8erZRef0vd+Hr439u8zmd5Htcz3lQC8VNMwU84cmuAYBSyr+N0gUXqkE35dWSI9cJ3Uf6nu8okpqvKvvz5Tq8qKkrN4sc/xm5/f7de+99+qxxx5TZ2enzMG9GnJgu66YVKHP1Q3SSK/J+PXgwI9O0pDADp07bqx2VMU5b8haKVb46sfIaN+mpfGipt5ZFlpCTBhiCj2bmM5zJtsgmO0sW6rvCaXytw/HJQpY5VPkTj8VZKtU9nacfKk0aNjA4z3BgX8PsYqM43VvMhWZF8Lnqpg5XDjv+9mp+s+hv1HL/7lan/vc51Qd+Ic6DnXrZy8GNe/XB3TnCwf1wd6OgT9/ksL7lja/qg/slCQt2NOu6p6e2ONIIVxJkpXV8BN+pMoj2gZc1qdRaazXM0kynuw/NGZ7ok2q18vm9+tEgXqqvfMohkeGyidgcXYeSlm6Z+AF9sQ+3v/vIdabUMPVMUKXEp+VmGx8uWgUG+ONcczzt2jRZz+iRy6p1BcmV6m6UjoQtLqnLajzf31At/3xLe3atSvu9/f/uZpXb9Z2e1Tv19XWhmarUhHnerZyj6rHrOwTsgb0xop7ZmhP9sE/jbMm07qed6Rzv99MPzD3fx5KyT808eE8N8qg+XH5LBGy3p5/+ZhWT3Qf5TKtn8nyd7Z/D+nUdBVqeT7RzyhJe7dqV0ePfvVSUCteDSoQlIKqVHDUqfrKlZ/TlYFf6NieHbG/P/wYTVy0SudXPKNP1vxKS0cdkbCIvVe/5cK4q4fdQ7T/9Ztjn0WYy9ezXNRgOf37zqRmL91x9T7H4y13UgyfsRIq2WGJUCqOveqKSbJPH/mYVk90H26e1nf6k9sfbkj9E3bkvvdu1YClvnT+HqKXS22cZbHA7sNvUoWYAUg0ax1+PRg1tELfPGOwfv+5Ifp8vVe7K4/We/9o149/do8uuPtNLV17UDv29cS93bE1Xj3ac5b+z8gxqYUraUCairt66Dmg/56v2L2xkr2eZfMcy3Y2MR/bViWaZYv3s6fzPOzz+pHmGJBcrl8TXDI7Vj6d3CN/3OUwo5FrqXQ1z0dH+mR/pG7siO/0tjUvrYh/Knz/gDHgU6NVb7uFEcdl/vcQr9u2dPjvLZXxOS1RF/Dwz7lz5Y06xn6gQPUxeq/+Up1/UkDHv/YrrXq5XQe6rH7zSlAPvxbU+SdW6ur6QfIdUdHnjXXh3FotXrlRnZWBuPsKZmP5i8t7e2P1kej1zInnWLbbTuV626o5N8eeATnh3Pg/ezrPw3g1btH3lezDSLnMoGcil68JLtoarHwClpT7P/pykSw8vbQi/huuk2+qmfyRFrrmzungmegTXyQIJFzqsJkvKyVdQtHhN5dCbBEU70048sY4+VKd8euhvQtNF1Q8ozuOeEBDPh7U16YM0W9fCeqXG4LqCFr97tUuPfZ6lz5T69U1i76q8CKjGut92rBnjR58r19fMIfs7NgZ/8J4r2eJtkaKfF+xixcwE/19pfM8TPQ6kcqJBC56k3elXL4muGi7ufJZIkRf2UyhJgo2kReWeJx8U020TJBtoW6uOP3JLdH3zbk5taWOTO47lduVDr/xxToTMZ0tUzKRwlJV9Jl5365coSEm1LpiRLXR/OmD9PgVQ3T96YNUU20UtB49vO0oXfyf9+k73/mOtmzZolVvr9Lvt90umdzUso4YPCL9b0q0NVKsZfJ0XgtcsvQiKfZZvSksC/cRbyYq0etEKicSUBifWC5Ldlx0Qlt5zWC5Vb6nklP9dBVvXIk+fSSaWne65i3ZDEWiywrF6U9u8W7POzL0u/rRxOQNHzO572RLKNLhx7vPbEOk9iscSFL9ZJ/p30iSWevIEt+/dP9FPvPBgMuPGGx0df0gXXJKlR7aMlz3bp+gffv26ZFHHtEf//hHfTDuAw06c5AGHTOo7zdGTh5KsT1DPO2de1R33yRVSOoxRmOGjumzpU5MiZZsgwHpsW8cblBqKiSZwycpJPp9FMOsTArLwik9j+bcHNoIPdasZCp/Ly56k3elXJbsFGrGPAZmsAqtEMXYqXy6SjSuRJ8+Er2AOF3ommiGIh+Ftplw+pNbvNv7zI8S12dle9/J3ii8I6VKb+hNKvp0+BHHacCbVrJP9sn+RrKYVWms9+mXp72rHw36RcIsNHyw0VW1HXr88cf1rW99SyNHjlR3d7d2tO7Qe3e9px2/3aGDOw/vQShjsg5X0bfTE76tHR07tOS5JVr19qr43xNzxjBKsONwoLI9A88Ajff7KIZZmWR/X6n2spt8aagdSaYngrh1Bt1NctVX0EUntDGDVWiFWC9O5dNVonFF6nViffqIV5Mz4rjc/DyJZihSqbnL9+xhLj65VXoP/668I0PhavKlh4NNPLkqbPeOlLoCmRca9/+dHOpI/Mae5azKaW/9P0kHk15PI8Zp6NCh+vznP6+mpiY98sgjuu6H1ymwP6D9L+3XgU0H5D3Rq6POPkqDxwxO6b4z0dndGb/4XTr8c8fbRicVsX5PxTAr4+Tf13m3SeNnZHZbyWbXkTsuOqGNgFVohXjRSmUKNdm44oWXYnphKdSSh1MnW8TqJdMV9e9Ez6Gmn8Xv/RO9SW90YIs252Zp5Vdi33asWbNUC41j/U7i2bvNmQ8oqfyt9XsODxkyJNQRfmq1/vePv6X3n9ql7o5udbzSoQOvH9CQjw7RyHNGqtpXndoY0pSw+F06/LOnvCdgP6Yi9LuIfgxdtPSSkJMnM2V6Wy56ky9LLjmhjSXCQivEVHIqU6iZjsutS3OxFMOSRyLJxp+oo3a8cPXIdX0DUmC31PLV2J3Xq4amN95UCo1Tqe2KGDHOmQ8o8R4n41Gy5/BFJ1+kL33+yzr+m8dr1PmjVFlTKVmp47UObfvZNm3/5XZ1bu1MfSwpGj10dPIrxfpbNCm+5McqiHfR0ktRKJWttZAxAlahFeJFK5UQlM24iuWFpRiWPBJJNv5EtTixtrFZ+ZWBG0BLsfcolKTKOMtg8d7EI4XGiZ576e5j5z0y9uXxjscS77l+0U8HPodj1Hs9ve1pVVRVaETDCB3/zeN1zIXHqPLIcNDa3KFtv9gm/31+Bd4LyImdM6o91VowbUFqV+7/tzj9S6nfUf8PG8X04QlwgZwvERpjPi1puSSPpJ9ba5fl+j6LSqGmkpNNoZbDFHexLHnEk2z8kd9V9JKfFPp3dCuNVJaRYgWfePsZ2p5QQEnQfyru88h7ZPLCfOnwG/sfbkh+3WRSfa7HWVLeOW5U71UqKit0xLQjNGzyMO1/eb92r9mtrr1dOvDGAQW2BOQd7w0tHR5fLZNKEXw4kKV1FmEi590mffim9M5fUrt+rL0pS+k1AMihnAYsY4xH0o8l/YukbZL+box51Fr7ai7vt+i49UXLreNKV7xC9lTqxVIpgi9Ux+ZUxh858aB/aEnU8T6WWKEzbsA77nDTx3Qfk+5Uis2jTpiIF/LiHY8nled6nCXZ0d1WOzx9w1JFZYWOmHqEhtcN1/5X9+vDP32orvYuHXjrgALvBVTtq9bIOSPlnegdGLTCoaqmp0eLPtyjeZVHObtf6u63U79usXzYAFwo1zNYp0t601r7tiQZY34j6UJJBCzkR6rb+sTbLDrZ9xayN1CqMy/ZLoUaT+yl4UQBL5Nw/tKK0BmDicSqFczXLGScx2vBh7u1ZIxPnd0Da62Mx2h43XANO2WY9m/ar91/3q3gnqACWwLafu92DR47WOPmj5PpF9A2bon+mTIoUk8k3WVYABnJdcDySYp+pdgm6ePRVzDGzJc0X5LGjx+f4+Gg7CQ7yyxREEjlDLVCb8uQSpBJFkKSdWOvHhH7PpxeRk52ckH/thLxAlmugkGcx3Fe5UjpzCVa/uJy7ejYoQpToZ5+G2Abj9HwSaGg1fFahz588kN1bOrQoQ8OKfBuQEM+MqT3ugMq2JwOi4ma0w4aWrolAUCeFbxNg7X2bkl3S1JDQ0Nu9ptA+cpm9iaV7y2GQvl0O973l2i5zcll5HTaSsRqUSHFbyvhhASP47yPzOtTF1V3f13MmzAVRsNOGabuQLc6t3VKlZLH6zl8BWt1yT/3Dbj9vPwcuXrcgDKV67MI/VLvvqiSNC58DMiPbNpgpPK9xdCxOeWO93Gkc0ZeNuLdT9XQgW/88do5DIpxXaekcBZdS5tfM5etkbXxC9h7unq0+6ndqhhUoeGThh9uSmqlse3H65rdVbK5PEuPswGBvMj1DNbfJZ1gjJmoULC6XNL/yvF9Aodl0/g0le8tlsaqqXS8/9HE1M7gy7dY7SAKNXOY4HFsafNr8cqNCgS7NXjwx1V15PMxd8vZ//J+de/rlqkyGnn2yMMXGMmOC2rdBQm6tDulVE5gAVwspzNY1touSV+TtFrSJkkrrLWv5PI+gT6y+bSeyveW0myAU2fk5eP+XThz2Lx6swLB0NY0B99vVHDPDFl7eN9nSbJdVrv/vFu2y2rIR4eoemzfTu8p7TUIoCjkvAbLWvu4pMdzfT9AXNl8Wk/le0tlNqDQfcHSuf8czhy2tPnVvHqztrcHNLbGq4Vza9VY70v6fdvb+y5ZHny/UQffb1TlEW0aOuZB9VR0a98r+9S1r2vg7FWUpHsNAigKdHIHEFLorVDSuf8czRxGlvn87QFZSf72gBav3Kj/bNmomcvWaOKiVZq5bI1a2gaWko6tid013+6bpo4d/6pBndWh2auglXeCV0eOjl/blnSvQQCuV/CzCAG4RKG796d7/zmYOYxe5osIBLv1wPPvKbLSFwldkvrMbC2cW9tbgxXhrfLo4uk+PbSuQlvXWB16f49kjEac8C3dOO0C3fXWl7SjY8eAcaS01yAAVyNgATis0MudBb7//st8Ef37xwSC3WpevblPwIr8O9byYv24EfrKF74n292lIyecotuuPV+N9T5VjVigJc8t6dOkNK29BotdoXZBAPKAgAUAYWNrvPLHCVn9xQpjjfW+mPVaNXvf0MQhB1UxYZTuvvuHqq8P1ZVF6qyWv7hcOzt2avTQ0dntNVhMCrkLApAHBCwACIu1zGc0cAZLil9z1V9PT4/++7//W4cOHdJpp52mqVOn9rm8f5PSslHoXRCAHKPIHQDCGut9WtpUJ1+NV0aSr8arK2aMl7fK0+d63iqPFs6tTek2n376aW3dulXV1dX66le/OnBz53JVDLsgAFlgBgsAosRa5ms4fmRGrRsis1cHDx5UfX29pk+fnqthF59CtwUBcowZLADIkWeffVZbtmzR4MGDdd111zF7Fa3QbUGAHGMGCwASiN4CR4rfpqE/a23v7FVdXZ1OP/30vIy3aBS6LQiQYwQsAEggXm+s/m0a+nvuuef01ltvMXuVSKHbggA5xBIhACQQrzdWvONSaPbqxz/+sQ4ePKiPfexjOvPMM3M1PAAuRcACgATitWOoMCbmljmS9Le//U2vv/46s1dAGSNgAUACC+fWDmjTIEnd1mrxyo0DQlZ07dXEiRM1c+bMfA0VgIsQsAAggUhvLE+MWahILVa0devW6bXXXuudvaqo4GUWKEf85QNAEo31PvXYWP3c+9ZiRWqvOjs7NX78eM2ePTtfQwTgMgQsAEhBvFqs6OPr16/Xyy+/zOwVAAIWAKRiwlEDA1b0ljnRs1c+n0+f/OQn8zxCAG5CwAKAJP6zZaOefWv3gOPTxo/o7YW1ceNGbdiwQYMHD9bXvvY1eTwDC+MBlA8CFgAk8f+9EGPPPEnPv72n99933XWXOjs7NXr0aJ199tn5GhoAl6KTOwDE0dLmV/PqzeqOU+Deba0mLlqlEZ07tX/tXzWiepC+9rWvqbKSl1ag3PEqAKDsRILT9vaAxtZ4tXBu7YBtb/rvQRiPlfTW0w/r0K69GnHiR/SpT30qhyMHUCwIWADKSqqbN8fagzCW4AdbdWj7ZtmKSh346BxmrwBIogYLQJlJtHlztER7DUbbt/4Pst2H5PGOUOCYSY6NE0BxI2ABKCupbt4cr+9VdEf34IfbdND/quSp0rD6z8h31HDnBgqgqBGwAJSVVBqGSon3IIzYt+GPsl1BebzDddRJH+/tiQUABCwAZSVWcIpuGBoR2YPQV+OVkQbsRRjcs10Ht74sVXh0/MxGLbtk2oBCeQDli2pMAGUlEoKSnUUYuW7k+MRFq/pctm99ePZqyAht+MUiDR48OPeDB1A0CFgAyk50cErV2Bqv/OE6rWD7Th18b6NU4dH4mY2EKwADsEQIACmIXlrcv2G1bNchVXmH6ftfv7LAIwPgRsxgAYCSNx+N/Pv7v3laO7a0qWpQla79t/m6dMZHCzVkAC5GwAJQ9lJtPtpY79OLD29ScJRXI0eO1PcXfKkg4wXgfiwRAih7qTYf3bp1q/74xz/K4/Hommuu0ZAhQ/I5TABFhIAFoOyl2nz03nvvVWdnp4YPH66LLrooH0MDUKQIWADKXirNR/1+v37/+9/L4/Ho6quv1tChQ/M1PABFiBosAEUtUXF6ssL1iIVza/vUYEmHm4/u3btXCxYsUGVlpTo7OzVixAhdfPHFefv5ABQnAhaAopWoOF1SSoXr0V/HCmMbN27Uxo0bVVFRIY/Ho2nTpqmigsl/AIkRsAAUrWTF6fEuS9a1PdrOnTvl8Xh07LHHKhAIaPXq1WpqatInPvEJB38SAKWGgAWgaKVanJ7qZbFs3bpV3d3d6u7u1u7du/WFL3xBs2bNSus2AJQf5rkBFK14xekVxqhmSFVa3xPPW2+9paqqKr3//vtqamrS4sWLWSIEkBSvEgCKVvT2NdG6rdX+zi5VeUyf45HC9XS888472rNnj84//3zdeOONhCsAKeGVAkDRaqz3aWlTnTzGDLgs2GM1dFClfDVeGUm+Gq+WNtWlvcnz7t27dfnll+s73/kO4QpAyoy1ttBj6NXQ0GBbW1sLPQwARWbiolWK9UpmJL2zbF5Wt71jxw6NGjVKlZWUrALoyxizzlrbEOsyPo4BKHqpNArN1JgxYwhXANJGwAJQ9GLVYmVSbwUATuFjGYCil6hRKAAUAgELQEmI1ygUAAqBJUIAAACHEbAAAAAcRsACAABwGDVYAFytpc1P8TqAokPAAuBaLW1+LV65UYFgtyTJ3x7Q4pUbJYmQBcDVWCIE4FrNqzf3hquIQLBbzas3F2hEAJAaAhYA19reHkjrOAC4BQELgGvlcgscAMglAhYA12ILHADFiiJ3AK7lpi1wOJsRQDoIWABczQ1b4HA2I4B0sUQIAElwNiOAdDGDBaDo5Hu5jrMZAaSLGSwARSWyXOdvD8jq8HJdS5s/Z/fJ2YwA0kXAAlBUCrFcx9mMANLFEiGAolKI5To3nc0IoDgQsAAUlbE1XvljhKlcL9e54WxGAMWDJUIARYXlOgDFgIAFoKg01vu0tKlONd6q3mPVVbyUAXAXXpUAFKWDXT29/95zIJjzMwkBIB0ELABFh8afANyOgAWg6ND4E4DbEbAAFB0afwJwOwIWgKLDmYQA3I6ABcD1Wtr8mrlsjSYuWqWZy9ZIkpY21clX45WR5KvxamlTHX2qALgGjUYBFFSyjZsjew9Gitojew8ubarTs4vOKdSwASAhZrAAFEwqGzdzxiCAYkTAAlAwqYQnzhgEUIwIWAAKJpXwxBmDAIoRAQtAwaQSnjhjEEAxImABKJhUwlNk70HOGARQTDiLEEDBREJSorMII9cjUAEoJgQsAAVFeAJQilgiBAAAcBgBCwAAwGEsEQJwnWTd3QHA7bIKWMaYZknnSzok6S1JX7LWtocvWyzpy5K6JX3dWrs6u6ECKFXRgapmSJX2d3Yp2GMlhbq7f/O369X67m79oLGuwCMFgNRku0T4pKRJ1trJkl6XtFiSjDGnSLpc0qmSPi3pLmOMJ+6tAChb/bfL2XMg2BuuIqykB55/r88WOgDgZlkFLGvtE9barvCXz0saF/73hZJ+Y609aK19R9Kbkk7P5r4AlKZY2+XEYsPXBYBi4GSR+9WS/hD+t0/S1qjLtoWPAUAf6ewpyP6DAIpF0hosY8yfJI2OcdFN1tpHwte5SVKXpAfSHYAxZr6k+ZI0fvz4dL8dQJEbW+OVP8XgxP6DAIpF0oBlrf1UosuNMVdJOk/SHGttpHDCL+m4qKuNCx+Ldft3S7pbkhoaGmys6wAoXQvn1mrxyo19lgk9FUbd/eqw2H8QQDHJ9izCT0v6tqRPWGsPRF30qKRfG2NukzRW0gmS/pbNfQEoTfG2y4l1jFYNAIqFOTzplME3G/OmpMGSPgwfet5ae234spsUqsvqkvQNa+0fYt/KYQ0NDba1tTXj8QAAAOSLMWadtbYh1mVZzWBZaz+W4LIfSvphNrcPAABQjNgqBwAAwGEELAAAAIcRsAAAABxGwAIAAHAYAQsAAMBhWZ1FCADZamnz0+8KQMkhYAEomJY2f58u7v72gBav3ChJhCwARY0lQgAF07x6c58tciQpEOxW8+rNBRoRADiDgAWgYLbH2eQ53nEAKBYELAAFM7bGm9ZxACgWBCwABbNwbq28VZ4+x7xVnt7NngGgWFHkDqBgIoXsnEUIoNQQsAAUVGO9j0AFoOSwRAgAAOAwAhYAAIDDCFgAAAAOI2ABAAA4jIAFAADgMAIWAACAwwhYAAAADiNgAQAAOIyABQAA4DACFgAAgMMIWAAAAA4jYAEAADiMgAUAAOAwAhYAAIDDCFgAAAAOI2ABAAA4jIAFAADgMAIWAACAwwhYAAAADiNgAQAAOIyABQAA4DACFgAAgMMIWAAAAA4jYAEAADiMgAUAAOAwAhYAAIDDCFgAAAAOI2ABAAA4jIAFAADgMAIWAACAwwhYAAAADiNgAQAAOIyABQAA4DACFgAAgMMIWAAAAA4jYAEAADiMgAUAAOAwAhYAAIDDCFgAAAAOI2ABAAA4jIAFAADgMAIWAACAwwhYAAAADiNgAQAAOIyABQAA4DACFgAAgMMIWAAAAA4jYAEAADiMgAUAAOAwAhYAAIDDCFgAAAAOI2ABAAA4jIAFAADgMAIWAACAwwhYAAAADqss9AAA5F9Lm1/Nqzdre3tAY2u8Wji3Vo31vkIPCwBKBgELKDMtbX4tXrlRgWC3JMnfHtDilRsliZAFAA5hiRAoM82rN/eGq4hAsFvNqzcXaEQAUHoIWECZ2d4eSOs4ACB9BCygzIyt8aZ1HACQPgIWUGYWzq2Vt8rT55i3yqOFc2sLNCIAKD0UuQNlJlLIzlmEAJA7BCygDDXW+whUAJBDLBECAAA4jIAFAADgMAIWAACAwwhYAAAADiNgAQAAOIyABQAA4DBHApYx5j+MMdYYc3T4a2OMudMY86Yx5iVjzDQn7gcAAKAYZB2wjDHHSTpX0ntRhz8j6YTwf/Ml/STb+wEAACgWTsxg3S7p25Js1LELJf3ShjwvqcYYM8aB+wIAAHC9rAKWMeZCSX5r7YZ+F/kkbY36elv4GAAAQMlLulWOMeZPkkbHuOgmSTcqtDyYMWPMfIWWETV+/PhsbgoAAMAVkgYsa+2nYh03xtRJmihpgzFGksZJetEYc7okv6Tjoq4+Lnws1u3fLeluSWpoaLCxrgMAAFBMMl4itNZutNYeY62dYK2doNAy4DRr7U5Jj0q6Mnw24QxJe621O5wZMgAAgLslncHK0OOSPivpTUkHJH0pR/cDoJ+WNr+aV2/W9vaAxtZ4tXBurRrrKYEEgHxyLGCFZ7Ei/7aSrnPqtgGkpqXNr8UrNyoQ7JYk+dsDWrxyoyQRsgAgj+jkDpSQ5tWbe8NVRCDYrebVmws0IgAoTwQsoIRsbw+kdRwAkBsELKCEjK3xpnUcAJAbBCyghCycWytvlafPMW+VRwvn1hZoRABQnnJ1FiGAAogUsnMWIQAUFgELKDGN9T4CFQAUGEuEAAAADiNgAQAAOIyABQAA4DACFgAAgMMIWAAAAA4jYAEAADiMgAUAAOAwAhYAAIDDCFgAAAAOI2ABAAA4jIAFAADgMAIWAACAwwhYAAAADiNgAQAAOIyABQAA4DACFgAAgMMIWAAAAA4jYAEAADjMWGsLPYZexphdkt4t9DhSdLSkDwo9iDLC450/PNb5w2OdXzze+VMuj/Xx1tpRsS5wVcAqJsaYVmttQ6HHUS54vPOHxzp/eKzzi8c7f3isWSIEAABwHAELAADAYQSszN1d6AGUGR7v/OGxzh8e6/zi8c6fsn+sqcECAABwGDNYAAAADiNgZcAYc70x5jVjzCvGmP+KOr7YGPOmMWazMWZuIcdYSowx/2GMscaYo8NfG2PMneHH+iVjzLRCj7EUGGOaw8/rl4wxDxtjaqIu47ntMGPMp8OP55vGmEWFHk8pMcYcZ4x5yhjzavh1ekH4+EhjzJPGmDfC/z+y0GMtFcYYjzGmzRjz+/DXE40xL4Sf3781xgwq9BjzjYCVJmPM2ZIulDTFWnuqpFvDx0+RdLmkUyV9WtJdxhhPwQZaIowxx0k6V9J7UYc/I+mE8H/zJf2kAEMrRU9KmmStnSzpdUmLJZ7buRB+/H6s0HP5FEmfCz/OcEaXpP+w1p4iaYak68KP7yJJf7bWniDpz+Gv4YwFkjZFff0jSbdbaz8maY+kLxdkVAVEwErfv0taZq09KEnW2n+Ej18o6TfW2oPW2nckvSnp9AKNsZTcLunbkqKLBS+U9Esb8rykGmPMmIKMroRYa5+w1naFv3xe0rjwv3luO+90SW9aa9+21h6S9BuFHmc4wFq7w1r7Yvjf+xR64/cp9BjfH77a/ZIaCzLAEmOMGSdpnqSfh782ks6R9GD4KmX5WBOw0neipFnhqc+/GGNOCx/3Sdoadb1t4WPIkDHmQkl+a+2GfhfxWOfe1ZL+EP43j7fzeEzzxBgzQVK9pBckHWut3RG+aKekYws1rhJzh0IfhHvCXx8lqT3qA1tZPr8rCz0ANzLG/EnS6BgX3aTQYzZSoWnn0yStMMZ8JI/DKylJHusbFVoehEMSPd7W2kfC17lJoSWWB/I5NsBpxphhkh6S9A1r7T9DEysh1lprjOE0+iwZY86T9A9r7TpjzCcLPBxXIWDFYK39VLzLjDH/LmmlDfW3+JsxpkehPZf8ko6Luuq48DEkEO+xNsbUSZooaUP4RXGcpBeNMaeLxzpjiZ7bkmSMuUrSeZLm2MM9XHi8ncdjmmPGmCqFwtUD1tqV4cPvG2PGWGt3hMsK/hH/FpCimZIuMMZ8VlK1pCMkLVeodKMyPItVls9vlgjT1yLpbEkyxpwoaZBCG1o+KulyY8xgY8xEhQqw/1aoQRY7a+1Ga+0x1toJ1toJCk0xT7PW7lTosb4yfDbhDEl7o6b9kSFjzKcVmua/wFp7IOointvO+7ukE8JnWg1S6CSCRws8ppIRrgH6haRN1trboi56VNIXw//+oqRH8j22UmOtXWytHRd+nb5c0hpr7RWSnpL0r+GrleVjzQxW+u6RdI8x5mVJhyR9MfxJ/xVjzApJryq0vHKdtba7gOMsZY9L+qxCxdYHJH2psMMpGf8tabCkJ8Ozhs9ba6+11vLcdpi1tssY8zVJqyV5JN1jrX2lwMMqJTMlfUHSRmPM+vCxGyUtU6is48uS3pV0aWGGVxZukPQbY8wPJLUpFHjLCp3cAQAAHMYSIQAAgMMIWAAAAA4jYAEAADiMgAUAAOAwAhYAAIDDCFgAAAAOI2ABAAA4jIAFAADgsP8f9c4ksq1NrcwAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 720x720 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def correlation(X):\n",
" \"\"\"\n",
" :param X: [N,d], 行是数据,列是变量\n",
" return: [d,d],变量之间的协方差\n",
" \"\"\"\n",
" # 每个变量的均值\n",
" feature_mean = np.mean(X, axis=0)\n",
" X_bar = X - feature_mean\n",
" feature_norm = np.sum(X_bar*X_bar, axis=0)**0.5\n",
" # 每个变量进行归一化\n",
" X_norm = X_bar / feature_norm\n",
" cov = X_norm.T.dot(X_norm)\n",
" return cov\n",
"\n",
"# a = np.array([[1,3,5], [5,4,1], [3,8,6]])\n",
"x = np.random.normal(2, 10, (100,1))\n",
"y = np.random.normal(5, 4, (100,1)) + 2*x\n",
"a = np.concatenate((x,y), axis=1)\n",
"print(\"相关系数:\\n\", correlation(a))\n",
"\n",
"plt.figure(figsize=(10, 10)) \n",
"plt.gca().set_aspect(1)\n",
"\n",
"def whiten(X):\n",
" plt.scatter(X[:,0], X[:,1], label=\"origin\")\n",
" corr = correlation(X)\n",
" cov = np.cov(X, rowvar=False)\n",
" eigVals, eigVecs = np.linalg.eig(cov)\n",
" idx = np.argsort(eigVals)[::-1]\n",
" eigVals = eigVals[idx]\n",
" eigVecs = eigVecs[:,idx]\n",
" print(\"特征值:\\n\", eigVals)\n",
" print(\"特征向量:\\n\", eigVecs)\n",
" for i in range(eigVecs.shape[1]):\n",
" # 每列是一个特征向量\n",
" vec = eigVecs[:,i]\n",
" plt.arrow(0,0, vec[0]*10, vec[1]*10, width=0.3, color='black', alpha=0.7)\n",
"\n",
" decorrelated = X.dot(eigVecs)\n",
" print(decorrelated.shape)\n",
" whitened = decorrelated / np.sqrt(eigVals + 1e-5)\n",
" plt.scatter(decorrelated[:,0], decorrelated[:,1], label=\"decorrelated\")\n",
" plt.scatter(whitened[:,0], whitened[:,1], label=\"whiten\")\n",
"\n",
"whiten(a)\n",
"plt.legend()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "jupyter",
"language": "python",
"name": "jupyter"
},
"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.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment