數值分析案例：Newton插值預測2019城市(Asian)溫度、Crout求解城市等溫性的因素係數

# -*- coding:utf-8 -*-
import numpy as np
import pandas as pd
import csv
import matplotlib.pyplot as plt
%matplotlib inline

path = 'DATA.csv'
data = pd.read_csv(path, index_col=0)
data.head(15)

Step 2：ui

數據可視化，可以直觀感覺城市溫度的對比，繪製城市溫度折線圖以下：spa

代碼以下：.net

# 將csv文件格式轉換爲np.array格式，數據格式爲float類型
Data = np.array(data.iloc[:,0:12], dtype=float) # 將數據轉爲float型
print(Data)

# 取出每一個城市的12個月的溫度統計數據
Beijing = Data[0,:]
Chongqing = Data[1,:]
Taipei = Data[2,:]
Tokyo = Data[3,:]
Sapporo = Data[4,:]
Seoul = Data[5,:]
Dikson = Data[6,:]
Vladivostok = Data[7,:]
ChiangMai = Data[8,:]
HatYai = Data[9,:]
DaNang = Data[10,:]
Hanoi = Data[11,:]
Mumbai = Data[12,:]
Erzurum = Data[13,:]

# 以折線圖的方式可視化城市的年度溫度變化數據
x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]

plt.figure(figsize=(13, 13))
plt.rcParams['font.sans-serif']=['SimHei'] #用來正常顯示中文標籤
plt.rcParams['axes.unicode_minus']=False #用來正常顯示負號
plt.plot(x, Beijing,  ms=2, label="北京")
plt.plot(x, Chongqing,ms=2, label="重慶")
plt.plot(x, Taipei,  ms=2, label="臺北")
plt.plot(x, Tokyo, ms=2, label="東京")

plt.plot(x, Sapporo, ms=2, label="札幌")
plt.plot(x, Seoul, ms=2, label="首爾")
plt.plot(x, Dikson, ms=2, label="Dikson")

plt.plot(x, Vladivostok, ms=2, label="海參崴")
plt.plot(x, ChiangMai,  ms=2, label="清邁")
plt.plot(x, HatYai,  ms=2, label="合艾")
plt.plot(x, DaNang, ms=2, label="峴港")

plt.plot(x, Hanoi,  ms=2, label="河內")
plt.plot(x, Mumbai,  ms=2, label="孟買")
plt.plot(x, Erzurum, ms=2, label="埃爾祖魯姆")
plt.xticks(rotation=45)
plt.xlabel("Month")
plt.ylabel("Temperature/degree")
plt.title("2019 temperature change in Asian cities")
plt.legend(loc="upper left")


# 在折線圖上顯示具體數值, ha參數控制水平對齊方式, va控制垂直對齊方式
for y in [Beijing, Chongqing, Taipei, Tokyo,Sapporo,Seoul,Dikson,Vladivostok,
          ChiangMai,HatYai,DaNang,Hanoi,Mumbai,Erzurum]:
    for x1, yy in zip(x, y):
        plt.text(x1, yy + 1, str(yy), ha='center', va='bottom', fontsize=12, rotation=0)
plt.savefig("a.jpg")
plt.show()

Step 3：3d

數據劃分，將城市溫度的數據劃分爲訓練集和測試集。爲了保證數據的分層採樣，訓練數據集：選取一、三、五、七、九、11月的溫度；測試數據集：選取二、四、六、八、十、12月的溫度。

代碼以下：

# train_data
Beijing_train = Data[0,[0,2,4,6,8,10]]
Chongqing_train = Data[1,[0,2,4,6,8,10]]
Taipei_train = Data[2,[0,2,4,6,8,10]]
Tokyo_train = Data[3,[0,2,4,6,8,10]]
Sapporo_train = Data[4,[0,2,4,6,8,10]]
Seoul_train = Data[5,[0,2,4,6,8,10]]
Dikson_train = Data[6,[0,2,4,6,8,10]]
Vladivostok_train = Data[7,[0,2,4,6,8,10]]
ChiangMai_train = Data[8,[0,2,4,6,8,10]]
HatYai_train = Data[9,[0,2,4,6,8,10]]
DaNang_train = Data[10,[0,2,4,6,8,10]]
Hanoi_train = Data[11,[0,2,4,6,8,10]]
Mumbai_train = Data[12,[0,2,4,6,8,10]]
Erzurum_train = Data[13,[0,2,4,6,8,10]]
# Test_data
Beijing_test = Data[0,[1,3,5,7,9,11]]
Chongqing_test = Data[1,[1,3,5,7,9,11]]
Taipei_test = Data[2,[1,3,5,7,9,11]]
Tokyo_test = Data[3,[1,3,5,7,9,11]]
Sapporo_test = Data[4,[1,3,5,7,9,11]]
Seoul_test = Data[5,[1,3,5,7,9,11]]
Dikson_test = Data[6,[1,3,5,7,9,11]]
Vladivostok_test = Data[7,[1,3,5,7,9,11]]
ChiangMai_test = Data[8,[1,3,5,7,9,11]]
HatYai_test = Data[9,[1,3,5,7,9,11]]
DaNang_test = Data[10,[1,3,5,7,9,11]]
Hanoi_test= Data[11,[1,3,5,7,9,11]]
Mumbai_test = Data[12,[1,3,5,7,9,11]]
Erzurum_test = Data[13,[1,3,5,7,9,11]] # 列表下標是從0開始的，下標0,2,4,6,8,10分別對應一、三、五、七、九、11月,而下標1,3,5,7,9,11分別對應二、四、六、八、十、12月

Step 4：

差值方法的選擇，考慮到Lagrange 插值的龍格現象，而Hermite插值須要導數的參與，分段插值的擬合性較弱。綜合考量，選擇Newton差商法進行差值預測。

以北京市爲例，預測二、四、六、八、十、12月的溫度以下：

代碼以下：

#Newton差商法
import numpy as np


def Qda_Table():
    for i in range(n+1):
        for j in range(n):
            if initial[i][j] is not None or i >= j + 2:
                continue
            else:
                initial[i][j] = (initial[i - 1][j] - initial[i - 1][j - 1]) / (
                        initial[0][j] - initial[0][j - i + 1])
            # print("%d，%d的當前值爲：" % (i, j) + str(initial[i][j]))


def calculate_factors(no, vari):  # 計算每次的相乘的自變量因子多項式
    factor = 1
    if no == 1:
        return factor
    else:
        for i in range(no-1):
            factor = factor * (vari-x[i])
    return factor


def final_calculate(k):
    result = 0
    for time in range(n):
        result = result + initial[time+1][time] * calculate_factors(time+1, v[k])
    return result



x = [1, 3, 5, 7, 9, 11]# 輸入的月份
y = Beijing_train.copy()
n = len(x)
sum = 0
v = [2, 4, 6, 8, 10, 12]# 須要預測的月份
y_predict = [] # 用來存放預測月份的溫度值

initial = []  # 初始化表格Table
for p in range(n+1):
    initial.append([])
    for q in range(n):
        initial[p].append(None)  # 生成n*n的二維列表，同時列表元素初始化爲0
for k in range(n):
    initial[0][k] = x[k]
    initial[1][k] = y[k]
# print("初始化後的表格以下：")
# print(initial) # 將x和y的值填入表格，完成Table的初始化。

Qda_Table()  # 計算每一個位置的差商填表
QDA_table = np.array(initial).T  # 轉置後輸出結果
# print("最終的差商表格以下：")
# print(QDA_table)

for i in range(n):
    y_predict.append(np.around(final_calculate(i),decimals=4))
print("預測月份的溫度信息以下：")
print(y_predict)

第二部分：「Crout求解分析城市的等溫性影響因素係數」

Step 1：

選擇第一部分城市所在國家(Korea, South、Japan、Turkey、India、Vietnam、Thailand、China、Russia)的環境因素數據，但因爲環境因素數據包含衆多的數據項，以及爲了構建後面的線性方程組求解因素係數並分析其相關性，這裏經查閱資料和篩選後，選取如下8項環境因素做爲城市等溫性因素分析的探索指標：*

elevation（海拔）

cropland_cover（農田覆蓋率）

tree_canopy_cover（植被覆蓋率）

rain_mean_annual（平均年降水量）

rain_seasonailty（季節性降水量）

temp_annual_range（年溫度浮動範圍）

cloudiness（雲層量）

temp_mean_annual（年平均溫度）

詳細數據以下：

代碼以下：

# 數據可視化
path = 'ENVI.csv'
source = pd.read_csv(path, index_col=0)
source.head(15)

Step 2：

構建線性方程組，假設可能性因素與等溫性存在線性關係。根據計算獲得的係數，初步分析其與等溫性的相關性狀況：假如|係數|近似爲0，則相關性較低；|係數|大於等於0.5，將其斷定爲較高的相關性。方程求解採用Crout分解法。

自變量：elevation、cropland_cover、tree_canopy_cover、rain_mean_annual、rain_seasonailty、temp_annual_range、cloudiness、temp_mean_annual

因變量：isothermality

結果以下：

代碼以下：

# Crout分解法求解
import numpy as np
import copy

def UAndL_Figure(k):
    # 按行計算U矩陣中Lj的值
    for j in range(mu):
        if j <= k:
            continue
        else:
            uj = A[k][j]/A[k][k]
            U[k][j] = uj
    # 根據U_Figure計算獲得的li係數改變A矩陣
    for j in range(mu):
        j = j+k
        for i in range(nu):
            if j < nu-1:
                A[i][j+1] = A[i][j+1] - U[k][j+1] * A[i][k]
            else:
                continue


mu, nu = 8, 8 # 方程和自變量個數爲8
source = np.array(source.iloc[:,0:9], dtype=float) # 將數據轉爲float型
A, L, U = [], [], []  # 初始化ALU矩陣
for p in range(nu):
    A.append([]), L.append([]), U.append([])
    for q in range(mu):
        x_in = source[p][q]
        A[p].append(x_in)
        if p == q:
            L[p].append(None), U[p].append(1)
        else:
            L[p].append(None), U[p].append(0)
A, L, U = np.array(A), np.array(L), np.array(U)


b_Const = source[:,[8]]  # 因變量y矩陣


for i in range(nu):
    UAndL_Figure(i)
L = copy.deepcopy(A)

y = np.dot(np.linalg.inv(L), b_Const)
x = np.dot(np.linalg.inv(U), y)
x = np.around(x,decimals=3)
# print("x的計算結果以下：")
# print(x)


t=[i for item in x for i in item]#將二維的X數據一維化
print("可能性因素的係數的結果以下：")
print(t)

Step 3：

通過Step 2的線性相關係數求解，咱們獲得與等溫性相關程度較高的因素：Cropland_cover(農田覆蓋率)、tree_canopy_cover(植被覆蓋率)、temp_annual_range(年溫度浮動範圍)。可視化全部城市在這些指標項上的數據狀況：

以cropland_over因素爲例,代碼以下：

# 將csv文件格式轉換爲np.array格式，數據格式爲float類型
source = np.array(source.iloc[:,0:9], dtype=float) # 將數據轉爲float型
print(source)

# 設置matplotlib正常顯示中文和負號
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False     

# 生成畫布
plt.figure(figsize=(20, 8), dpi=70)
# 橫座標城市名字
country_name = ['Korea, South','Japan','Turkey','India','Vietnam','Thailand','China','Russia']
# 縱座標海拔
y = source[:,1]
x=range(len(country_name))

plt.bar(x,y,color=['b','r','g','y','c','m','y','k','c','g','g'])
plt.xticks(x, country_name)
plt.xlabel("Country")
plt.ylabel("Cropland_cover")
plt.title("2019 Cropland_cover of the countries in Asian")

# 在折線圖上顯示具體數值, ha參數控制水平對齊方式, va控制垂直對齊方式
for x1, yy in zip(x, y):
    plt.text(x1, yy + 1, str(yy), ha='center', va='bottom', fontsize=12, rotation=0)

plt.show()

Step 4：

經過熱力圖分析潛在的因素，由於所選的國家個數以及選擇的環境數據項侷限性，並不可以很好的直接創建城市(Asian)等溫性與其餘環境因素之間的關係，因此*根據皮爾森相關係數(Pearson correlation coefficient)，用來反映兩個變量線性相關程度的統計量。

該熱力圖的解釋以下：

橫軸指標分別爲：Elevation、 rain_mean_annual、 rain_seasonailty、cloudiness、temp_mean_annual、Isothermality；縱軸數據分別爲：Korea, South、Japan、Turkey、India、Vietnam、Thailand、China、Russia各個國家。

該熱力圖用於探索不一樣國家在相同的環境因素指標上是否具備相關性。能夠發現Elevation和rain _mean_annual在不一樣國家的數據統計上具備潛在的相關性。

該熱力圖不包括Cropland_cover(農田覆蓋率)、tree_canopy_cover(植被覆蓋率)、temp_ annual_range(年溫度浮動範圍)因素分析。

代碼以下：

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
f, ax = plt.subplots(figsize = (6,6))

#顯示熱力圖
sns.heatmap(source[:,[0,3,4,6,7,8]], annot=True,fmt='.2f',cbar=True)
ax.set_title('City Statistics ')
ax.set_xlabel('Index Of Environmental Attributes')
ax.set_ylabel('City name')

3、實驗結果與分析

「採用Newton插值預測2019城市(Asian)溫度」，該部分採用了Newton插值法進行對2019城市(Asian)溫度的插值預測，經檢驗，插值的結果是比較好的，與實際統計的偏差較小。實際的意義：當監測月份的溫度數據缺失時，可採用Newton插值法進行數據的補全以供使用。

「Crout求解分析城市的等溫性影響因素係數」，該部分爲探索性實驗，其結論能夠做爲研究和探索方向，但未經大量驗證是不可以斷章取義的。與等溫性相關的可能性因素計算結果爲：Cropland_cover(農田覆蓋率)、tree_canopy_cover(植被覆蓋率)、temp_annual_range(年溫度浮動範圍),也是符合地理常識的認知。

最後結合熱力圖分析，以皮爾森相關係數做爲理論依據，分析潛在因素有：Elevation（海拔）和rain_mean_annual（年平均降水量），也是實驗能夠繼續的因素方向。以及能夠考慮非線性關係的構建，利用神經網絡的易擬合性進行進一步的強化分析*。