プログラミング学習：モンテカルロ法で円周率を求める（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()

x
y

0
0.790832
0.774845
1
0.568365
0.950229
2
0.835321
0.895938
3
0.967860
0.286567
4
0.624064
0.281333

x,yを列名とし、それぞれの配列をデータフレームに変換しました。

・√x²+y²を計算する。

montecarlo["hypot"]=np.hypot(montecarlo["x"],montecarlo["y"])
montecarlo.head()

x
y
hypot

0
0.790832
0.774845
1.107159
1
0.568365
0.950229
1.107238
2
0.835321
0.895938
1.224935
3
0.967860
0.286567
1.009393
4
0.624064
0.281333
0.684547

行ごとにnp.hypotで√x²+y²を計算し、新しい列 hypotに値を格納します。

・半径1の円内に点が入るかどうかの判定。

montecarlo_small=montecarlo[montecarlo["hypot"]<=1]
montecarlo_large=montecarlo[montecarlo["hypot"]>1]
montecarlo_small.head()

x
y
hypot

4
0.624064
0.281333
0.684547
5
0.411707
0.040501
0.413694
6
0.278791
0.077373
0.289328
7
0.415050
0.306739
0.516096
9
0.799246
0.122930
0.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

プログラミング学習:numpyをつかって1から50までの自然数の和を計算する。(python GCI chapter2）

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

1から50までの自然数の和を計算するプログラムを書いて、最後の計算結果を表示させるプログラムを書いてください。ただし、Numpyを使ってください。

まずはnumpyをインストール。以後、npとして呼び出しします。

import numpy as np

np.arangeで配列を作成

sample_array = np.arange(50)
sample_array

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
       34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49])

np.arangeでは1つだけ数字を指定すると、0から指定した数字の1つ前までの連番配列を作ります。50なので、0から49までの配列ができました。

配列に1を足す

sample_array=sample_array+1
sample_array

array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
       35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50])

sample_array+1で作った配列に1を足すと、全ての要素に1が足され、1から50の配列ができます。

あとはsum関数で要素全ての和を計算すればできあがり。

sum(sample_array)

1275

forループを使うよりすっきり簡単。

変数を使わなくても、printの中で計算することもできます。これだと1行で計算できる。

print(sum(np.arange(50)+1))

1275

最初から1から50の配列を作りたい場合はこのようになります。

np.arange(1,51)

array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
       35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50])

1つめの因数が配列の始めの数字、2つ目が最後の数字＋1になります。pythonでは、0から始まって指定した数の1つ前（51と指定すると50）で終わることに注意。これ、なかなか慣れずによく失敗します。

参考：
東京大学グローバル消費インテリジェンス寄付講座
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
プログラミング学習:10までの素数を表示させる。(python GCI chapter1） https://eneprog.blogspot.com/2018/04/python10gci-chapter1.html

プログラミング学習:10までの素数を表示させる。(python GCI chapter1）

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

（1）10までの素数を表示させるプログラムを書いてください。なお、素数とは、1とそれ自身の数以外は約数をもたない正の整数のことをいいます。

k=10

for t in range(2,k+1):    
    for i in range(2,t+1):
        if t%i ==0:
            break

    if i==t:
        print(t)

2
3
5
7

特に工夫も何もなく、2から順に割り算を繰り返し、割り切れたらbreak、そのものの数まで割り切れなかったら素数とみなしています。

ついでに素数判定。
こちらも同様で、2から順に割っていって、割り切れる数があれば素数。そのものの数まで割り切れなかったら素数ではないと表示。

t=211

for i in range(2,t+1):
    if t%i ==0:
        break

if i==t:
    print(str(t)+"は素数")
else:
    print(str(t)+"は素数ではない")

211は素数

（2）上記をさらに一般化して、Nを自然数として、Nまでの素数を表示する関数を書いてください。

def sosu(N):

    for t in range(2,N+1):    
        for i in range(2,t+1):
            if t%i ==0:
                break

        if i==t:
            print(t)

課題1の素数表示をdefで関数化。関数呼び出しは、sosu()のかっこの中に数字を入れる。

sosu(10)

2
3
5
7

参考：
東京大学グローバル消費インテリジェンス寄付講座
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

追伸：bloggerでコードを表示させるのに、どうすればいいかと悩んで、jupyter notebookのhtmlをそのまま貼り付けしましたが、結構面倒です。もっといい方法を探す予定。

プログラミング学習：東京大学のデータサイエンティスト育成講座（GCI）の演習コンテンツ（無料）

プログラミング、データ解析を勉強しようとオンラインコースを探していたところ、下記の記事を見つけました。

新米データサイエンティスト、研修プログラムを考える。
https://qiita.com/s_yaginuma/items/a9c97f62d9bec129c33b

この中で、東大松尾研データサイエンス演習コンテンツがすばらしかったのでご紹介します。

東京大学松尾研　GCIデータサイエンティスト育成講座演習コンテンツ 公開ページ
http://weblab.t.u-tokyo.ac.jp/gci_contents/

これは、東京大学グローバル消費インテリジェンス寄付講座（Chair for Global Consumer Interigence:GCI)のデータ解析パートのコンテンツを無料で公開したものとのこと。GCIでは、2017年度に社会人一般向け Data Science Online Courseを開講していたのですが、残念ながら2018年度は東京大学在学生のみが対象のようです。

GCI;第2期 Data Science Online Course
http://gci.t.u-tokyo.ac.jp/dsonline/

コースの概要

• 大量のデータを自由自在に解析・分析し、隠れた関係性を発見する。そんなスキルを身につけた「データサイエンティスト」に対する需要は、工学分野だけならず、医療・経済・経営・ライフサイエンスなど非常に多くの分野で高まる一方です。
• 本コースでは、あらゆる分野で武器になるデータの解析・分析スキルのコアとなる機械学習およびビッグデータを扱う技術、分析結果を効果的に可視化する技術の基盤を網羅的に身につけ、一人前のデータサイエンティストとして活躍する入り口に立つことを目指します。

上記ページにカリキュラムがありますが、Pythonの基礎から始まり、データ加工、可視化、機械学習まで学ぶことができるようになっています。なお、この講座はPythonでの演習になっていて、コンテンツはJupyter notebookの.ipynb形式で公開されています。(Python、Jupyter notebookがインストールされていないとコンテンツを見ることができませんので、ご注意を)

参考：AnacondaからPythonをインストール（windows、Jupyter notebook）
http://eneprog.blogspot.com/2018/04/anacondapythonwindowsjupyter-notebook.html

データサイエンティスト講座　Chapter 1


コンテンツの中身をJupyter notebookで見たところ。とても詳しい解説が記載されており、さらにJupyter notebook上でそのままプログラミングを動かして演習することができます。



Chapter 1では、Pythonの基礎、Jupyter notebookの使い方のほか、ポイントとしてプログラミングのコツ（？）などもあり、とても丁寧な教材となっています。これからPythonをやってみようという人はChapter 1だけでも一通り見ておくとよさそうです。

この講座では各Chapterの最後に課題があります。これから、この課題をやっていってブログで紹介していこうと思います。目指せデータサイエンティスト！


参考：
新米データサイエンティスト、研修プログラムを考える。
https://qiita.com/s_yaginuma/items/a9c97f62d9bec129c33b
東京大学グローバル消費インテリジェンス寄付講座
http://gci.t.u-tokyo.ac.jp/
東京大学松尾研　GCIデータサイエンティスト育成講演習コンテンツ 公開ページ
http://weblab.t.u-tokyo.ac.jp/gci_contents/

AnacondaからPythonをインストール（windows、Jupyter notebook）
http://eneprog.blogspot.com/2018/04/anacondapythonwindowsjupyter-notebook.html