情報基礎A 「Cプログラミング」(ステップ6・「総和型」のアルゴリズム)

このステップの目標

1. Cによるデータ処理の基礎

統計解析ソフト

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

しかし、こうした出来合いのソフトは、中身がブラックボックスであるし、融通が効きにくい。 このステップと次のステップにかけて,C言語を使った基本的なデータ処理について紹介する。

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

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

例題6(ex6.c)
shakyou

#include <stdio.h>
main( )
{
     int n,i ;
     float sum,data,average ;

     printf("データの件数を入れてください:") ;
     scanf("%d",&n) ;
     sum = 0 ;
     for (i=1 ; i<=n ; i=i+1) {
        printf("%d 番目のデータ:", i) ;
        scanf("%f",&data) ;
        sum = sum + data ;
     }
     printf("データの総和は %f です\n",sum) ;
     average = sum / n ;
     printf("平均値は %f です\n",average) ;
}

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

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

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

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} $$ と計算できる。

けれども,分散は『二乗の平均値から 平均の二乗を引いた値』に等しいことを簡単に示せるので、具体的には,以下の手順を踏めばよい:

  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}$

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

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

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


ファイルに格納されたデータの利用

データの件数が多いと、数値をキーボードから入力するのは大変な手間となる。 TurtleEditを使わず、Linux付属のソフトウェアで作業している場合は、まず、エディタ(gedit等)で、

10
3.4 5.1 5.4 3.9 6.5 5.6 4.9 6.8 7.1 5.8

のようなデータファイルをあらかじめ作成しておいて(仮に、そのファイル名をsample.datとすると)、端末エミュレータ(Terminal)で、コマンド

a.out < sample.dat

標準出入力
リダイレクション

を実行すると、ファイルの中身が自動的に「入力」される。ここで < sample.dat は、ファイル(sample.dat) の内容を(あたかもキーボードからタイプしたかのように)プログラムに「食べさせなさい」、という意味の指示である。 このデータの例では、最初の10はデータの件数、その後、10個の数値が並んでいる。 (この方法は、UNIXだけでなく、Windowsのコマンドプロンプトでも使うことができる。 さらに詳しくは「標準入出力, リダイレクション」等のキーワードで検索すれば調べることができる。)


2.級数和の計算プログラム

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

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

けれども,こうした関数も,加減乗除の組み合わせだけで求めることができる。数学(解析学)を学ぶと、指数関数などの値はテイラー展開で求められることがわかるはずだ。例えば指数関数 $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という風に急激に小さくなっているので、先にすすめば進むほど項の値は急激に小さくなっていくであろう(蛇足ながら0!は1である)。

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

この数式をCのプログラムで(近似的に)表すと,次のようになるだろう(ex7.c):

例題7 (ex7.c)
shakyou

#include <stdio.h>
main()
{
  int n ;
  float sum, x, p, fact ;

  printf("xの値:") ;
  scanf("%f",&x) ;
  p=1 ;
  fact=1 ;
  sum = 0 ;
  for (n=0; n<10; n=n+1) {
    if (n>0) {
      fact = fact * n ;
      p = p * x ;
    }
    sum = sum + 1/fact * p ;
  }
  printf("EXP( %f ) = %f \n",x,sum) ;

}

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

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

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

icon-teacher 公式の選択は大切

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

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

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

ex7b.c

#include <stdio.h>
main()
{
  int k,n ;
  float e,x ;
  printf("x: ") ;
  scanf("%f",&x) ;
  printf("nの値(大きな数):") ;
  scanf("%d",&n) ;
  e=1 ;
  for (k=0; k<n; k=k+1) {
     e = e * (1.0+x/n) ;
  }
  printf("exp(%f) = %f \n",x,e) ;
}

nを大きくすればするほど、よい近似値が得られるはずだが、果たして現実はどうだろうか?

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

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

$$ \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 $$

ポイント