tak0kadaの何でもノート

発声練習、生存確認用。

医学関連は 医学ノート

降水量について調べてみる - 「ThinkStats」第1章

どこかからデータを拾ってきて「統計学的に有意」かどうかはともかく計算してみる練習をする。今回は手間だがnumpyなどの関数は使わずに自作する。

  • Q1: 一年を通じて雨が多い月はいつか
  • Q2: スキー部の活動で例年5月にBBQをするのだが雨が振りやすい曜日、週はあるのだろうか

について考えてみる。

準備

気象庁ページから10年分、滋賀県大津市の降水量データを拾ってくる。

データはcsvファイルで下のような感じ

ダウンロードした時刻:2014/04/18 18:25:46

,,大津,大津,大津
年月日,曜日,降水量の合計(mm),降水量の合計(mm),降水量の合計(mm)
,,,品質情報,均質番号
1995/4/17,月,4,8,1
1995/4/18,火,1,8,1
.. to be continued

データになぜか欠損している日があることと途中から降水量の測定精度が小数点1桁まで上がることが気になったが今回は無視する。品質に関してはどれも同じレベルだった。

Q1: 一年を通じて雨が多い月はいつか

import csv

def ave(L):
    return sum(L)/len(L)
    
def std(L):
    a = ave(L)
    return (sum([(i-a)*2 for i in L])/len(L)) * (1/2)

#make an array R of 24rows each represent earle(<16th) and late(>15th) months
#assuming that all months have about 30 days
R = [[str(i+1)+j] for i in range(12) for j in ('early','late')]

with open('shiga-rain.csv','r') as csvfile:
    csvreader = csv.reader(csvfile)
    for ncol, row in enumerate(csvreader):
        """remove the headers and skip the date if nodata"""
        if ncol < 5 or row[2] == '': continue

        month, date = int(row[0].split('/')[1]),\
                      int(row[0].split('/')[2])

        for i in range(12):
            if month == i+1:
                if date < 16:
                    R[2*i].append(float(row[2]))
                else:
                    R[2*i+1].append(float(row[2]))

for i in range(24):
    Ave = ave(R[i][1:])
    count = len(R[i][1:])
    print(str(i//2+1)+'月', Ave, count)

結果は

1月 1.7017543859649122 285
1月 1.764406779661017 295
2月 2.1543859649122807 285
2月 3.367063492063492 252
3月 3.743859649122807 285
3月 3.9654605263157894 304
4月 3.823321554770318 283
4月 3.535211267605634 284
5月 5.910526315789474 285
5月 5.368421052631579 304
6月 3.9701754385964914 285
6月 10.257894736842106 285
7月 8.764912280701754 285
7月 5.401315789473684 304
8月 3.8473684210526318 285
8月 5.296052631578948 304
9月 5.395017793594306 281
9月 6.128070175438596 285
10月 4.626315789473685 285
10月 4.243421052631579 304
11月 2.966666666666667 285
11月 2.542105263157895 285
12月 1.919298245614035 285
12月 1.824013157894737 304

だった。6月の後半まではひどく雨が降ることはないのか。

Q2: スキー部の活動で例年5月にBBQをするのだが雨が振りやすい曜日、週はあるのだろうか

days = {'日': 0,
        '月': 1,
        '火': 2,
        '水': 3,
        '木': 4,
        '金': 5,
        '土': 6}
#'guess' if there is a day inclined to rain
R1 = [[str(i)] for i in range(7)]

with open('shiga-rain.csv','r') as csvfile:
    csvreader = csv.reader(csvfile)
    for ncol, row in enumerate(csvreader):
        if ncol < 5 or row[2] == '': # remove the headers and skip the date if nodata
            continue

        month = int(row[0].split('/')[1])
        day = days[row[1]] # convert Japanese to numbers

        if month == 5:
            R1[day].append(float(row[2]))

for i in range(7):
    Ave = ave(R1[i][1:])
    Std = std(R1[i][1:])
    count = len(R1[i][1:])
    print(i, Ave, Std, count)

結果は

0 7.506024096385542 -1.7870577842159148e-15 83
1 5.505952380952381 1.0150610510858574e-15 84
2 4.720238095238095 3.383536836952858e-16 84
3 5.647058823529412 -4.702121045471251e-16 85
4 6.458823529411765 -2.716781048494501e-16 85
5 5.273809523809524 2.389622891097956e-15 84
6 4.315476190476191 4.440892098500626e-16 84

だった。日曜日は雨が多い、、のか?

今回の分析では標本の抽出的な意味では大津市の天気でBBQ場の天気が近いのか、くらいだろうか。
使ったデータは生データで実際のデータベースには他の(気温etc)のフィールドがあって本来は自分の使うフィールドからだけ取ってくる作業があるはずだが気象庁のサイトの容量制限に引っかかったのでパス。後は抽出する情報が平均で良かったのかという検討やどう要約された統計量を解釈するかが必要。