統計プログラム集

　　サンプルデータ
　　マウスで右クリックし，「対象を保存」を選んでください。
　　　ex1.txt　　 data1.txt　　 data2.txt　　 data3.txt　　 data4.txt
　　（乱数を用いて適当に生成したもので，特に由来はありません。）
　　　Linux，Macでは行末コードを修正してからご利用ください。

　　プログラムをデータファイルと同じフォルダに保存して実行してください。
　　(プログラムの保存名は，“NoName”を含まない名前に変更してください。)

　　　記述統計

ファイルから数値データを読み込んで件数，平均，標準偏差を計算する。

100 OPEN #1: NAME "ex1.txt"
110 LET n=0   ! 件数
120 LET s=0   ! 和
130 LET s2=0  ! 平方の和 
140 DO
150    INPUT #1,IF MISSING THEN EXIT DO : x
160    LET n=n+1
170    LET s=s+x
180    LET s2=s2+x^2
190 LOOP
200 PRINT "件数",n
210 PRINT "平均",s/n
220 PRINT "標準偏差",SQR(s2/n-(s/n)^2) 
230 CLOSE #1
240 END

度数分布を調べる～10点満点の試験結果の集計

120 DIM f(0 TO 10)
130 MAT f=ZER
140 OPEN #1: NAME "data2.txt"
150 DO
160    INPUT #1,IF MISSING THEN EXIT DO : x
170    LET f(x)=f(x)+1
180 LOOP
190 CLOSE #1
200 FOR x=0 TO 10
210    PRINT x;"点",f(x);"人"
220 NEXT x
230 END

度数分布～100点満点の試験を10点刻みで集計

110 DIM f(0 TO 10)
120 MAT f=ZER
130 OPEN #1: NAME "data3.txt"
140 DO
150    INPUT #1,IF MISSING THEN EXIT DO : x
160    LET i=INT(x/10)
170    LET f(i)=f(i)+1
180 LOOP
190 CLOSE #1
200 FOR i=0 TO 10
210    PRINT i*10;"点台",f(i);"人"
220 NEXT i
230 END

2変量データの集計～平均，標準偏差，相関係数

100 OPEN #1:NAME "data4.txt"
110 LET n=0   ! 件数
120 LET mx=0  ! xの和
130 LET my=0  ! yの和 
140 LET xx=0  ! x^2の和
150 LET yy=0  ! y^2の和
160 LET xy=0  ! xyの和
170 DO
180    INPUT #1,IF MISSING THEN EXIT DO:x,y
190    LET n=n+1
200    LET mx=mx+x
210    LET my=my+y
220    LET xx=xx+x^2
230    LET yy=yy+y^2
240    LET xy=xy+x*y
250 LOOP
260 LET mx=mx/n            ! xの平均
270 LET my=my/n            ! yの平均
280 LET sx=SQR(xx/n-mx^2)  ! xの標準偏差
290 LET sy=SQR(yy/n-my^2)  ! yの標準偏差
300 LET sxy=xy/n-mx*my
310 PRINT sxy/(sx*sy)      ! 相関係数
320 CLOSE #1
330 END

散布図を描く

110 OPEN #1:NAME "data4.txt"
120 SET WINDOW 100,200,0,100
130 DRAW grid WITH SCALE(10,10)
140 DO
150    INPUT #1,IF MISSING THEN EXIT DO: x,y
160    PLOT POINTS:x,y
170 LOOP
180 CLOSE #1
190 END

2つの回帰直線（xからy，yからx）を描く。

100 OPEN #1:NAME "data4.txt"
110 ! 準備
120 LET n=0   ! 件数
130 LET mx=0  ! xの和
140 LET my=0  ! yの和 
150 LET xx=0  ! x^2の和
160 LET yy=0  ! y^2の和
170 LET xy=0  ! xyの和
180 ! 座標設定
190 LET left=100
200 LET right=200
210 LET bottom=0
220 LET top=100
230 SET WINDOW left,right,bottom,top
240 DRAW grid WITH SCALE(10,10) 
250 ! データ読込み 
260 DO
270    INPUT #1,IF MISSING THEN EXIT DO: x,y
280    LET n=n+1
290    LET mx=mx+x
300    LET my=my+y
310    LET xx=xx+x^2
320    LET yy=yy+y^2
330    LET xy=xy+x*y
340    PLOT POINTS:x,y
350 LOOP
360 ! 計算
370 LET mx=mx/n            ! xの平均
380 LET my=my/n            ! yの平均
390 LET sx=SQR(xx/n-mx^2)  ! xの標準偏差
400 LET sy=SQR(yy/n-my^2)  ! yの標準偏差
410 LET sxy=xy/n-mx*my
420 LET r=sxy/(sx*sy)      ! 相関係数
430 ! 回帰直線を描く
440 DEF f(x)=a*(x-mx)+my
450 SET LINE STYLE 2
460 LET a=r*sy/sx
470 PLOT LINES:left,f(left);right,f(right)
480 SET LINE STYLE 3
490 LET a=1/r*sy/sx
500 PLOT LINES:left,f(left);right,f(right)
510 
520 CLOSE #1
530 END

確率論

さいころをn回投げるとき1の目の出る回数の確率分布の図示

110 OPTION ARITHMETIC NATIVE 
120 SET WINDOW -1,49,-0.01,0.49
130 DRAW axes0
140 INPUT n
150 FOR k=0 TO n
160    LET p=Comb(n,k)/6^n*5^(n-k)
170    PLOT LINES: k-0.5,0 ; k-0.5,p ; k+0.5,p ; k+0.5,0
180 NEXT k
190 END

標本平均の分布のシミュレーション

100 DIM a(1000)  ! 母集団
110 LET n=10     ! 標本の大きさ
120 DIM s(n)     ! 標本
130 OPEN #1: NAME "data3.txt"
140 FOR i=1 TO 1000
150    INPUT #1:a(i)
160 NEXT i
170 CLOSE #1
180 DIM d(0 TO 100)
190 MAT d=ZER 
200 RANDOMIZE
210 LET times=2000  ! 実験回数
220 FOR t=1 TO times
230    FOR i=1 TO n
240       LET s(i)=a(INT(RND*1000)+1)
250    NEXT i
260    LET m=0
270    FOR i=1 TO n
280       LET m=m+s(i)
290    NEXT i
300    LET m=m/n
310    LET i=INT(m)
320    LET d(i)=d(i)+1
330 NEXT t
340 SET WINDOW 0,100,0,0.4
350 FOR i=0 TO 100
360    PLOT AREA:i,0;i+1,0;i+1,d(i)/times;i,d(i)/times
370 NEXT i
380 END

さいころの目の数の平均の分布とそれを近似する正規分布

100 DIM f(10,60)
110 MAT f=ZER
120 FOR x=1 TO 6
130    LET f(1,x)=1
140 NEXT x
150 FOR k=2 TO 10
160    FOR x=k TO 6*k
170       FOR y=x-6 TO x-1
180          IF k-1<=y AND y<=6*(k-1) THEN LET f(k,x)=f(k,x)+f(k-1,y)
190       NEXT y
200    NEXT x
210 NEXT k
220 SET WINDOW -1, 7, -0.03, 1
230 FOR k=1 TO 10
240    CLEAR
250    SET LINE COLOR 1
260    DRAW axes
270    LET w=1/k
280    FOR x=k TO 6*k
290       LET h=f(k,x)/6^k/w
300       PLOT LINES: x/k-w/2,0; x/k-w/2,h; x/k+w/2,h; x/k+w/2,0
310    NEXT x
320    WAIT DELAY 1
330    LET m=7/2
340    LET s2=35/12/k
350    LET s=SQR(s2)
360    SET LINE COLOR 4
370    FOR x=0 TO 7 STEP 0.01
380       PLOT LINES:x,1/(SQR(2*PI)*s)*EXP(-(x-m)^2/(2*s2));
390    NEXT x
400    PLOT LINES
410    WAIT DELAY 1
420 NEXT k     
430 END

2項分布B(n,p)の累積確率分布

100 OPTION ARITHMETIC NATIVE
110 DECLARE EXTERNAL FUNCTION C
120 LET n=180
130 LET p=1/6
140 LET t=0
150 FOR k=0 TO n
160    LET t=t+C(n,k)*p^k*(1-p)^(n-k)
170    PRINT USING "###  #.######":k,t
180 NEXT k
190 END
200 ! 組合せの数(2項係数)を計算する外部関数定義
210 EXTERNAL FUNCTION C(n,r)
220 OPTION ARITHMETIC NATIVE
230 IF r=0 THEN LET C=1 ELSE LET C=C(n-1,r-1)*n/r
240 END FUNCTION

このプログラムが使えるのはn=1020まで。それよりnが大きいと2項係数の計算で桁あふれが生じる。

2項分布B(n,p)の累積確率分布(改良版)

100 OPTION ARITHMETIC NATIVE
110 LET n=1200
120 LET p=1/6
130 LET t = (1-p)^n       ! t 累積確率
140 LET u = 0             ! u 2項係数の対数             
150 FOR k = 1 TO n
160    LET u = u + LOG2(n-k+1)-LOG2(k)
170    LET t = t + 2^( u + k*LOG2(p) + (n-k)*LOG2(1-p))
180    PRINT USING "#####  #.######":k,t
190 NEXT k
200 END

対数を利用して桁あふれを抑えた改良版。大きなnに対しても使える。

その他，いくつかのプログラム例をこちらで公開しています。

Ver. 7.2.2より，確率分布を計算するためのライブラリを追加しました。
PROBDIST.LIBに，正規分布，t分布，χ²分布，F分布を，
DiscDist.LIBに，二項分布，超幾何分布，ポアソン分布，負の二項分布を収録しています。
これらのファイルは，BASICをインストールしたフォルダのLibraryサブフォルダにあります。

使用例　3000個の不良品を含む18000個の製品の山から6個非復元抽出したときの不良品の個数Xの累積確率分布Pr(X≦i)

10 DECLARE EXTERNAL FUNCTION HyperGeomLCum
20 DECLARE NUMERIC NN,M,n,i
30 LET NN=18000
40 LET M=3000
50 LET n=6
60 FOR i=0 TO n 
70    PRINT i,HyperGeomLCum(NN,M,n,i)
80 NEXT i
90 END
100 MERGE "discdist.lib"