中学・高校数学でオイラーのγの近似値を求めるプログラミング：数学×Pythonプログラミング入門（2/2 ページ）
調和数列を使ってオイラーのγ（ガンマ）の近似値を求める問題を通して、Σを使った総和の計算をプログラミングする。また、Pythonのリストや繰り返し処理についてまとめ、身近な事例を使って、繰り返し処理の制御変数とリストのインデックスをどう対応させるかという問題について考える。
4. 可処分所得の前年度比を求めるサンプルプログラム
ここまで、総和（累計）や個数を求めるためのプログラミングについてざっと見てきました。リストと繰り返し処理は「相性」がいいので、多くの面倒な計算を自動化するのに役立ちます。ただ、意のままにプログラミングするには繰り返しの制御変数やリストのインデックスの操作に慣れる必要があります。
- 制御変数： for i in range(0, n)のiにあたる変数
- リストのインデックス： data[i]のiにあたる変数
なおPythonでは、インデックスを使わなくてもfor x in dataという形式で、dataの個々の値を変数xに代入しながら繰り返し処理したり、インデックスが必要な場合は、for i, x in enumerate(data)のような形式で、インデックスiと個々の値xを取得しながら繰り返し処理したりできます。しかし今回はインデックスの操作について学ぶために、for i in range(0, n)のようにしてインデックスだけを取得しながら繰り返し処理するコードを書いてみることにします。
さて、これまでのサンプルプログラムは繰り返し処理だけでできるものばかりで、インデックスを使ってリストを操作する例はありませんでした（頻度表の作成で少し触れましたが）。プログラミングの入門書などには、合計を求めたり、平均を求めたりするプログラムを例に、インデックスを使ってリストを操作する例がよくあります。以下のようなサンプルプログラムをご覧になった方も多いでしょう。
def average(x):
sum = 0
for i in range(0, len(x)):
sum += x[i]
return sum / len(x)
data = [67, 85, 78, 72]
average(data)
# 出力例：75.5
「全ての和の総和を求めて、個数で割る」という、数学での平均値の求め方に沿って素直に書くとこのようなコードになる。ただし、上でも記したように個々の値を取得する書き方（リスト3でも紹介した書き方）の方がPythonでは一般的。また、平均値を求めるだけなら、合計を求めるsum関数とリストの長さを求めるlen関数を使えば、sum(x)/len(x)という1行のコードで済む。実際のプログラミングでは、既に作られているものがあればそれを利用するのが得策だが、ここでは、インデックスの働きをおさらいするために、あえて上のような（冗長な）コードを示した。なお、ここではsumを変数名として使っているが、sum関数と同じ名前なので、同じ関数定義などの中でそれらを同時に使うとエラーになることに注意。
しかし。しかし、です。この程度のプログラムなら理解もできるし、簡単に作ることもできるのに、別の問題を自分でプログラミングしようとすると途方に暮れてしまう、という人も多いようです。一体なぜでしょう。それは、恐らく制御変数とインデックスの対応が複雑になった場合にあまり慣れていないからだと思います。リスト9のプログラムだと制御変数iがそのままリストのインデックスとして使えたので、それらの対応について深く考えなくてもよかったからです。
そこで、制御変数とリストのインデックスの対応を考える例を通して、それらの操作に慣れることにしましょう。いかにも数学的な例ばかりが続いているので、ちょっと気分を変えるために、日常生活に関わりの深いデータを使ってみます。
以下のデータは、2007年から2020年までの勤労者世帯の可処分所得の推移です（出典は総務省統計局)。
図2 勤労者世帯の可処分所得の推移
2018年〜2019年になってようやく2008年の水準に戻り、それ以降は上昇の傾向が見られる。しかし、データが「全世帯」ではなく「勤労者世帯」であることや、ここには示されていないが消費支出が減少していることなどから、一概に景気がよくなっているとは判断できないことに注意。
このデータが以下のようなリストとして表されているものとします。
income = [402116, 402932, 383960, 389848, 380863, 383851, 380966, 381929, 381193, 376576, 382434, 400964, 416980, 431992]
例えば、income[0]が2007年の可処分所得の値（つまり402116）です。
このデータを基に、対前年度比（％）を最近の値から順に求めたリストincrを作ってみましょう（関数にしなくても構いません）。つまり、元のデータと順序が逆になるわけです。先に結果を示しておきます。
incr
# 出力例：[103.6, 104.0, 104.8, 101.6, 98.8, 99.8, 100.3, 99.2, 100.8, 97.7, 101.5, 95.3, 100.2]
最近の値からのリストなので、incr[0]は2020年の値/2019年の値*100となる。つまり、incr[0] = income[13]/income[12]*100となる。実際に計算してみると、431992÷416980×100＝103.6となる。
対前年度比（％）は以下の数式で表されます。数式そのものは簡単ですね。
ある年の値 ÷ その前の年の値 × 100 （%）
では、プログラミングに取り組みましょう。ここでは、これまでとちょっとやり方を変えて、まず、incrというリストをあらかじめ作っておき、先頭から順に結果を入れていくようにしてみましょう（appendメソッドは使いません。またreverseメソッドを使うとリストを逆順にできますが、練習のためなのでそれも使わないことにします）。
というわけで、以下のコードの続きを考えてみてください。なお、小数点以下1桁で四捨五入するには、Pythonの標準組み込み関数であるround関数を使って、round(数値, 1)と書きます。
income = [402116, 402932, 383960, 389848, 380863,
383851, 380966, 381929, 381193, 376576,
382434, 400964, 416980, 431992]
incr = [0]*(len(income)-1) # incomeよりも1つ短いリストを作っておく
n = len(incr)
for i in range(n):
# for文の内容を考える
incr # 結果の表示
for文の内容を考えてみよう。incomeの長さは14だが、incrの長さは13。iの値は0〜n-1までとなる。
いかがでしょう。対前年度比を求める数式は簡単なのに、for文の内容を考えると急にややこしくなった感じがしますね。そのような場合には、以下のポイントに留意しながら進めるといいでしょう。
- 具体的な値（インデックス）を幾つか使って考え、一般化する
- 端から「はみ出す」場合に注意する
では、具体的に考えてみます。以下の図をご覧ください。
図3 制御変数とインデックスの対応を考える
iの値は0からn-1まで1ずつ増えていく。incrには順に結果を入れていくので、incr[0]〜incr[n-1]までを処理していけばいい。つまり、制御変数とインデックスは一致している。一方、incomeについては、iが0のときはincome[13]とincome[12]を使って計算する。そのため、制御変数からインデックスを求めるルールを見付ける必要がある。
制御変数iとincrのインデックスは一致させればいいので、これに関しては難しくありませんね（incomeのインデックスと一致させて処理する方法もありますが）。問題はincomeのインデックスです。iの値が0のとき、income[13]とincome[12]を使って計算するので、0を13と12に対応させてやる必要があります。また、iの値が1のとき、income[12]とincome[11]を使って計算するので、1を12と11に対応させてやる必要があります。このようにして、いくつか具体例を考えてみると、一般的なルールを導き出す手がかりが得られます。つまり、以下のような例からルールを導き出せばいいということが分かります。
（1） iの値が増えれば、incomeのインデックスが減る
（2） iの0番を使って、incomeの13番にあたる値を取り出す
（3） iの0番を使って、incomeの12番にあたる値を取り出す
（1）から、制御変数の順序を逆転させてやれば、インデックスの順序に対応させられるということが分かります。そのためには、iの値にマイナスを付けてやればいいですね。
（2）からは、13を足せばいいということが分かります。この13というのはnの値です。従って、-i+n、つまりn-iをインデックスとして使えばいいということが分かります。
（3）からは、-iにn-1を足せばいいということが分かります。従って、n-i-1をインデックスとして使います。念のため図でも見ておきましょう。
図4 制御変数の逆順に処理する方法の例
制御変数の逆の順序でリストのインデックスを利用するには、制御変数にマイナスを付けて順序を逆転させ、先頭を合わせるためにリストの個数（＝n）から引くというのが基本的。先頭がずれる場合は上の（3）のようにずらした数だけ引く（あるいは足す）とよい。
以下の動画は、図4の動きを解説したものです(プログラミングではなく、図解の動画です)。ぜひ、参考にしてみてください。
動画4 制御変数とインデックスの取り扱い
というわけで、
incr[i] ← income[n-i]とincome[n-i-1]を使って計算した結果
とを求めればいいということが分かります。iの値は増えていきますがn未満（＝12まで）なので、n-i-1の最小値は13−12−1＝0です。従って、incomeの端からはみ出すこともありません。というわけで、プログラムを完成させましょう。
income = [402116, 402932, 383960, 389848, 380863,
383851, 380966, 381929, 381193, 376576,
382434, 400964, 416980, 431992]
incr = [0]*(len(income)-1)
n = len(incr)
for i in range(n):
incr[i] = round(income[n-i] / income[n-i-1] * 100, 1) # これが答え
incr
# 出力例：[103.6, 104.0, 104.8, 101.6, 98.8, 99.8, 100.3, 99.2, 100.8, 97.7, 101.5, 95.3, 100.2]
インデックスの扱いさえ分かれば、数式に沿ってコードを書くだけ。小数点以下1桁まで求めるので、round関数を使って四捨五入している。
今回取り扱ったリストはインデックスが1つだけの（1次元の）例ばかりでした。今後、行列などの計算を行う場合には、リストが二次元以上になることもあるので、制御変数とインデックスの対応を2つ以上考えなければなりません。しかも、ある制御変数がリストの行を表したり列を表したりすることがあるので、頭が付いていかなくなることがあるかもしれません。その場合でも、上で述べたような「具体例から一般化する」「はみ出す場合に注意する」といったアプローチが助けになるのではないかと思います。ですが、そういった事例については、また回を改めて。
というわけで、まだまだ基礎固めのレベルですが、もう少し制御変数やインデックスの取り扱いに慣れるため、ぜひ以下の練習問題に取り組んでみてください。多少難しいな、と感じても一歩ずつ丁寧に考えていけば必ずできます。
練習問題
（1）ネイピア数(e)の近似値を求める
ネイピア数（自然対数の底e＝2.71828...）は、さまざまな形で表されますが、以下の定義もその一つです。nの値を指定して、eの近似値を求める関数myexpを作成してみてください。
（ヒント） n!はnの階乗を表します。例えば、4!＝4×3×2×1＝24です。階乗はmathモジュールに含まれているfactorial関数を使えば求められます。
import math
math.factorial(4)
# 出力例：24
myexp関数を作成し、例えば、10を与えて呼び出すと以下のような結果になります。
myexp(10)
# 出力例：2.7182815255731922
ちなみに、n!はnが少し大きくなるだけで、結果がとても大きくなります。例えば、10!＝3628800で、20!＝2432902008176640000です。従って、1/n!の値はすぐにとても小さな値になります。というわけで、関数myexpにはそれほど大きな値を与えなくても（今回のように10でも）良い近似値が得られます。
（2）共分散を求める
共分散を求める関数covarを作成してみましょう。共分散は以下の式で求められます。
ただし、X＝{xi},Y＝{yi}とし、
はそれぞれX,Yの平均、nはデータの個数とします（ここでは、Xの個数とYの個数は同じものとするので、Xの個数を使えばよい）。
共分散はXとYの関係を表す値で、Xが増えればYも増える場合は正の大きな値になります。逆にXが増えればYが減る場合には負の大きな値になります。例えば、aを気温のデータとし、bをビールの売上本数とし、
a = [12.1, 15.3, 18.6, 21.7, 26.1, 32.1]
b = [45, 520, 2864, 6874, 25487, 102870]
というデータが与えられたとき、関数covarの実行結果は以下のようになります（気温が上がるとビールの売り上げも増える）。
a = [12.1, 15.3, 18.6, 21.7, 26.1, 32.1]
b = [45, 520, 2864, 6874, 25487, 102870]
covar(a, b)
# 出力例：211454.23333333337
（ヒント） コードは以下のように途中まで書かれているので、続きを書いてみてください。なお、平均値はリスト9のaverage関数を使って求めています。
def covar(x, y):
n = len(x) # データの個数
avex = average(x) # xの平均
avey = average(y) # yの平均
sumxy = 0 # 合計
# この部分を書く
（3）株価のn日移動平均を求める
以下のデータは2021年3月22日〜8月11日までのファイザー社の株価（終値、単位ドル）の一覧です。これらのデータを基に、n日移動平均を求め、そのリストを返す関数movingaverageを作成してください。
図5 ファイザー社の株価（終値）の5日移動平均と15日移動平均
5日移動平均とは、直前の5日間の平均のこと。例えば、3月26日の35.11という値は、3月22日〜3月26日の5日間の平均。また3月29日の35.23は3月23日〜3月29日の5日間の平均（土日などを挟むので日付は連続していないが5日分）。このように、5日間の値を1日ずつずらしながら順に平均を求めていく。短期移動平均が長期移動平均を下から上に抜ける箇所はゴールデンクロスと呼ばれ、株価上昇のシグナルとなる。逆に、短期移動平均が長期移動平均を上から下に抜ける箇所はデッドクロス（またはデスクロス）と呼ばれ、株価下落のシグナルとなる。
（ヒント） 以下のコードを利用すれば、株価データを取得し、dataというリストに入れることができます。pandas_datareaderは、インターネット上で公開されているデータをデータフレームと呼ばれるデータ構造に読み込むためのライブラリです。これを使った株価データ取得の詳細については、今回のテーマから外れるのでここでは省略しますが、興味のある方は「「Python」と「Google Colaboratory」で株価データ分析に挑戦：「Python」×「株価データ」で学ぶデータ分析のいろは（1） - ＠IT」を参照してください。
!pip install pandas-datareader
import pandas_datareader.data as pdr
df = pdr.DataReader("PFE.US","stooq").sort_index()
data = df.loc['2021-3-22':'2021-08-11', 'Close' ].tolist()
data
# 出力例：[35.329, 34.701, 34.946, ..., 46.31]
Google Colabを使っている場合はpandas_datareaderがあらかじめインストールされているので、1行目はなくても構わない（書いても構わないが、必要なものは既にあるというメッセージが表示される）。ここでは、pandas_datareaderのDataReaderメソッドを使って、ファイザーの株価（PFE.US）をStooqと呼ばれるサイトから取得している。取得されたデータは日付の新しい順になっているので、sort_indexメソッドを使って日付の古い順に並べ替え、2021年3月22日〜8月11までの終値（Close）を取り出し、リストにしている。
これでデータが用意できました。次に、移動平均を求めるためのデータを順に取り出します。このとき、リストのインデックスに[下限:上限]を指定すると、下限以上、上限未満のインデックスを持つ要素をリストとして取り出せます（この機能をスライスと呼びます）。
例えば、x = [10, 20, 30, 40]のとき、x[0:3]とすると、[10, 20, 30]というリストが返されます。0:3は0以上3未満という意味になります（0番から3個という意味ではありません）。この機能を使って5つ値を取り出し、5日間の平均を求めるといった計算ができます。取り出す位置（下限と上限）を順に変えていけば、次々と移動平均が求められます。
ということで、以下のコードの空欄になっている箇所（_____の部分）を埋めて完成させてください。
def movingaverage(data, period): # 移動平均の期間はperiodとする
n = len(data)
result = []
for i in range(0, _____):
result.append(round(average(data[_____]), 2))
return result
movingaverage(data, 5)
# 出力例：[120.78, 120.25, 120.66, ... , 146.52]
仮引数のdataには株価のデータをリストとして与え、periodには移動平均の期間（nに当たる値）を渡すようにする
練習問題の解答例
練習問題の解答例は以下の動画でも紹介しています。ぜひ参考にしてみてください。
動画5 練習問題の解答
（1）ネイピア数(e)の近似値を求める関数myexpの作成例
これは数式をそのままプログラムとして表現しただけなので難しくないですね。
import math
def myexp(n):
sum = 0
for i in range(n):
sum += 1 / math.factorial(i)
return sum
（2）共分散を求める関数covarの作成例
これも数式をそのままプログラムとして表現しただけです。制御変数iがそのままxやyのインデックスになっているので簡単です。
def covar(x, y):
n = len(x)
avex = average(x)
avey = average(y)
sumxy = 0
for i in range(0, n):
sumxy += (x[i] - avex) * (y[i] - avey)
return sumxy/n
数式に従って素直にコードを書くと上のようになりますが、Pythonではfor文の部分を以下のように書いた方が簡潔です。
for x_i, y_i in zip(x, y):
sumxy += (x_i - avex) * (y_i - avey)
zip関数は複数の値をまとめ、繰り返しに使える値を返す関数。この例なら(x[0], y[0])に当たる値をそれぞれx_i, y_iに代入してfor文の中の文を実行し、次に(x[1], y[1])に当たる値をそれぞれx_i, y_iに代入してfor文の中の文を実行する、という具合に繰り返し処理が進められる。どちらか少ない方の要素がなくなったらzip関数が値を返さなくなるので、繰り返しが終わる。
（3）n日移動平均を求める関数movingaverageの作成例
平均を求めるためには、リスト9のavarage関数を呼び出すだけなので簡単ですが、プログラムを完成させるには、制御変数とインデックスの対応を考える必要があります。
例えば、5日移動平均であれば、繰り返しの0回目は0〜4番の平均を求めます。繰り返しの1回目は1〜5番の平均です。データの件数が100件あれば、最後のインデックスは99番なので、最後の繰り返しは95回目となり、95番〜99番の平均を求めます。
従って、繰り返しは0回〜96回未満までとなります。n=100、period=5なので、96はn-period+1で求められます。つまり、range(0, n-period+1)を指定すればいいですね。
スライスについては、i番からi+period番未満を取り出せばいいので、data[i:i+period]とします。
def movingaverage(data, period): # 移動平均の期間はperiodとする
n = len(data)
result = []
for i in range(0, n-period+1):
result.append(round(average(data[i:i+period]), 2))
return result
コラム 共分散から相関係数を求めるには
共分散は元の値の大小によって、結果の値の範囲が大きく変わります。そこで、結果をそれぞれの標準偏差の積で割ると、値を-1〜1の範囲に調整できます。その値が相関係数です。標準偏差は以下の式で求められます。
この式を基に標準偏差を求める関数stdevpを作成し、相関係数を求めてみると、以下のようになります。
import math
def stdevp(x):
n = len(x)
avex = sum(x)/n
sqsum = 0
for i in range(0, n):
sqsum += (x[i] - avex)**2
return math.sqrt(sqsum/n)
covar(a, b)/(stdevp(a)*stdevp(b))
# 出力例：0.8633475827899139
今回は、オイラーのγを求めるというトピックを通じて、数学に関しては数列や数列の和について見ました。そして、数式で表された総和などをプログラムとして表すことに取り組みました。また、繰り返しの制御変数とリストのインデックスを対応させる方法について、身近な例で考えてみました。
前回、今回と計算が続いたので、次回はちょっと気分を変えて、図形やグラフの描画をテーマとし、数式や数式で表されるさまざまな現象を視覚化する方法を見ていきたいと思います。
Copyright© Digital Advantage Corp. All Rights Reserved.