{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "798076ab-a8f3-47d4-a221-00b465568f41",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statistics import mean\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import explained_variance_score,r2_score,median_absolute_error,mean_squared_error,mean_absolute_error\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "plt.rcParams[\"font.sans-serif\"] = [\"SimHei\"]  # 设置字体\n",
    "plt.rcParams[\"font.size\"] = 16\n",
    "plt.rcParams[\"axes.unicode_minus\"] = False  # 正常显示负号"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "697bdeac-4ca7-4cf9-ad01-2ee00b4c7c72",
   "metadata": {},
   "outputs": [],
   "source": [
    "def scatter_out_1(x,y): ## x,y为两个需要做对比分析的两个量。\n",
    "    # ==========计算评价指标==========\n",
    "    BIAS = mean(x - y)\n",
    "    MSE = mean_squared_error(x, y)\n",
    "    RMSE = np.power(MSE, 0.5)\n",
    "    R2 = r2_score(x, y)\n",
    "    MAE = mean_absolute_error(x, y)\n",
    "    EV = explained_variance_score(x, y)\n",
    "    print('==========算法评价指标==========')\n",
    "    print('BIAS:', '%.3f' % (BIAS))\n",
    "    print('Explained Variance(EV):', '%.3f' % (EV))\n",
    "    print('Mean Absolute Error(MAE):', '%.3f' % (MAE))\n",
    "    print('Mean squared error(MSE):', '%.3f' % (MSE))\n",
    "    print('Root Mean Squard Error(RMSE):', '%.3f' % (RMSE))\n",
    "    print('R_squared:', '%.3f' % (R2))\n",
    "    # ===========Calculate the point density==========\n",
    "    xy = np.vstack([x, y])\n",
    "    z = stats.gaussian_kde(xy)(xy)\n",
    "    # ===========Sort the points by density, so that the densest points are plotted last===========\n",
    "    idx = z.argsort()\n",
    "    x, y, z = x[idx], y[idx], z[idx]\n",
    "    def best_fit_slope_and_intercept(xs, ys):\n",
    "        m = (((mean(xs) * mean(ys)) - mean(xs * ys)) / ((mean(xs) * mean(xs)) - mean(xs * xs)))\n",
    "        b = mean(ys) - m * mean(xs)\n",
    "        return m, b\n",
    "    m, b = best_fit_slope_and_intercept(x, y)\n",
    "    regression_line = []\n",
    "    for a in x:\n",
    "        regression_line.append((m * a) + b)\n",
    "    fig,ax=plt.subplots(figsize=(12,9),dpi=600)\n",
    "    scatter=ax.scatter(x,y,marker='o',c=z, edgecolors='b',s=15,label='LST',cmap='Spectral_r')\n",
    "    cbar=plt.colorbar(scatter,shrink=1,orientation='vertical',extend='both',pad=0.015,aspect=30,label='frequency', )\n",
    "    plt.plot([0,35],[0,35],'black',lw=1.5)  # 画的1:1线，线的颜色为black，线宽为0.8\n",
    "    plt.plot(x,regression_line,'red',lw=1.5)      # 预测与实测数据之间的回归线\n",
    "    plt.axis([0,35,0,35])  # 设置线的范围\n",
    "    plt.title(\"总孔体积拟合结果 $10^2 cm^3$\", fontdict={\"fontsize\":16})\n",
    "    plt.xlabel('预测值', fontdict={\"fontsize\":16})\n",
    "    plt.ylabel('真实值', fontdict={\"fontsize\":16})\n",
    "    plt.text(0.5,34, '$N=%.f$' % len(y), fontdict={\"fontsize\":16}) # text的位置需要根据x,y的大小范围进行调整。\n",
    "    plt.text(0.5,33, '$R^2=%.3f$' % R2, fontdict={\"fontsize\":16})\n",
    "    plt.text(0.5,32, '$BIAS=%.4f$' % BIAS, fontdict={\"fontsize\":16})\n",
    "    plt.text(0.5,31, '$RMSE=%.3f$' % RMSE, fontdict={\"fontsize\":16})\n",
    "    plt.xlim(0,35)                                  # 设置x坐标轴的显示范围\n",
    "    plt.ylim(0,35)                                  # 设置y坐标轴的显示范围\n",
    "    plt.savefig('./总孔体积.png',dpi=300, bbox_inches='tight',pad_inches=0)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "9e1279b5-4b18-4f57-bdc0-197ba84cda33",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "2d5bcdf5-cc45-40b7-8000-a0e41b18ef93",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('./rst/总孔体积_比表.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "8908a8b3-6bd5-4d71-b012-cbdd690b6a31",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>真实值</th>\n",
       "      <th>预测值</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>184.000000</td>\n",
       "      <td>184.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>267.871020</td>\n",
       "      <td>272.776239</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>696.475264</td>\n",
       "      <td>693.395059</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>0.060000</td>\n",
       "      <td>0.069085</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>0.539250</td>\n",
       "      <td>0.570501</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>0.877000</td>\n",
       "      <td>0.889113</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>1.673250</td>\n",
       "      <td>1.551479</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>3322.000000</td>\n",
       "      <td>3225.575700</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               真实值          预测值\n",
       "count   184.000000   184.000000\n",
       "mean    267.871020   272.776239\n",
       "std     696.475264   693.395059\n",
       "min       0.060000     0.069085\n",
       "25%       0.539250     0.570501\n",
       "50%       0.877000     0.889113\n",
       "75%       1.673250     1.551479\n",
       "max    3322.000000  3225.575700"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "d6bfbaae-a201-4963-be4e-4731a3fe8a37",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==========算法评价指标==========\n",
      "BIAS: 0.049\n",
      "Explained Variance(EV): 0.921\n",
      "Mean Absolute Error(MAE): 0.622\n",
      "Mean squared error(MSE): 3.797\n",
      "Root Mean Squard Error(RMSE): 1.949\n",
      "R_squared: 0.921\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 849340 characters]",
      "text/plain": [
       "<Figure size 7200x5400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "scatter_out_1(df['预测值'].values/100, df['真实值'].values/100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "6e95d4f4-eb89-4263-8f81-f1694ad7e97f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 128548 characters]",
      "text/plain": [
       "<Figure size 1200x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "samples = np.random.choice(df.index.values, 50, replace=False)\n",
    "plt.figure(figsize=(12, 9))\n",
    "plt.plot(range(len(samples)), df.iloc[samples]['预测值'].values, 'o-', label='预测值')\n",
    "plt.plot(range(len(samples)), df.iloc[samples]['真实值'].values, '*-', label='真实值')\n",
    "plt.xlabel('预测值 $(10^2 cm^3/g)$', fontdict={\"fontsize\":16})\n",
    "plt.ylabel('真实值 $(10^2 cm^3/g)$', fontdict={\"fontsize\":16})\n",
    "plt.title('氮气吸附量拟合结果')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37c6c196-acb0-43c3-a234-c05a2f44dda4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
