• 投稿 2010/09/12
  • 体験
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言語」との関係はだんだん薄くなってしまった。

 

だが、これも私にとっては貴重な良い経験であった。

 

(おそらく、読者の皆さんのなかには、この記事をごらんになって、自分だったらもっとスマートにコーディングできるのにと思われる方もおられるとおもう。なにしろ、あくまでも、当時の私の経験したことなので、ご寛容を願いたい。)

 

(2010-9-12)

 

 

(追記)

 

もっとスマートにコーディングするとしたら、こんな感じになるでしょうか。

 

(2010-9-12)