05 散布図と相関係数

２次元のデータを扱います。

題材は今まで同様，新体力テストのデータとします。

散布図を描くものを少し汎用に作って，関数として定義しましょう。引数はいろいろ考えられますが，ひとまずはシンプルに２つのデータのリストと，全体のタイトル，とします。また，データの最初は項目名があるものとします。スプレッドシートでいうと，１行目が項目名の表で，列ごと選択することに相当します。

では，すこしずつ作っていきます。

関数の名前と引数，エラーチェック

scatterplot(data1,data2,title):=(

areasize=20; // 表示領域のサイズ

textsize(16);

if(length(data1)!=length(data2),

drawtext([0,5],"２つのデータの個数が異なります");

,

areasize は，描画面で散布図を表示する領域の大きさです。あとで変えやすいように最初に書いておきます。

textsize(16) もあとで変えやすいようにはじめに書いておきます。

その次がエラーチェックです。２つのデータの個数が合わないときはメッセージを出すようにしています。

本体は，コンマのあとに書いていきます。

項目名の分離

引数のリストの先頭には項目名が入っている前提です。これを name , name2とし，-- 演算子で削除します。

また，データの個数 datan を求めておきます。

name1=data1_1; // 項目名

name2=data2_1;

data1=data1--[name1];

data2=data2--[name2];

datan=length(data1);

データ範囲などの設定

散布図の作成で一番面倒なのはここです。汎用に作ろうとしているので，与えられたデータからその範囲を求めて，プロットする位置を決めるのですが，新体力テストのデータを見ればわかる通り，小数もあれば２桁，３桁のデータもあります。ある程度は試行錯誤しながら作っていくことになります。いくつか想定してみましょう。表示する範囲を０〜20 としたとき，どのように配置するかです。

・50m走 小数で 6.3 から 8 まで ： 6.0 から 8.0 まで，範囲２とするとちょうどよい

・立ち幅跳び 173から 282 まで ： 170 から 290 とするか 160 から 300 とするか

・反復横跳び 45から74 ： 45 から 75 とするか

人間なら最大値と最小値を見て適当に決め，具合が悪ければ修正することができますが，これを自動でやろうとしているわけです。

まず，データの桁数を見てから判断しましょう。立ち幅跳びの場合，172から 284 などとすることはまずなくて，170 から 290 のように10の倍数とするでしょう。そのためには桁数を調べて端数を切ることになります。

ちょっと本体から外れて，次の関数を定義しておきます。

// 指数表記にしたときの指数部を求め，実際より±１増減する

exponent(n):=(

regional(k);

k=0;

if(n>=10,

while(n>=10,n=n/10;k=k+1);

k=k-1;

);

if(n<1,

while(n<1,n=n*10;k=k-1);

k=k+1;

);

k;

);

±１増減する意味は，反復横跳びのように２桁ならそのままで，３桁になったら１０の位を切るためです。

本体に戻ります。データの下限 min1,min2 と 上限 max1 , max2 を決め，範囲 range1 , range2 を求めます。

この結果，幅 20 の目盛にぴったり入るとは限らないので（範囲が９となるなど）方眼に合わせて目盛をとるのではなく，範囲に合わせて目盛をとることにします。それが unit1 ，unitscale1 です。目盛を打つときの幅が scalex,scaley (x軸方向，y軸方向)です。

// 目盛幅を調整するため指数表記で整数化して戻す

e=exponent(data1_1);

max1=ceil(max(data1)/10^e)*10^e;

min1=floor(min(data1)/10^e)*10^e;

range1=max1-min1;

if(min1>10,

unit1=ceil(range1/10);

,

unit1=range1/10;

);

unitscale1=unit1*areasize/range1; // 1目盛の値

e=exponent(data2_1);

max2=ceil(max(data2)/10^e)*10^e;

min2=floor(min(data2)/10^e)*10^e;

range2=max2-min2;

if(min2>10,

unit2=ceil(range2/10);

,

unit2=range2/10;

);

unitscale2=unit2*areasize/range2;

scalex=areasize/range1; // １目盛の長さ

scaley=areasize/range2;

軸と目盛を描く

ここまでで基本的な設定はできたので，あとはその値を使って表示していくだけです。

// 目盛

apply(1..10,drawtext([#*unitscale1,-0.8],min1+#*unit1,align->"center"));

apply(1..10,drawtext([-1.2,#*unitscale2-0.2],min2+#*unit2));

draw([[0,0],[areasize+1,0]],color->[0,0,0]);

draw([[0,0],[0,areasize+1]],color->[0,0,0]);

forall(1..10,draw([[#*unitscale1,0],[#*unitscale1,0.2]],color->[0,0,0]));

forall(1..10,draw([[0,#*unitscale2],[0.2,#*unitscale2]],color->[0,0,0]));

グリッドを表示するかどうかは迷うところ。ひとまずスクリプトは書いておいて，先頭に // を書いてコメント化しておきましょう。

// グリッド

// forall(1..10,draw([[#*unitscale1,0],[#*unitscale1,areasize]],dashtype->3,color->[0,0,0]));

// forall(1..10,draw([[0,#*unitscale2],[areasize,#*unitscale2]],dashtype->3,color->[0,0,0]));

// タイトルと項目名

drawtext([5,20],title,size->20);

drawtext([21,-1],name1);

drawtext([-3,21],name2);

プロットする

// 点を打つ

repeat(datan,

draw([(data1_#-min1)*scalex,(data2_#-min2)*scaley]);

);

相関係数の算出

公式に従って相関係数を算出します。平均値にラインを引くかどうかも迷うところ。これもスクリプトは書いておいてコメント化。また，回帰直線も求めておいて，ひとまずはコメント化しておきましょう。

// 相関係数算出

mean1=sum(data1)/datan;

mean2=sum(data2)/datan;

sx=0;

sy=0;

sxy=0;

repeat(datan,

sx=sx+(data1_#-mean1)^2;

sy=sy+(data2_#-mean2)^2;

sxy=sxy+(data1_#-mean1)*(data2_#-mean2);

);

correlation=sxy/(sqrt(sx)*sqrt(sy));

// 平均値にラインを引く

// draw([[(mean1-min1)*scalex,0],[(mean1-min1)*scalex,areasize]],size->1,dashtype->1,color->[0,1,0]);

// draw([[0,(mean2-min2)*scaley],[areasize,(mean2-min2)*scaley]],size->1,dashtype->1,color->[0,1,0]);

drawtext([15,19],"相関係数："+correlation);

// 回帰直線 y=sxy/sx*(x-mean1)+mean2

yy1= (sxy/sx*(min1-mean1)+mean2-min2)*scaley;

yy2= (sxy/sx*(max1-mean1)+mean2-min2)*scaley;

// draw([[0,yy1],[areasize,yy2]],size->2);

); // 全体を関数定義とした最初のif文のとじかっこ

戻り値

最後に戻り値です。普通は散布図がかけてしまえば戻り値はいらないのですが，いまは最小値と範囲を戻り値としておきます。すると，これを利用して，散布図上にプロットされている点をクリックしてそのデータを表示することができるのです。

// 戻り値

[min1,range1,min2,range2];

);

ローカル変数の宣言

最後に，はじめに戻って，ここまでに使った変数をローカル変数として宣言しておきます。

regional(・・・)はローカル変数の宣言ですが，これは

regional(datan,e,areasize,max1,min1,range1,max2,min2,range2,scale1,scale2,

mean1,mean2,sx,sy,sxy,yy1,yy2);

scatterplot() を使う

Drawスロットに戻って，引数のリストを作って散布図を表示します。

データはこれまでと同様です。

if(Orgdata_(-1)==unicode("000a"),

Orgdata=substring(Orgdata,0,length(Orgdata)-1);

);

datalist=tokenize(Orgdata,[unicode("000a"),","]);

でデータリストを作るのでした。

どの種目をとるか，やりやすいように変数にしておきます。列データは column() でリストにできます。

Datano1=2;

Datano2=9;

data1=column(datalist,Datano1);

data2=column(datalist,Datano2);

scatterplot(data1,data2,"新体力テスト");

これで散布図が描けます。

クリッカブル散布図もしくは検索可能な散布図

プロットされた点をクリックすると，そのデータが示されるようにします。

たとえば，このデータは新体力テストの結果ですが，上図で左下にちょっと外れた値があるのがわかります。これはどの生徒なのかを知りたいとします。普通なら，もとの表を目視で調べていくことになるでしょう。Excelなどの表計算ソフトでは１つの値の検索は簡単にできますが，２つの値はなかなか面倒です。

原理は簡単です。クリックした位置からデータ範囲を求め，リスト処理で該当するものを選択すればよいのです。そのために，戻り値を使うのです。

まず，Mouse Click スロットに，次の２行を書きます。

Mpx=mouse().x;

Mpy=mouse().y;

マウスをクリックすると，その座標が Mpx,Mpy に入ります。

Draw スロットに戻ります。

ret=scatterplot(data1,data2,"新体力テスト");

min1=ret_1;

range1=ret_2;

min2=ret_3;

range2=ret_4;

areasize=20;

xa=min1+Mpx/areasize*range1-range1/areasize/2;

xb=min1+Mpx/areasize*range1+range1/areasize/2;

ya=min2+Mpy/areasize*range2-range2/areasize/2;

yb=min2+Mpy/areasize*range2+range2/areasize/2;

plist=select(datalist,xa<=#_Datano1 & #_Datano1<=xb & ya<=#_Datano2 & #_Datano2<=yb);

stable=[[datalist_1_1,datalist_1_Datano1,datalist_1_Datano2]];

stable=stable++apply(plist,[#_1,#_Datano1,#_Datano2]);

drawtable([21,0],stable,width->80);

ret が戻り値です。ここから，最小値と範囲を取得します。

xa,xb,ya,ub は検索する範囲です。最後の /2 は±0.5 の目盛の範囲としています。ここは適当に。

範囲がきまったら，select() で該当する要素を選択します。

stable が結果の表ですが，まず，項目名を入れ，それから生徒番号と，２つのデータを入れています。

最後に，drawtable() で，できたリストを表示します。

さきほどの，左下にちょっと外れた点をクリックすると次のように表示されます。

中央付近の，ちょっと密集したところをクリックすると，３つの値が表示されます。

このページの最後に，scattarplot() の全体を再掲しておきます。

// 散布図を描く データリストの先頭には項目名があるものとする

scatterplot(data1,data2,title):=(

regional(datan,e,areasize,max1,min1,range1,max2,min2,range2,

scale1,scale2,mean1,mean2,sx,sy,sxy,yy1,yy2);

areasize=20; // 表示領域のサイズ

textsize(16);

if(length(data1)!=length(data2),

drawtext([0,5],"２つのデータの個数が異なります");

,

// 以下が処理

name1=data1_1; // 項目名

name2=data2_1;

data1=data1--[name1];

data2=data2--[name2];

datan=length(data1);

// 目盛幅を調整するため指数表記で整数化して戻す

e=exponent(data1_1);

max1=ceil(max(data1)/10^e)*10^e;

min1=floor(min(data1)/10^e)*10^e;

range1=max1-min1;

if(min1>10,

unit1=ceil(range1/10);

,

unit1=range1/10;

);

unitscale1=unit1*areasize/range1; // 1目盛の値

e=exponent(data2_1);

max2=ceil(max(data2)/10^e)*10^e;

min2=floor(min(data2)/10^e)*10^e;

range2=max2-min2;

if(min2>10,

unit2=ceil(range2/10);

,

unit2=range2/10;

);

unitscale2=unit2*areasize/range2;

scalex=areasize/range1; // １目盛の長さ

scaley=areasize/range2;

// 目盛

apply(1..10,drawtext([#*unitscale1,-0.8],min1+#*unit1,align->"center"));

apply(1..10,drawtext([-1.2,#*unitscale2-0.2],min2+#*unit2));

draw([[0,0],[areasize+1,0]],color->[0,0,0]);

draw([[0,0],[0,areasize+1]],color->[0,0,0]);

forall(1..10,draw([[#*unitscale1,0],[#*unitscale1,0.2]],color->[0,0,0]));

forall(1..10,draw([[0,#*unitscale2],[0.2,#*unitscale2]],color->[0,0,0]));

// グリッド

// forall(1..10,draw([[#*unitscale1,0],[#*unitscale1,areasize]],dashtype->3,color->[0,0,0]));

// forall(1..10,draw([[0,#*unitscale2],[areasize,#*unitscale2]],dashtype->3,color->[0,0,0]));

// タイトルと項目名

drawtext([5,20],title,size->20);

drawtext([21,-1],name1);

drawtext([-3,21],name2);

// 点を打つ

repeat(datan,

draw([(data1_#-min1)*scalex,(data2_#-min2)*scaley]);

);

// 相関係数算出

mean1=sum(data1)/datan;

mean2=sum(data2)/datan;

sx=0;

sy=0;

sxy=0;

repeat(datan,

sx=sx+(data1_#-mean1)^2;

sy=sy+(data2_#-mean2)^2;

sxy=sxy+(data1_#-mean1)*(data2_#-mean2);

);

correlation=sxy/(sqrt(sx)*sqrt(sy));

// 平均値にラインを引く

// draw([[(mean1-min1)*scalex,0],[(mean1-min1)*scalex,areasize]],size->1,dashtype->1,color->[0,1,0]);

// draw([[0,(mean2-min2)*scaley],[areasize,(mean2-min2)*scaley]],size->1,dashtype->1,color->[0,1,0]);

drawtext([15,19],"相関係数："+correlation);

// 回帰直線 y=sxy/sx*(x-mean1)+mean2

yy1= (sxy/sx*(min1-mean1)+mean2-min2)*scaley;

yy2= (sxy/sx*(max1-mean1)+mean2-min2)*scaley;

// draw([[0,yy1],[areasize,yy2]],size->2);

); // 全体を関数定義とした最初のif文の綴じかっこ

// 戻り値

[min1,range1,min2,range2];

);

// 指数表記にしたときの指数部を求め，実際より±１増減する

exponent(n):=(

regional(k);

k=0;

if(n>=10,

while(n>=10,n=n/10;k=k+1);

k=k-1;

);

if(n<1,

while(n<1,n=n*10;k=k-1);

k=k+1;

);

k;

);

参考：スピアマンの順位相関係数

GeoGebraでは，スピアマンの順位相関係数も算出しているので，ちょっと紹介しておきましょう。

実データの代わりに順位を用いて相関係数を算出するのです。ただし，順位はよくある順位（同点は同順位）と異なり，同点の場合は同点の数だけ均等割にします。

たとえば，よく使われる順位は，1,2,3,3,5,6,6,6,9,・・ですが，1,2,3.5,3.5,5,7,7,7,9,・・ となっていきます。スピアマンの相関係数の式も２通りあるようですが，GeoGebraでは順位をそのままピアソンの相関係数と同じ式で求めています。

ちょっとやってみましょう。

データリストから，上記のような順付けをします。通常の順付けと異なるので若干手間がかかります。

まずスピアマンの順位相関係数のための順位付けをしてそのリストを返す関数を定義します。

makeord(list):=(

regional(dn,dt,d,c,od,g,ret);

dn=length(list);

// データに番号を付けて昇順にソート

dt=apply(1..dn,[#,list_#]);

dt=sort(dt,#_2);

// 順位のリストを作る。同順位のものをカウントし，順位と個数のリストを作る

d=1;

c=1;

od=[];

repeat(dn-1,s,

if(dt_s_2==dt_(s+1)_2,

c=c+1;

,

od=append(od,[d,c]);

c=1;

d=s+1;

);

);

od=append(od,[d,c]);

// 同順位のものを均等割の数に変換

g=[];

repeat(length(od),s,

g=append(g,[sum(od_s_1..(od_s_1+od_s_2-1))/od_s_2,od_s_2]);

);

// できた順位をデータの数だけリストにする

od=[];

repeat(length(g),s,

od=od++apply(1..(g_s_2),g_s_1);

);

// 各データに順位を追加して番号順に戻す

ret=apply(1..dn,dt_#++od_#);

ret=sort(ret,#_1);

// 戻り値は順位だけのリスト

apply(1..dn,ret_#_3);

);

これを用いて，実データを順位に変換し，

d1=makeord(data1);

d2=makeord(data2);

このあとは相関係数を出す手順と同じ方法を用います。

データを順位に変換したあと，スピアマンの順位相関係数はもう一つ式がありますので，そちらでやってみるのもよいでしょう。