การวิเคราะห์ข้อมูลด้วย pandas

ทำไมต้องรู้

ข้อมูลจากการทดสอบหลายสิบเที่ยวบินมีหลายตัวแปรพร้อมกัน ทั้งยี่ห้อแบตเตอรี่ น้ำหนักบรรทุก และลำดับการบิน คำถามที่ฝ่ายวางแผนภารกิจต้องการคือ “บรรทุกเพิ่มหนึ่งกิโลกรัม เวลาบินลดลงเท่าใด” pandas ที่เรียนใน UAT 104 ช่วยจัดข้อมูล ส่วนสหสัมพันธ์และการถดถอยช่วยตอบคำถามเชิงปริมาณ แต่ตัวเลขสรุปอาจหลอกได้ถ้าไม่ดูกราฟ

ตรวจข้อมูลก่อนวิเคราะห์

ไฟล์ endurance.csv เป็นข้อมูลสังเคราะห์ 24 เที่ยวบิน ทดสอบแบตเตอรี่สองยี่ห้อ น้ำหนักบรรทุกสามระดับ ระดับละสี่เที่ยว บินตาม run_order ที่สุ่มไว้

import pandas as pd

df = pd.read_csv("endurance.csv")
print(df.dtypes.to_string())
print("missing values:", int(df.isna().sum().sum()), "  rows:", len(df))
print(pd.crosstab(df["battery"], df["payload_kg"]))

flight_id            str
run_order          int64
battery              str
payload_kg       float64
endurance_min    float64
missing values: 0   rows: 24
payload_kg  0.0  0.5  1.0
battery
A             4    4    4
B             4    4    4

ตรวจสามเรื่องก่อนคำนวณอะไร คือชนิดข้อมูลถูกต้องหรือไม่ (ตัวเลขไม่ถูกอ่านเป็นข้อความ) มีค่าหายหรือไม่ และทุกกลุ่มมีจำนวนเที่ยวบินครบตามแผนหรือไม่ ตารางไขว้ยืนยันว่าทุกคู่มีสี่เที่ยว

summary = df.groupby(["payload_kg", "battery"])["endurance_min"].agg(["count", "mean", "std"]).round(2)
print(summary)
print(df.pivot_table(index="payload_kg", columns="battery", values="endurance_min", aggfunc="mean").round(2))

                    count   mean   std
payload_kg battery
0.0        A            4  25.62  0.79
           B            4  25.46  1.03
0.5        A            4  24.08  1.13
           B            4  24.78  0.79
1.0        A            4  22.56  0.94
           B            4  22.49  0.93
battery         A      B
payload_kg
0.0         25.62  25.46
0.5         24.08  24.78
1.0         22.56  22.49

ค่าเฉลี่ยลดลงชัดเจนเมื่อน้ำหนักบรรทุกเพิ่ม ส่วนความต่างระหว่างยี่ห้อเล็กกว่าการกระจายภายในกลุ่ม ซึ่งสอดคล้องกับผลการทดสอบในโมดูล 3

สหสัมพันธ์

สัมประสิทธิ์สหสัมพันธ์ของเพียร์สัน (Pearson correlation, $r$) วัดความสัมพันธ์ เชิงเส้น ระหว่างสองตัวแปร มีค่าตั้งแต่ −1 ถึง +1 ค่าลบหมายถึงตัวหนึ่งเพิ่มอีกตัวลด

from scipy import stats

r = stats.pearsonr(df["payload_kg"], df["endurance_min"])
ci = r.confidence_interval(confidence_level=0.95)
print(f"r {r.statistic:.3f}   p {r.pvalue:.2e}   95% CI [{ci.low:.3f}, {ci.high:.3f}]")

r -0.819   p 9.98e-07   95% CI [-0.919, -0.620]

การถดถอยเชิงเส้น

การถดถอยเชิงเส้น (linear regression) หาเส้นตรง $\hat{y}=b_0+b_{1}x$ ที่ผลรวมกำลังสองของ ค่าคงเหลือ (residual) $e_i=y_i-{\hat{y}}_i$ น้อยที่สุด ความชัน $b_1$ ตอบคำถามว่าบรรทุกเพิ่มหนึ่งกิโลกรัม เวลาบินเปลี่ยนไปเท่าใด

fit = stats.linregress(df["payload_kg"], df["endurance_min"])
print(f"endurance = {fit.intercept:.2f} + ({fit.slope:.2f}) x payload   R^2 {fit.rvalue ** 2:.3f}")
print(f"slope standard error {fit.stderr:.3f} min/kg")
df["predicted"] = fit.intercept + fit.slope * df["payload_kg"]
df["residual"] = df["endurance_min"] - df["predicted"]
print(df.groupby("payload_kg")["residual"].mean().round(3).to_string())
print(f"residual sd {df['residual'].std():.2f} min   prediction at 0.75 kg {fit.intercept + fit.slope * 0.75:.2f} min")

endurance = 25.68 + (-3.02) x payload   R^2 0.671
slope standard error 0.451 min/kg
payload_kg
0.0   -0.131
0.5    0.263
1.0   -0.131
residual sd 0.88 min   prediction at 0.75 kg 23.41 min

ความชัน −3.02 นาทีต่อกิโลกรัม และ $R^2=0.67$ บอกว่าน้ำหนักบรรทุกอธิบายความแปรปรวนของเวลาบินได้ประมาณสองในสาม ค่าคงเหลือเฉลี่ยของทั้งสามระดับอยู่ใกล้ศูนย์ แต่ที่ 0.5 kg เป็นบวกและที่ 0 กับ 1.0 kg เป็นลบเล็กน้อย รูปแบบนี้อาจเป็นสัญญาณว่าความสัมพันธ์โค้งเล็กน้อย ควรตรวจด้วยข้อมูลมากขึ้นก่อนใช้แบบจำลองเส้นตรง

ข้อมูลสี่ชุดของ Anscombe

ปี 1973 นักสถิติ F. J. Anscombe สร้างข้อมูลสี่ชุดที่มีค่าเฉลี่ย สหสัมพันธ์ และเส้นถดถอยเกือบเท่ากันทุกประการ แต่เมื่อวาดกราฟกลับต่างกันสิ้นเชิง

import statistics

x123 = [10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5]
anscombe = {
    "I": (x123, [8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68]),
    "II": (x123, [9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74]),
    "III": (x123, [7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73]),
    "IV": ([8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8], [6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89]),
}
for name, (x, y) in anscombe.items():
    line = statistics.linear_regression(x, y)
    print(f"{name:<3} mean x {statistics.mean(x):.2f}  mean y {statistics.mean(y):.2f}  "
          f"r {statistics.correlation(x, y):.3f}  y = {line.intercept:.2f} + {line.slope:.3f}x")

I   mean x 9.00  mean y 7.50  r 0.816  y = 3.00 + 0.500x
II  mean x 9.00  mean y 7.50  r 0.816  y = 3.00 + 0.500x
III mean x 9.00  mean y 7.50  r 0.816  y = 3.00 + 0.500x
IV  mean x 9.00  mean y 7.50  r 0.817  y = 3.00 + 0.500x

ชุด I เป็นความสัมพันธ์เชิงเส้นปกติ ชุด II เป็นเส้นโค้งที่เส้นตรงอธิบายไม่ได้ ชุด III เส้นตรงถูกดึงด้วยจุดผิดปกติจุดเดียว และชุด IV ความชันทั้งหมดมาจากจุดเดียว ข้อสรุปคือ วาดกราฟก่อนเชื่อตัวเลขสรุปเสมอ ทั้งกราฟกระจายของข้อมูลและกราฟค่าคงเหลือ

จาก log ของ PX4 สู่ตาราง

log ของ PX4 เป็นไฟล์ ULog (.ulg) โมดูล pyulog อ่านไฟล์นี้ได้ ทั้งผ่านคำสั่ง ulog2csv ที่แปลงแต่ละหัวข้อ (topic) เป็นไฟล์ CSV และผ่าน Python โดยตรง pyulog ใช้กับ ULog ของ PX4 เท่านั้น ไม่ใช่ DataFlash log ของ ArduPilot โค้ดนี้ต้องมีไฟล์ log จริง จึงไม่ถูกรันในการตรวจอัตโนมัติ

import pandas as pd
from pyulog import ULog

ulog = ULog("flight.ulg", message_name_filter_list=["vehicle_local_position", "battery_status"])
pos = pd.DataFrame(ulog.get_dataset("vehicle_local_position").data)
bat = pd.DataFrame(ulog.get_dataset("battery_status").data)
pos["time_s"] = pos["timestamp"] / 1e6
print(pos[["time_s", "x", "y", "z"]].describe())
print(bat[["timestamp", "voltage_v", "remaining"]].tail())

ค่า timestamp ของ ULog มีหน่วยไมโครวินาทีนับจากเปิดเครื่อง และแกน z ของ vehicle_local_position ชี้ลงตามกรอบ NED ชื่อฟิลด์อาจเปลี่ยนตามรุ่นเฟิร์มแวร์ จึงต้องตรวจด้วย ulog_info ก่อนเขียนโค้ด

ปฏิบัติการประจำโมดูล

ข้อผิดพลาดที่พบบ่อย

สรุป

• ตรวจชนิดข้อมูล ค่าหาย และจำนวนต่อกลุ่มก่อนวิเคราะห์ แล้วใช้ groupby กับ pivot_table สรุปตามกลุ่ม
• $r$ วัดความสัมพันธ์เชิงเส้นเท่านั้น และไม่บอกเหตุและผลถ้าไม่ได้มาจากการทดลองที่ควบคุม
• ความชันของการถดถอยตอบคำถามเชิงปริมาณ ต้องตรวจค่าคงเหลือและไม่ประมาณนอกช่วงข้อมูล
• ข้อมูล Anscombe สอนว่าต้องวาดกราฟเสมอ และ pyulog แปลง ULog ของ PX4 เป็นตารางสำหรับวิเคราะห์ต่อได้

แบบฝึกตรวจความเข้าใจ

1. $r=-0.82$ บอกอะไรเกี่ยวกับน้ำหนักบรรทุกกับเวลาบิน
2. จากสมการ $\hat{y}=25.68-3.02x$ บรรทุก 0.4 kg ทำนายเวลาบินได้เท่าใด
3. เที่ยวบินหนึ่งบรรทุก 1.0 kg บินได้ 23.0 นาที ค่าคงเหลือเท่าใด
4. $R^2=0.67$ มีความหมายอย่างไร
5. ทำไมชุด IV ของ Anscombe จึงไม่ควรใช้เส้นถดถอยสรุป