C言語の習い始めと「最小自乗法」
これも、いまでは表計算ソフトのExcelやCalcなどで簡単にできるとおもうのだが、およそ25年ほど前はまだそれほど便利ではなかった。ちょっとした実験データを整理するのに使った「最小自乗法」のコードは、自分でコードを書いて使っていた。
その当時は、パソコンも8ビットから16ビットへの移行の時期で、プログラミング言語はもちろんBASICがメインであった。
当時、作成したコードの一部をご紹介しよう。
測定するデータとして、x と y の組を何個か(n組としよう)を準備し、いわゆる y=ax+b の係数 a と b を求めるものであり、直線への回帰の度合いを r^2 という指標で表すものである。
これをBASIC言語で記述すると、次のような感じになる。
=====(ここから)=====
1 ' file"LSQL"
10 ' least square method (linear)
20 'SCREEN 0:CLS:CLEAR:KEY OFF
30 DIM X(30),Y(30),Z(30),W(30),V(30)
40 PRINT"formula Y = A X + B":PRINT
45 PRINT
50 INPUT"quantity of set ( N ) = ";N
60 PRINT:PRINT
70 FOR I=1 TO N
75 PRINT I;
80 INPUT" ( X , Y ) = ";X(I),Y(I)
85 IF X(I)=-1 THEN I=I-1:GOTO 75
90 NEXT I
100 FOR I=1 TO N
110 Z(I)=X(I)*Y(I):W(I)=(X(I))^2:V(I)=(Y(I))^2
120 NEXT I
130 Z1=0:X1=0:Y1=0:W1=0:V1=0
140 FOR I=1 TO N
150 Z1=Z1+Z(I):X1=X1+X(I):Y1=Y1+Y(I):W1=W1+W(I):V1=V1+V(I)
160 NEXT I
170 A=(Z1-(X1*Y1)⁄N)⁄(W1-(X1)^2/N)
180 B=Y1/N-A*(X1⁄N)
190 R1=(Z1-(X1*Y1)⁄N)^2
200 R2=W1-(X1)^2⁄N
210 R3=V1-(Y1)^2⁄N
220 R=R1⁄(R2*R3)
225 PRINT:PRINT
230 PRINT" A = ";A
240 PRINT
250 PRINT" B = ";B
260 PRINT
270 PRINT" R^2 = ";R
280 PRINT
290 END
=====(ここまで)=====
BASICはインタプリターとよばれるものであって、一行ずつ解釈しながら実行する。そのため、スピードは遅い。これに対して、コンパイラとよばれるものは、一度機械語に変換してから実行するので速い。
C言語のコンパイラが16ビットパソコンで使えるときき、なんとかしてBASICで書いたコードをCに変換できないかを考えたものであった。
C言語のコンパイラは、当時、LSI-C86試食版ver3.30 というものが、たしか、雑誌の付録で無償提供されていたのであった。
そんなとき、街中の本屋さんで偶然みつけたのが、科学技術関係のフリーソフトを納めたCD-ROMであり、その中におもしろいものがいくつかあった。
そのひとつが、「B to C」という変換ソフトで、BASIC言語のコードをC言語に強制的に変換するものであった。変換後の多少の手直しは必要だが、そのときは、とても便利であった。たとえば、上述の「最小自乗法」のコードは、次のようになる。
=====(ここから)=====
⁄**** Least Square Method by Osamu Furukawa ****⁄
#include "stdio.h"
#include "math.h"
#include "string.h"
#include "stdlib.h"
unsigned char _s0[256],_s1[256],_s2[256];
float a;
float b;
int i;
int n;
float r;
float r1;
float r2;
float r3;
float v[31];
float v1;
float w[31];
float w1;
float x[31];
float x1;
float y[31];
float y1;
float z[31];
float z1;
main()
{
extern float a;
extern float b;
extern int i;
extern int n;
extern float r;
extern float r1;
extern float r2;
extern float r3;
extern float v[31];
extern float v1;
extern float w[31];
extern float w1;
extern float x[31];
extern float x1;
extern float y[31];
extern float y1;
extern float z[31];
extern float z1;
putchar('\n');
printf("Least Square Method \n");
putchar('\n');
printf(" formula Y = A X + B\n");
putchar('\n');
printf(" quantity of set ( N ) = "" ? ");
scanf(" %d",&n);
putchar('\n');
for(i=1;i<=n;i+=1){
l_75: ;
printf(" %d ",i);
printf(" ( X , Y ) = "" ? ");
scanf(" %f, %f",&x[i],&y[i]);
if(x[i]==-1){
i=i-1;
goto l_75;}
}
for(i=1;i<=n;i+=1){
z[i]=x[i]*y[i];
w[i]=pow((x[i]),2);
v[i]=pow((y[i]),2);
}
z1=0;
x1=0;
y1=0;
w1=0;
v1=0;
for(i=1;i<=n;i+=1){
z1=z1+z[i];
x1=x1+x[i];
y1=y1+y[i];
w1=w1+w[i];
v1=v1+v[i];
}
a=(z1-(x1*y1)⁄n)⁄(w1-pow((x1),2)⁄n);
b=y1⁄n-a*(x1⁄n);
r1=pow((z1-(x1*y1)⁄n),2);
r2=w1-pow((x1),2)⁄n;
r3=v1-pow((y1),2)⁄n;
r=r1⁄(r2*r3);
putchar('\n');
printf(" A = %g\n",a);
putchar('\n');
printf(" B = %g\n",b);
putchar('\n');
printf(" R^2 = %g\n",r);
putchar('\n');
}
=====(ここまで)=====
と、こんな感じになるのであった。もちろん、むだの部分も多いとおもうが、計算処理は速くできたのであった。
このようなきっかけが、私のC言語の習い始めであった。ただ、その後、表計算ソフトが発達し、DOSベースでのLotus1-2-3や、Windows3.1 ベースでExcelなどが簡単に使えるようになってからは、この「C言語」との関係はだんだん薄くなってしまった。
だが、これも私にとっては貴重な良い経験であった。
(おそらく、読者の皆さんのなかには、この記事をごらんになって、自分だったらもっとスマートにコーディングできるのにと思われる方もおられるとおもう。なにしろ、あくまでも、当時の私の経験したことなので、ご寛容を願いたい。)
(追記)
もっとスマートにコーディングするとしたら、こんな感じになるでしょうか。
=====(ここから)=====
#include <stdio.h>
⁄* 最小自乗法 *⁄
main()
{
⁄* 変数の宣言 *⁄
double x[50],y[50];
double a,b;
double sum_x2,sum_x,sum_xy,sum_y;
double delta;
int n;
int i;
⁄* データ入力 *⁄
printf("Input number of data sets = ");
scanf("%d",&n);
for(i=0;i<n;i++){
printf("Input X[%2d]=",i);
scanf("%lf",&x[i]);
printf("Input Y[%2d]=",i);
scanf("%lf",&y[i]);
}
⁄* 計算 *⁄
sum_x2=0.0;
sum_x=0.0;
sum_xy=0.0;
sum_y=0.0;
for(i=0;i<n;i++){
sum_x2 += x[i]*x[i];
sum_x += x[i];
sum_y += y[i];
sum_xy += x[i]*y[i];
}
delta=n*sum_x2-sum_x*sum_x;
a=(n*sum_xy-sum_x*sum_y)/delta;
b=(-sum_x*sum_xy+sum_x2*sum_y)/delta;
⁄* 出力 *⁄
printf("a = %10.3lf\n",a);
printf("b = %10.3lf\n",b);
⁄* 関数値 *⁄
return(1);
}
⁄* end of main *⁄
=====(ここまで)=====
(2010-9-12)
■ 以前の記事
→ 「スクリプト」に関する以前の記事は、こちらをごらん下さい。
・ZIPファイル解凍時の文字化けと対処法(Linux)
・アクセスカウンタの設置方法
