2018年5月21日月曜日

サンフランシスコのレンタサイクルのデータを見てみる その2。よく使われているステーションを調べる(python、pandas、groupby、pivot_table)

edX UCバークレー データサイエンス基礎講座(Foundation of Data Science)の演習、前回の続きです。

サンフランシスコ ベイエリアのレンタサイクル "Ford Go Bike"の利用状況を把握する。
  • 利用時間の傾向、分布
  • よく使われているステーション
  • ステーションの地図の可視化

今回はこのうち、よく使われているステーションを調べてみます。

元データ
サンフランシスコ ベイエリア レンタサイクル Ford Go Bike
https://www.fordgobike.com/
システムデータ
https://s3.amazonaws.com/fordgobike-data/index.html
2017年のトリップデータ(このデータを使います)
https://s3.amazonaws.com/fordgobike-data/2017-fordgobike-tripdata.csv

import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt

trips = pd.read_csv('tripdata.csv')

duration= trips[["duration_sec","start_station_name","end_station_name"]]
duration=duration.rename(columns={'duration_sec': 'Duration','start_station_name': 'Start Station','end_station_name':"End Station"})

データ読み込み、列名変更。ここまで前回と同じ。

commute = duration[duration["Duration"] <= 2000]
commute.head()

前回見たとおり、利用時間が2000秒以上は少ないので今回は2000秒以下のみを対象にします。

common_start=commute.groupby(["Start Station"]).count().sort_values(by=["Duration"],ascending=False).reset_index().drop("End Station",axis=1)
common_start=common_start.rename(columns={'Duration':'start_count'})
common_start.head()
DurationEnd Station
Start Station
San Francisco Ferry Building (Harry Bridges Plaza)1333813338
San Francisco Caltrain (Townsend St at 4th St)1232012320
San Francisco Caltrain Station 2 (Townsend St at 4th St)1189011890
The Embarcadero at Sansome St1186111861
Market St at 10th St1157711577

借りたステーション(Start Station)でグループ化(groupby)し、回数を集計(count)して並べ替え(sort)。列名が階層化されています。

common_start=common_start.reset_index().drop("End Station",axis=1).rename(columns={'Duration':'start_count'})
common_start.head()
Start Stationstart_count
0San Francisco Ferry Building (Harry Bridges Pl...13338
1San Francisco Caltrain (Townsend St at 4th St)12320
2San Francisco Caltrain Station 2 (Townsend St...11890
3The Embarcadero at Sansome St11861
4Market St at 10th St11577


reset_index()でインデックスを振り直して列名の階層を解除、dropで"End Station"列を削除、renameで列名を変更しました。

同様に返却ステーション(End Station)でも利用回数を計算しました。
End Stationend_count
0San Francisco Caltrain (Townsend St at 4th St)17225
1San Francisco Ferry Building (Harry Bridges Pl...15770
2San Francisco Caltrain Station 2 (Townsend St...13532
3Montgomery St BART Station (Market St at 2nd St)13115
4The Embarcadero at Sansome St13038


借りたステーション、返却したステーション、それぞれ利用回数が多い順に並んでいます。全体に返却ステーションの方が回数が多いようです。あちこちから目的地に向かって集まっているというイメージでしょうか。

common_list=pd.concat([common_start,common_end],axis=1)
common_list.head() 
Start Stationstart_countEnd Stationend_count
0San Francisco Ferry Building (Harry Bridges Pl...13338San Francisco Caltrain (Townsend St at 4th St)17225
1San Francisco Caltrain (Townsend St at 4th St)12320San Francisco Ferry Building (Harry Bridges Pl...15770
2San Francisco Caltrain Station 2 (Townsend St...11890San Francisco Caltrain Station 2 (Townsend St...13532
3The Embarcadero at Sansome St11861Montgomery St BART Station (Market St at 2nd St)13115
4Market St at 10th St11577The Embarcadero at Sansome St13038


借りたステーションと返却したステーションの情報を結合します。pd.concat、axis=1で列の横方向結合です。こちらも借りたステーション、返却したステーション、それぞれ利用回数が多い順に並んでいます。

common_station=common_start.set_index('Start Station').join(common_end.set_index('End Station')).reset_index()
common_station=common_station.rename(columns={'Start Station':'Station'})
common_station.head()
Stationstart_countend_count
0San Francisco Ferry Building (Harry Bridges Pl...1333815770
1San Francisco Caltrain (Townsend St at 4th St)1232017225
2San Francisco Caltrain Station 2 (Townsend St...1189013532
3The Embarcadero at Sansome St1186113038
4Market St at 10th St1157711022


借りたステーションと返却ステーションの利用回数情報をステーション名で結合(join)してみました。joinはindexでの結合になるので、set_indexでStart Station、End Stationをそれぞれindexに設定してから結合し、借りたステーションの利用回数順に並べ替えています。

これを見ると、借りたステーションと返したステーションで利用回数にかなり違いがあります。返却が多いステーションから借りる方が多いステーションへ自転車を移動する必要がありそうです。

commute_pivot=pd.pivot_table(commute,index='Start Station', columns='End Station',values='Duration',aggfunc="count")
commute_pivot.loc[common_start["Start Station"],common_end["End Station"]].iloc[0:10,0:10]
End StationSan Francisco Caltrain (Townsend St at 4th St)San Francisco Ferry Building (Harry Bridges Plaza)San Francisco Caltrain Station 2 (Townsend St at 4th St)Montgomery St BART Station (Market St at 2nd St)The Embarcadero at Sansome StMarket St at 10th StPowell St BART Station (Market St at 4th St)Berry St at 4th StSteuart St at Market StPowell St BART Station (Market St at 5th St)
Start Station
San Francisco Ferry Building (Harry Bridges Plaza)553.0304.0104.0329.02842.0145.0249.01381.028.0216.0
San Francisco Caltrain (Townsend St at 4th St)52.0734.011.0459.0361.0149.0193.010.0268.0221.0
San Francisco Caltrain Station 2 (Townsend St at 4th St)11.0580.048.0548.0191.0618.0328.010.0226.0241.0
The Embarcadero at Sansome St453.01639.068.0296.0486.063.0358.0323.01764.0279.0
Market St at 10th St200.0273.01146.0726.069.0144.0825.095.0179.0661.0
Montgomery St BART Station (Market St at 2nd St)985.0588.0295.091.0346.0434.0123.0363.0213.0117.0
Berry St at 4th St16.01671.010.0430.0459.0122.0161.0178.0642.0233.0
Howard St at Beale St1181.0174.0232.066.0581.0267.0139.0381.084.093.0
Powell St BART Station (Market St at 4th St)664.0356.0315.0170.0339.0607.0146.0187.0230.0134.0
Steuart St at Market St707.027.0149.0121.01290.054.0145.0839.086.083.0


pivot_tableで借りたステーションと返却したステーションの対応を見てみます。行(index)と列(columns)をそれぞれ選択、aggfuncで集計方法をします。.locで行、列の並び順を利用回数の多い順に並んだシリーズで指定し、ilocでトップ10のみを表示させました。

San Francisco Ferry Buildingで借りた人がThe Embarcadero at Sansome Stに返した回数は2842回、Berry St at 4th Stに返した回数は1381回です。けっこうバラつきがあり、自転車を過不足なく各ステーションに配置するのは大変そうです。

つづく

参考:
edx UC Berkeley Foundations of Data Science(UCバークレー データサイエンス基礎講座)
ttps://www.edx.org/professional-certificate/berkeleyx-foundations-of-data-science
Computational Thinking with Python(pythonによるプログラミング的思考)
https://www.edx.org/course/foundations-data-science-computational-uc-berkeleyx-data8-1x

サンフランシスコ ベイエリア レンタサイクル Ford Go Bike
https://www.fordgobike.com/
システムデータ
https://s3.amazonaws.com/fordgobike-data/index.html

サンフランシスコのレンタサイクルのデータを見てみる(基本統計量、python、pandas、ヒストグラム)
https://eneprog.blogspot.com/2018/05/pythonpandas_17.html

プログラミング学習:edx UCバークレー データサイエンス基礎講座の紹介(python)
https://eneprog.blogspot.com/2018/04/edx-uc-python.html
プログラミング学習:pandasでweb上の表を取得する。(python)
https://eneprog.blogspot.com/2018/04/pandaswebpython.html プログラミング学習:各国の女性の平均就学期間をグラフ化する。(python,pandas)
https://eneprog.blogspot.com/2018/05/pythonpandas.html

2018年5月17日木曜日

サンフランシスコのレンタサイクルのデータを見てみる(基本統計量、python、pandas、ヒストグラム)

edX UCバークレー データサイエンス基礎講座(Foundation of Data Science)の演習です。

サンフランシスコ ベイエリアのレンタサイクル "Ford Go Bike"の利用状況を把握する。
  • 利用時間の傾向、分布
  • よく使われているステーション
  • ステーションの地図の可視化

今回はこのうち、利用時間の傾向・分布を見て見ます。

サンフランシスコ ベイエリア レンタサイクル Ford Go Bike
システムデータ

使用するデータは上記システムデータの2017年tripdata(https://s3.amazonaws.com/fordgobike-data/2017-fordgobike-tripdata.csv)にします。

・インポート
import numpy as np
import pandas as pd

import matplotlib
import matplotlib.pyplot as plt

データ読み込み(データが118MBと大きいため、一旦ダウンロードしてから読み込んでいます)
trips = pd.read_csv('tripdata.csv')
trips.head()
duration_secstart_timeend_timestart_station_idstart_station_namestart_station_latitudestart_station_longitudeend_station_idend_station_nameend_station_latitudeend_station_longitudebike_iduser_typemember_birth_yearmember_gender
0801102017-12-31 16:57:39.65402018-01-01 15:12:50.245074Laguna St at Hayes St37.776435-122.42624443San Francisco Public Library (Grove St at Hyde...37.778768-122.41592996Customer1987.0Male
1788002017-12-31 15:56:34.84202018-01-01 13:49:55.6170284Yerba Buena Center for the Arts (Howard St at ...37.784872-122.40087696Dolores St at 15th St37.766210-122.42661488Customer1965.0Female
2457682017-12-31 22:45:48.41102018-01-01 11:28:36.8830245Downtown Berkeley BART37.870348-122.267764245Downtown Berkeley BART37.870348-122.2677641094CustomerNaNNaN
3621722017-12-31 17:31:10.63602018-01-01 10:47:23.5310608th St at Ringold St37.774520-122.4094495Powell St BART Station (Market St at 5th St)37.783899-122.4084452831CustomerNaNNaN
4436032017-12-31 14:23:14.00102018-01-01 02:29:57.5710239Bancroft Way at Telegraph Ave37.868813-122.258764247Fulton St at Bancroft Way37.867789-122.2658963167Subscriber1997.0Female

duration_secで利用時間(秒)、借りたステーション、返したステーションなどの情報が入っています。今回はこのうち、duration_sec、start_station_name、end_station_nameだけを使います。

・必要な列の選択、列名の更新
duration= trips[["duration_sec","start_station_name","end_station_name"]]
duration=duration.rename(columns={'duration_sec': 'Duration','start_station_name': 'Start Station','end_station_name':"End Station"})
duration.head() 
DurationStart StationEnd Station
080110Laguna St at Hayes StSan Francisco Public Library (Grove St at Hyde...
178800Yerba Buena Center for the Arts (Howard St at ...Dolores St at 15th St
245768Downtown Berkeley BARTDowntown Berkeley BART
3621728th St at Ringold StPowell St BART Station (Market St at 5th St)
443603Bancroft Way at Telegraph AveFulton St at Bancroft Way

必要な3つの列のみを取り出し、renameで列名(columns)を変更しました。

・並べ替え
duration.sort_values(by=["Duration"],ascending=False).head()
DurationStart StationEnd Station
13886286369San Pedro SquareSan Pedro Square
22350886355Folsom St at 19th StCentral Ave at Fell St
5488886325Mosswood ParkMosswood Park
5488986281Mosswood ParkMosswood Park
12042386252Downtown Berkeley BARTTelegraph Ave at Alcatraz Ave

利用時間が大きい順に並べ替えしてみました。sort_valuesで値で並べ替え、byで並べ替える項目を、ascending=Falseで降順(値の大きい順)を指定しています。ascendingをTrueまたは省略すると、昇順(値の小さい順)での並べ替えになります。

利用時間の最大値は86369秒、約1日となっています。

・利用時間のヒストグラムを描く
duration.hist("Duration")


単純に利用時間の全データでヒストグラムを描いてみると、ほぼ10000秒(約170分)未満だけのグラフになってしまいました。

・基本統計量を見る。
duration.describe()
Duration
count519700.000000
mean1099.009521
std3444.146451
min61.000000
25%382.000000
50%596.000000
75%938.000000
max86369.000000

describeで基本統計量をまとめて表示できます。
count:データ数 519700件
mean:平均値 約1099秒
std:標準偏差 3444
min:最低値 61秒
25%:25%タイル値 382秒
50%:50%タイル値=中央値 596秒
75%:75%タイル値 938秒
max:最大値 86369秒
最大利用時間が86369秒(約1日)なのに対し、中央値は596秒(約10分)。75%タイル値も938秒(約15分)でほとんどが短時間の利用になっているようです。

・パーセンタイル
duration["Duration"].quantile(.90)
1502.0

quantileでタイル値を表示します。.90で90%タイル値です。1502秒未満の利用が90%を占めています。

・範囲指定でヒストグラムを描く
duration.hist("Duration",range=(0,2000))


約1500秒未満の利用が90%を占めているので、2000秒までの分布を見てみます。rangeで範囲を0以上、2000未満に設定してヒストグラムを描きました。利用時間は500秒あたりにピークがあって、その後は減っていってるようです。

・分割数(bin)を指定してヒストグラムを描く
duration.hist("Duration",range=(0,2000),bins=60)


binsで分割数を60に指定してヒストグラムを描いてみました。約400秒の利用が最も多く、その後は徐々に減っていく、右のすそが長い分布となりました。

レンタサイクル "Ford to go"の料金プランを見てみると、1回利用が30分以内で$2(追加の15分は$3)、1日(24時間)利用が$10になっています。

Ford to go 料金プラン
https://www.fordgobike.com/pricing

利用時間の分布から見ると、ほとんどが1回利用で30分(1800秒)以内に返却していて、1日利用は少ないようです。10分以下の利用が多くなっていて、近場での移動に使っているようです。1回利用で30分を過ぎてしまった利用者もそれなりにいるように見えます。

つづく

参考:
edx UC Berkeley Foundations of Data Science(UCバークレー データサイエンス基礎講座)
ttps://www.edx.org/professional-certificate/berkeleyx-foundations-of-data-science
Computational Thinking with Python(pythonによるプログラミング的思考)
https://www.edx.org/course/foundations-data-science-computational-uc-berkeleyx-data8-1x

サンフランシスコ ベイエリア レンタサイクル Ford Go Bike
https://www.fordgobike.com/
システムデータ
https://s3.amazonaws.com/fordgobike-data/index.html

プログラミング学習:edx UCバークレー データサイエンス基礎講座の紹介(python)
https://eneprog.blogspot.com/2018/04/edx-uc-python.html
プログラミング学習:pandasでweb上の表を取得する。(python)
https://eneprog.blogspot.com/2018/04/pandaswebpython.html プログラミング学習:各国の女性の平均就学期間をグラフ化する。(python,pandas) https://eneprog.blogspot.com/2018/05/pythonpandas.html


追伸:edX UC Berkley Foundations of Data Scienceの最初の講座(Computational Thinking with Python)は終了し、2番目のInferential Thinking by Resamplingが2018年5月22日から開始になります。

2018年5月10日木曜日

プログラミング学習:モンテカルロ法で円周率を求める(GCI chapter2 python、pandas)

東大 GCIデータサイエンティスト育成講座演習コンテンツ Chapter 2の課題です。

乱数を発生させる方法を使って、円周率を求めるプログラムを作成する。
  • 一様分布に従う乱数を2組発生させて、それぞれ10,000個の一様乱数を作る。
  • x−y軸を使った中心(0,0)、半径1の円と、長さ1の正方形を考える。ここで先ほどのxとyの組み合わせの乱数10,000個のうち、円の内部に入る点は何組あるか。
  • 半径1の1/4の円の面積と長さ1の長方形の面積の比は、π/4:1となる。これを利用して、円周率を求めよ。

モンテカルロ法については下記ページが参考になりました。

 モンテカルロ法と円周率の近似計算

1x1の正方形にランダムに点をN個打ち、そのうち半径1の円(の4分の1)内に入った点の数をXとすると、4X/Nがπの近似値になる。図を見るとイメージがつきやすいかと思います。

これをpandasを使ってやってみます。

・モジュールのインポート(numpy,pandas)
import numpy as np
import numpy.random as random
import pandas as pd

・乱数の発生
x=np.random.uniform(0, 1, 10000)
y=np.random.uniform(0, 1, 10000)
print(x,y)
[0.15469463 0.35501088 0.31071841 ... 0.67647336 0.10477222 0.77432113] [0.56009742 0.42374295 0.5874093  ... 0.62031945 0.66022409 0.92282591]

np.random.uniform(0, 1, 10000)で、0以上1未満の一様乱数を10,000個発生させ、それぞれの配列をx、yに格納します。

・x、yの配列データフレームに変換
montecarlo = pd.DataFrame({'x':x,'y':y})
montecarlo.head()
xy
00.7908320.774845
10.5683650.950229
20.8353210.895938
30.9678600.286567
40.6240640.281333
x,yを列名とし、それぞれの配列をデータフレームに変換しました。

・√x2+y2を計算する。
montecarlo["hypot"]=np.hypot(montecarlo["x"],montecarlo["y"])
montecarlo.head()
xyhypot
00.7908320.7748451.107159
10.5683650.9502291.107238
20.8353210.8959381.224935
30.9678600.2865671.009393
40.6240640.2813330.684547
行ごとにnp.hypotで√x2+y2を計算し、新しい列 hypotに値を格納します。

・半径1の円内に点が入るかどうかの判定。
montecarlo_small=montecarlo[montecarlo["hypot"]<=1]
montecarlo_large=montecarlo[montecarlo["hypot"]>1]
montecarlo_small.head()
xyhypot
40.6240640.2813330.684547
50.4117070.0405010.413694
60.2787910.0773730.289328
70.4150500.3067390.516096
90.7992460.1229300.808644
hypotの値が1以下なら円内、1より大きければ円外とし、それぞれmontecarlo_small、montecarlo_largeのdataframeに格納しました。

・グラフ描画
large=montecarlo_large.plot.scatter(x="x",y="y",figsize=(8,8),color="red")
montecarlo_small.plot.scatter(x="x",y="y",ax=large,color="blue")

円外の点を赤、円内に入る点を青でそれぞれ散布図で描画します。ポイントはfigsizeで図の大きさを正方形にすること。

1番目の図を変数largeに格納し、2番目の図でax=largeと指定することで1番目の図に上書きして2番目の図が描画されます。

・出力


円内の点が青、円外の点が赤で描かれました。所々空白がありますが、ちゃんと円が浮き上がって見えます。

・近似円周率の計算
(len(montecarlo_small)*4)/10000
3.1724

半径1の内に入っている点の数を数え、4倍して全体の点の数(10000)で割ると円周率の近似値となります。lenでmontecarlo_smallの行数を数えて、これを円内に入っている点の数としました。

円周率近似値の結果は、3.1724。もっと点の数を多くすると精度は上がりそうです。

グラフ描画などもっとスマートな方法がありそうな気がしますが、ひとまずこれで。

参考:
東京大学グローバル消費インテリジェンス寄付講座
http://gci.t.u-tokyo.ac.jp/
東京大学松尾研 GCIデータサイエンティスト育成講演習コンテンツ 公開ページ
http://weblab.t.u-tokyo.ac.jp/gci_contents/
プログラミング学習:東京大学のデータサイエンティスト育成講座(GCI)の演習コンテンツ
https://eneprog.blogspot.com/2018/04/python10gci-chapter1.html

モンテカルロ法と円周率の近似計算
https://mathtrain.jp/montecarlo

プログラミング学習:10までの素数を表示させる。(python GCI chapter1)
https://eneprog.blogspot.com/2018/04/python10gci-chapter1.html
プログラミング学習:numpyをつかって1から50までの自然数の和を計算する。(python GCI chapter2)
https://eneprog.blogspot.com/2018/04/numpy150python-gci-chapter2.html
//SyntaxHighlighter CDNより https://cdnjs.com/libraries/SyntaxHighlighter // 対応言語