cormort/DataAnalysis
0
1import os2 3# Set Streamlit's home directory to a writable location4# This ensures Streamlit internal files are written to /tmp5os.environ['STREAMLIT_HOME'] = '/tmp' 6 7# Existing environment variables for Streamlit8os.environ['STREAMLIT_SERVER_HEADLESS'] = 'true'9os.environ['STREAMLIT_BROWSER_GATHERUSAGETSTATS'] = 'false'10os.environ['STREAMLIT_SERVER_ENABLE_STATIC_SERVING'] = 'true' 11os.environ['STREAMLIT_SERVER_FOLDER_WATCHER_TYPE'] = 'inotify' 12 13import streamlit as st14import pandas as pd15import numpy as np16import re17from scipy.stats import gmean18from sklearn.linear_model import LinearRegression19from sklearn.metrics import mean_absolute_error, mean_squared_error20import plotly.graph_objects as go21import matplotlib.pyplot as plt22 23# --- New Imports for Advanced Features ---24from statsmodels.tsa.seasonal import STL25from statsmodels.tsa.stattools import adfuller26from statsmodels.graphics.tsaplots import plot_acf27from pmdarima import auto_arima 28 29 30# --- 1. 定義預測方法及其描述 ---31METHOD_DEFINITIONS = {32 '平均 (近3年)': {'definition': '計算最近 3 年數據的平均值', 'formula': '$mean(x_{t-2}, x_{t-1}, x_t)$'},33 '平均 (近4年)': {'definition': '計算最近 4 年數據的平均值', 'formula': '$mean(x_{t-3}, ..., x_t)$'},34 '平均 (近5年)': {'definition': '計算最近 5 年數據的平均值', 'formula': '$mean(x_{t-4}, ..., x_t)$'},35 '平均 (值)': {'definition': '計算所有數據的平均值', 'formula': '$mean(x_1, ..., x_t)$'},36 '平均 (去極值)': {'definition': '移除最大最小值後的平均值', 'formula': '$mean(x_1, ..., x_t \\setminus \\{max, min\\})$'},37 '幾何平均 (近3年)': {'definition': '最近 3 年數據的幾何平均值', 'formula': '$(x_{t-2} \\cdot x_{t-1} \\cdot x_t)^{1/3}$'},38 '幾何平均 (近4年)': {'definition': '最近 4 年數據的幾何平均值', 'formula': '$(x_{t-3} \\cdot ... \\cdot x_t)^{1/4}$'},39 '幾何平均 (近5年)': {'definition': '最近 5 年數據的幾何平均值', 'formula': '$(x_{t-4} \\cdot ... \\cdot x_t)^{1/5}$'},40 '幾何平均 (值)': {'definition': '所有數據的幾何平均值（負值或零視為 1）', 'formula': '$(x_1 \\cdot ... \\cdot x_t)^{1/t}$'},41 '幾何平均 (去極值)': {'definition': '移除最大最小值後的幾何平均值', 'formula': '$(x_1 \\cdot ... \\cdot x_t \\setminus \\{max, min\\})^{1/(t-2)}$'},42 '最大值': {'definition': '所有數據的最大值', 'formula': '$max(x_1, ..., x_t)$'},43 '最小值': {'definition': '所有數據的最小值', 'formula': '$min(x_1, ..., x_t)$'},44 '指數平滑 (α=0.3)': {'definition': '簡單指數平滑，α=0.3', 'formula': '$S_{t+1} = αx_t + (1-α)S_t, S_1=x_1$'},45 '指數平滑 (α=0.4)': {'definition': '簡單指數平滑，α=0.4', 'formula': '$S_{t+1} = αx_t + (1-α)S_t, S_1=x_1$'},46 '指數平滑 (α=0.5)': {'definition': '簡單指數平滑，α=0.5', 'formula': '$S_{t+1} = αx_t + (1-α)S_t, S_1=x_1$'},47 '線性趨勢預測（OLS）': {'definition': '基於最小平方法(OLS)計算的線性趨勢線，預測下一期', 'formula': '$Y_{t+1} = a + b(t+1)$'},48 'Holt 線性趨勢法': {'definition': 'Holt線性趨勢法 (雙指數平滑 α=0.5, β=0.3)', 'formula': '$L_t=αY_t+(1-α)(L_{t-1}+B_{t-1}); B_t=β(L_t-L_{t-1})+(1-β)B_{t-1}; F_{t+1}=L_t+B_t$'},49 '簡化 ARIMA（差分 + 平均變動）': {'definition': '一階差分序列的平均值，作為下一期的預期變動量', 'formula': '$Y_{t+1} = Y_t + mean(Y_i - Y_{i-1})$'},50 'Auto ARIMA': {'definition': '自動搜尋最佳 (p,d,q) 參數的 ARIMA 模型', 'formula': 'Auto-selected ARIMA(p,d,q)'}51}52 53# --- 2. 數據解析函數 (Python) ---54def parse_data(raw_input):55 processed_input = raw_input.replace('\t', ',')56 raw_rows = [r.strip() for r in re.split(r'[\n\r;]+', processed_input) if r.strip()]57 paired_data_list, value_only_list = [], []58 has_paired_data, has_value_only_data = False, False59 60 if not raw_rows: return pd.DataFrame()61 62 for row_string in raw_rows:63 parts = [p.strip() for p in row_string.split(',') if p.strip()]64 if len(parts) == 2:65 try:66 year_num = int(parts[0]); value_num = float(parts[1])67 if not np.isnan(value_num):68 paired_data_list.append({'year': year_num, 'value': value_num}); has_paired_data = True; continue69 except (ValueError, IndexError): pass70 try:71 comma_values_in_row = [float(v) for v in parts]72 if all(not np.isnan(v) for v in comma_values_in_row):73 value_only_list.extend(comma_values_in_row); has_value_only_data = True74 except ValueError: pass75 76 if has_paired_data and has_value_only_data:77 st.warning("偵測到混合數據格式。將優先處理 `年份,數值` 的成對數據。")78 return pd.DataFrame(paired_data_list).sort_values(by='year').reset_index(drop=True)79 if has_paired_data:80 return pd.DataFrame(paired_data_list).sort_values(by='year').reset_index(drop=True)81 if has_value_only_data:82 if st.session_state.input_start_year:83 start_year = st.session_state.input_start_year84 data_with_years = [{'year': start_year + i, 'value': val} for i, val in enumerate(value_only_list)]85 return pd.DataFrame(data_with_years)86 else:87 st.info("偵測到純數值序列。請在側邊欄輸入【起始年份】後，再點擊'讀取與分析'。")88 return pd.DataFrame()89 return pd.DataFrame()90 91# --- 3. 統計與預測方法函數 (Python) ---92# Helper functions93def get_values_no_extremes(values):94 return sorted(values)[1:-1] if len(values) > 2 else []95def geometric_mean_python(values):96 positive_values = [v for v in values if v > 0]97 return gmean(positive_values) if positive_values else 0.098 99# Prediction methods100def exponential_smoothing(values, alpha):101 if not values: return np.nan102 smoothed = values[0]103 for i in range(1, len(values)): smoothed = alpha * values[i] + (1 - alpha) * smoothed104 return smoothed105def linear_trend_forecast(values):106 if len(values) < 2: return np.nan107 X = np.arange(1, len(values) + 1).reshape(-1, 1); y = np.array(values)108 model = LinearRegression().fit(X, y)109 return model.predict(np.array([[len(values) + 1]]))[0]110def holt_linear_trend_forecast(values, alpha=0.5, beta=0.3):111 if len(values) < 2: return np.nan112 level, trend = values[0], values[1] - values[0]113 for i in range(1, len(values)):114 last_level = level115 level = alpha * values[i] + (1 - alpha) * (last_level + trend)116 trend = beta * (level - last_level) + (1 - beta) * trend117 return level + trend118def simple_arima_forecast(values):119 if len(values) < 2: return np.nan120 return values[-1] + np.mean(np.diff(values))121def auto_arima_forecast(values):122 if len(values) < 4: return np.nan123 try:124 model = auto_arima(values, seasonal=False, suppress_warnings=True, error_action='ignore')125 return model.predict(n_periods=1)[0]126 except:127 return np.nan128 129# --- 4. 主計算函數 ---130ALL_METHODS = {131 '平均 (近3年)': lambda v: np.mean(v[-3:]) if len(v) >= 3 else np.nan,132 '平均 (近4年)': lambda v: np.mean(v[-4:]) if len(v) >= 4 else np.nan,133 '平均 (近5年)': lambda v: np.mean(v[-5:]) if len(v) >= 5 else np.nan,134 '平均 (值)': lambda v: np.mean(v),135 '平均 (去極值)': lambda v: np.mean(get_values_no_extremes(v)),136 '幾何平均 (近3年)': lambda v: geometric_mean_python(v[-3:]) if len(v) >= 3 else np.nan,137 '幾何平均 (近4年)': lambda v: geometric_mean_python(v[-4:]) if len(v) >= 4 else np.nan,138 '幾何平均 (近5年)': lambda v: geometric_mean_python(v[-5:]) if len(v) >= 5 else np.nan,139 '幾何平均 (值)': geometric_mean_python,140 '幾何平均 (去極值)': lambda v: geometric_mean_python(get_values_no_extremes(v)),141 '最大值': lambda v: np.max(v),142 '最小值': lambda v: np.min(v),143 '指數平滑 (α=0.3)': lambda v: exponential_smoothing(v, 0.3),144 '指數平滑 (α=0.4)': lambda v: exponential_smoothing(v, 0.4),145 '指數平滑 (α=0.5)': lambda v: exponential_smoothing(v, 0.5),146 '線性趨勢預測（OLS）': linear_trend_forecast,147 'Holt 線性趨勢法': holt_linear_trend_forecast,148 '簡化 ARIMA（差分 + 平均變動）': simple_arima_forecast,149 'Auto ARIMA': auto_arima_forecast150}151 152def calculate_all_predictions(values):153 predictions = []154 for name, func in ALL_METHODS.items():155 try:156 pred_val = func(values)157 predictions.append({'method': name, 'value': pred_val})158 except:159 predictions.append({'method': name, 'value': np.nan})160 return pd.DataFrame(predictions)161 162def run_model_evaluation(values, test_size):163 train_data = values[:-test_size]164 test_data = values[-test_size:]165 evaluation_results = []166 167 for name, func in ALL_METHODS.items():168 try:169 predictions = []170 current_train = list(train_data)171 for _ in range(test_size):172 pred = func(current_train)173 if pd.isna(pred): # If a model fails, stop its evaluation174 break175 predictions.append(pred)176 current_train.append(test_data[len(predictions)-1])177 178 if len(predictions) == test_size:179 mae = mean_absolute_error(test_data, predictions)180 rmse = np.sqrt(mean_squared_error(test_data, predictions))181 evaluation_results.append({'方法': name, 'MAE': mae, 'RMSE': rmse})182 except:183 continue184 return pd.DataFrame(evaluation_results).sort_values(by='RMSE')185 186# --- 5. 繪圖函數 ---187@st.cache_data188def plot_decomposition(df, period):189 df_for_plot = df.set_index('year')190 stl = STL(df_for_plot['value'], period=period, robust=True)191 result = stl.fit()192 fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(10, 8), sharex=True)193 result.observed.plot(ax=ax1, legend=False); ax1.set_ylabel('Observed')194 result.trend.plot(ax=ax2, legend=False); ax2.set_ylabel('Trend')195 result.seasonal.plot(ax=ax3, legend=False); ax3.set_ylabel('Seasonal')196 result.resid.plot(ax=ax4, legend=False); ax4.set_ylabel('Residual')197 plt.tight_layout(); return fig198@st.cache_data199def plot_acf_chart(values, lags):200 fig, ax = plt.subplots(figsize=(10, 4))201 plot_acf(values, lags=lags, ax=ax)202 plt.tight_layout(); return fig203def plot_forecast_chart(df, method_name, predicted_value):204 fig = go.Figure()205 fig.add_trace(go.Bar(x=df['year'], y=df['value'], name='歷史數據', marker_color='royalblue'))206 next_year = df['year'].max() + 1207 fig.add_trace(go.Bar(x=[next_year], y=[predicted_value], name='預測值', marker_color='crimson'))208 fig.update_layout(209 title=f"圖表: {method_name}",210 xaxis_title='年份', yaxis_title='數值', legend_title="數據類型",211 xaxis=dict(212 tickmode='array',213 tickvals=df['year'].tolist() + [next_year],214 ticktext=[str(y) for y in df['year']] + [f"<b>{next_year}<br>(預測)</b>"]215 )216 )217 return fig218 219# --- Streamlit UI ---220st.set_page_config(layout="wide", page_title="時間序列分析與預測")221st.title("📊 時間序列分析與預測工具")222 223# --- Session State ---224if 'data' not in st.session_state: st.session_state.data = pd.DataFrame()225if 'predictions' not in st.session_state: st.session_state.predictions = pd.DataFrame()226if 'evaluation' not in st.session_state: st.session_state.evaluation = pd.DataFrame()227 228# --- Sidebar ---229with st.sidebar:230 st.header("Step 1: 輸入數據")231 data_input = st.text_area(232 "貼上數據 (Excel, CSV, etc.):\n- 格式一: `年份,數值`\n- 格式二: `數值,數值...`",233 "2015,100\n2016,120\n2017,130\n2018,110\n2019,140\n2020,150\n2021,160\n2022,180\n2023,170",234 height=200, key="data_input_area"235 )236 237 st.header("Step 2: 數據選項")238 st.session_state.input_start_year = st.number_input(239 "若無年份，請指定起始年份:",240 min_value=1, max_value=3000, value=2015, format="%d"241 )242 243 if st.button("讀取與分析資料", type="primary"):244 st.session_state.data = parse_data(data_input)245 if not st.session_state.data.empty:246 st.session_state.predictions = calculate_all_predictions(st.session_state.data['value'].tolist())247 st.session_state.evaluation = pd.DataFrame() # Clear old evaluation248 else:249 st.session_state.predictions = pd.DataFrame()250 251# --- Main Panel ---252if st.session_state.data.empty:253 st.info("⬅️ 請在左側側邊欄輸入數據，然後點擊 **'讀取與分析資料'**。")254else:255 tab1, tab2, tab3 = st.tabs(["📈 預測與分析", "🔬 數據探索", "⚖️ 模型評估"])256 257 with tab1:258 st.header("預測結果總覽")259 if not st.session_state.predictions.empty:260 df_pred = st.session_state.predictions.copy()261 df_pred['定義'] = df_pred['method'].map(lambda m: METHOD_DEFINITIONS.get(m, {}).get('definition', ''))262 df_pred['公式'] = df_pred['method'].map(lambda m: METHOD_DEFINITIONS.get(m, {}).get('formula', ''))263 st.dataframe(264 df_pred.rename(columns={'method': '方法', 'value': '預測值'})[['方法', '預測值', '定義', '公式']],265 hide_index=True, use_container_width=True,266 column_config={"預測值": st.column_config.NumberColumn(format="%.2f")}267 )268 else:269 st.info("點擊側邊欄按鈕以產生預測。")270 271 st.header("預測圖表")272 if not st.session_state.predictions.empty:273 chart_methods = ['(不顯示)'] + st.session_state.predictions['method'].tolist()274 selected_method = st.selectbox("選擇要顯示圖表的預測方法:", options=chart_methods)275 276 if selected_method != '(不顯示)':277 pred_row = st.session_state.predictions[st.session_state.predictions['method'] == selected_method].iloc[0]278 if pd.notna(pred_row['value']):279 fig = plot_forecast_chart(st.session_state.data, pred_row['method'], pred_row['value'])280 st.plotly_chart(fig, use_container_width=True)281 else:282 st.warning(f"方法 '{selected_method}' 無法計算預測值，因此無法繪製圖表。")283 284 with tab2:285 st.header("數據探索分析")286 st.info("此區塊幫助您深入了解數據的內在特性，如趨勢、季節性和相關性。")287 288 st.subheader("1. 時間序列分解")289 period = st.number_input("請輸入季節性週期", min_value=2, value=max(2, len(st.session_state.data)//4))290 if st.button("顯示分解圖"):291 if len(st.session_state.data) >= 2 * period:292 fig = plot_decomposition(st.session_state.data.copy(), period)293 st.pyplot(fig)294 else:295 st.warning(f"數據點不足，無法進行週期為 {period} 的分解。至少需要 {2 * period} 個數據點。")296 297 st.subheader("2. 平穩性檢驗 (ADF Test)")298 if st.button("執行 ADF 檢驗"):299 values = st.session_state.data['value'].dropna()300 result = adfuller(values)301 st.write(f"ADF 統計量: **{result[0]:.4f}**")302 st.write(f"P-value: **{result[1]:.4f}**")303 if result[1] <= 0.05:304 st.success("結論: 數據很可能是**平穩的** (因為 P-value ≤ 0.05)。")305 else:306 st.error("結論: 數據很可能是**非平穩的** (因為 P-value > 0.05)。")307 308 st.subheader("3. 自相關圖 (ACF)")309 max_lags = len(st.session_state.data) // 2 - 1310 if max_lags > 0:311 lags = st.slider("選擇最大延遲期數 (Lags)", 1, max_lags, min(10, max_lags))312 if st.button("顯示 ACF 圖"):313 fig = plot_acf_chart(st.session_state.data['value'], lags)314 st.pyplot(fig)315 else:316 st.warning("數據點不足，無法繪製 ACF 圖。")317 318 with tab3:319 st.header("模型評估")320 st.info("此功能會將您的數據分為訓練集和測試集，透過在測試集上的表現來評估各個模型的預測準確度。RMSE/MAE 越小，代表模型表現越好。")321 max_test_size = len(st.session_state.data) // 2322 if max_test_size > 0:323 test_size = st.slider("選擇測試集大小 (最後 N 個數據點)", 1, max_test_size, max(1, len(st.session_state.data)//5))324 325 if st.button("評估所有模型", type="primary"):326 with st.spinner("正在執行滾動預測與評估，請稍候..."):327 st.session_state.evaluation = run_model_evaluation(st.session_state.data['value'].tolist(), test_size)328 329 if not st.session_state.evaluation.empty:330 st.subheader("評估結果")331 st.dataframe(332 st.session_state.evaluation.style.highlight_min(color='#D4EDDA', subset=['MAE', 'RMSE']),333 use_container_width=True, hide_index=True334 )335 else:336 st.info("點擊上方按鈕開始評估。")337 else:338 st.warning("數據點太少，無法進行模型評估。請提供更多數據。")339