Pythonプログラミング（ステップ６・「総和型」のアルゴリズム）

このステップの目標

1. 統計量の計算

統計解析ソフト

文系・理系を問わず、実験や調査で得られたデータを分析する際に統計処理は欠かせない。 件数の多いデータを処理したり、推定や検定などを行う際にコンピュータを活用するのは、今や常識と言える。 統計計算には、所謂、統計解析ソフトを利用するケースもあるし（例えば、フリーソフトとしては R が有名)、 Excel等の表計算ソフトもこうした用途によく使われる。

しかし、出来合いのソフトは、中身がブラックボックスであるし、融通が効きにくい。 このステップと次のステップにかけて，件数の大きなデータ処理を想定し、 Pythonを使った基本的な統計計算の方法を紹介する。

特に、この種の計算の多くは「総和のパターン」の応用として捉えることができるので、「総和」が自在に使いこなせるようになっておきたい。

平均値を計算するプログラム

例題６(ex6.py)
shakyou 

# coding: utf-8
import math

n=int(input("データの件数を入れてください:"))

sum=0
for i in range(1,n+1,1):
   x=float(input("{0}番目のデータ:".format(i)))
   sum=sum+x

print("データの総和は",sum,"です")
average=sum/n
print("平均値は",average,"です")

† input()の中の"{0}番目のデータ:".format(i)の箇所に動作については、こちらのページの説明を参照のこと。

このプログラムを実行すると，最初にデータの件数を尋ねられるので，何か数値（例えば５）を入れる。すると，それぞれのデータの値の入力を求められる。データの入力が終わると，総和と平均値が表示される。

icon-pc　練習：標準偏差の計算

例題６を発展させて，平均値だけでなく標準偏差も計算して，（両方を）表示するプログラムを書いてみなさい。

icon-hint　ヒント

$n$ 個のサンプル（標本とも言う） $x_1, x_2, \cdots, x_{n}$ の標本平均は，まずデータをすべて合計し，それを $n$ で割ることによって得られる。 これを式で表すと、 $$ \bar{x}=\frac{1}{n} \sum_{i=1}^{n} x_i $$ である。

標準偏差の計算手順

標準偏差は『データのばらつきの程度』を表す重要な統計量である。 それを得るには、まず標本分散$V$を求め、その平方根を取ればよい。 標本分散は『平均値からの差を二乗した値の合計を求め，それを件数で割った量』、すなわち、 式で書くと、 $$ V = \frac{1}{n} \sum_{i=1}^{n} \left( x_i - \bar{x}\right)^2 $$ である。ここで、$\bar{x}$はサンプルの平均。これを使って、サンプルの標準偏差 $s$ は $$ s = \sqrt{V} $$ と計算できる。

このやり方では、標本平均を求める反復計算と、平均からの差の二乗和を求める反復計算の２つが必要となる。ここで、 $$ \begin{eqnarray} \frac{1}{n} \sum_{i=1}^{n} \left( x_i - \bar{x}\right)^2 = \frac{1}{n} \sum_{i=1}^n \left[ {x_i}^2 - 2 {x_i} \bar{x} + \bar{x}^2 \right] \\ =\frac{1}{n} \sum_{i=1}^n {x_i}^2 - 2 \bar{x} \frac{1}{n} \sum_{i=1}^n {x_i} + \bar{x}^2 \frac{1}{n} \sum_{i=1}^n 1 \\ = \frac{1}{n} \sum_{i=1}^n {x_i}^2 - \bar{x}^2 \end{eqnarray} $$ すなわち、標本分散は『二乗の平均値から　平均の二乗を引いた値』に等しいから、以下の手順でも計算できることがわかる：

  1. データの件数を $n$ とし、$x_i$ を $i$ 番目のデータとする$(i=1,2, \cdots , n)$。
  2. データの総和 $S$ を計算する：\[ S= \sum_{i=1}^{n} x_i \]
  3. 二乗の総和 $U$ を計算する \[U = \sum_{i=1}^{n} \left( x_i \right)^2\]
  4. 標本平均 $E=\left. {S} \middle/ {n} \right.$ を求める
  5. 二乗の平均を $W = \left. {U} \middle/ {n} \right.$ によって求める
  6. 標本分散を $V = W - E^2$ で計算する
  7. そうすると標準偏差は $s = \sqrt{V}$

ここで、２と３の総和は以下のように、ひとつのループでまとめて計算できることに注意しよう:

1: SとUを0にセットする
2: データ入力が終わるまで以下を反復する:
3:    データを入力する
4:    S ← S + (入力データ)
5:    U ← U + (入力データの2乗)
6: 反復終了

平方根の計算には、これまでにも何度か登場した math.sqrt( ) を用いる。 その際に、プログラム冒頭のimport mathを忘れぬように。


２．級数和の計算プログラム

総和のパターンで超越関数を計算するプログラム

三角関数や指数関数など，数学ではおなじみの関数は，もちろんプログラミングでもよく登場するので，sin(), exp( ) といった名前で簡単に数式中で使うことができる（例題５などを参照）。

けれども，こうした関数も，加減乗除の組み合わせだけで求めることができる。数学（解析学）を学ぶと、指数関数などの値はテイラー展開で求められることがわかるはずだ。例えば指数関数 $e^x$ は

$$ e^x = \sum_{n=0}^{\infty} \frac{1}{n!} x^n $$

を使って評価できるので、これも総和のパターンの一種と言える。最初の数項を書いてみると、

$$ e^x = 1 + x + \frac{1}{2} x^2 + \frac{1}{6} x^3 + \frac{1}{24} x^4 + \cdots $$

となり、係数の部分が1, 1/2, 1/6, 1/24という風に急激に小さくなっているので、先にすすめば進むほど項の値は急激に小さくなっていくであろう（蛇足ながら０！は１である）。

本来ならば（数学的には）無限項の総和を行わなければならないところだけれども、この様子から、級数を途中で打ち切っても結果はそれほど違わないと予想できる。そもそも、コンピュータが扱える実数（float）には有限の精度しかないのだから、十分「先」で計算を打ち切っても、結果には全く影響が無いはずである。

この数式をPythonのプログラムで(近似的に)表すと，次のようになるだろう（ex7.py）：

例題７ (ex7.py)
shakyou 

# coding: utf-8
import math

x=float(input("xの値:"))

fact=1
p=1
sum=0
for n in range(0,10,1):
    if n>0:
        fact=fact*n
        p=p*x
    sum=sum+1/fact*p

print("EXP(",x,")=",sum)

変数factが階乗の箇所，pは $x^n$ の箇所をそれぞれ計算しながら，級数の各項をsumに足し合わせている様子がわかるだろうか。factとpの初期値は 1 にセットしてあるのに注意すること。このプログラムでは、級数和を10項目までで打ち切っている。

たとえば，$x$ を１にしたら　$e^1=2.71828\cdots$ にちゃんとなるだろうか？

こちらのページにもう少し詳しい説明を書いたので参照のこと。

icon-teacher　公式の選択は大切

例題７ではテイラー級数を使って指数関数を計算したが、数学の教科書を見ると、異なる（あるいは本来の）定義として

$$ e^x = \lim_{n\to\infty} \left(1 + \frac{x}{n}\right)^{n} $$

のように書いてあるから、こちらを使って（例題７とは異なるアルゴリズムで）計算したほうが簡単そうだ。 そこで、こちらの公式に沿って次のプログラムを作ってみた（ex7b.py）。

ex7b.py

# coding: utf-8

x=float(input("xの値:"))
n=int(input("nの値(大きな数):"))

e=1
for k in range(0,n,1):
    e=e*(1+x/n)

print("EXP(",x,")=",e)

よい近似値を得るにはnをどれくらい大きくしなければならないだろうか？

icon-pc　練習：正弦関数の計算

例題７を元にして，正弦関数を加減乗除の組み合わせによって計算するプログラムを作成しなさい。正弦関数は以下のようなテイラー級数を使って計算できることが知られている：

$$ \sin(x) = \sum_{n=0}^{\infty} \frac{(-1)^n}{(2n+1)!} x^{2n+1} $$

具体的に最初の数項のみ書き下すと：

$$ \sin(x)=x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \frac{x^9}{9!} - \cdots $$

ポイント