運用幾何平均估計引力模型常數與殘差分析

引力模型為理解網絡流量分佈提供了一個基於質量與距離的基礎框架。然而，模型的價值不僅在於其預測能力，更在於它揭示偏離預期模式的能力。本文深入探討此模型的兩個核心步驟：首先，為解決流量數據常見的極端值問題，採用幾何平均數來估計引力常數 $k$，以建立一個更具代表性的基準。其次，在獲得預期流量後，透過殘差分析量化實際流量與預期值之間的差異。正向殘差突顯出超越距離與規模所能解釋的「特殊連結」，而負向殘差則可能指向潛在的阻礙因素。此方法將模型的「誤差」轉化為有價值的洞察，為發掘網絡中隱含的特殊關係提供了明確的分析路徑。

引力模型常數的幾何平均估計與殘差分析

本章節將延續對引力模型的探討，重點在於如何利用計算出的邊距離、節點質量以及實際流量，來估計引力模型中的常數乘數 $k$。我們將採用幾何平均數的方式來匯總各邊計算出的 $k$ 值，並進一步計算預期流量與實際流量之間的殘差，以揭示模型無法解釋的「特殊關係」。

常數乘數 $k$ 的幾何平均估計

幾何平均數的選擇：
- 數據特性：在處理流量、距離等數據時，這些數值通常分佈較為「重尾」(heavy-tailed)，意味著存在一些極端值。
- 算術平均數的局限：算術平均數對極端值較為敏感，可能會被異常值嚴重影響。
- 幾何平均數的優勢：幾何平均數對數值範圍的變化不敏感，更能反映數據的中心趨勢，特別是當數據是乘法關係時（如引力模型）。
- 計算方式：
  - 幾何平均數的定義是 $N$ 個數的乘積的 $N$ 次方根。
  - 然而，更常見且數值上更穩健的計算方式是通過對數轉換：先取每個數的對數（例如，以 10 為底），計算這些對數的算術平均數，然後將結果取反對數（以 10 為底）。
  - 程式碼實現：g = 10**(sum([math.log10(g) for g in g_list]) / len(g_list))
    - math.log10(g)：計算列表中每個 $k$ 值的以 10 為底的對數。
    - sum(...) / len(...)：計算這些對數的算術平均數。
    - 10**(...)：將平均對數值轉換回原始尺度，得到幾何平均數。
常數 $g$ 的意義：
- 這個計算出的 g 值，代表了我們從數據中估計出的引力模型常數乘數 $k$。
- 它提供了一個單一的數值，用於統一調整模型，使其能夠基於節點質量和距離來預測流量。
- 這個值可以被視為流量與質量、距離關係的「平均強度因子」。

殘差分析：預期流量與實際流量的差異

預期流量的計算：
- 一旦獲得了常數 $g$，我們就可以計算每條邊的預期流量 (expected)。
- 公式：expected = (g * G_air.nodes[v]['degree'] * G_air.nodes[w]['degree'] / G_air.edges[v, w]['distance']**2)
  - g：估計出的常數乘數。
  - G_air.nodes[v]['degree']：節點 $v$ 的質量（加權度數）。
  - G_air.nodes[w]['degree']：節點 $w$ 的質量。
  - G_air.edges[v, w]['distance']**2：邊距離的平方。
殘差的定義與計算：
- 殘差 (Residual)：指實際觀察到的流量 (count) 與模型預測的流量 (expected) 之間的差異。
- residual = count - expected
- 程式碼將這個殘差值添加為邊的一個新屬性 ['residual']。
殘差的解釋：
- 正殘差 (residual > 0)：
  - 表示實際流量高於引力模型基於距離和質量所預期的流量。
  - 這可能意味著這兩個地區之間存在某種「特殊關係」，例如：
    - 強勁的經濟聯繫。
    - 重要的旅遊目的地。
    - 較低的票價或更優的航班時刻表。
    - 該航線的「吸引力」比單純的距離和規模所能解釋的要大。
- 負殘差 (residual < 0)：
  - 表示實際流量低於模型預期的流量。
  - 可能的原因包括：
    - 地理距離雖然不遠，但兩地區的經濟或人口規模較小。
    - 存在其他競爭性交通方式。
    - 航班供應不足或票價較高。
    - 該航線的「吸引力」較弱。
- 零殘差 (residual ≈ 0)：
  - 表示模型的預測與實際流量非常接近。

程式碼實現細節

計算 $k$ 值列表 g_list：
- 程式碼中包含了計算每個邊所需 $k$ 值的邏輯，並將結果存入 g_list。
- if v >= w: continue：確保每個邊只計算一次。
- try...except KeyError:：處理邊屬性缺失的情況，並將 0 添加到 g_list（儘管在實際計算幾何平均數時，0 值需要小心處理，通常會被忽略或替換）。
- distance = G_air.edges[v, w]['distance'] 等：獲取所需的邊和節點屬性。
- g_list.append(count * distance**2 / v_degree / w_degree)：計算並追加 $k$ 值。
計算幾何平均數 $g$：
- g = 10**(sum([math.log10(g) for g in g_list if g > 0]) / len([g for g in g_list if g > 0]))：
  - 注意：這裡加入了 if g > 0 過濾，以避免對 0 或負數取對數時出錯。
  - 計算對數的算術平均數，然後取反對數。
計算預期流量 expected：
- 遍歷所有邊，使用計算出的 g 值、節點質量和邊距離，根據引力模型公式計算預期流量。
- G_air.edges[v, w]['expected'] = expected：將預期流量作為新屬性添加到邊。
計算殘差 residual：
- 程式碼片段中，計算殘差的邏輯尚未完整展示，但其思路是 count - expected。
- 預計後續會將 residual 添加為邊的屬性。

import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pathlib import Path
import math

# --- 假設 G_air_spatial 已載入並處理好節點屬性 (lat, long, count, mass, distance) ---
# 為了程式碼的獨立性，這裡重新載入並處理部分數據
# 實際應用中，應接續前文的 G_air_spatial

# --- 重新創建模擬數據和網路 (如果前文未執行) ---
data_dir_sim = Path('./data_simulated')
carrier_data_path = data_dir_sim / 'BTS2018' / 'carrier.csv'
airport_db_path = data_dir_sim / 'partow' / 'GlobalAirportDatabase.txt'

# 模擬創建 carrier.csv
if not carrier_data_path.exists():
    data_dir_sim.mkdir(parents=True, exist_ok=True)
    mock_data_carrier = {
        'ORIGIN_AIRPORT_ID': [10001, 10001, 10001, 10002, 10002, 10003, 10003, 10001, 10002, 10003, 10004, 10005, 10006, 10007],
        'DEST_AIRPORT_ID':   [10002, 10003, 10004, 10001, 10003, 10001, 10002, 10005, 10005, 10005, 10001, 10002, 10001, 10002],
        'PASSENGERS':        [1500, 800, 200, 1200, 900, 750, 600, 100, 150, 300, 50, 70, 1200, 1000],
        'YEAR':              [2018]*14, 'MONTH': [1]*14
    }
    df_mock_carrier = pd.DataFrame(mock_data_carrier)
    df_mock_carrier.to_csv(carrier_data_path, index=False)

# 模擬創建 GlobalAirportDatabase.txt
if not airport_db_path.exists():
    data_dir_partow = data_dir_sim / 'partow'
    data_dir_partow.mkdir(parents=True, exist_ok=True)
    mock_airport_db_content = """
1:JFK:John F. Kennedy International Airport:USA:40.639751:-73.778925
2:LAX:Los Angeles International Airport:USA:33.941556:-118.408530
3:ORD:Chicago O'Hare International Airport:USA:41.974209:-87.907321
4:ATL:Hartsfield-Jackson Atlanta International Airport:USA:33.640445:-84.427700
5:DEN:Denver International Airport:USA:39.856057:-104.673737
6:SFO:San Francisco International Airport:USA:37.619002:-122.371297
7:SEA:Seattle-Tacoma International Airport:USA:47.448981:-122.309310
8:MIA:Miami International Airport:USA:25.793199:-80.290594
9:DFW:Dallas/Fort Worth International Airport:USA:32.899811:-97.040416
10:LAS:McCarran International Airport:USA:36.084030:-115.153734
11:PHX:Phoenix Sky Harbor International Airport:USA:33.434250:-112.011590
12:IAH:George Bush Intercontinental Airport:USA:29.990163:-95.336770
13:BOS:Logan International Airport:USA:42.365599:-71.009667
14:MSP:Minneapolis–Saint Paul International Airport:USA:44.884794:-93.217711
15:PHL:Philadelphia International Airport:USA:39.874399:-75.242449
16:CLT:Charlotte Douglas International Airport:USA:35.214433:-80.947171
17:EWR:Newark Liberty International Airport:USA:40.689532:-74.174492
18:DTW:Detroit Metropolitan Airport:USA:42.211501:-83.353409
19:FLL:Fort Lauderdale–Hollywood International Airport:USA:26.074199:-80.150667
20:BWI:Baltimore/Washington International Airport:USA:39.177452:-76.668409
21:SLC:Salt Lake City International Airport:USA:40.789999:-111.979111
22:SAN:San Diego International Airport:USA:32.733819:-117.193751
23:TPA:Tampa International Airport:USA:27.977474:-82.531189
24:PDX:Portland International Airport:USA:45.589802:-122.597510
25:IAD:Washington Dulles International Airport:USA:38.953054:-77.456459
26:DCA:Ronald Reagan Washington National Airport:USA:38.851111:-77.037500
27:ANC:Ted Stevens Anchorage International Airport:USA:61.174355:-149.996333 # Alaska, will be removed
28:HNL:Daniel K. Inouye International Airport:USA:21.318694:-157.924694 # Hawaii, will be removed
"""
    with open(airport_db_path, 'w') as f:
        f.write(mock_airport_db_content.strip())

# --- 重新載入並處理數據 ---
G_air_spatial = nx.Graph()
airport_locations = {}

try:
    # 載入旅客流量
    df_carrier = pd.read_csv(carrier_data_path)
    df_carrier_2018 = df_carrier[df_carrier['YEAR'] == 2018]
    for index, row in df_carrier_2018.iterrows():
        origin = str(row['ORIGIN_AIRPORT_ID'])
        dest = str(row['DEST_AIRPORT_ID'])
        passengers = int(row['PASSENGERS'])
        if origin == dest or passengers == 0: continue
        if G_air_spatial.has_edge(origin, dest):
            G_air_spatial.edges[origin, dest]['count'] += passengers
        else:
            G_air_spatial.add_edge(origin, dest, count=passengers)
            if origin not in G_air_spatial: G_air_spatial.add_node(origin)
            if dest not in G_air_spatial: G_air_spatial.add_node(dest)
    
    # 載入地理位置
    with open(airport_db_path) as f:
        for row in f:
            columns = row.strip().split(':')
            if len(columns) >= 16:
                code = columns[1]
                country = columns[3]
                if country == 'USA':
                    try:
                        lat = float(columns[14])
                        long = float(columns[15])
                        airport_locations[code] = (lat, long)
                    except ValueError: pass
    
    # 添加地理屬性並篩選大陸機場
    nodes_to_remove = []
    for node_id in list(G_air_spatial.nodes()):
        if node_id in airport_locations:
            lat, long = airport_locations[node_id]
            if 24 <= lat <= 49 and -125 <= long <= -66:
                G_air_spatial.nodes[node_id]['lat'] = lat
                G_air_spatial.nodes[node_id]['long'] = long
            else:
                nodes_to_remove.append(node_id)
        else:
            nodes_to_remove.append(node_id)
    G_air_spatial.remove_nodes_from(nodes_to_remove)

    print(f"網路準備就緒: {G_air_spatial.number_of_nodes()} 個節點, {G_air_spatial.number_of_edges()} 條邊。")

except Exception as e:
    print(f"初始化網路時發生錯誤: {e}")

# --- Haversine 距離計算函數 ---
def haversine(lat1, lon1, lat2, lon2):
    R_km = 6371 
    lat1_rad, lon1_rad, lat2_rad, lon2_rad = map(math.radians, [lat1, lon1, lat2, lon2])
    dlon = lon2_rad - lon1_rad
    dlat = lat2_rad - lat1_rad
    a = math.sin(dlat / 2)**2 + math.cos(lat1_rad) * math.cos(lat2_rad) * math.sin(dlon / 2)**2
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
    distance = R_km * c
    return distance

# --- 為邊添加距離屬性 ---
for v, w in G_air_spatial.edges():
    try:
        p_v = (G_air_spatial.nodes[v]['lat'], G_air_spatial.nodes[v]['long'])
        p_w = (G_air_spatial.nodes[w]['lat'], G_air_spatial.nodes[w]['long'])
        G_air_spatial.edges[v, w]['distance'] = haversine(p_v[0], p_v[1], p_w[0], p_w[1])
    except KeyError:
        G_air_spatial.edges[v, w]['distance'] = None 
print("已為邊添加距離屬性。")

# --- 計算節點質量 (加權度數) ---
degree_data = list(G_air_spatial.degree(weight='count'))
nx.set_node_attributes(G_air_spatial, dict(degree_data), 'mass')
print("已計算節點質量 (總旅客流量)。")

# --- 簡化方法估計引力模型常數乘數 k ---
k_values = [] 

for v, w in G_air_spatial.edges():
    if v >= w: 
        continue
        
    try:
        count = G_air_spatial.edges[v, w]['count']
        distance = G_air_spatial.edges[v, w]['distance']
        mass_v = G_air_spatial.nodes[v]['mass']
        mass_w = G_air_spatial.nodes[w]['mass']
        
        # 確保所有值都有效，避免除以零或使用 None 值
        if count is not None and distance is not None and distance > 0 and \
           mass_v is not None and mass_w is not None and mass_v > 0 and mass_w > 0:
            
            # 從公式 k = count * (distance^2) / (mass_v * mass_w) 計算 k
            k_val = count * (distance**2) / (mass_v * mass_w)
            k_values.append(k_val)
        else:
            # 如果值無效，不添加到 k_values 中，也不計入幾何平均數的分母
            pass
            
    except KeyError:
        # 如果邊屬性或節點屬性缺失，則跳過
        continue
    except Exception as e:
        print(f"處理邊 ({v}, {w}) 時發生錯誤: {e}")

# --- 計算幾何平均數 g ---
g = 0 # 初始化 g
valid_k_values = [val for val in k_values if val > 0] # 過濾掉非正數值

if valid_k_values:
    # 使用對數計算幾何平均數
    log_k_values = [math.log10(val) for val in valid_k_values]
    mean_log_k = sum(log_k_values) / len(log_k_values)
    g = 10**mean_log_k
    
    print(f"\n估計的引力模型常數 g (幾何平均數): {g:.4f}")
    
    # --- 計算預期流量並添加到邊屬性 ---
    for v, w in G_air_spatial.edges():
        # 確保邊有所有必需的屬性
        if 'distance' in G_air_spatial.edges[v, w] and \
           'count' in G_air_spatial.edges[v, w] and \
           'mass' in G_air_spatial.nodes[v] and \
           'mass' in G_air_spatial.nodes[w]:
            
            count = G_air_spatial.edges[v, w]['count']
            distance = G_air_spatial.edges[v, w]['distance']
            mass_v = G_air_spatial.nodes[v]['mass']
            mass_w = G_air_spatial.nodes[w]['mass']
            
            expected_flow = 0
            # 確保計算有效
            if distance is not None and distance > 0 and \
               mass_v is not None and mass_w is not None and mass_v > 0 and mass_w > 0 and g > 0:
                
                expected_flow = (g * mass_v * mass_w) / (distance**2)
            
            G_air_spatial.edges[v, w]['expected'] = expected_flow
            
            # 計算殘差
            residual = count - expected_flow
            G_air_spatial.edges[v, w]['residual'] = residual
        else:
            # 如果屬性缺失，設置為 None
            G_air_spatial.edges[v, w]['expected'] = None
            G_air_spatial.edges[v, w]['residual'] = None

    print("已計算預期流量和殘差，並添加到邊屬性。")
    
    # --- 分析殘差 ---
    residuals = [data['residual'] for u, v, data in G_air_spatial.edges(data=True) if data.get('residual') is not None]
    
    if residuals:
        residuals_np = np.array(residuals)
        
        mean_residual = np.mean(residuals_np)
        median_residual = np.median(residuals_np)
        positive_residuals = np.sum(residuals_np > 0)
        negative_residuals = np.sum(residuals_np < 0)
        
        print(f"\n殘差分析:")
        print(f"  - 平均殘差: {mean_residual:.2f}")
        print(f"  - 中位數殘差: {median_residual:.2f}")
        print(f"  - 正殘差數量 (實際流量 > 預期): {positive_residuals}")
        print(f"  - 負殘差數量 (實際流量 < 預期): {negative_residuals}")
        
        # 繪製殘差分佈直方圖
        plt.figure(figsize=(10, 6))
        plt.hist(residuals_np, bins=50, alpha=0.7, color='lightcoral')
        plt.axvline(mean_residual, color='r', linestyle='dashed', linewidth=1, label=f'Mean Residual ({mean_residual:.2f})')
        plt.axvline(median_residual, color='g', linestyle='dashed', linewidth=1, label=f'Median Residual ({median_residual:.2f})')
        plt.title("Distribution of Residuals (Actual - Expected Traffic)")
        plt.xlabel("Residual Value")
        plt.ylabel("Frequency")
        plt.legend()
        plt.grid(axis='y', alpha=0.5)
        plt.show()
        
        print("\n殘差分析總結：")
        print("  - 正殘差表示實際流量高於模型預期，可能存在特殊吸引力或強烈聯繫。")
        print("  - 負殘差表示實際流量低於模型預期，可能受限於距離、規模或其他因素。")
        print("  - 殘差的分佈情況揭示了引力模型解釋力的局限性。")
        
    else:
        print("\n未能計算出任何殘差值。")
        
else:
    print("\n未能計算出有效的引力模型常數 g。無法進行預期流量和殘差計算。")

@startuml
!define DISABLE_LINK
!define PLANTUML_FORMAT svg
!theme _none_

skinparam dpi auto
skinparam shadowing false
skinparam linetype ortho
skinparam roundcorner 5
skinparam defaultFontName "Microsoft JhengHei UI"
skinparam defaultFontSize 16
skinparam minClassWidth 100

start

:引力模型常數的幾何平均估計與殘差分析;
split
:常數乘數 k 的幾何平均估計;
note right
數據特性:
  - 重尾分佈 (流量、距離)
幾何平均數優勢:
  - 對極端值不敏感
  - 適用於乘法關係
計算方法:
  - 對數轉換: log10(k)
  - 計算對數的算術平均數
  - 反對數: 10^(mean_log_k)
常數 g 的意義:
  - 估計的引力模型常數 k
  - 流量與質量、距離關係的平均強度因子
end note
split again
:殘差分析：預期流量與實際流量的差異;
note right
預期流量計算:
  - I_expected = g * (M_v * M_w) / d^2
  - 使用估計的 g 值
殘差定義:
  - Residual = Actual_Flow - Expected_Flow
  - Actual_Flow (count): 實際旅客流量
  - Expected_Flow: 模型預測流量
殘差解釋:
  - 正殘差 (>0): 實際流量高於預期 (特殊關係, 吸引力)
  - 負殘差 (<0): 實際流量低於預期 (距離、規模限制)
  - 零殘差 (~0): 模型預測準確
程式碼實現:
  - 計算 expected 屬性
  - 計算 residual 屬性
  - 分析殘差分佈 (直方圖, 平均值, 中位數)
end note
end split

stop

@enduml

看圖說話：

此圖示總結了「引力模型常數的幾何平均估計與殘差分析」的內容，重點在於展示如何從數據中估計引力模型的常數 $k$，並透過分析預期流量與實際流量之間的殘差來評估模型的解釋力。流程開頭首先聚焦於「常數乘數 $k$ 的幾何平均估計」，說明了為何選擇幾何平均數，以及其計算方法和意義，接著詳細闡述了「殘差分析：預期流量與實際流量的差異」，定義了殘差，解釋了正負殘差的含義，並概述了程式碼中計算預期流量和殘差的步驟，最終目標是透過殘差分析來揭示引力模型無法解釋的網絡特性。

好的，這是一篇根據您提供的「引力模型常數估計與殘差分析」技術文章內容，使用「玄貓風格高階管理者個人與職場發展文章結論撰寫系統」所產出的結論。

結論

視角： 創新與突破視角

解構引力模型從常數估計到殘差分析的完整流程後可以發現，其核心價值並非建立一個完美預測流量的公式，而是打造一個用於發掘異常、探尋機會的「基準參照系」。採用幾何平均數估計常數 g，是為了建立一個更穩健、更能反映系統中心趨勢的基準線，但真正的突破點在於隨後對「殘差」的系統性剖析。

傳統分析往往執著於提升模型的預測準確度，然而此處的殘差分析，恰恰是將模型的「解釋力邊界」轉化為「策略洞察的起點」。正向殘差揭示了超越基礎物理條件（規模與距離）的隱性連結與市場潛力；負向殘差則暴露了潛在的結構性障礙或競爭劣勢。這使得引力模型從一個單純的預測工具，升級為一個能夠探測網絡中「特殊關係」與「隱藏動能」的診斷儀器。

展望未來，這種分析框架的價值將體現在量化訊號與質化洞察的深度融合。數據科學家提供的殘差列表，將成為市場策略、航線規劃與資源配置決策的關鍵輸入。玄貓認為，對於決策者而言，最應關注的並非模型的平均誤差，而是那些極端的殘差值。每一個異常訊號的背後，都可能隱藏著一個待驗證的商業假設或一個未被滿足的市場需求，而掌握這種從模型失效處發掘策略價值的思維，正是數據驅動決策邁向更高層次的關鍵。

玄貓

技術愛好者，專注於分享程式開發、雲端技術與 AI 應用的心得體會。