first

README.md

@@ -1,2 +0,0 @@
-# URL
-

案例11 铁路旅客客流量预测/code/.idea/.gitignore

@@ -0,0 +1,8 @@
+# 默认忽略的文件
+/shelf/
+/workspace.xml
+# 基于编辑器的 HTTP 客户端请求
+/httpRequests/
+# Datasource local storage ignored files
+/dataSources/
+/dataSources.local.xml

案例11 铁路旅客客流量预测/code/.idea/code.iml

@@ -0,0 +1,8 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<module type="PYTHON_MODULE" version="4">
+  <component name="NewModuleRootManager">
+    <content url="file://$MODULE_DIR$" />
+    <orderEntry type="jdk" jdkName="Python 3.12 (sohu)" jdkType="Python SDK" />
+    <orderEntry type="sourceFolder" forTests="false" />
+  </component>
+</module>

案例11 铁路旅客客流量预测/code/.idea/inspectionProfiles/profiles_settings.xml

@@ -0,0 +1,6 @@
+<component name="InspectionProjectProfileManager">
+  <settings>
+    <option name="USE_PROJECT_PROFILE" value="false" />
+    <version value="1.0" />
+  </settings>
+</component>

案例11 铁路旅客客流量预测/code/.idea/misc.xml

@@ -0,0 +1,4 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.12 (sohu)" project-jdk-type="Python SDK" />
+</project>

案例11 铁路旅客客流量预测/code/.idea/modules.xml

@@ -0,0 +1,8 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="ProjectModuleManager">
+    <modules>
+      <module fileurl="file://$PROJECT_DIR$/.idea/code.iml" filepath="$PROJECT_DIR$/.idea/code.iml" />
+    </modules>
+  </component>
+</project>

案例11 铁路旅客客流量预测/code/.idea/other.xml

@@ -0,0 +1,6 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="PySciProjectComponent">
+    <option name="PY_INTERACTIVE_PLOTS_SUGGESTED" value="true" />
+  </component>
+</project>

案例11 铁路旅客客流量预测/code/9.2 预处理客流量数据.py

@@ -0,0 +1,166 @@
+# 代码9-1
+import numpy as np
+import os
+from itertools import chain
+import pandas as pd
+import re
+import matplotlib.pyplot as plt 
+
+
+path = '../data/201501-201603'  # 文件夹目录
+filename = os.listdir(path)  # 得到文件夹下的所有文件名称
+n_filename = len(filename)
+datalist = []
+date = []
+
+for i in range(n_filename):
+    name = '../data/201501-201603/' + filename[i]
+    datalist.append(['../data/201501-201603/' +
+                        filename[i] + '/' + j for j in os.listdir(name)])
+    date.append(len(datalist[i]))
+filedata = list(chain.from_iterable(datalist))  # 将二维列表转换为一维列表
+n_file = sum(date)
+SaveFile_Name = '../tmp/Station.csv'  # 数据合并后要保存的文件名（Station表合并）
+
+
+def Deal_Fun(b):
+    for m in range(b):
+        data = pd.read_excel(filedata[m])
+        row = data.shape[0]
+        Line =[]  # 存放始发日期和上车站信息
+        Head_d = []  # 存放始发日期
+        for i in range(row):
+            if '始发日期' in data.iloc[i, 0]:
+                Head_d.append(re.findall('[0-9\—]+',data.iloc[i, 0]))
+            if data.iloc[i,0] == '上车站':
+                 Line.append(re.findall('[A-Z]{2}[0-9]{2} ', data.iloc[i-1, 0]))
+        Line = pd.DataFrame(Line)  # 以数据框形式存放上车站
+        Head_d = pd.DataFrame(Head_d)
+        Line['Head'] = 0  # 新加存放始发日期的列
+        for i in range(len(Line)):
+            Line.iloc[i,1] = Head_d.iloc[0, 0]
+        # 提取列车站点信息
+        # 提取上车站点位置信息
+        on_station = [i for i, x in enumerate(data.iloc[:, 0]) if x == '上车站']
+        # 提取上车人数合计位置
+        on_count = [i for i, x in enumerate(data.iloc[:, 0]) if x == '上车人数合计']
+        Size = pd.DataFrame(np.zeros([len(on_station), 2]),
+                            columns=['on_station', 'on_count'])    
+        Size['on_station'] = on_station
+        Size['on_count'] = on_count
+        Size['off_count'] = 0
+        for h in range(len(Size.iloc[:, 0])):
+            Size.loc[h,'off_count'] = [i for i, x in enumerate(
+                    data.iloc[Size.iloc[h,0],:]) if x == '下车人数合计'][0]
+        # 提取上下车站点、人数和时间
+        # 下车站点
+        guodu = pd.DataFrame(data.iloc[:,0])
+        off_station = []
+        for j in range(len(Size.iloc[:,0])):
+            off_station.append(guodu.iloc[Size.iloc[j, 0]+2:Size.loc[j, 'on_count'], 0]) 
+        sum_station = 0
+        for i in range(len(off_station)):
+            sum_station = sum_station + len(off_station[i])        
+        Out_off = pd.DataFrame(np.zeros([sum_station,3]),
+                               columns = ['off_station', 'off_man', 'off_time'])
+        h = 0
+        for i in range(len(off_station)):
+            Out_off.iloc[h:h+len(off_station[i]), 0] = list(off_station[i])
+            h = h + len(off_station[i])
+        # 下车人数
+        off_man = []
+        for i in range(len(Size)):
+            data1 = pd.DataFrame(data.iloc[:, Size.loc[i, 'off_count']])
+            off_man.append(data1.iloc[Size.iloc[i,0]+2:Size.loc[i, 'on_count'], 0])
+        h = 0
+        for i in range(len(off_man)):
+            Out_off.iloc[h:h+len(off_man[i]),1] = list(off_man[i])
+            h = h + len(off_man[i])  
+        # 下车时间
+        off_time = []
+        for i in range(len(Size)):
+            data1 = pd.DataFrame(data.iloc[:,1])
+            off_time.append(data1.iloc[Size.iloc[i,0]+2:Size.loc[i, 'on_count'], 0])
+        h = 0        
+        for i in range(len(off_time)):
+             Out_off.iloc[h:h+len(off_time[i]), 2] = list(off_time[i])
+             h = h + len(off_time[i])  
+        # 上车信息
+        Come = pd.DataFrame(np.zeros([sum_station, 3]),
+                            columns = ['on_station', 'on_man', 'on_time'])
+        # 上车站点
+        on_station1 = []
+        for i in range(len(Size)):
+            data1 = pd.DataFrame(data.iloc[Size.loc[i,'on_station'], :])
+            on_station1.append(data1.iloc[2:Size.loc[i,'off_count'], 0])
+        h = 0
+        for i in range(len(on_station1)):
+            Come.iloc[h:h+len(on_station1[i]),0] = list(on_station1[i])
+            h = h + len(on_station1[i])  
+        # 上车人数
+        on_man = []
+        for i in range(len(Size)):
+            data1 = pd.DataFrame(data.iloc[Size.loc[i,'on_count'],:])
+            on_man.append(data1.iloc[2:Size.loc[i,'off_count'], 0])
+        h = 0
+        for i in range(len(on_man)):
+            Come.iloc[h:h+len(on_man[i]), 1] = list(on_man[i])
+            h = h + len(on_man[i])  
+        # 上车时间
+        on_time = []
+        for i in range(len(Size)):
+            data1 = pd.DataFrame(data.iloc[Size.loc[i,'on_station']+1, :])
+            on_time.append(data1.iloc[2:Size.loc[i,'off_count'], 0])
+        h = 0
+        for i in range(len(on_time)):
+            Come.iloc[h:h+len(on_time[i]),2] = list(on_time[i])
+            h = h + len(on_time[i]) 
+        Station = pd.DataFrame(np.zeros([len(Out_off),7]),
+                               columns = ['on_station', 'on_man', 'on_time',
+                                          'off_man', 'off_time', 'date', 'train'])
+        Station['on_station'] = list(Out_off.iloc[:, 0])
+        Station['off_man'] = list(Out_off['off_man'])
+        Station['off_time'] = list(Out_off['off_time'])
+        k = 0
+        for i in range(len(on_man)):
+                Station.loc[k:k+len(on_man[i])-1, 'on_man'] = list(on_man[i])
+                Station.loc[k:k+len(on_time[i])-1, 'on_time'] = list(on_time[i])
+                Station.loc[k+len(on_time[i]):k-1+len(off_time[i]), 'on_time'] = 0.1
+                Station.loc[k+len(on_man[i]):k-1+len(off_man[i]), 'on_man'] = 0.1
+                Station.loc[k:k-1+len(off_man[i]), 'date'] = Line.iloc[i, 1]
+                Station.loc[k:k-1+len(off_man[i]), 'train'] = Line.iloc[i, 0]
+                k = k + len(off_man[i])          
+        Station.to_csv(SaveFile_Name, encoding="utf_8",
+                       index=False, header=False, mode='a+')
+Deal_Fun(n_file)
+
+
+
+# 代码9-2
+Train_Station = pd.read_csv('../tmp/Station.csv',
+                            header=None, encoding='utf-8')
+Train_Station.columns = ['on_station', 'on_man', 'on_time',
+                         'off_man', 'off_time', 'date','Station']
+Train_Station.fillna(value=0, inplace=True)  # 处理nan值
+
+for i in range(len(Train_Station)):
+    for j in range(len(Train_Station.iloc[0, :])):
+        if Train_Station.iloc[i, j] == '0.1':
+            Train_Station.iloc[i, j] = 0
+            
+s_date = [re.findall('[0-9]+', i)[0][0:4] +
+          '-' + re.findall('[0-9]+', i)[0][4:6] +
+          '-' + re.findall('[0-9]+', i)[0][6:8] 
+          for i in Train_Station.loc[:, 'date'] ]
+
+Train_Station.loc[:,'date'] = s_date
+# 部分数据为空格，将其替换为0
+ind_on = [i for i in Train_Station.index if Train_Station.loc[i, 'on_man'] == ' ']
+ind_off = [i for i in Train_Station.index if Train_Station.loc[i, 'off_man'] == ' ']
+Train_Station.loc[ind_on, 'on_man'] = 0
+Train_Station.loc[ind_off, 'off_man'] = 0
+Train_Station['on_man'] = Train_Station['on_man'].astype(float)
+Train_Station['off_man'] = Train_Station['off_man'].astype(float)
+Train_Station.to_csv('../tmp/Train_Station.csv', encoding = 'utf-8')
+Train_Station = pd.read_csv('../tmp/Train_Station.csv', index_col=0,
+                            encoding = 'utf-8')

案例11 铁路旅客客流量预测/code/9.3 探索客流量数据.py

@@ -0,0 +1,148 @@
+# 代码9-3
+import pandas as pd 
+import matplotlib.pyplot as plt
+import re
+
+
+# 读取处理后的CSV文件
+Train_Station = pd.read_csv('../tmp/Train_Station.csv', index_col=0, encoding='utf-8')
+print(Train_Station.head())
+# 打印前几行数据
+
+# 提取出上车站点信息
+on = pd.DataFrame(Train_Station['on_station'])
+on = on.drop_duplicates()
+on['on_mean'] = 0
+on['off_mean'] = 0
+print(on.head())  # 打印前几行数据
+
+# 计算每个站点的平均上车人数和下车人数
+for i in range(len(on)):
+    # on数据框是上下车人数
+    data = Train_Station[Train_Station.iloc[:, 0] == on.iloc[i, 0]]
+    # 遍历on数据框中的每个站点，根据当前站点名称从Train_Station数据框中筛选出相关数据。
+    if not data.empty:
+        on.iloc[i, 1] = sum(data.iloc[:, 1]) / len(data)
+        on.iloc[i, 2] = sum(data.iloc[:, 3]) / len(data)
+        # 将计算得到的平均上车人数存储到on数据框的第i行
+    else:
+        print(f"Warning: No data for station {on.iloc[i, 0]}")
+print("每个站点的上车人数")
+print(on.head())  # 打印前几行数据
+
+# 随机采样20个站点
+on_sample = on.sample(20, random_state=44)
+# random_state=44确保每次运行时抽样的结果一致。
+on_sample.index = on_sample['on_station']
+print(on_sample.head())  # 打印前几行数据
+
+# 设置Matplotlib支持中文显示
+plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
+plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
+
+# 绘制条形图
+fig, ax = plt.subplots(figsize=(15, 10))  # 创建一个较大的图形区域
+on_sample.plot(kind='bar', title='各站点上下车人数分析', ax=ax)
+plt.xticks(rotation=45)  # 旋转x轴标签以便更好地显示
+plt.ylabel('上下车人数（人）')
+plt.xlabel('站点')
+plt.tight_layout()  # 自动调整子图参数，使之填充整个图像区域
+plt.show()  # 显示图表
+
+
+
+# 代码9-4
+# ST111-01站点上下车客流人数分析
+Train_ST111_01 = Train_Station[Train_Station.iloc[:, 0] == 'ST111-01']
+
+def process_data(df, columns):
+    df = df.iloc[:, columns]
+    df.index = range(len(df))
+    df.iloc[:, 1] = df.iloc[:, 1].astype(str).apply(lambda x: int(x[0:2]) if x != '0' else 0)
+    return df
+
+# 上车客流分析
+On_t = process_data(Train_ST111_01, [1, 2])
+on_mean_t = On_t.groupby('on_time')['on_man'].mean()
+
+# 下车客流分析
+Off_t = process_data(Train_ST111_01, [3, 4])
+off_mean_t = Off_t.groupby('off_time')['off_man'].mean()
+
+# 创建一个图形窗口，并在其中添加两个子图
+fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
+
+# 上车客流图
+on_mean_t.plot(kind='bar', ax=ax1, color='blue', title='每个时段上车客流量')
+ax1.set_xticklabels(on_mean_t.index)
+ax1.set_ylabel('上车人数（人）')
+ax1.set_xlabel('时间')
+
+# 下车客流图
+off_mean_t.plot(kind='bar', ax=ax2, color='red', title='每个时段下车客流量')
+ax2.set_xticklabels(off_mean_t.index)
+ax2.set_ylabel('下车人数（人）')
+ax2.set_xlabel('时间')
+
+# 调整子图之间的间距
+plt.tight_layout()
+
+# 显示图形
+plt.show()
+
+
+
+
+# 代码9-5
+# ST111-01站点节假日客流变化
+holiday = open('../data/2015-2016节假日.csv', encoding='utf-8')
+holiday = pd.read_csv(holiday)
+for i in range(len(holiday)):
+    s = re.findall('[0-9]+', holiday.iloc[i, 0])
+    l = s[0]
+    if len(s[1]) < 2:
+        l = l + '-' + '0' +s[1]
+    else:
+        l = l + '-' +s[1]
+    if len(s[2]) < 2:
+        l = l + '-' + '0' + s[2]
+    else:
+        l = l + '-' + s[2]
+        # 格式化日期
+    holiday.iloc[i, 0] = l
+    #将格式化后的日期保存回 holiday DataFrame 中。
+holiday.to_csv('../tmp/holiday.csv', encoding='utf-8')
+Train_ST111_01 = Train_Station[Train_Station.iloc[:,0] == 'ST111-01']
+on_h = Train_ST111_01.groupby('date')['on_man'].sum()
+#  按日期分组，计算每天的上车人数总和。
+on_h = pd.DataFrame(on_h)
+on_h['date'] = 0
+on_h['holiday'] = 0
+
+# 添加日期和类型（工作日或者小长假）
+#遍历 holiday DataFrame，将对应的日期和节假日类型添加到 on_h DataFrame 中。
+for i in range(len(holiday)):
+    for j in range(len(on_h)):
+        if holiday.iloc[i, 0] == on_h.index[j]:
+            on_h.loc[on_h.index[j], 'holiday'] = holiday.iloc[i, 1]
+            on_h.loc[on_h.index[j], 'date'] = holiday.iloc[i, 0]
+# 节假日影响图                  
+fig = plt.figure(figsize=(20, 12))  # 设置画布
+ax = fig.add_subplot(1, 1, 1)
+for i in range(len(on_h) - 2):
+    for j in range(i+1, len(on_h)-1):
+        if on_h.iloc[i, 2] == on_h.iloc[j,2] and on_h.iloc[i, 2] != on_h.iloc[j+1, 2]:
+            if on_h.iloc[i, 2] == '小长假':
+                ax.scatter(on_h.iloc[i, 1], on_h.iloc[i, 0], color='red', linewidth=5)
+                ax.plot(on_h.iloc[i:j+1, 1], on_h.iloc[i:j+1, 0], color='black')
+            else:
+                ax.plot(on_h.iloc[i:j+1, 1], on_h.iloc[i:j+1, 0], color='black')
+on_h.to_csv('../tmp/on_h.csv', encoding='utf-8')
+plt.xlabel('日期')
+plt.ylabel('上车人数（人）')
+plt.legend(['点加连线-节假日', '黑线-全部日期'],loc = 8)
+plt.title('节假日客流量变化')
+plt.xticks((0, 50, 100, 150, 200, 250, 300, 350, 400), rotation=30)
+plt.tight_layout()
+plt.rcParams.update({'font.size': 33})
+plt.show()

案例11 铁路旅客客流量预测/code/9.4 构建模型并预测节假日客流量.py

@@ -0,0 +1,318 @@
+# 代码9-6
+#ARIMA模型
+import pandas as pd 
+import matplotlib.pyplot as plt
+import warnings
+import math
+import itertools
+import numpy as np
+import statsmodels.api as sm
+from statsmodels.graphics.tsaplots import plot_acf
+from statsmodels.stats.diagnostic import acorr_ljungbox
+from statsmodels.tsa.stattools import adfuller as ADF
+from statsmodels.tsa.arima_model import ARIMA
+import matplotlib.ticker as ticker 
+
+
+on_h = pd.read_csv('../tmp/on_h.csv', index_col=0
+                            ,encoding='utf-8')
+train = pd.DataFrame(on_h.iloc[0:426, 0])
+test = pd.DataFrame(on_h.iloc[426:, 0])
+#train.plot(title = '训练集时序图')  # 画训练集时序图
+x = train.index
+y = train.on_man
+fig, ax = plt.subplots(1,1)
+ax.plot(x, y)
+ticker_spacing = 70
+ax.xaxis.set_major_locator(ticker.MultipleLocator(ticker_spacing))
+plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
+plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
+plt.xticks(rotation = 45)
+plt.xlabel('日期')
+plt.ylabel('客流量')
+plt.title('训练集时序图')
+plt.tight_layout()  # tight_layout()方法可以保证图像的完整度
+plt.show()
+
+plot_acf(train,lags=400)
+plt.xlabel('日期索引')
+plt.ylabel('自相关系数')
+plt.title('训练集自相关图')
+plt.tight_layout()
+plt.show()  # 画训练集自相关图
+
+print('原始序列的ADF检验结果为：', ADF(train['on_man']))  # 检验训练集平稳性
+print('白噪声检验结果为：', acorr_ljungbox(train['on_man'], lags=1))
+
+
+
+# 代码9-7
+import statsmodels.api as sm
+train['on_man'] = train['on_man'].astype(float)
+# 定阶
+bic_matrix = []  # bic矩阵，选择bic矩阵中最小值对应的行(p),列(q)
+# 存在部分报错，所以用try来跳过报错
+for p in range(11):
+    tmp = []
+    for q in range(5):
+        try:
+            tmp.append(sm.tsa.ARIMA(train, order=(p, 1, q)).fit().bic)  
+        except:
+            tmp.append(None)
+    bic_matrix.append(tmp)
+bic_matrix = pd.DataFrame(bic_matrix)  # 从中可以找出最小值
+p, q = bic_matrix.stack().idxmin()  # 先用stack展平，然后用idxmin找出最小值位置
+print('BIC最小的p值和q值为：%s、%s' % (p, q))
+
+model = sm.tsa.ARIMA(train, order=(p, 1, q)).fit()  # 建立ARIMA(p, 1, 1)模型
+summary = model.summary()  # 给出一份模型报告
+forecast = model.forecast(10)
+print('10天的预测结果、标准误差和置信区间分别为：\n', forecast)
+
+# pre = pd.DataFrame(forecast, columns = ['predict'])
+pre = pd.DataFrame({'predict':forecast})
+pre.index = test.index #如果索引不同将pre加到test中时会出错
+test['pre'] = pre
+plt.plot(test.index, test.on_man)
+plt.plot(test.index, test.pre, linestyle=':')
+plt.xticks(rotation = 45)
+plt.legend(['真实', '预测'])
+plt.title('预测结果和实际结果的对比')
+plt.xlabel('日期')
+plt.ylabel('上车人数')
+plt.tight_layout()  # tight_layout()方法可以保证图像的完整度
+plt.show()
+
+
+
+# 代码9-8
+# 剔除节假日之后的时序图
+on_nh = on_h[on_h.iloc[:,2] != '小长假']  # 剔除假期
+x = on_nh.index
+y = on_nh.on_man
+fig, ax = plt.subplots(1,1)
+ax.plot(x, y)
+plt.xticks(range(len(x)), on_nh['date.1'])
+ticker_spacing = 70
+ax.xaxis.set_major_locator(ticker.MultipleLocator(ticker_spacing))
+plt.xticks(rotation = 45)
+plt.title('剔除节假日后ST111-01站点客流量')
+plt.xlabel('日期')
+plt.ylabel('客流量（人）')
+plt.tight_layout()  # tight_layout()方法可以保证图像的完整度
+plt.show()
+
+
+
+# 代码9-9
+# 剔除节假日及前一天后的时序图
+new_h = on_h #存放节假日更新的数据
+new_h.index = range(len(new_h))
+for i in range(1, len(new_h)):
+    if new_h.iloc[i, 2] == '小长假' :
+        new_h.iloc[i-1, 2] = '小长假'
+new_nh = new_h[new_h.iloc[:,2] != '小长假']  # 剔除节假日及其前一天
+
+plt.plot(new_nh['date.1'], new_nh['on_man'],color='green')
+plt.gca().xaxis.set_major_locator(ticker.MultipleLocator(70))
+plt.xticks(rotation = 45)
+plt.title('剔除节假日及其前一天后ST111-01站点客流量')
+plt.legend(['on_man'], loc=2)
+plt.xlabel('日期')
+plt.ylabel('客流量（人）')
+plt.tight_layout()  # tight_layout()方法可以保证图像的完整度
+plt.savefig('../tmp/9-10.png', dpi=1080)
+plt.show()
+
+
+# 代码9-10
+new_nh.index = range(len(new_nh))
+train1 = pd.DataFrame(new_nh.iloc[0:378, 0])
+test1 = pd.DataFrame(new_nh.iloc[378:, 0])  
+p = q = range(4)
+d = range(2)
+pdq = list(itertools.product(p, d, q))
+seasonal_pdq = [(x[0], x[1], x[2], 7)for x in list(itertools.product(p, d, q))]
+print('Examples of parameter combinations for Seasonal ARIMA...')
+print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[1]))
+print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[2]))
+print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[3]))
+print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[4]))
+warnings.filterwarnings("ignore")  # specify to ignore warning messages
+sa =[]
+for param in pdq:
+    for param_seasonal in seasonal_pdq:
+        try:
+            mod = sm.tsa.statespace.SARIMAX(train1,
+                                            order=param,
+                                            seasonal_order=param_seasonal,
+                                            enforce_stationarity=False,
+                                            enforce_invertibility=False)
+            results = mod.fit() 
+            print('ARIMA{}x{}7 - AIC:{}'.format(param, param_seasonal, results.aic))
+            sa.append(param)
+            sa.append(param_seasonal)
+            sa.append(results.aic)
+        except:
+            continue
+AIC = [i for i in sa if type(i) == np.float64]
+AIC_min = min(AIC)
+for i in np.arange(2,len(sa),3):
+    if sa[i] == min(AIC):
+        param = sa[i-2]
+        param_seasonal = sa[i-1]
+mod = sm.tsa.statespace.SARIMAX(train1,
+                                order=(param),
+                                seasonal_order=(param_seasonal),
+                                enforce_stationarity=False,
+                                enforce_invertibility=False) 
+print('模型最终定阶为：', (param, param_seasonal))
+results = mod.fit() 
+print(results.summary().tables[1])
+fig = plt.figure(figsize=(15, 12))
+results.plot_diagnostics(figsize=(15, 12), fig=fig)
+plt.show()
+pre_10 = results.predict(start=378, end=387,dynamic=True)
+out_pre = pd.DataFrame(np.zeros([10,3]),columns = ['real', 'pre', 'error'])
+out_pre['real'] = list(test1['on_man'])
+out_pre['pre'] = list(pre_10)
+# 计算相对误差
+error_seasonal = (out_pre.loc[:, 'pre']-out_pre.loc[:,'real'])/out_pre.loc[:,'real']
+# 平均相对误差
+error_mean = abs(error_seasonal).mean()
+print('预测平均相对误差为：', error_mean)
+
+
+
+# 代码9-11
+# 节假日客流规律
+# 由于之前步骤对on_h有修改，所以此处重新载入
+holiday = pd.read_csv('../tmp/holiday.csv', index_col=0
+                            ,encoding = 'utf-8')
+Train_Station = pd.read_csv('../tmp/Train_Station.csv', index_col=0
+                            ,encoding='utf-8')
+Train_ST111_01 = Train_Station[Train_Station.iloc[:, 0] == 'ST111-01']
+on_h = Train_ST111_01.groupby('date')['on_man'].sum()
+on_h = pd.DataFrame(on_h)
+on_h['date'] = 0
+on_h['holiday'] = 0
+# 添加日期和类型（工作日或者小长假）
+for i in range(len(holiday)):
+    for j in range(len(on_h)):
+        if holiday.iloc[i,0] == on_h.index[j]:
+            on_h.loc[on_h.index[j], 'holiday'] = holiday.iloc[i,1]
+            on_h.loc[on_h.index[j], 'date'] = holiday.iloc[i,0]
+# 2015春节
+fig = plt.figure(figsize=(12, 6))  # 设置画布
+ax = fig.add_subplot(1, 1, 1)
+ax.plot(on_h.loc['2015-01-19':'2015-02-18', 'on_man'], color = 'blue')
+ax.plot(on_h.loc['2015-02-18':'2015-02-25', 'on_man'], color = 'red', linestyle=':')
+ax.plot(on_h.loc['2015-02-25':'2015-03-01', 'on_man'], color = 'blue')
+plt.xlabel('日期')
+plt.ylabel('上车人数')
+plt.title('2015春节客流量')
+plt.legend(['工作日','节假日'])
+plt.xticks(rotation = 45)
+plt.show() 
+
+# 2015劳动节
+fig1 = plt.figure(figsize=(12, 6))  # 设置画布
+ax1 = fig1.add_subplot(1, 1, 1)
+ax1.plot(on_h.loc['2015-04-27':'2015-05-01', 'on_man'], color = 'blue')
+ax1.plot(on_h.loc['2015-05-01':'2015-05-04', 'on_man'], color = 'red', linestyle=':')
+ax1.plot(on_h.loc['2015-05-04':'2015-05-11', 'on_man'], color = 'blue')
+plt.xlabel('日期')
+plt.ylabel('上车人数')
+plt.title('五一客流量')
+plt.legend(['工作日','节假日'])
+plt.xticks(rotation = 45)
+plt.show() 
+
+# 国庆和中秋
+fig2 = plt.figure(figsize=(12, 6))  # 设置画布
+ax2 = fig2.add_subplot(1, 1, 1)
+ax2.plot(on_h.loc['2015-09-21':'2015-09-26', 'on_man'],color = 'blue')
+ax2.plot(on_h.loc['2015-09-26':'2015-09-28', 'on_man'],color = 'red', linestyle=':')
+ax2.plot(on_h.loc['2015-09-28':'2015-09-30', 'on_man'],color = 'blue')
+ax2.plot(on_h.loc['2015-09-30':'2015-10-08', 'on_man'],color = 'red', linestyle=':')
+ax2.plot(on_h.loc['2015-10-08':'2015-10-12', 'on_man'],color = 'blue')
+plt.xlabel('日期')
+plt.ylabel('上车人数')
+plt.title('中秋和国庆客流量')
+plt.legend(['工作日', '节假日'])
+plt.xticks(rotation = 45)
+plt.show()
+
+
+
+# 代码9-12
+# 2015和2016春节客流量比较
+compare = pd.DataFrame(on_h.loc['2015-02-05':'2015-02-26', 'on_man'])
+compare['2016'] = list(on_h.loc['2016-01-25':'2016-02-15', 'on_man'])
+compare.columns = ['2015', '2016']
+compare.index = range(len(compare))
+plt.plot(compare.index, compare['2015'], linestyle=':')
+plt.plot(compare.index, compare['2016'])
+plt.legend(['2015', '2016'])
+plt.xlabel('日期')
+plt.ylabel('客流量')
+plt.title('2015和2016春节客流量比较')
+
+
+
+# 代码9-13
+# 2015春节节假日波动系数
+M = on_h.loc['2015-01-01':'2015-02-17']
+M = M[M.loc[:,'holiday'] != '小长假']
+M1 =  on_h.loc['2015-01-01':'2015-02-07']
+M1 = M1[M1.loc[:,'holiday'] != '小长假']
+B_coef = []
+# 春节前10天
+for i in on_h.loc['2015-02-08':'2015-02-17',:].index:
+    B_coef.append('%.2f' % (on_h.loc[i, 'on_man']/math.ceil(M1.iloc[-30:].mean())))
+# 春节及春节后两天
+for i in on_h.loc['2015-02-18':'2015-02-26',:].index:
+    B_coef.append('%.2f' % (on_h.loc[i, 'on_man']/math.ceil(M.iloc[-30:].mean())))
+B_coef = pd.DataFrame(B_coef)
+B_coef.columns = ['on_man']
+B_coef[u'on_man'] = B_coef[u'on_man'].astype(float)
+fig3 = plt.figure(figsize=(8,6))  # 设置画布
+ax3 = fig3.add_subplot(1, 1, 1)
+ax3.plot(B_coef.iloc[:,0], color='blue')
+plt.xlabel('Index')
+plt.ylabel('系数')
+plt.title('2015春节客流量波动系数')
+# 设置数字标签
+for a, b in zip(B_coef.index, B_coef.iloc[:, 0]):
+    plt.text(a, b, b, ha='center', va='bottom', fontsize=20) 
+plt.legend()
+plt.show()
+
+
+# 预测2016年春节客流量
+MM = on_h.loc['2015-01-01':'2016-02-06',:]
+MM = MM[MM.loc[:, 'holiday'] != '小长假']
+MM_mean = math.ceil(MM.iloc[-30:].mean())
+MM1 = on_h.loc['2015-01-01':'2016-01-27',:]
+MM1 = MM1[MM1.loc[:, 'holiday'] != '小长假']
+MM1_mean = math.ceil(MM1.iloc[-30:].mean())
+pre_2016_b = B_coef.iloc[0:10, 0] * MM1_mean
+pre_2016_a = B_coef.iloc[10:, 0] * MM_mean
+pre_2016 = pd.DataFrame(on_h.loc['2016-01-28':'2016-02-15', 'on_man'])
+pre_2016['pre'] = 0
+pre_2016.loc[0:10,'pre'] = list(pre_2016_b)
+pre_2016.loc[10:,'pre'] = list(pre_2016_a)
+pre_2016.columns=['real', 'pre']
+plt.plot(pre_2016.index, pre_2016.real)
+plt.plot(pre_2016.index, pre_2016.pre, linestyle=':')
+plt.xticks(rotation=45)
+plt.legend(['real', 'pre'])
+plt.xlabel('日期')
+plt.ylabel('客流量')
+plt.title('预测2016年春节客流量')
+
+# 计算相对误差
+error_pre = (pre_2016.loc[:, 'pre'] - pre_2016.loc[:, 'real'])/pre_2016.loc[:, 'real']
+# 平均相对误差
+error_pre_mean = abs(error_pre).mean()
+print('预测的平均相对误差为：', error_pre_mean)

案例11 铁路旅客客流量预测/data/.Rhistory

@@ -0,0 +1,512 @@
+f<-function(x){x^3-x-1}
+fe(f,1,2,e)
+hujiao<-c( 0,1,2,3,4,5,6)
+pinshu<-c(7,10,12,8,3,2,0)
+lambda<-sum(hujiao*pinshu)/sum(pinshu)
+pinshu1<-c(7,10,12,8,5)
+p<-matrix(0,5,1)
+for(i in 1:4){
+p[i,]<-(lambda^(i-1)*exp(-lambda))/factorial(i-1)
+}
+p
+p[5,]<-1-sum(p[1:4,])
+p
+sum(p[1:4,])
+Newtons=function(fun,x,ep=1e-5,it_max=100){
+index=0;k=1
+while(k<=it_max){
+x1=x;obj=fun(x);
+x=x-solve(obj$J,obj$f);
+norm=sqrt((x-x1)%*%(x-x1))
+if(norm<ep){
+index=1;break
+}
+k=k+1
+}
+obj=fun(x)
+list(root=x,it=k,index=index,Funval=obj$f)
+}
+#求解方程的函数
+funs=function(x){
+f=c(x[1]^2+x[2]^2-5,(x[1]+1)*x[2]-(3*x[1]+1))
+J=matrix(c(2*x[1],2*x[2],x[2]-3,x[1]+1),nrow=2,byrow=T)
+list(f=f,J=J)
+}
+Newtons(funs,c(0,1))
+fe<-function(f,a,b,eps){
+if(f(a)*f(b)>0)
+break
+if(f(a)*f(b)<0){
+repeat{
+if(f(a)*f((a+b)/2)<0) b=(a+b)/2 else a=(a+b)/2
+if(abs(a-b)<eps)
+return((a+b)/2)
+}
+}
+}
+e<-10^(-6)
+f<-function(x){x^3-x-1}
+fe(f,1,2,e)
+x<-c(2300,1700,2000,2500,1700,1400,1500,1200,1200,1500,1600,1700,2000,1400,1600,1200,1900,1100,1800,1300)
+mean(x)
+x<-c(2300,1700,2000,1500,1700,1400,1500,1200,1200,1500,1600,1700,2000,1400,1600,1200,1900,1100,1800,1300)
+mean(x)
+y<-x-mean(x)
+y^2/19
+y^2
+sum(y^2)/19
+((1-1/10)/20)*sum(y^2)/19
+z<-c(200,150,170,150,160,130,140,100,110,140,150,160,180,130,150,100,180,100,170,120)
+mean(z)
+f<-mean(z)
+sum(y*f)/19
+y*f
+f<-z-mean(z)
+sum(y*f)/19
+((1-1/10)/20)*(826.0256-2*0.0915*8831.579-0.0915^2*99578.95)
+((1-1/10)/20)*(826.0256-2*0.0915*8831.579+0.0915^2*99578.95)
+v<-((1-1/10)/20)*(826.0256-2*0.0915*8831.579+0.0915^2*99578.95)
+146.3-1.96*sqrt(v)
+146.3+1.96*sqrt(v)
+sum(y*f)/19
+(sum(y*f)/19)/(sqrt(99578.95*826.0156))
+yu<-c(95,97,87,120,110,115,103,102,92,105)
+xi<-c(150,155,140,180,175,185,165,160,150,170)
+ym<-yu-mean(yu)
+xm<-xi-mean(xi)
+mean(yu)
+mean(xi)
+sum(y^2)/9
+sum(xm^2)/9
+((1-1/12)/10)*sum(xm^2)/9
+((1-1/12)/10)*sum(ym^2)/9
+sqrt(((1-1/12)/10)*sum(ym^2)/9)
+sum(ym^2)/9
+sqrt(sum(ym^2)/9)
+sum(ym*m)/9
+sum(ym*xm)/9
+sqrt(sum(xm^2)/9)
+146.333/(14.56785*10.34086)
+0.971*(14.56785/10.34086)
+163+1.368*(100-102.6)
+((1-1/12)/10)*(212.222+1.368^2*106.933-2*1.368*146.333)
+((1-1/12)/10)*212.222*(1-1.368^2)
+((1-1/12)/10)*212.222*(1-1.368^2)
+1.368^2
+((1-1/12)/10)*212.222*(1-0.971^2)
+wh<-c(8.2,6.5,13.7,5.6,11.8,11.6,17,9.8,9.8,7)
+yh<-c(89,56,102,76,97,79,83,52,36,52)
+y<-sum(wh*yh)
+y
+y<-sum(wh*yh/100)
+y
+wh<-c(8.2,6.5,13.7,5.6,11.8,11.6,17,9.8,8.8,7)
+yh<-c(89,56,102,76,97,79,83,52,36,52)
+y<-sum(wh*yh/100)
+y
+nh<-c(16,13,27,11,24,23,34,20,18,14)
+sum(nh)
+sh<-c(105,74,186,97,106,89,112,73,44,65)
+vy<-sum(wh^2*sh^2/nh)
+vy
+vy<-sum((wh/100)^2*sh^2/nh)
+vy
+y
+y-1.96*sqrt(vy)
+y+1.96*sqrt(vy)
+wh1<-c(0.2,0.3,0.5)
+ph<-c(0.1,0.2,0.4)
+sum(wh1*ph)
+0.28*(1-0.28)/100
+s<-c(0.09,0.16,0.24)
+sum(wh1*s)/(0.002016)
+s1<-(1/43)*(1-1/43)
+s2<-(2/57)*(1-2/57)
+(0.7*s1+0.3*s2)/100+(0.3*s1+0.7*s2)/10000
+s3<-(0.7*s1+0.3*s2)/100+(0.3*s1+0.7*s2)/10000
+sqrt(s3)
+rs<-c(6.342,5.925,6.476,7.017,6.891,6.216,6.602,6.185,6.6)
+rc<-c(6.387,6.188,6.457,6.947,6.875,6.439,6.243,6.227,6.6)
+rs_m<-mean(rs)
+rc_m<-mean(rc)
+rs_m
+rc_m
+sum((rs-rs_m)^2)/9
+sum((rs-rs_m)^2)/8
+sum((rc-rc_m)^2)/8
+sum((rs-rs_m)^2)/8+0.027^2
+sum((rc-rc_m)^2)/8+0.015^2
+x<-c(51,62,49,73,101,48,65,49,73,61,58,52,65,49,55)
+z<-c(42,53,40,45,63,31,38,30,54,45,51,29,46,37,42)
+y<-x/z
+y
+y<-z/x
+y
+M<-sum(x)/15
+M
+sum(y)/(15*M)
+22/800
+sum(z)/(15*M)
+10/726
+10/726-1
+10/33-1
+1- 15/87
+sum(z)/(15)
+(1 - 15/87)*sum((z-sum(z)/(15))^2)/(15*M^2*14)
+sqrt((1 - 15/87)*sum((z-sum(z)/(15))^2)/(15*M^2*14))
+0.0064*M^2/(1 - 15/87)*sum((z-sum(z)/(15))^2)
+(1 - 15/87)*sum((z-sum(z)/(15))^2)/0.0064*M^2
+M
+((1 - 15/87)*sum((z-sum(z)/(15))^2))/(0.0064*M^2)
+6.2*5.2
+sqrt(44.97624)
+x<c(42,51,49,55,47,58,43,59,48,41,60,52,61,49,57,63,45,46,62,58)
+x<-c(42,51,49,55,47,58,43,59,48,41,60,52,61,49,57,63,45,46,62,58)
+yi<-c(6.2,5.8,6.7,4.9,5.2,6.9,4.3,5.2,5.7,6.1,6.3,6.7,5.9,6.1,6,4.9,5.3,6.7,6.1,7)
+m<-sum(x)/20
+m
+mean(yi)/m
+sum(x*yi)/(20*m)
+sum(x*yi)/(20)
+sum(yi)/20
+((1 - 20/386)*sum((x*yi-sum(x*yi)/(20))^2))/(20*19*m^2)
+sqrt(((1 - 20/386)*sum((x*yi-sum(x*yi)/(20))^2))/(20*19*m^2))*1.96+5.91
+5.91-sqrt(((1 - 20/386)*sum((x*yi-sum(x*yi)/(20))^2))/(20*19*m^2))*1.96
+326^2-188^2/6
+sqrt(326^2-188^2/6)
+(188/sqrt(326^2-188^2/6))*sqrt(10)
+1650/500
+Mi<-c(32,45,36,54)
+mi<-c(4,5,4,6)
+y1<-c(4,2,3,6)
+y2<-c(2,2,4,3,6)
+y3<-c(3,2,5,8)
+y4<-c(4,3,6,2,4,6)
+Y<-c(sum(y1),sum(y2),sum(y3),sum(y4))
+Y
+Y_y<-sum(Y)/4
+Y_y
+Y_u<-sum(Y)/4
+S<-c(sum((y1-mean(y1))^2)/3,sum((y2-mean(y2))^2)/4,sum((y3-mean(y3))^2)/3,sum((y4-mean(y4))^2)/5)
+S
+f<-c(7/8,8/9,8/9,8/9)
+((100*0.6sum((Y-Y_u)^2))/12+2.5*sum((Mi^2*fS)/mi))/(500^2)
+((100*0.6*sum((Y-Y_u)^2))/12+2.5*sum((Mi^2*fS)/mi))/(500^2)
+((100*0.6*sum((Y-Y_u)^2))/12+2.5*sum((Mi^2*f*S)/mi))/(500^2)
+MI<-Mi^2
+MI*f*S)/mi
+(MI*f*S)/mi
+((100*0.6*sum((Y-Y_u)^2))/12+2.5*sum((MI*f*S)/mi))/(500^2)
+a<-(100*0.6*sum((Y-Y_u)^2))/12
+((100*0.6*sum((Y-Y_u)^2))/12+2.5*sum((MI*f*S)/mi))/(500^2))
+b<-2.5*sum((MI*f*S)/mi)
+c<-500^2
+d<-a+b
+d
+d/c
+e<-d/c
+e^2
+f<-c(8/9,8/9,8/9,8/9)
+((100*0.6*sum((Y-Y_u)^2))/12+2.5*sum((MI*f*S)/mi))/(500^2)
+MI
+32^2
+Mi<-c(32,45,36,54)
+MI<-Mi^2
+mi<-c(4,5,4,6)
+y1<-c(4,2,3,6)
+y2<-c(2,2,4,3,6)
+y3<-c(3,2,5,8)
+y4<-c(4,3,6,2,4,6)
+Y<-c(sum(y1),sum(y2),sum(y3),sum(y4))
+Y
+Y_u<-sum(Y)/4
+S<-c(sum((y1-mean(y1))^2)/3,sum((y2-mean(y2))^2)/4,sum((y3-mean(y3))^2)/3,sum((y4-mean(y4))^2)/5)
+f<-c(7/8,8/9,8/9,8/9)
+((100*0.6*sum((Y-Y_u)^2))/12+2.5*sum((MI*f*S)/mi))/(500^2)
+(MI*f*S)/mi
+sqrt(0.04899633)
+500*sum(Y)/(sum(Mi))
+sum(Y)/(sum(Mi))
+(500*sum(Y)/(sum(Mi)))/(sum(Mi))
+Y1<-c(mean(y1),mean(y2),mean(y3),mean(y4))
+Y<-Y1*Mi
+Y_u<-sum(Y)/4
+S<-c(sum((y1-mean(y1))^2)/3,sum((y2-mean(y2))^2)/4,sum((y3-mean(y3))^2)/3,sum((y4-mean(y4))^2)/5)
+f<-c(7/8,8/9,8/9,8/9)
+((100*0.6*sum((Y-Y_u)^2))/12+2.5*sum((MI*f*S)/mi))/(500^2)
+q<-((100*0.6*sum((Y-Y_u)^2))/12+2.5*sum((MI*f*S)/mi))/(500^2)
+sqrt(q)
+500*sum(Y)/(sum(Mi))
+sum(Y)/(sum(Mi))
+((100*0.6*sum(MI*((Y1-3.95)^2)))/12+2.5*sum((MI*f*S)/mi))/(500^2)
+sqrt((100*0.6*sum(MI*((Y1-3.95)^2)))/12+2.5*sum((MI*f*S)/mi))/(500^2)
+sqrt(((100*0.6*sum(MI*((Y1-3.95)^2)))/12+2.5*sum((MI*f*S)/mi))/(500^2))
+0.404254/3.3
+0.2674/3.95
+x<-c(51,62,49,73,101,48,65,49,73,61,58,52,65,49,55)
+z<-c(42,53,40,45,63,31,38,30,54,45,51,29,46,37,42)
+y<-z/x
+y
+M<-sum(x)/15
+M
+sum(z)/(15)
+sum((z-mean(z))^2)
+s<-sum((z-mean(z))^2)
+s1<-s/(mean(x)^2)
+s2<-s1/87
+0.0064*(6.2^2)+(s2-0.0064)*6.2
+s1
+0.0064*(7^2)+(s2-0.0064)*7
+s/(0.0064*mean(z)^2)
+sqrt(s/(0.0064*mean(z)^2))
+s/(15*mean(x)^2)
+ss<-s/(15*mean(x)^2)
+sss<-ss/0.0064
+(6.2^2)+(sss-1)*6.2
+15*sss
+(6.3^2)+(sss-1)*6.3
+2
+(6.2^2)+(sss-1)*6.2
+ss<-s/(87*mean(x)^2)
+sss<-ss/0.0064
+(6.2^2)+(sss-1)*6.2
+87*sss
+mean(z)
+mean(x)
+(6.22^2)+(sss-1)*6.22
+ss<-s/(15*mean(x)^2)
+sss<-ss/0.0064
+(6.22^2)+(sss-1)*6.22
+15*sss
+(6.21^2)+(sss-1)*6.21
+(6.201^2)+(sss-1)*6.201
+(6.2001^2)+(sss-1)*6.2001
+(6.200001^2)+(sss-1)*6.200001
+15*sss
+(6.1800001^2)+(sss-1)*6.1800001
+(6.1600001^2)+(sss-1)*6.1600001
+(6.1700001^2)+(sss-1)*6.1700001
+(6.1781^2)+(sss-1)*6.1781
+15*sss
+(6.1761^2)+(sss-1)*6.1761
+zc<-c(200,160,170,100,120,180,150,100,170,130)
+zz<-c(2300,1700,2000,1200,1300,2000,1600,1200,1800,1400)
+mean(zc)
+sum(zc)
+y<-mean(zc)
+sum((zc-y)^2)
+sum((zc-y)^2)/9
+10560/9
+ss<-sum((zc-y)^2)/9
+(1173.3*(1-1/15))/10
+sqrt((1173.3*(1-1/15))/10)
+148-10.46
+148+10.46
+sum(zz)
+1600*148/1650
+x<-mean(zz)
+(sum((zc-y)*(zz-x)))/9
+sum((zc-y)*(zz-x))
+114000/9
+sum((zz-x)^2)/9
+148/1650
+(1-1/15)*(1173.3-2*0.09*12666.67+(0.09^2)*142777.8)/10
+(1173.3*1.96^2)/(0.01*148^2)
+21/150
+21/(1+21/150)
+19/0.8
+28/30
+26/30
+w<-c(0.26,0.28,0.24,0.22)
+p<-c(0.9,0.867,0.933,0.9)
+sum(w*p)
+pp<-c(0.9,0.9333,0.9,0.8667,0.9333,0.9667)
+ww<-c(0.18,0.21,0.14,0.09,0.16,0.22)
+n<-c(27,28,27,26,28,29)
+1-pp
+pw<-1-pp
+nn<-n-1
+nn
+sum(ww^2*pp*pw/nn)
+sum(((ww^2)*pp*pw)/nn)
+1-pp
+ww^2
+pp*pw
+sum(((ww^2)*pp*pw)/nn)
+sqrt(sum(((ww^2)*pp*pw)/nn))
+www<-ww^2
+wq<-www*pp
+wq<-www*pp*pw
+q<-wq/nn
+q
+q<-sum(wq/nn)
+q
+sum((w^2)*p*(1-p)/(n-1))
+sum(((w^2)*p*(1-p))/(n-1))
+w<-c(0.26,0.28,0.24,0.22)
+p<-c(0.9,0.867,0.933,0.9)
+n<-c(27,26,28,27)
+sum(((w^2)*p*(1-p))/(n-1))
+sqrt(0.0008965)
+sqrt(sum(((w^2)*p*(1-p))/(n-1)))
+x<-c(4,2,4)
+y<-c(2,1,2)
+x/y
+sqrt(sum(((w^2)*p*(1-p))/(n-1)))
+sqrt(0.0008965)
+pp*(1-PP)
+pp<-c(0.9,0.9333,0.9,0.8667,0.9333,0.9667)
+pp*(1-PP)
+pp*(1-pp)
+P<-pp*(1-pp)
+P<-pp*(1-pp)
+s<-sum(P*ww^2)
+q<-(0.1*0.924/1.96)^2
+a<-sum(ww*P)/1650000
+s/(a+q)
+s
+a
+q
+q<-(0.01*0.924/1.96)^2
+a<-sum(ww*P)/1650000
+s/(a+q)
+s/a
+s<-(sum(sqrt(P)*ww))^2
+q<-(0.01*0.924/1.96)^2
+a<-sum(ww*P)/1650000
+s/(a+q)
+a<-sum(ww*P)/180
+s/(a+q)
+1650000
+1650000
+a<-sum(ww*P)/1650000
+s/(a+q)
+1-pp
+P
+P
+pp<-c(0.9,0.9333,0.9,0.8667,0.9333,0.9667)
+ww<-c(0.18,0.21,0.14,0.09,0.16,0.22)
+P<-pp*(1-pp)
+s<-(sum(sqrt(P)*ww))^2
+q<-(0.01*0.924/1.96)^2
+a<-sum(ww*P)/1650000
+s/(a+q)
+pp<-c(0.9,0.9333,0.9,0.8667,0.9333,0.9667)
+ww<-c(0.18,0.21,0.14,0.09,0.16,0.22)
+P<-pp*(1-pp)
+sh<-pp*(1-pp)
+s<-(sum(sqrt(sh)*ww))^2
+q<-(0.01*0.924/1.96)^2
+a<-sum(ww*P)/1650000
+s/(a+q)
+pp<-c(0.9,0.9333,0.9,0.8667,0.9333,0.9667)
+ww<-c(0.18,0.21,0.14,0.09,0.16,0.22)
+sh<-pp*(1-pp)
+s<-(sum(sqrt(sh)*ww))^2
+q<-(0.01*0.924/1.96)^2
+a<-sum(ww*sh)/1650000
+s/(a+q)
+2999*0.18
+160/(1/5)
+s<-124/(1/6)
+s
+s*(5/6)+160
+s<-150/(1/15)
+s*(1-1/6-1/5)+160+124
+150*15*19/30
+ss<-150*15*19/30
+ss+124+160
+(800+780+1709)/3
+t<-c(800,780,1709)
+yy<-1096.333
+sum((t-yy)^2)/6
+t-yy
+(t-yy)^2
+sum((t-yy)^2)/6
+x1<-c(6,7,8,12)
+x2<-c(6,8,9,12)
+x3<-c(7,9,10,15)
+y1<-c(6,6,7)
+y2<-c(7,8,9)
+y3<-c(8,9,10)
+y4<-c(12,12,15)
+mean(x1)
+mean(x2)
+mean(x3)
+mean(y1)
+mean(y2)
+mean(y3)
+mean(y4)
+z<-c(6.33,8,9,13)
+mean(z)
+sum((x1-mean(x1))^2)/3
+sum((x2-mean(x2))^2)/3
+sum((x3-mean(x3))^2)/3
+sum((y1-mean(y1))^2)/2
+sum((y2-mean(y2))^2)/2
+sum((y3-mean(y3))^2)/2
+sum((y4-mean(y4))^2)/2
+(sum(x1)+sum(x2)+sum(x3))/12
+xx<-c(11,12,13,14,15)
+6*sum((xx-13)^2)/29
+(sum((x1-9.08)^2)+sum((x2-9.08)^2)+sum((x3-9.08)^2))/11
+(6.92+6.25+11.58)/3
+(0.33+5)/4
+29/30*2.07-6*4*2.5/30
+29/30*2.07
+11/12*7.54-3*3*8.25/12
+11/12*7.54-3*3*1.33/12
+(1-1/3)*7.54/4
+(1-1/4)*7.54/3
+(12-3-1)/18+1/4
+1/4-(12-3-1)/18
+(12-3-1)/18
+1/4-(4-3-1)/18
+1/4+(4-3-1)/18
+(6+8+9+12)/4
+w<-c(0.26,0.28,0.24,0.22)
+p<-c(0.9,0.867,0.933,0.9)
+sum(w*p)
+p*(1-p)
+s<-p*(1-p)
+sqrt(s)
+q<-sqrt(s)
+s[1]
+s[2]
+w[1]*q[1]/sum(w*q)
+w[2]*q[2]/sum(w*q)
+w[3]*q[3]/sum(w*q)
+w[4]*q[4]/sum(w*q)
+pp<-c(0.9,0.9333,0.9,0.8667,0.9333,0.9667)
+ww<-c(0.18,0.21,0.14,0.09,0.16,0.22)
+sh<-pp*(1-pp)
+s<-(sum(sqrt(sh)*ww))^2
+q<-(0.01/1.96)^2
+a<-sum(ww*sh)/1650000
+s/(a+q)
+s/(a)
+pp<-c(0.9,0.9333,0.9,0.8667,0.9333,0.9667)
+ww<-c(0.18,0.21,0.14,0.09,0.16,0.22)
+sh<-pp*(1-pp)
+s<-(sum(sqrt(sh)*ww))^2
+q<-(0.01/1.96)^2
+a<-sum(ww*sh)/1650000
+s/(a)
+s/a
+s/q
+w<-c(0.26,0.28,0.24,0.22)
+p<-c(0.9,0.867,0.933,0.9)
+sum(w*p)
+s<-p*(1-p)
+q<-sqrt(s)
+w[4]*q[4]/sum(w*q)
+(sum(w*q)^2)/((0.05/1.96)^2)
+138*0.26
+138*0.32
+138*0.2
+138*0.22
+536+520+417+304+396+392
+148-1.96*10.46
+148+1.96*10.46
+setwd("F:\\交通案例库\\铁路客流量预测")
+data.ST250 <- Train_Station[grep("ST250",Train_Station$on.station),]  #取出ST250站点的数据

案例11 铁路旅客客流量预测/data/2015-2016节假日.csv

@@ -0,0 +1,446 @@
+date,holiday
+2015/1/1,小长假
+2015/1/2,小长假
+2015/1/3,小长假
+2015/1/4,工作日
+2015/1/5,工作日
+2015/1/6,工作日
+2015/1/7,工作日
+2015/1/8,工作日
+2015/1/9,工作日
+2015/1/10,工作日
+2015/1/11,工作日
+2015/1/12,工作日
+2015/1/13,工作日
+2015/1/14,工作日
+2015/1/15,工作日
+2015/1/16,工作日
+2015/1/17,工作日
+2015/1/18,工作日
+2015/1/19,工作日
+2015/1/20,工作日
+2015/1/21,工作日
+2015/1/22,工作日
+2015/1/23,工作日
+2015/1/24,工作日
+2015/1/25,工作日
+2015/1/26,工作日
+2015/1/27,工作日
+2015/1/28,工作日
+2015/1/29,工作日
+2015/1/30,工作日
+2015/1/31,工作日
+2015/2/1,工作日
+2015/2/2,工作日
+2015/2/3,工作日
+2015/2/4,工作日
+2015/2/5,工作日
+2015/2/6,工作日
+2015/2/7,工作日
+2015/2/8,工作日
+2015/2/9,工作日
+2015/2/10,工作日
+2015/2/11,工作日
+2015/2/12,工作日
+2015/2/13,工作日
+2015/2/14,工作日
+2015/2/15,工作日
+2015/2/16,工作日
+2015/2/17,工作日
+2015/2/18,小长假
+2015/2/19,小长假
+2015/2/20,小长假
+2015/2/21,小长假
+2015/2/22,小长假
+2015/2/23,小长假
+2015/2/24,小长假
+2015/2/25,工作日
+2015/2/26,工作日
+2015/2/27,工作日
+2015/2/28,工作日
+2015/3/1,工作日
+2015/3/2,工作日
+2015/3/3,工作日
+2015/3/4,工作日
+2015/3/5,工作日
+2015/3/6,工作日
+2015/3/7,工作日
+2015/3/8,工作日
+2015/3/9,工作日
+2015/3/10,工作日
+2015/3/11,工作日
+2015/3/12,工作日
+2015/3/13,工作日
+2015/3/14,工作日
+2015/3/15,工作日
+2015/3/16,工作日
+2015/3/17,工作日
+2015/3/18,工作日
+2015/3/19,工作日
+2015/3/20,工作日
+2015/3/21,工作日
+2015/3/22,工作日
+2015/3/23,工作日
+2015/3/24,工作日
+2015/3/25,工作日
+2015/3/26,工作日
+2015/3/27,工作日
+2015/3/28,工作日
+2015/3/29,工作日
+2015/3/30,工作日
+2015/3/31,工作日
+2015/4/1,工作日
+2015/4/2,工作日
+2015/4/3,工作日
+2015/4/4,小长假
+2015/4/5,小长假
+2015/4/6,小长假
+2015/4/7,工作日
+2015/4/8,工作日
+2015/4/9,工作日
+2015/4/10,工作日
+2015/4/11,工作日
+2015/4/12,工作日
+2015/4/13,工作日
+2015/4/14,工作日
+2015/4/15,工作日
+2015/4/16,工作日
+2015/4/17,工作日
+2015/4/18,工作日
+2015/4/19,工作日
+2015/4/20,工作日
+2015/4/21,工作日
+2015/4/22,工作日
+2015/4/23,工作日
+2015/4/24,工作日
+2015/4/25,工作日
+2015/4/26,工作日
+2015/4/27,工作日
+2015/4/28,工作日
+2015/4/29,工作日
+2015/4/30,小长假
+2015/5/1,小长假
+2015/5/2,小长假
+2015/5/3,小长假
+2015/5/4,工作日
+2015/5/5,工作日
+2015/5/6,工作日
+2015/5/7,工作日
+2015/5/8,工作日
+2015/5/9,工作日
+2015/5/10,工作日
+2015/5/11,工作日
+2015/5/12,工作日
+2015/5/13,工作日
+2015/5/14,工作日
+2015/5/15,工作日
+2015/5/16,工作日
+2015/5/17,工作日
+2015/5/18,工作日
+2015/5/19,工作日
+2015/5/20,工作日
+2015/5/21,工作日
+2015/5/22,工作日
+2015/5/23,工作日
+2015/5/24,工作日
+2015/5/25,工作日
+2015/5/26,工作日
+2015/5/27,工作日
+2015/5/28,工作日
+2015/5/29,工作日
+2015/5/30,工作日
+2015/5/31,工作日
+2015/6/1,工作日
+2015/6/2,工作日
+2015/6/3,工作日
+2015/6/4,工作日
+2015/6/5,工作日
+2015/6/6,工作日
+2015/6/7,工作日
+2015/6/8,工作日
+2015/6/9,工作日
+2015/6/10,工作日
+2015/6/11,工作日
+2015/6/12,工作日
+2015/6/13,工作日
+2015/6/14,工作日
+2015/6/15,工作日
+2015/6/16,工作日
+2015/6/17,工作日
+2015/6/18,工作日
+2015/6/19,小长假
+2015/6/20,小长假
+2015/6/21,小长假
+2015/6/22,小长假
+2015/6/23,工作日
+2015/6/24,工作日
+2015/6/25,工作日
+2015/6/26,工作日
+2015/6/27,工作日
+2015/6/28,工作日
+2015/6/29,工作日
+2015/6/30,工作日
+2015/7/1,工作日
+2015/7/2,工作日
+2015/7/3,工作日
+2015/7/4,工作日
+2015/7/5,工作日
+2015/7/6,工作日
+2015/7/7,工作日
+2015/7/8,工作日
+2015/7/9,工作日
+2015/7/10,工作日
+2015/7/11,工作日
+2015/7/12,工作日
+2015/7/13,工作日
+2015/7/14,工作日
+2015/7/15,工作日
+2015/7/16,工作日
+2015/7/17,工作日
+2015/7/18,工作日
+2015/7/19,工作日
+2015/7/20,工作日
+2015/7/21,工作日
+2015/7/22,工作日
+2015/7/23,工作日
+2015/7/24,工作日
+2015/7/25,工作日
+2015/7/26,工作日
+2015/7/27,工作日
+2015/7/28,工作日
+2015/7/29,工作日
+2015/7/30,工作日
+2015/7/31,工作日
+2015/8/1,工作日
+2015/8/2,工作日
+2015/8/3,工作日
+2015/8/4,工作日
+2015/8/5,工作日
+2015/8/6,工作日
+2015/8/7,工作日
+2015/8/8,工作日
+2015/8/9,工作日
+2015/8/10,工作日
+2015/8/11,工作日
+2015/8/12,工作日
+2015/8/13,工作日
+2015/8/14,工作日
+2015/8/15,工作日
+2015/8/16,工作日
+2015/8/17,工作日
+2015/8/18,工作日
+2015/8/19,工作日
+2015/8/20,工作日
+2015/8/21,工作日
+2015/8/22,工作日
+2015/8/23,工作日
+2015/8/24,工作日
+2015/8/25,工作日
+2015/8/26,工作日
+2015/8/27,工作日
+2015/8/28,工作日
+2015/8/29,工作日
+2015/8/30,工作日
+2015/8/31,工作日
+2015/9/1,工作日
+2015/9/2,小长假
+2015/9/3,小长假
+2015/9/4,小长假
+2015/9/5,小长假
+2015/9/6,工作日
+2015/9/7,工作日
+2015/9/8,工作日
+2015/9/9,工作日
+2015/9/10,工作日
+2015/9/11,工作日
+2015/9/12,工作日
+2015/9/13,工作日
+2015/9/14,工作日
+2015/9/15,工作日
+2015/9/16,工作日
+2015/9/17,工作日
+2015/9/18,工作日
+2015/9/19,工作日
+2015/9/20,工作日
+2015/9/21,工作日
+2015/9/22,工作日
+2015/9/23,工作日
+2015/9/24,工作日
+2015/9/25,工作日
+2015/9/26,小长假
+2015/9/27,小长假
+2015/9/28,工作日
+2015/9/29,工作日
+2015/9/30,小长假
+2015/10/1,小长假
+2015/10/2,小长假
+2015/10/3,小长假
+2015/10/4,小长假
+2015/10/5,小长假
+2015/10/6,小长假
+2015/10/7,小长假
+2015/10/8,工作日
+2015/10/9,工作日
+2015/10/10,工作日
+2015/10/11,工作日
+2015/10/12,工作日
+2015/10/13,工作日
+2015/10/14,工作日
+2015/10/15,工作日
+2015/10/16,工作日
+2015/10/17,工作日
+2015/10/18,工作日
+2015/10/19,工作日
+2015/10/20,工作日
+2015/10/21,工作日
+2015/10/22,工作日
+2015/10/23,工作日
+2015/10/24,工作日
+2015/10/25,工作日
+2015/10/26,工作日
+2015/10/27,工作日
+2015/10/28,工作日
+2015/10/29,工作日
+2015/10/30,工作日
+2015/10/31,工作日
+2015/11/1,工作日
+2015/11/2,工作日
+2015/11/3,工作日
+2015/11/4,工作日
+2015/11/5,工作日
+2015/11/6,工作日
+2015/11/7,工作日
+2015/11/8,工作日
+2015/11/9,工作日
+2015/11/10,工作日
+2015/11/11,工作日
+2015/11/12,工作日
+2015/11/13,工作日
+2015/11/14,工作日
+2015/11/15,工作日
+2015/11/16,工作日
+2015/11/17,工作日
+2015/11/18,工作日
+2015/11/19,工作日
+2015/11/20,工作日
+2015/11/21,工作日
+2015/11/22,工作日
+2015/11/23,工作日
+2015/11/24,工作日
+2015/11/25,工作日
+2015/11/26,工作日
+2015/11/27,工作日
+2015/11/28,工作日
+2015/11/29,工作日
+2015/11/30,工作日
+2015/12/1,工作日
+2015/12/2,工作日
+2015/12/3,工作日
+2015/12/4,工作日
+2015/12/5,工作日
+2015/12/6,工作日
+2015/12/7,工作日
+2015/12/8,工作日
+2015/12/9,工作日
+2015/12/10,工作日
+2015/12/11,工作日
+2015/12/12,工作日
+2015/12/13,工作日
+2015/12/14,工作日
+2015/12/15,工作日
+2015/12/16,工作日
+2015/12/17,工作日
+2015/12/18,工作日
+2015/12/19,工作日
+2015/12/20,工作日
+2015/12/21,工作日
+2015/12/22,工作日
+2015/12/23,工作日
+2015/12/24,工作日
+2015/12/25,工作日
+2015/12/26,工作日
+2015/12/27,工作日
+2015/12/28,工作日
+2015/12/29,工作日
+2015/12/30,工作日
+2015/12/31,小长假
+2016/1/1,小长假
+2016/1/2,小长假
+2016/1/3,小长假
+2016/1/4,工作日
+2016/1/5,工作日
+2016/1/6,工作日
+2016/1/7,工作日
+2016/1/8,工作日
+2016/1/9,工作日
+2016/1/10,工作日
+2016/1/11,工作日
+2016/1/12,工作日
+2016/1/13,工作日
+2016/1/14,工作日
+2016/1/15,工作日
+2016/1/16,工作日
+2016/1/17,工作日
+2016/1/18,工作日
+2016/1/19,工作日
+2016/1/20,工作日
+2016/1/21,工作日
+2016/1/22,工作日
+2016/1/23,工作日
+2016/1/24,工作日
+2016/1/25,工作日
+2016/1/26,工作日
+2016/1/27,工作日
+2016/1/28,工作日
+2016/1/29,工作日
+2016/1/30,工作日
+2016/1/31,工作日
+2016/2/1,工作日
+2016/2/2,工作日
+2016/2/3,工作日
+2016/2/4,工作日
+2016/2/5,工作日
+2016/2/6,工作日
+2016/2/7,小长假
+2016/2/8,小长假
+2016/2/9,小长假
+2016/2/10,小长假
+2016/2/11,小长假
+2016/2/12,小长假
+2016/2/13,小长假
+2016/2/14,工作日
+2016/2/15,工作日
+2016/2/16,工作日
+2016/2/17,工作日
+2016/2/18,工作日
+2016/2/19,工作日
+2016/2/20,工作日
+2016/2/21,工作日
+2016/2/22,工作日
+2016/2/23,工作日
+2016/2/24,工作日
+2016/2/25,工作日
+2016/2/26,工作日
+2016/2/27,工作日
+2016/2/28,工作日
+2016/2/29,工作日
+2016/3/1,工作日
+2016/3/2,工作日
+2016/3/3,工作日
+2016/3/4,工作日
+2016/3/5,工作日
+2016/3/6,工作日
+2016/3/7,工作日
+2016/3/8,工作日
+2016/3/9,工作日
+2016/3/10,工作日
+2016/3/11,工作日
+2016/3/12,工作日
+2016/3/13,工作日
+2016/3/14,工作日
+2016/3/15,工作日
+2016/3/16,工作日
+2016/3/17,工作日
+2016/3/18,工作日
+2016/3/19,工作日
+2016/3/20,工作日