軌跡にはパンくずを

戻れるように置くのです

退屈な実験AはPythonにやらせよう

はじめに

UEC Advent Calender 2019 22日目の記事になります。

昨日はid:penguin22231さんでした。

penguin-uec.hatenablog.com

音楽はいいですよね。
あえて好きなアーティストを1つ挙げるならゆずです。
本日12月22日は冬至の日。
冬至と言ったらゆずを浮かべたお風呂でゆずを聴く。
最高ですね。
永遠に語れそうなので、これくらいで本題に移りましょうか。

この記事はなに?

1年次に基礎科学実験A(以下実験A)という科目があります。
12日目のid:baki01さんの言葉を借りるなら

基礎科学実験AはレポートはTeXやWordを使っていい代わりに不確かさの計算をしてきっちり値を求める必要があるため計算がしんどく、対して基礎科学実験Bはレポートが手書きのため手がしんどいという特徴があります。 これらの科目に関してはこれから嫌というほど先輩たちから噂を聞くと思うのでこのくらいに...

電気通信大学に合格した人へ - uec19bのブログ

要は計算が面倒なのです。
さて、その計算を何で行うかにも色々な選択肢があるでしょう。
筆算、電卓、Excelシート、これらの組み合わせ。
今回はこの選択肢にPythonを入れてあげようって話になります。

Pythonとは?

科学計算によく使われるプログラミング言語です。
コードの書き方は1年生の基礎プロ前半で学ぶRubyに近いため、比較的習得しやすいかと。

Pythonを使う方法は色々とありますが、環境設定などもあり、面倒です。
お試し程度なら無料のクラウドサービスをお勧めします。
実験Aの計算もそれで十分です。
例えば、Google Colaboratory
Googleアカウントさえあれば使えます。

環境

Python 3.7.4
ライブラリ
matplotlib 3.1.2
numpy 1.17.4
pandas 0.25.3
scikit-learn 0.21.2
scipy 1.4.1


では、試しにPythonを電卓がわりに使ってみましょう。

>>> 1 / 2 + 1 * 3 - 1  # 四則演算
2.5

# 変数を保持して使い回せるのがいいところ
>>> R_c = 10 * 1.0e3  # Ω
>>> L_c = 203.8 * 1.0e-3  # H
>>> r_c = R_c / (2 * L_c)
>>> r_c
24533.85672227674

>>> round(r_c, 2) # 偶数への丸めなので注意
24533.86

>>> round(r_c, -3)
25000.0

端数処理(偶数への丸め) - Wikipedia

数学関数はPythonライブラリのmathを使えばok。

>>> import math

>>> math.sin(1)  # 三角関数
0.8414709848078965

>>> math.log(2)  # 自然対数
0.6931471805599453

>>> math.sqrt(3)  # 平方根
1.7320508075688772

これだけで便利な電卓として使えそうですね!

同じ処理を何回も行うのが億劫なら一気に計算させましょう。
まずはそのための導入です。

>>> [1, 2, 3]  # リストの作成
[1, 2, 3]

>>> a = [1, 2, 3] 
>>> b = [4, 5, 6]
>>> a + b  # 要素同士の加算ではなくリストの結合
[1, 2, 3, 4, 5, 6]

>>> a[0] # 要素の取り出し
1

>>> sum([1, 2, 3]) # 総和
6

繰り返す処理はfor文などを使ってもいいのですが、あえて紹介しません。
今回は数値計算ライブラリのnumpyを使っていきます。

>>> import numpy as np
>>> a = np.array([1, 2, 3])  # リストからndarrayの作成
>>> b = np.array([4, 5, 6])
>>> a + b  # 要素ごとの四則演算が可能
array([5, 7, 9])

>>> a + 1  # 要素全てに対して四則演算も可能(ブロードキャスト)
array([2, 3, 4])

>>> a[0]  # 要素の取り出し
1

>>> a[0:2]  # 複数取り出す
array([1, 2])

>>> np.sum([1, 2, 3])  # 総和
6

>>> np.diff([1, 2, 3])  # 差分
array([1, 1])

>>> np.cumsum([1, 2, 3])  # 累積和
array([1, 3, 6])

>>> np.sqrt([1, 2, 3])  # 平方根
array([1.        , 1.41421356, 1.73205081])

>>> np.sin([1, 2, 3])  # 三角関数
array([0.84147098, 0.90929743, 0.14112001])

>>> np.log([1, 2, 3])  # 自然対数
array([0.        , 0.69314718, 1.09861229])

>>> np.mean([1, 2, 3])  # 平均
2.0

>>> np.var([1, 2, 3])  # 標本分散
0.6666666666666666

>>> np.var([1, 2, 3], ddof=1)  # 不偏分散
1.0

>>> np.std([1, 2, 3])  # 標準偏差
0.816496580927726

>>> np.std([1, 2, 3], ddof=1)  # 不偏標準偏差
1.0

平均や標準偏差は実験Aでたくさん出てくるのでどんどん使いましょう。

最小二乗法

最小二乗法は実験Aでよく出くわす計算ですね。
一次関数y = ax+bとすれば、既知はxy、未知なのはabです。
こんな時は機械学習ライブラリscikit learnを使いましょう。
試しにy = 2xを学習させてみます。

>>> from sklearn.linear_model import LinearRegression
>>> lr = LinearRegression()
>>> x = [[1], [2], [3], [4]]  # 2次元配列やリストのリストにする必要あり
>>> y = [2, 4, 6, 8]
>>> lr.fit(x, y)  # 学習

まだ学習させてない5, 6, 7, 8を代入させて、確認してみます。
それぞれの2倍が出力されれば良いわけです。

>>> x2 = [[5], [6], [7], [8]]
>>> pred = lr.predict(x2)  # 予測
>>> pred
array([10., 12., 14., 16.])

>>> lr.coef_[0]  # 傾き
2.000000000000001
>>> lr.intercept_  # 切片
0.0

ほぼy = 2xが学習できてますね。
レポートのまとめ方はy=ax+bを式変形して
xyを代入したらabが出てきたよ、と
いかにも計算しました感を出すのがポイントです。

さて最小二乗法に付き物なのは、求めた一次関数のグラフを方眼紙に描くことです。
こんなときに参考にできるグラフが予めあれば描きやすいことでしょう。
そんなわけで、グラフ描画ライブラリmatplotlibを使って関数の可視化を行います。

>>> import matplotlib.pyplot as plt
>>> plt.plot(x2, pred)  # 線の描写
>>> plt.plot(x2, pred, 'bo')  # 点の描写
>>> plt.show()

f:id:puman03:20191222190450p:plain

タイトルつけたり、ラベルをつけたりと色々と出来ますが、シンプルな表示にしました。

CSV読み込み

実験ノートだけでなく、表計算ソフトにメモすることもあるかもしれません。
その際、表計算ソフトの表の値を使ってPythonで計算させたいときの話になります。
そこ、Excelシート使えとか言わない。
仮にNumbersで表を作ったとしたら、ファイル>書き出す>CSV で可能です。
Excelはもう知りません。
f:id:puman03:20191221233341p:plain:w400

csvを読み込む方法は色々あります。
ですが、その後に計算が待っているのでデータ解析用ライブラリpandasを使っていきます。
例えば、作成したcsvファイルをcsvディレクトリ下においたとしましょう。

f:id:puman03:20191222153006p:plain

hinechu.ipynbからは以下のような感じでcsvを読み込めます。

>>> import pandas as pd
>>> df = pd.read_csv("csv/hinechu.csv")  # csv読み込み
>>> df

f:id:puman03:20191221235438p:plain

NaNはNot a Numberの略でセルに数値が入ってないときに割り当てられます。

Pandasと戯れる

csvをpandasを使って読み込んだのなら、pandasの使い方を知らないと不便でしょう。
なので少しだけお勉強です。

>>> a = pd.DataFrame([1, 2, 3], columns=["a"])  # リストからDataframe作成
>>> a

f:id:puman03:20191222152251p:plain

>>> bc = pd.DataFrame([[4, 7], [5, 8], [6, 9]], columns=["b", "c"])
>>> bc

f:id:puman03:20191222152335p:plain

>>> abc = pd.concat([a, bc], axis=1)  # 横に結合
>>> abc

f:id:puman03:20191222152351p:plain

>>> abc.columns = ["b", "a", "c"]  # 全ての列名を変更
>>> abc

f:id:puman03:20191222152410p:plain

abc = abc.rename(columns={"a":"b","b":"a"})  # 指定して列名を変更
abc

f:id:puman03:20191222152420p:plain

>>> abc[["a", "c"]]  # 列の取り出し

f:id:puman03:20191222152436p:plain

>>> abc[1:3]  # 行の取り出し

f:id:puman03:20191222152452p:plain

>>> abc.add(abc)  # 加算(今回はabc + abcでも可能)

f:id:puman03:20191222155251p:plain

加算addの他に減算sub、乗算mul、除算divで四則演算が可能です。

>>> abc.mean()  # 縦に平均
b    2.0
a    5.0
c    8.0
dtype: float64

>>> abc.mean(axis=1)  # 横に平均
0    4.0
1    5.0
2    6.0
dtype: float64

>>> abc.sum()  # 縦に総和
b     6
a    15
c    24
dtype: int64

>>> abc.std()  # 縦に不偏標準偏差
b    1.0
a    1.0
c    1.0
dtype: float64

>>> abc.std(ddof=0)  # 縦に標準偏差
b    0.816497
a    0.816497
c    0.816497
dtype: float64

>>> abc.index  # index取得
RangeIndex(start=0, stop=3, step=1)

これでもう十分戦えるかと。

実際に放射線をやってみる

Pythonが便利なのはわかるけど、実際に実験Aの計算をPythonでどうやるの?って疑問もあるかと思います。
なので、実験Aのテーマの一つ「放射線」を題材にして計算していきますよ。
選定理由は頻度、棒グラフ、ポアソン分布などやることが豊富だからです。


さて、表計算ソフトで実験結果をメモったとして、csvファイルを読み込むところから始めます。
csvを置く位置は先ほどと同じです。

import pandas as pd
>>> df = pd.read_csv("csv/housyasen.csv")
>>> df

f:id:puman03:20191222083829p:plain:w400

(今回は前半のみをやるので、後半の金属の厚さなどは一部省略して表示)

では出現確率を出してみようと思います。
\beta線源の出現回数(頻度)を求めてから、それぞれの総和で割ると求めることが出来るでしょう。
頻度はpd.value_countsです。

>>> frequency = df[["ß線源(100回)", "ß線源(300回)"]].apply(pd.value_counts)
>>> incidence = frequency / frequency.sum()
>>> incidence

f:id:puman03:20191222085019p:plain

次にポアソン分布を算出しましょう。

\displaystyle{
P(N, \overline{N}) = \frac{\overline{N}^{N}}{N!}e^{-\overline{N}} 
}

Nは計測値、\overline{N}は期待値(平均値)としました。
ということは計測値と平均値が必要なようです。

>>> N = incidence.index
>>> N
Float64Index([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0], dtype='float64')

>>> N_mean = df[["ß線源(100回)", "ß線源(300回)"]].mean()
>>> N_mean
ß線源(100回)    3.00
ß線源(300回)    3.15
dtype: float64

ポアソン分布は数値解析ライブラリのscipyを使います。
poisson.pmfってやつです。
np.stackは算出した二つのポアソン分布を結合させるために用いています。

>>> import numpy as np
>>> from scipy.stats import poisson
>>> poisson_100 =  poisson.pmf(N, N_mean["ß線源(100回)"])
>>> poisson_300 =  poisson.pmf(N, N_mean["ß線源(300回)"])
# 後々のためにテーブル状にしておきます
>>> poisson_pmf = pd.DataFrame(np.stack([poisson_100, poisson_300], axis=1), columns=["ß線源(100回)", "ß線源(300回)"])
>>> poisson_pmf

f:id:puman03:20191222091353p:plain

ここでグラフの可視化をしておきましょう。棒グラフの場合はplt.barです。

>>> import matplotlib.pyplot as plt
>>> plt.bar(N, incidence["ß線源(100回)"], color="lightgrey", label="incidence")
>>> plt.plot(N, poisson_pmf["ß線源(100回)"], 'bo', ms=8,  color="dimgrey", label='poisson pmf')
>>> plt.legend(loc='upper right')  # 凡例の位置を右上に
>>> plt.xlabel("count")  # x軸のラベル付け
>>> plt.ylabel("incidence")  # y軸のラベル付け
>>> plt.show()

f:id:puman03:20191222135206p:plain

\beta線源(300回)も同じように可視化できるので省略しました。

今度は計測値と出現回数の積を求めます。

>>> nN_mul = frequency.mul(N, axis=0)
>>> nN_mul

f:id:puman03:20191222094452p:plain

次は2 乗誤差と出現回数の積です。\beta線源(100回)と\beta線源(300回)の二つを同時に計算したかったのでnp.tileを使って二つ分作成しています。

>>> nN_mse = pd.DataFrame(np.tile(N, (2, 1)).T, columns=["ß線源(100回)", "ß線源(300回)"])
>>> nN_mse = (nN_mse.sub(N_mean, axis=1)**2).mul(frequency, axis=1)
>>> nN_mse

f:id:puman03:20191222094728p:plain

最後は全てを結合させて、合計値を付けたら完成です。
\beta線源(100回)のみを載せておきます。

>>> beta100 = pd.concat([frequency["ß線源(100回)"], incidence["ß線源(100回)"], nN_mul["ß線源(100回)"], nN_mse["ß線源(100回)"], poisson_pmf["ß線源(100回)"]], axis=1)
beta100.columns = ["出現回数", "出現確率", "計測値とnNの積", "2乗誤差とnNの積", "ポアソン分布"]
>>> beta100 = beta100.append(beta100.sum(), ignore_index=True)
>>> beta100

f:id:puman03:20191222095138p:plain

私はlatex派なので、to_latexを使いlatexに変換します。

>>> print(beta100.to_latex(na_rep="-", column_format="cccccc").replace("{}", "計測値").replace("11", "合計値"))

\begin{tabular}{cccccc}
\toprule
計測値 &   出現回数 &  出現確率 &  計測値とnNの積 &  2乗誤差とnNの積 &    ポアソン分布 \\
\midrule
0  &    2.0 &  0.02 &       0.0 &       18.0 &  0.049787 \\
1  &   14.0 &  0.14 &      14.0 &       56.0 &  0.149361 \\
2  &   32.0 &  0.32 &      64.0 &       32.0 &  0.224042 \\
3  &   18.0 &  0.18 &      54.0 &        0.0 &  0.224042 \\
4  &   15.0 &  0.15 &      60.0 &       15.0 &  0.168031 \\
5  &   12.0 &  0.12 &      60.0 &       48.0 &  0.100819 \\
6  &    5.0 &  0.05 &      30.0 &       45.0 &  0.050409 \\
7  &      - &     - &         - &          - &  0.021604 \\
8  &    1.0 &  0.01 &       8.0 &       25.0 &  0.008102 \\
9  &      - &     - &         - &          - &  0.002701 \\
10 &    1.0 &  0.01 &      10.0 &       49.0 &  0.000810 \\
合計値 &  100.0 &  1.00 &     300.0 &      288.0 &  0.999708 \\
\bottomrule
\end{tabular}

NaNは-に置き換えて、センタリングもしておきました。
計測値と合計値のところも置き換えてあります 。
pdfでは以下のような感じで表示されます。

f:id:puman03:20191222103522p:plain

表示に用いたtexファイルはこんな感じです。

\documentclass{jsarticle}
\usepackage{booktabs}

\begin{document}

\begin{tabular}{cccccc}
\toprule
計測値 &   出現回数 &  出現確率 &  計測値とnNの積 &  2乗誤差とnNの積 &    ポアソン分布 \\
\midrule
0  &    2.0 &  0.02 &       0.0 &       18.0 &  0.049787 \\
1  &   14.0 &  0.14 &      14.0 &       56.0 &  0.149361 \\
2  &   32.0 &  0.32 &      64.0 &       32.0 &  0.224042 \\
3  &   18.0 &  0.18 &      54.0 &        0.0 &  0.224042 \\
4  &   15.0 &  0.15 &      60.0 &       15.0 &  0.168031 \\
5  &   12.0 &  0.12 &      60.0 &       48.0 &  0.100819 \\
6  &    5.0 &  0.05 &      30.0 &       45.0 &  0.050409 \\
7  &      - &     - &         - &          - &  0.021604 \\
8  &    1.0 &  0.01 &       8.0 &       25.0 &  0.008102 \\
9  &      - &     - &         - &          - &  0.002701 \\
10 &    1.0 &  0.01 &      10.0 &       49.0 &  0.000810 \\
合計値 &  100.0 &  1.00 &     300.0 &      288.0 &  0.999708 \\
\bottomrule
\end{tabular}

\end{document}

良さげですね。
有効数字やタイトルなど、細かいところを調整する必要があるのですが、こんな感じの流れでやっていきます。
不備もありますが、とりあえず前半戦は終了。
後半戦は自分で試してみてくださいね。

最後に

記事を書きながら、なんで休日に実験Aをやっているのだろうと悲しい気持ちに襲われていました。
少しでもPythonで計算してみようかなっと思っていただけたら幸いです。
それでは、素敵な実験A生活を!


明日はUECの学生の音楽事情を取り上げてくれるそうですよ。

wktk037.hatenablog.com

前後で音楽に囲まれているので、この記事も実質音楽の話題ですね!

2次元配列の複数の部分を別の2次元配列の同じ位置に代入したい[Numpy]

備忘録

結論から言うと、Numpyのix_関数を使えばできる。

numpy.ix_ — NumPy v1.16 Manual

求める条件

  • 単純なスライスでは対応できない飛び飛びの位置でも可能
  • Pythonのfor文は出来るだけ避けたい

前提条件

  • 代入したい行と列のインデックスをリストで持っている
# 例
>>> row_idx = [1, 2]
>>> column_idx = [0, 1]

この例だと (1, 0), (1, 1), (2, 0), (2, 1) の位置にある値をどうにかして別の配列の同じ位置に代入したいということ。

確認用に3×3の配列を作っておく。

>>> import numpy as np

>>> before_arr = np.arange(1, 10).reshape(3, 3)
>>> before_arr
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])
>>> after_arr = np.zeros_like(before_arr)
>>> after_arr
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])

イメージとしてはこんな感じ。

f:id:puman03:20190322121930p:plain
3×3の配列の特定の部分を別の配列に代入するときのイメージ

idxの太字(row_idxとcolumn_idx)の位置から同じ位置の灰色部分にそれぞれ代入したい。

解決案

Numpyのix_関数を使えばいいらしい。

>>> ixgrid = np.ix_(row_idx, column_idx)
>>> ixgrid
(array([[1],
        [2]]), array([[0, 1]]))
>>> after_arr[ixgrid] = before_arr[ixgrid]
>>> after_arr
array([[0, 0, 0],
       [4, 5, 0],
       [7, 8, 0]])

Cは嫌だ?でも速度が欲しい?なら「Nim」を学ぶのです

これはUEC Advent Calendar 2018 1日目の記事です。

お陰様で1日目にしてカレンダーが全て埋っています。参加表明してくださった方々、ありがとうございます。
今日から毎日いろんな人の記事が読めるみたいですよ。楽しみですね。

今日は記念すべき1日目なので、Nimを知らずに人生を損している皆様のためにも、Nimの紹介をしようと思います。


Nimとは

Nimは今後流行ればいいなと個人的に思っているプログラミング言語のひとつです。

Nimのおすすめな点は

Rubyのようなわかりやすい構文で書けて
C言語に迫る実行速度を実現できる
夢のようなコンパイル言語

と言えます。 とにかくこのことだけでも覚えてもらえれば僕はもう満足です。

続きを読む