teacup. [ 掲示板 ] [ 掲示板作成 ] [ 有料掲示板 ] [ ブログ ]

 投稿者
  題名
  内容 入力補助画像・ファイル
    
 URL
[ ケータイで使う ] [ BBSティッカー ] [ 書込み通知 ] [ 検索 ]


2×2分割表再訪:12 その3 片側・両側確率の明示など

 投稿者:加藤 厚  投稿日:2018年 5月25日(金)07時20分56秒 em119-72-193-10.pool.e-mobile.ne.jp
   11のfunction 正確()は順不問の観察度数入力から可能な全パターンと確率の全体像を示します。その推移の中に観察パターンの位置と片側及び両側の確率が明示されれば一層便利でしょう。

 前節の画像にも示されているように、偏りの大きいパターンほど確率は小さいため、値は低く始まって中間部分で高まりその後低下します。従って、最も偏りが大きいa=0から観察度数パターンまでの確率の和が片側確率、その後観察度数パターンのそれを下回った確率からaの最大値の確率までの和が逆側の片側確率、両者の和が両側確率です。そこで、並べ替え後の観察度数パターンとその確率をまず示し、続いてパターンと確率の推移を「観察度数パターンまでの片側確率」のラベルを含めて示し、最後に逆側の片側確率と両側確率を提示するよう加筆したfunction 正確()は以下の通りです:

  if(d==m){d=a;x=b;b=c;c=x;a=m}(この行までは11と同じ)
  終=a+b;if(c<b)終=a+c;出力=出力+"\n\n観察度数パターン\n"  //見出し
  e=a+b;f=c+d;g=a+c;h=b+d;n=g+h  //入替を反映
  定数=階(e)*階(f)*階(g)*階(h)/階(n);観察=r8(定数/階(a)/階(b)/階(c)/階(d))
  出力=出力+"("+a+","+b+","+c+","+d+")\tp="+観察+"\n\n可能な全パターン\n"
  両=0;片=0;l="■片側確率="  //0は初期値 lはラベル ↑見出し
  for(j=0;j<=終;j++){
    確率=r8(定数/階(j)/階(e-j)/階(g-j)/階(f-g+j))
    if(確率<=観察){両=両+確率;片=片+確率}  //観察度数の確率以下なら加算
    出力=出力+"("+j+","+(e-j)+","+(g-j)+","+(f-g+j)+")\tp="+確率+"\n"
    if(a==j){出力=出力+l+r8(片)+"\n";片=0}
  }  //最初の片側確率の追記と片の初期化↑
  出力=出力+l+r8(片)+"\n□両側確率="+r8(両)
}

 まず観察度数パターンという見出しと入替後の(a,b,c,d)及びその確率=観察を出力に追加します。続いてforループ内で観察度数の確率以下の確率を変数:両と片に加算し、最初にa==jつまり可能パターンが観察度数と一致した時点でラベル(l)■片側確率=とそれまでの確率の和(片)を追記して変数:片を0に戻します。ループが終了した時点で変数:片に加算されている逆側の片側確率と両側確率をそれぞれのラベル付きで出力に追記して入出力欄に代入します。

 変更を含む終=a+b;‥以下の12行を差し替えて2×2分割表12.htmなどとして保存・実行すると下の画像の結果が表示されます。これで、2×2分割表において観察されたある度数パターンが、周辺度数同一という条件の下で可能な全パターン中のどこに位置し、またそれ以上に偏ったパターンも含めた発生確率の合計はどれほどかを、両側確率及び逆側の片側確率も含めて明確に把握できるようになりました。

 結局、同一周辺度数で可能な全12パターン中、観察された「成果なし1人」とより偏った「成果なし0人」の確率は約.024と約.002で計約.026、つまり危険率5%で有意(p<.05)でした。08で示したように「オッズ比の信頼区間ではp<.05、修正なしのχ2乗ではp<.05、修正χ2乗でもp<.05」で、今回の例の場合はどの方法でも同様の結果が得られています。また、09で示したようにピアソンのχ2乗の片側確率は.02274の半分で約.011、連続修正後のそれは.05927の半分で約.030ですから、直接確率(約.026)に最も近い値を示す近似方法は確かに「連続修正したχ2乗」でした。

 
 

2×2分割表再訪:11 その2 可能な全度数パターンとその確率

 投稿者:加藤 厚  投稿日:2018年 5月18日(金)07時42分39秒 em119-72-197-116.pool.e-mobile.ne.jp
編集済
   10のfunctionで正しい値を得るにはa~d中の最小値をbの位置に入力する必要があり面倒です。また、観察パターンを含む可能な全パターンとその確率の一覧は「全体像を踏まえた理解」を促します。そこで順不問の観察度数入力を処理して可能な全パターンとその確率を一覧的に示すscriptを工夫してみます。

 入力順については「a:bとc:dの各対の反転によるaの最小化」で対処します。それは以下の4類型の対応により実現できます:
1. 元々の入力においてaが最小ならそのまま
2. 〃 bが最小なら両対の左右を入れ替える:b⇔aとd⇔c
3. 〃 cが最小なら両対の上下を入れ替える:c⇔aとd⇔b
4. 〃 dが最小なら両対の上下と左右を入れ替える:d⇔aとc⇔b

 最小値をaに置く理由は、度数パターン(a,b,c,d)の先頭のaを0から可能な最大値※まで変化させることで「度数の推移」を最も明瞭に示せ、推移の把握が容易そうだからです(下の画像参照)。
※aの可能な最大値はa+bとa+cの小さい方と一致します。なぜならaがその最大値をとる時、b・cのいずれかが0になるからです。

 以上の処理を追記したfunction 正確()などは以下の通りです。

function 正確(){
  m=a;if(b<m)m=b;if(c<m)m=c;if(d<m)m=d
  if(b==m){b=a;x=d;d=c;c=x;a=m}  //左右
  if(c==m){c=a;x=d;d=b;b=x;a=m}  //上下
  if(d==m){d=a;x=b;b=c;c=x;a=m}
  終=a+b;if(c<b)終=a+c;e=a+b;f=c+d;g=a+c;h=b+d;n=g+h  //入替を反映
  出力=出力+"\n\n可能な全パターン\n";定数=階(e)*階(f)*階(g)*階(h)/階(n)
  for(j=0;j<=終;j++){
    確率=r8(定数/階(j)/階(e-j)/階(g-j)/階(f-g+j))
    出力=出力+"("+j+","+(e-j)+","+(g-j)+","+(f-g+j)+")\tp="+確率+"\n"
  }
}
function r8(x){return Math.round(100000000*x)/100000000}
function 階(x){r=1;for(i=1;i<=x;i++){r*=i}return r}

 まずaを最小値mに初期値として代入し、以後b,c,dと比べてより小さい値があればそれをmに代入します。次にmがb,c,dのいずれかであったら順に左右・上下・左右&上下を入れ替えて新しいa~dなどを定めます。また、aが取りうる最大値(a+bとa+cの小さい方)を変数:終に代入します。

 パターン発生確率の計算式:e!*f!*g!*h!/(n!*a!*b!*c!*d!)中のe!*f!*g!*h!/n!は周辺度数が同一の場合は一定なのでまずその値を求めて定数とし、aの値が0から最大値に至る各パターン※の確率をforループで順次求めてパターン(a,b,c,d)と共に出力に追記します。
※2×2の場合、同一周辺度数ならaが決まればb,c,dも決まります(自由度1)。従って、j=aに対してa,b,c,dはその関数(j,e-j,g-j,f-g+j)になります。

 前節の正確()と置き換え、小数第9位を四捨五入して8位に丸めるfunction r8()を追加して2×2分割表11.htm などとして保存・実行すると下の画像の結果が表示されます(前節から入出力欄を21 行で打ち止めとしたため表示には要スクロール)。

 例の度数パターン:10,1,10,10が1,10,10,10に並べ替えられ、最小値の1をaとして同一周辺度数で可能なaの範囲0~11の全パターンとその確率※が表示されています。パターンは10(前節)の「全度数パターン」と一致し、例はaが最も少ない側から2番目のパターンであることが確認できます。
※1e-8は1÷(10の8乗)=0.00000001

 

2×2分割表再訪:10 正確検定:その1 観察度数パターンなどの確率

 投稿者:加藤 厚  投稿日:2018年 5月11日(金)07時18分24秒 em119-72-193-233.pool.e-mobile.ne.jp
編集済
   さて、χ2乗に関する07~09の内容につきまとってきた「スッキリしない感じ」を、08で予告した“本道”つまり直接確率法で晴らす時が来ました(‥07の「コラム」の方法も直接確率法です)。

 フィッシャーの正確検定は以下の3段階で「観察された度数以上の偏りの発生確率」を直接求める方法です。
1. 観察された度数パターンの発生確率を求める。
2. 同一周辺度数(行と列の小計)において1.のパターンと連続する「確率がより小さいパターン」の確率の総和を1.の値と合計する(これが「片側確率」)。
3. 同一周辺度数で起こりうる度数パターンの確率を全て求め、1.の値以下の値を総和する(これが「両側確率」、その値-2.の値が逆側の「片側確率」)。

 まず1.、つまり例のa=10:b=1対c=10:d=10のパターン(下の画像:上側の黄色)の確率を求めましょう。新しい方法で指導された計11人から10人を選ぶ組合せ(同緑色)は(1人を選ぶ組合せと同じ)11通り、従来の方法で指導された計20人から10人を選ぶ組合せ(20C10:同紫色)は20!/(20-10)!/10!で184756通り(←nCr=n!/(n-r)!/r!)。前者の組合せと後者のそれは独立なので両者の組合せは結局11×184756=2032316通りあります。

 他方、31人から(a+cの)20人を選ぶ全ての組合せ(同灰色)は(11人を選ぶ組合せと同じ)31C20なので31!/(31-20)!/20!=84672315通り。

 度数パターンの発生確率は「該当する全ての組合せ÷起こりうる全ての組合せ」なので(11C10)×(20C10)÷(31C20)=11×184756÷84672315=2032316÷84672315≒0.024で約2.4%。そして、この値は計算式:e!*f!*g!*h!/(n!*a!*b!*c!*d!)でも求められます。

 次に2.、つまり観察されたパターンに連続する確率がより小さい=より偏ったパターンの全確率を求めます。同一周辺度数において1.のパターンから連続するより偏ったパターンは11,0,9,11の1つのみ(同水色)でその確率は11!*20!*20!*11!/31!/11!/0!/9!/11!≒0.00198で約0.2%。従って、片側確率は約2.6%(←約2.4+約0.2)となり、結局、オッズ比の信頼区間の.01<p<.05やピアソンのχ2乗及び連続修正χ2乗の.01<p<.03(両側確率の半分の値)と一貫する値が得られました。

 ここまでの過程を行うfunctionの一例※は以下の通りです:
※a~d中bが最小を前提としたscriptなので他のパターンには非対応!

function 正確(){
  定数=階(e)*階(f)*階(g)*階(h)/階(n);出力=出力+"\n\n度数パターン\n"
  for(j=b;0<=j;j--){
    p=Math.round(100000000*定数/階(e-j)/階(j)/階(g-e+j)/階(h-j))/100000000
    出力=出力+"("+(e-j)+","+j+","+(g-e+j)+","+(h-j)+")\tp="+p+"\n"
  }
}
function 階(x){r=1;for(i=1;i<=x;i++){r*=i}return r}

 まず関数funciton 階(x)で階乗を定義し、function 正確()では計算式の定数部分を計算→“度数パターン”と改行を出力に追加→b=1と0の場合の確率pを計算→(度数パターン)と共に出力に追加という処理を行っています。function 処理()のchi2();の後ろに正確();を追記の上2×2分割表10.htmなどとして保存・実行すると下の画像:下側の結果※が表示されます。
※出力を下にスクロールした状態

 

2×2分割表再訪:09 確率的有意性と近似確率の自動表示

 投稿者:加藤 厚  投稿日:2018年 5月 4日(金)09時41分13秒 em1-115-198-223.pool.e-mobile.ne.jp
編集済
   自由度が1のχ2乗分布の場合、上側確率10,5,2,1%に対応する値は順に約2.71,3.84,5.41,6.64です。また、偏りの一方を無視する片側検定の場合、具体的には対立仮説の内容が≠ではなくて>や<(例:新しい方法は従来の方法より有効)の場合、各値は半分の確率(順に5,2.5,1,0.5%)に対応します。

 毎回これらの値を参照するのは面倒なので、上側確率pをn.s.,p<.10,p<.05,p<.01などと自動表示させたいならfunction chi2()に以下の3行を追加し、変数:出力の定義の(CHI2)の後ろに+p、(YCC2)の後ろに+yと追記すればOKです。

s=new Array(" n.s."," p<.10"," p<.05"," p<.02"," p<.01")
x=CHI2;p=s[0];if(2.71<x)p=s[1];if(3.84<x)p=s[2];if(5.41<x)p=s[3];if(6.64<x)p=s[4]
x=YCC2;y=s[0];if(2.71<x)y=s[1];if(3.84<x)y=s[2];if(5.41<x)y=s[3];if(6.64<x)y=s[4]

 また、確率の近似値は、標準正規分布の累積確率密度の近似計算式(例:Box-Muller法)にχ2乗値の正・負の平方根を代入して得られる値の和(=両側確率)※として求められます(自由度1のχ2乗分布=標準正規分布の2乗)。
※実際は正の平方根を代入して得た片側確率×2=両側(上側)確率でOK。

 Box-Muller法の精度は以下のscriptの実行で確認できます。

<script>
a=1.64485;b=1.95996;c=2.57583;s=" "  //標準正規分布の上側確率 5,2.5,0.5%値
document.write(a+s+nd(a)+"<BR※>"+b+s+nd(b)+"<BR※>"+c+s+nd(c))
function nd(x){
  p=0.2316419;b1=0.31938153;b2=-0.356563782;b3=1.781477937;b4=-1.821255978
  b5=1.330274429;t=1/(1+p*Math.abs(x));Z=Math.exp(-x*x/2)/Math.sqrt(2*Math.PI)
  return Z*((((b5*t+b4)*t+b3)*t+b2)*t+b1)*t
}
</script>
※BRは本来は半角(BRでもbrでもOK)

 起動画面には下の画像:上側の各値が表示されます。左側の1.64‥~2.57‥の3値は数表から求めた標準正規分布の上側確率5,2.5,0.5%の値(05で信頼区間の計算にも使用)、それらをfunction nd(x)にxとして代入した戻り値(returnで返る値)が右側の0.0500‥、0.0250‥、0.00500‥の3値で、小数第6位まで0が続く高い精度を示しています。

 以上の自動表示機能を組み込んだfunctionは以下の通りです:

function chi2(){tmp=Math.abs(a*d-b*c);CHI2=n*tmp*tmp/(a+b)/(c+d)/(a+c)/(b+d)
  YCC2=n*(tmp-n/2)*(tmp-n/2)/(a+b)/(c+d)/(a+c)/(b+d);if(tmp<=n/4)YCC2=0
  me=e*g/n;if(e*h/n<me)me=e*h/n;if(f*g/n<me)me=f*g/n;if(f*h/n<me)me=f*h/n
  s=new Array(" n.s."," p<.10"," p<.05"," p<.02"," p<.01");l="\n上(両)側確率p="
  x=CHI2;p=s[0];if(2.71<x)p=s[1];if(3.84<x)p=s[2];if(5.41<x)p=s[3];if(6.64<x)p=s[4]
  x=YCC2;y=s[0];if(2.71<x)y=s[1];if(3.84<x)y=s[2];if(5.41<x)y=s[3];if(6.64<x)y=s[4]
  出力=出力+"\n\nピアソンのχ2乗="+r(CHI2)+p+l+r(2*nd(Math.sqrt(CHI2)))
  出力=出力+"\nχ2乗(連続修正)="+r(YCC2)+y+l+r(2*nd(Math.sqrt(YCC2)))
  出力=出力+"\n最小期待度数="+r(me)
}
function nd(x){  //Box-Muller法
  p=0.2316419;b1=0.31938153;b2=-0.356563782;b3=1.781477937;b4=-1.821255978
  b5=1.330274429;t=1/(1+p*Math.abs(x));Z=Math.exp(-x*x/2)/Math.sqrt(2*Math.PI)
  return Z*((((b5*t+b4)*t+b3)*t+b2)*t+b1)*t
}

textareaを rows="21"にして2×2分割表09.htmなどとして保存・実行すると下の画像:下側の結果が表示され、p<‥※と一貫する確率の近似値※が得られています。
※共にχ2乗の上側(=正規分布の両側)確率なので、対立仮説に方向性があれば値は半分#にできます。
#両側のp<.02/.10=片側のp<.01/.05

 

2×2分割表再訪:08 独立性のχ2乗検定と連続修正

 投稿者:加藤 厚  投稿日:2018年 4月27日(金)07時41分37秒 em119-72-196-94.pool.e-mobile.ne.jp
編集済
   07の適合度検定を「2つのカテゴリ変数は独立」と仮定した場合に拡張し、仮定が真の場合の期待値と実際の度数の差からχ2乗を求めるのが独立性のχ2乗検定です。実際の度数a,b,c,d(a+b=e,c+d=f,a+c=g,b+d=h,a+b+c+d=n)の期待値は、仮定が真なら総数nの均等配分値なので順にe*g/n[←n*(e/n)*(g/n)]、e*h/n、f*g/n、f*h/nです。そこで、以下のfunction chi2()を追加し、function 処理()の推定();の次にchi2();を追記の上、rows="17"にして例の10,1,10,10を処理すると「ピアソンのχ2乗=5.18802」と表示されます。

function chi2(){ea=e*g/n;eb=e*h/n;ec=f*g/n;ed=f*h/n
  CHI2=(a-ea)*(a-ea)/ea+(b-eb)*(b-eb)/eb+(c-ec)*(c-ec)/ec+(d-ed)*(d-ed)/ed
  出力=出力+"\n\nピアソンのχ2乗="+r(CHI2)
}

 例で検討したいのは新しい方法の有効性なので逆(1:10)側の偏りは無視すると、上側確率10%の約2.71と2%の約5.41(約2.326の2乗)が片側確率5%と1%の値となり、約5.19は2.71<5.19<5.41のため危険率5%で有意です。他方、a=10,b=1,c=10,d=10の期待値は順に、7.1,3.9,12.9,7.1でbのそれが5未満です。従って近似に関する2条件※中の「期待値5未満のセルが全体の20%以上あってはいけない」は満たされていません。
※これらが課される理由は、離散変数のχ2乗を連続分布で近似する証明にスターリングの公式(階乗の指数関数近似の公式)を使っているためだそうです。

 このような場合に「イェーツの連続修正」を勧める意見があります。逆に「修正は常に不要」という意見もあります。連続修正でフィッシャーの正確検定に「近づける」なら、正確検定で直接確率を求めるのが“本道”とも言えます(10 参照)。ここではその議論はさておいて処理を進めます。上のfunction chi2()の計算式は「考え方をそのまま式にした」定義式でした。理解しやすい点が長所ですが、式中の全値が揃っていないと複数段階の処理が必要、除算が多いと丸め誤差が大きくなるなどの短所があります。そこで、簡潔さと精度を重視して「最少の段階と除算で定義式と同じ答を得られるように工夫した式」が計算式です。χ2乗のそれは下のfunction chi2()の1行目、「差が0.5より大きければ0.5を引いてから2乗する」(差が0.5以下なら0とする)が定義のイェーツの連続修正を加えたχ2乗(YCC2)のそれは2行目です。最小期待度数(me)を求めて出力に追記するscriptを含む以下のfunction chi2()を以前のそれと置き換え、textareaをrows="19"にして2×2分割表08.htmなどとして保存します。

function chi2(){tmp=Math.abs(a*d-b*c);CHI2=n*tmp*tmp/(a+b)/(c+d)/(a+c)/(b+d)
  YCC2=n*(tmp-n/2)*(tmp-n/2)/(a+b)/(c+d)/(a+c)/(b+d);if(tmp<=n/4)YCC2=0
  me=e*g/n;if(e*h/n<me)me=e*h/n;if(f*g/n<me)me=f*g/n;if(f*h/n<me)me=f*h/n
  出力=出力+"\n\nピアソンのχ2乗="+r(CHI2)+"\nχ2乗(連続修正)="+r(YCC2)
  出力=出力+"\n最小期待度数="+r(me)
}

 起動して実行をクリックすると下の画像の結果が表示されます。連続修正したχ2乗の値は約3.55※で 2.71<3.55<5.41なのでやはり5%水準で有意です。
※連続修正はχ2乗の値を減らします。

 オッズ比の信頼区間ではp<.05(05参照)、修正なしのχ2乗ではp<.05、修正χ2乗でもp<.05で、今回の例については連続修正は「大まかな判断」を変えませんでした。

 

2×2分割表再訪:07 さて、やっとχ2乗ですが‥

 投稿者:加藤 厚  投稿日:2018年 4月20日(金)07時36分46秒 em1-115-195-178.pool.e-mobile.ne.jp
編集済
   統計量としてのχ2乗は(2×2以上の)分割表一般でそれを求められる点がオッズ比より優れています(χ2乗から求める連関の指標のクラメールのVは2×2分割表の場合にはφ係数と同値)。他方、引き出せる意味は、p値しか与えないχ2乗に対し(p値も判断できる※)信頼区間が示せるオッズ比の方が「豊か」です。
※95%信頼区間が1を含まない(例:1.2~4.8、0.6~0.8)ならp<.05で有意。なぜなら「母オッズ比=1(2条件に優劣なし)の確率が5%未満」だから。

 加えて、χ2乗検定には「期待度数が1未満のセルがあってはいけない」とか「期待度数が5未満のセルが全体の20%以上あってはいけない」(以上は近似関連)とか「イェーツの修正の使用が望ましい」とか逆に「使用は不要」(以上は正確検定との関連)とか様々な条件・意見が付随しており「スッキリしない感じ」がつきまといます。が、(順序として)まずその来歴と求め方などを紹介しましょう。

 1897年にゴールトン(ダーウィンの従兄弟で「指紋の独自性」の発見者としても有名)からロンドンの生物測定研究所を引き継いだカール・ピアソンの功績の1つが「適合度検定」(当てはまりのよさの検定)の考案です。人間を含む生物について様々な測定を行っていた彼らには、収集できた資料(標本)がその全体(母集団)の実態からあまり偏っていないことを示す必要がありました。そこでまず「理論から期待される分布と実際の標本における分布とのズレ」を評価する指標としてのχ2乗を発案したワケです。

 簡単な例として「コイン投げ」を取り上げましょう。10回あるいは10個投げて、表と裏の割合の偏りがどの程度までなら(表裏半々の意味で)「正しいコイン」と言えるでしょう‥5:5なら期待通り、4:6や6:4もアリガチ、では3:7や2:8ならどうでしょうか?

 カテゴリ数が2のこの問題は組合せで容易に解けます(下のコラム参照)が、3以上の場合にも簡単に近似値を得るための方法※が、理論から期待される度数(期待度数)と標本における実際の度数(観察度数)の差の2乗を期待度数で割った値の総和(=χ2乗)を指標とする適合度検定です。
※様々な場合(自由度)のχ2乗分布の確率数表が整備されています。

 今回のχ2乗は5:5なら0、4:6なら(-1*-1)/5+(1*1)/5=0.4、3:7なら(-2*-2)/5+(2*2)/5=1.6、2:8なら(-3*-3)/5+(3*3)/5=3.6、1:9なら(-4*-4)/5+(4*4)/5=6.4‥です。他方、カテゴリ数が2=自由度1のχ2乗分布の5%値は標準正規分布の上・下側確率2.5%値約±1.96(05の「信頼限界」参照)の2乗の約3.84です。従って、偏りが2:8では3.6<3.84@両側検定のため危険率5%では「正しくない」とは言えません。

 ここで要注意なのは-∞~+∞に分布する標準正規分布(下の画像の左側)を2乗した値は0~+∞に分布(同右側)するため、その上側5%には標準正規分布の上・下側が2.5%ずつ含まれる点です。従って、前段落の判断は8:2側の偏りも含んでおり、2:8側の偏りだけに注目するなら上側10%値の約2.71(約1.645の2乗)が5%値に当たるため、2.71<3.6で「危険率5%で正しくない」と言えます(片側検定)。他方、この判断は組合せに基づく結果(p=.055:コラム参照)とは異なります(これが「近似の悪さ」の問題)。

■コラム:確率=組合せ数の割合
 10回中r回表が出る組合せの数は10Cr※なのでr=0から順に1,10,45,120,210,252,210,120‥、組合せ総数1024(2の10乗)で割った確率は0から順に約.001,.010,.044,.117,.205,.246,.205,.117‥。従ってそれ以上の偏りが偶然起こる確率は2:8や8:2で約.055(=.044+.010+.001)、1:9や9:1で約.011です。
※nCr=n!/(n-r)!/r!

 

2×2分割表再訪:06 「実例」で再確認

 投稿者:加藤 厚  投稿日:2018年 4月13日(金)07時20分23秒 em117-55-68-171.emobile.ad.jp
   これまでの内容を実例で再確認してみましょう。ビタミンが発見されていない明治時代、脚気は原因不明の病気でした。そんな中、海軍医官高木兼寛は食事に注目し、白米に大麦を混ぜるなどの変更を加えて1884年に287日間の遠洋航海を行い、前年に同様の航海を行った場合と発病状況を比較しました。その内訳は以下の分割表の通りです(杉 2013 pp.96-97)※。
※杉 晴夫 2013 「栄養学を拓いた巨人たち」 講談社 ブルーバックス B-1811

      非発病  発病  小計
  ───────────────
  改善食  319   14   333
  従来食  207   169   376
  ───────────────
  小 計  526   183   709

 これまでのパターンに従って上の行が望ましい条件(改善食)、左側が望ましい結果(非発病)とすると、入力順は改善食&非発病,改善食&発病,従来食&非発病,従来食&発病で各数値は319,14,207,169です。例の10,1,10,10をこれらの値に置き換えて実行すると下の画像:上側の結果が得られます。

 リスク比(14/333÷169/376)の0.09は改善食により脚気発病の危険度が1/10以下に減少することを、オッズ比(319/14÷207/169=(319*169)÷(14*207))の18.6は従来食の約19倍の「非発病」効果が改善食に見込めることを示すものです。ユールの連関係数Qは.90と高く、効果量としてのφ係数も.46で「中程度の連関」が認められます。オッズ比の信頼区間は全て1.0(優劣なし)を大きく上回り、オッズ比の真の値(母オッズ比)が10.5~33.0の区間内にあることが95%の確率で推定できる結果が得られました※。
※a~dの値が大きいと標準誤差は小さくなり、より狭い(=特定的な)信頼区間が得られます。

 入出力欄にポインタ(矢印)を重ねると簡単な入力ガイドを表示する機能(title="表示内容")を追加したscriptを以下に提示します:

<form><textarea rows="15" cols="30" title="条件&結果:1&1,1&0,0&1,0&0の4値を入力">
10,1,10,10</textarea></form>
<input value="実行" type="button" onClick="処理()">
<script>
白=new Array(""," ","  ","   ","    ","     ","      ")
function s(x){return 白[7-(x+"").length]+x}  //spacing
function r(x){return Math.round(100000*x)/100000}  //rounding
function 作表(){
  v=document.forms[0].elements[0].value.split(",")
  a=v[0]*1;b=v[1]*1;c=v[2]*1;d=v[3]*1;e=a+b;f=c+d;g=a+c;h=b+d;n=g+h
  出力=s(a)+s(b)+s(e)+"\n"+s(c)+s(d)+s(f)+"\n"+s(g)+s(h)+s(n)+"\n\n"
}
function 計算(){
  OR=(a/b)/(c/d);l=""
  if(0==a*b*c*d){OR=((a+.5)/(b+.5))/((c+.5)/(d+.5));l="(+.5修正)"}
  出力=出力+"リスク比="+r((b/e)/(d/f))+"\n"+"オッズ比="+r(OR)+l
  出力=出力+"\n\nユールのQ="+r((a*d-b*c)/(a*d+b*c))
  出力=出力+"\nφ係数="+r((a*d-b*c)/Math.sqrt(e*f*g*h))
}
function 推定(){
  LOR=Math.log(OR);v=1/a+1/b+1/c+1/d
  if(0==a*b*c*d)v=1/(a+.5)+1/(b+.5)+1/(c+.5)+1/(d+.5)
  p=new Array(90,95,99);z=new Array(1.64485,1.95996,2.57583)
  出力=出力+"\n\nオッズ比の信頼区間(CI)"+l
  for(i=0;i<3;i++){  //p(%)とz(値)のi番要素を代入して3回実行
    幅=z[i]*Math.sqrt(v);出力=出力+"\n"+p[i]+"%:"
    出力=出力+r(Math.exp(LOR-幅))+"~"+r(Math.exp(LOR+幅))}
}
function 処理(){作表();計算();推定();document.forms[0].elements[0].value=出力}
</script>

 このscriptを2×2分割表06.htmなどとして保存し、入出力欄にポインタ(矢印)を重ねると「入力ガイド」が下の画像:下側のようにポップアップします。

 

2×2分割表再訪:05 標本の結果から母数のある範囲を推定する

 投稿者:加藤 厚  投稿日:2018年 4月 6日(金)07時25分58秒 em119-72-192-26.pool.e-mobile.ne.jp
編集済
   これまで、2×2分割表の集計結果から意味を引き出す以下の数値指標の特徴などを検討して来ました:
①リスク比(Risk Ratio:RR)は「危険率がn倍」と表現できて分かりやすいが、事後に対照(control)群を設定する症例対照研究の場合「リスクなし群の件数」に影響されて値が変わってしまう。
②オッズ比(Odds Ratio:OR)は事後の群設定でも事前のそれと同一の値が得られるが、0~∞に分布するため値の意味が分かりにくい。
③オッズ比の意味は「値の範囲」を-1~+1に変換したユールの連関係数Qで補足できる。また相関係数の一種のφ係数=四分点相関係数はその「筋の良さ」により効果量として使える(03参照)。

 しかし、いずれを使うにせよ、その値の意味(例:「新しい方法」の有望さ)は今回の「標本」での結果にすぎません。“一般的な意味”を実証的に示すには追試(replication:同様の実験や調査の反復実施)が確実※ですが、それには大きな手間がかかります。
※例:追試のオッズ比が順に8,11,9,12,10,11‥なら「真の値」は約10。

 1つの標本の結果から“一般的な意味”をある確からしさで合理的に求める方法があったらとても便利なことでしょう。そんな(魔法のような?)方法が区間推定(Interval Estimation)です。

 区間推定では、1つの標本から得た統計量(例:標本平均)とその分布の推測された標準偏差(例:√不偏分散)から「真の値(例:母平均)を一定の確率(例:95%)で含む区間」(信頼区間=Confidence Interval=CI)を推定します。そして上記の指標の中ではオッズ比の区間推定が最も広く利用されているようです(‥というのは疫学や社会調査などの分野の話で、教育心理学ではまず見ませんが)。

 03でも触れたように、標本オッズ比は「1を中心として0以上∞未満に上側に広く裾を引いて分布」しますが、その値を自然対数に変換した対数標本オッズ比の分布は「分散(v)が1/a+1/b+1/c+1/dの正規分布」で近似できます。従って、対数標本オッズ比(Log OR)±1.96×√v※の値を指数変換して真数に戻せば、その2値(95%信頼限界)を両端とする範囲がオッズ比の95%信頼区間になります。
※母平均μを区間推定する場合の「標本平均±1.96×不偏標準偏差」と同様です。

 以上の処理※を行う以下のfunction 推定()を追加し、function 処理(){‥}内の計算();の後ろに推定();を追記し、<form>のtextareaのrows="10"をrows="15"にしたら2×2分割表05.htmなどとして保存します。起動→実行クリックで下の画像の結果が表示されます。
※95%信頼区間に加え90%と99%のそれらも求めています。

function 推定(){
  LOR=Math.log(OR);v=1/a+1/b+1/c+1/d
  if(0==a*b*c*d)v=1/(a+.5)+1/(b+.5)+1/(c+.5)+1/(d+.5)
  p=new Array(90,95,99);z=new Array(1.64485,1.95996,2.57583)
  出力=出力+"\n\nオッズ比の信頼区間(CI)"+l
  for(i=0;i<3;i++){  //p(%)とz(値)のi番要素を代入して3回実行
    幅=z[i]*Math.sqrt(v);出力=出力+"\n"+p[i]+"%:"
    出力=出力+r(Math.exp(LOR-幅))+"~"+r(Math.exp(LOR+幅))}
}
 オッズ比の真の値(母オッズ比)は指定した確率※で信頼区間内にあります。従って、信頼区間に1(同オッズ=優劣なし)を含まない「95%:1.07‥~93.43‥」という結果は「新しい方法は従来の方法よりも有効」という“一般的な意味”を95%以上の確率で示すものです。
※例:95%なら100回の試行中の95回。

 

2×2分割表再訪:04 桁揃えや四捨五入で体裁を整える

 投稿者:加藤 厚  投稿日:2018年 3月30日(金)07時37分31秒 em119-72-195-242.pool.e-mobile.ne.jp
編集済
   これまでは「結果」優先で、出力書式などには無頓着でした。しかし、表は数値右詰が読み易く、また無闇に多い小数点以下の桁数は無意味です。そこで03のscriptに表の数値右詰表示と小数第5位への丸め(第6位四捨五入)の機能(function=関数)を追加します。併せて、長くなってきた処理()を作表()と計算()にfunction化します。

<form><textarea rows="10" cols="30">10,1,10,10</textarea></form>
<input value="実行" type="button" onClick="処理()">
<script>
白=new Array(""," ","  ","   ","    ","     ","      ")
function s(x){return 白[7-(x+"").length]+x}  //spacing
function r(x){return Math.round(100000*x)/100000}  //rounding
function 作表(){
  v=document.forms[0].elements[0].value.split(",")
  a=v[0]*1;b=v[1]*1;c=v[2]*1;d=v[3]*1;e=a+b;f=c+d;g=a+c;h=b+d;n=g+h
  出力=s(a)+s(b)+s(e)+"\n"+s(c)+s(d)+s(f)+"\n"+s(g)+s(h)+s(n)+"\n\n"
}
function 計算(){
  OR=(a/b)/(c/d);l=""
  if(0==a*b*c*d){OR=((a+.5)/(b+.5))/((c+.5)/(d+.5));l="(+.5修正)"}
  出力=出力+"リスク比="+r((b/e)/(d/f))+"\n"+"オッズ比="+r(OR)+l
  出力=出力+"\n\nユールのQ="+r((a*d-b*c)/(a*d+b*c))
  出力=出力+"\nφ係数="+r((a*d-b*c)/Math.sqrt(e*f*g*h))
}
function 処理(){作表();計算();document.forms[0].elements[0].value=出力}
</script>

 右詰7桁表示化には、まず半角ブランク0~6個を要素とする配列:白を用意し、その(7-添え字xの文字としての長さ)番目の要素をxの前に付けて返す関数:function s(x)を追加します。この関数はs(x)で呼び出され、例えばaが10ならs(a)は文字10の長さ2を7から引いた5に対応する配列:白の5番目の要素”     ”+10を返します。

 小数第5位までの丸めは、添え字xを100000倍して四捨五入(Math.round)した値を100000で割って返す関数:function r()を追加してr(x)で呼び出せばOKです。例えばxが0.123456‥なら100000倍した12345.6‥を四捨五入して100000で割った0.12346を返します。

 03の処理()中の作表部分を作表()、計算部分を計算()にfunction化し、処理()の中にはそれぞれの名称:作表()と計算()及び出力指示行のみを書き並べます。<form>から</script>までの20行をコピーして2×2分割表04.htmといったファイル名で保存し、10,1,10,10のまま実行すると下の画像:上側の結果が、10,0,10,10にして実行すると下の画像:下側の結果が示されます。分割表の値は右詰、指標は小数第6位四捨五入で示され、読み易い出力が得られています。

【展開に関する補足】
この連載はオッズ比→χ2乗→直接確率‥と展開しますが、内容の「積み上げ的理解」には逆の展開が(実は)適切です。ということで、まずは個別の内容についてscriptで体験・確認し、15まで行ってから逆向き※に理解の「積み上げ」を試みることをオススメしておきます。
※「組合せの比率」に基づくマクネマー検定@15→「組合せの積」を使う直接確率@10→「標準正規分布とdf=1のχ2乗分布」に基づく“近似”としてのχ2乗検定@07→「標本オッズ比とその標準誤差」に基づく区間推定@05

 

2×2分割表再訪:03 値を常に定め分かりやすい範囲に収める

 投稿者:加藤 厚  投稿日:2018年 3月23日(金)07時27分48秒 em1-115-194-139.pool.e-mobile.ne.jp
   2群の設定方法※に影響されないオッズ比の欠点は2つ:
※事前(prospective)・事後(retrospective)・同時(cross-sectional)‥

1.(a/b)/(c/d)の式が示すように、b,c,dのどれかが0だと値が∞=不定になる点。この欠点はa~dのいずれか※が0ならa~dの全てに0.5を加えて近似値を求めることで補います(ウールフの修正)。
※a~dのいずれかが0だとオッズ比の標準誤差(05参照)が∞になるため。

2.値が、2つのオッズが等しい時の1を中心として0以上∞未満に上側に広く裾を引いて分布するため「値の意味が分かりにくい」点。この欠点は、オッズ比の分布を-1~+1に変換したユールの連関係数Q※や元々-1~+1に分布する相関係数の一種のφ(ファイ=四分点相関)係数#との併用で補います。
※(オッズ比-1)/(オッズ比+1)=Q=(a*d-b*c)/(a*d+b*c)
#φ=(a*d-b*c)/√((a+b)*(c+d)*(a+c)*(b+d))

 ウールフの修正とQ・φの計算を加えたfunction 処理()の一例:
function 処理(){
  v=document.forms[0].elements[0].value.split(",")
  a=v[0]*1;b=v[1]*1;c=v[2]*1;d=v[3]*1
  e=a+b;f=c+d;g=a+c;h=b+d;n=g+h
  出力=a+"\t"+b+"\t"+e+"\n"+c+"\t"+d+"\t"+f+"\n"+g+"\t"+h+"\t"+n
  OR=(a/b)/(c/d);l=""  //←初期値""で条件成立なら入替↓(l=label)
  if(0==a*b*c*d){OR=((a+.5)/(b+.5))/((c+.5)/(d+.5));l="(+.5修正)"}
  出力=出力+"\n\n"+"リスク比="+(b/e)/(d/f)+"\n"+"オッズ比="+OR+l
  出力=出力+"\n\nユールのQ="+(a*d-b*c)/(a*d+b*c)
  出力=出力+"\nφ係数="+(a*d-b*c)/Math.sqrt(e*f*g*h)
  document.forms[0].elements[0].value=出力
}

 例では「a~dのいずれかが0」を4値の積で判断した結果(=l)及びQとφの説明文字列・値を出力に追加しています。functionの内容を入れ替え、出力行数の増加に対応して<form>のrows="6"をrows="10"にしたら2×2分割表03.htmなどとして保存します。まず10,0,10,10にして実行すると、a~dに.5を加えた値に基づくオッズ比とその旨の注(+.5修正)が出力されます(下の画像:上側)。

 次に、リスクあり群の件数を半分にして5,0,5,10で実行しても(オッズ比が同じなので)Qの値は変わりません。他方、φの値はリスクあり群の件数及び周辺度数(e,f,g,h)に影響されて変わりますが、こちらには「効果量」としての価値※があります(下の画像:下側)。
※クラメールのVやピアソンの積率相関係数と値が一致します。

 リスクあり群を倍にした20,1,20,10、リスクなし群を倍にした10,2,10,20、また弱・中・強の連関に対応する10,5,5,10、10,3,3,10、10,1,1,10などで実行してみましょう(一般に|φ|<|Q|)※。
※「処理の繰り返し」には「再読み込み」(Windowsなら[F5]キー)が便利です。

 

2×2分割表再訪:02 成果あり:なしの度数に注目すると

 投稿者:加藤 厚  投稿日:2018年 3月16日(金)07時16分32秒 em117-55-68-183.emobile.ad.jp
   01の例では、新しい方法での成果の有無(a:b)は10:1、従来の方法のそれ(c:d)は10:10でした。これらの値から求められる指標の1つにリスク比(Risk Ratio:RR=相対危険)※があります。
※成果ではなくてリスク(危険)に注目する理由は、この指標が疾病などの発生条件の解明に取り組む「疫学」の分野で用いられてきたからです。

 そこで「指導を受けたのに成果なし」をリスクと見なすと、新しい方法での発生率は1/11、従来の方法でのそれは10/20で、リスク比(新しい方法のリスク÷従来の方法のリスク)は(1/11)/(10/20)≒0.18、つまり新しい方法は従来の方法と比べて「成果なし」に終わるリスク(危険)が5分の1未満であることが示されます。

 では「リスク比を指標にすればOKか」というと(例のように最初から件数の定まった2群に異なる条件=指導を与えた場合=介入研究ならOKなのですが)「リスク群がまずあり、その後に集めたリスクなし群との間で条件の有無を比べる」という資料の集め方=症例対照研究(case-control study)の場合にはNGです。なぜなら、リスク比は「リスクなし群の件数」によって変動してしまうからです。

 Scriptで計算して確認してみましょう。例ではa=10:b=1→e=a+b=11、c=10:d=10→f=20なのでリスク比(1/11)/(10/20)は(b/e)/(d/f)、そこでfunction 処理()の出力指示行の前=上に内容を追加する次の1行を書き込んで「2×2分割表02.htm」などとして保存します:
  出力=出力+"\n\n"+"リスク比="+(b/e)/(d/f)

 10,1,10,10で実行するとリスク比は0.18‥です。次に比率は同じままリスクなし群(aとc)の件数を半分の5と5に書き換えて5,1,5,10で実行するとリスク比は0.25=(1/6)/(10/15)に変わってしまいます。

 横断調査(断面研究)や縦断・追跡調査に加え問題(例:食中毒)の発生後に条件(例:汚染食材)を遡及的に探すことも多い疫学では、事前・事後いずれの群設定※であっても共に使える指標が望まれます。
※設定基準:条件(例:指導法)の差@事前 vs.問題(例:食中毒)の有無@事後

 そんな都合のいい指標がオッズ比(Odds Ratio:OR)です。リスクの計算で分母に行小計(例:20@従来の方法)を用いるのとは異なり、オッズの計算では分母にリスクなしの件数(例:10@従来の方法)を用います(つまり「あり÷なし」)。従って、新しい方法のオッズは10/1=10、従来の方法のそれは10/10=1でオッズ比は10/1=10、つまり新しい方法は従来の方法と比べて掛け率(=見込み)で「10倍有利」なのです。

 では、リスクなし群の件数の変化のオッズ比への影響をscriptで確認してみましょう。リスク比を求める行にオッズ比関連の指示を次行のように追記して同名(2×2分割表02.htmなど)で保存します:
  出力=出力+"\n\n"+"リスク比="+(b/e)/(d/f)+"\n"+"オッズ比="+(a/b)/(c/d)

 まず10,1,10,10で実行するとリスク比=0.18‥、オッズ比=10。次に5,1,5,10にして実行するとリスク比=0.25、オッズ比=10(下の画像参照)。フィッシャー(Ronald Fisher 1890-1962)が交差積比(cross-product ratio)と呼んだように(a/b)/(c/d)=a*d/(b*c)=(d/b)/(c/a)であるオッズ比は、aとcを(0倍以外の)何倍にしても“分子分母打ち消し”の結果「同値」なのです。

 

安曇野での統計のお勉強

 投稿者:唐澤俊英  投稿日:2018年 3月 9日(金)21時45分30秒 mo1-74-25-148.air.mopera.net
  加藤君

統計基礎についての教材を作成していただき感謝です。
とにかく、まずは私が勉強です。

今、安曇野です。
玉ねぎの肥料追加とシイタケ栽培の原木準備作業です。
机の上には、統計や調査法の参考書。
眠い目をこすりながらですね。

パソコンでは、中農研や「農林水産研究に関する国内の論文・情報が探せる
データベース(アグリナレッジ)」   http://agriknowledge.affrc.go.jp/
検索やっています。

推計学の始祖と言われる、フィッシャーは農事試験場で推計学始めたんだものね。
今では、調査・統計は実践の学として、農業ばかりではなく、社会生活を送る上でも避けて通れない。

とにかく、納得できるまで、お勉強します。

心から感謝感謝です。
 

2×2分割表再訪:01 いきなりχ2乗に跳ぶ前に

 投稿者:加藤 厚  投稿日:2018年 3月 9日(金)07時46分51秒 em1-115-194-139.pool.e-mobile.ne.jp
編集済
   ボクたちが学生だった1970年代、周り(教育心理学周辺)では、対応のないカテゴリ変数間の関連の検討はχ2乗検定で“決り”でした※。
※直接確率法=フィッシャーの正確検定は「知って」はいましたが、卒論などで実際に使った人はいたのでしょうか?

 しかし、今振り返ると「期待度数と観察度数の差」に注目するχ2乗※というのは随分と工夫を凝らした考え方の指標です。
※考案者はカール・ピアソン(Karl Pearson 1857-1936)。標準偏差や相関係数の考案者でもあります。名言:「数学は科学の言語、統計学は科学の文法。」

 もしボクが、心理統計の授業を受けずに新たな指導法などを試みて従来の方法を上回る成果を得たら、まずは成果(例:理解や納得)が認められた事例とそうでない事例の割合に注目することでしょう。

 例えば、従来の方法では20人に実施して10対10(五分五分)だった「成果の有無」が、新しい方法を11人に試してみたら10対1(9割以上に成果あり)だったなら、従来の方法と比べて新しい方法は「かなり有望そう」です。

 本連載では、上の例のような「2×2分割表」という“素朴な集計結果”から意味を引き出すいくつかの方法と指標について、その特徴(求めやすさとか分かりやすさとか汎用性とか発展性とか‥)を具体的な「求め方」と共に検討していきます。

 さて、まずは2組の成果の有無の度数a:bとc:d(上の例ならa=10:b=1とc=10:d=10)から下の画像のような小計・総計つきの2×2分割表を作成するscriptが必要でしょう※。
※各群における成果の有無の件数は「集計済」とします。

 下に示した<form>から</script>までの11行をコピーしてメモ帳などに貼り付け、「2×2分割表01.htm」※といったファイル名で保存します。Wクリックで起動したら[実行]をクリックしてみましょう(//以降はコメント)。起動時の例:10,1,10,10(‥有10=a:無1=b@新しい方法 vs.有10=c:無10=d@従来の方法)が集計され、下の画像の結果が表示されるハズです。
※.txtではなく.htmとして保存します(文字コードはUTF-8がオススメ)。
↓この一行(実行許可@IE)はInternet Explorer以外のブラウザでは不要。
<!-- saved from url=(0008)about:internet -->
<form><textarea rows="6" cols="30">10,1,10,10</textarea></form>
<input value="実行" type="button" onClick="処理()">
<script>
function 処理(){
  v=document.forms[0].elements[0].value.split(",")  //文字列を,で分割
  a=v[0]*1;b=v[1]*1;c=v[2]*1;d=v[3]*1  //v[0]~[3]を数値としてa~dに代入
  e=a+b;f=c+d;g=a+c;h=b+d;n=g+h  //行小計e,f、列小計g,h、総計nを計算
  出力=a+"\t"+b+"\t"+e+"\n"+c+"\t"+d+"\t"+f+"\n"+g+"\t"+h+"\t"+n
  document.forms[0].elements[0].value=出力  //値をタブと改行で整形して出力
}  //↑出力指示行
</script>

 <form>~</form>と<input ‥>で6行30列の入出力欄及びクリックでfunction 処理(){}を実行するボタンを定義し、function 処理()の{}内で文字列の分割→数値化(*1=1倍)と変数a~dへの代入→計算→簡易整形(+"\t"と+"\n")及び変数:出力の定義→入出力欄への代入を行っています。例の10,1,10,10を別の値に書き換えて[実行]をクリックし、結果を確認してみましょう。

 

9月からの連載(?)内容に加筆して‥

 投稿者:加藤 厚  投稿日:2017年12月20日(水)07時50分0秒 em119-72-199-101.pool.e-mobile.ne.jp
  「ブラウザに計算をさせてみよう!」を冊子↓※にまとめました。
http://mmua.html.xdomain.jp/kato/pdf/keisan@browser.pdf
※keisan@browser.pdf(2.98MB)

保存の上、SumatraPDF, Kinoppyなどで「全画面見開き閲読」してください。
 

ブラウザに計算をさせてみよう!:15 その5 算術幾何平均法@多倍長小数

 投稿者:加藤 厚  投稿日:2017年12月13日(水)08時00分4秒 em119-72-198-5.pool.e-mobile.ne.jp
   前回のスクリプトを任意精度演算ライブラリbig.js※用に書き直すと以下のようになります:
※12と同様、このスクリプトの作動にはbig.js(MikeMcl, 2017 22KB)が同一フォルダ内に必要です。
<script src="big.js"></script>
<script>
回数=7;試行=0;Big.DP=50
a=Big(1);b=Big(2).sqrt().pow(-1);s=b.pow(2);t=Big(1)
while(試行<回数){
  x=(a.plus(b)).div(2);y=(a.times(b)).sqrt()
  c=(a.minus(b)).div(2);t=t.times(2);s=s.plus(t.times(c.pow(2)))
  π=a.pow(2).times(2).div(Big(1).minus(s))
  document.write(試行+": "+π+"<BR※>")
  a=x;b=y;試行++
}
</script>
※BRは本来は半角(brでもBRでもOK)

 Big.DP=50で計算する小数位(Decimal Places)を50までに指定し※全ての初期値を定数ならBig(n)、平方根なら.sqrt()、べき乗なら.pow(n)で代入します。加算はx.plus(y)、減算はx.minus(y)、乗算はx.times(y)、除算はx.div(y)です。
※指定がない場合はBig.DP=20‥小数第20位まで計算します。

 このスクリプトをコピペして算術幾何平均法@多倍長小数.htmとして保存し、ブラウザにD&Dすると下の画像の結果が示されることでしょう。

 第3~5試行で小数第9~42位まで、そして第6試行では小数第49位まで正しい値が得られています(48位以降の真値は5105820‥)※。語呂合わせの「才子異国に婿さ、子は苦無く身ふさわし」の2倍以上の桁数に達しており、ライブラリの効果が示されています。
※末尾(今回は小数第50位)の値は全計算過程の「丸めの誤差」を含みます。

 さて、今回の紹介はここまで。ブラウザ、スクリプト、アルゴリズム、ライブラリといった素敵な道具を創ってくれた先人たちの肩に乗る※ことで、ボクもかなり遠くまで見晴らすことができました。Alan Key(1940-)がAltoを作った1973年からまもなく半世紀、Steve Jobs(1955-2011)がLisaを作った1983年からでも35年が過ぎました。彼らが目指した「個人の能力を拡張することでその活動を支援する道具」としてのコンピュータとネットの利用をこれからも追求していきたいと思います。
※If I have seen further it is by standing on the shoulders of Giants.(Isaac Newton 1642-1727)


 

先日(12/5)、「SQPの使い方」の改訂※について唐澤君にmailしたところ‥

 投稿者:加藤 厚  投稿日:2017年12月 7日(木)15時05分38秒 em1-115-196-216.pool.e-mobile.ne.jp
編集済
  卒論サポート@J女子大に関する以下の1.~3.の近況紹介(?)がありました。

関心のある人がいるかも、と考え、この BBSで情報提供しておきます。


1.「どういう時にどんな分析=分散・平均値比較など」が有効かを簡単に説明できることを目指しています。

 分析の目的、測定の水準、分布の特徴、変数・群の数と対応の有無などに基づく手法の使い分けの再確認と整理には、次行 URLの(率直?でマジメな)論文:「理学療法領域における統計解析法の選択」(小林 2005)が有益と考えます。
https://www.jstage.jst.go.jp/article/mpta/16/1/16_1_25/_pdf

 他方、「最小限の装備でまずまずの成果を目指す」(←SQP がそうです‥笑)観点からは、尺度は「ほぼ間隔」、分布は「まあ山型」となるように測定法を工夫してパラメトリックな手法に持ち込み、3群以上の比較は(マジメに多重比較したりはせず)平均値の図示で済ませる‥という態度・姿勢もアリ&現実的というのが「私見」です。

 適切な仮説検定は確かに「画竜点睛」ですが、本質的に重要なのは研究自体の内容と達成(←「効果量重視」と共通する姿勢)、「点の打ち方」に拘り過ぎて肝心の竜が蛇尾のようではそれこそ本末転倒でしょう。;-)
---------/---------/---------/---------/---------/---------/---------/

2.「因子分析」「重回帰」など学生が持ってきますが、私には、まだちょっと無理。私自身、「簡単にできる因子分析」などができないかと考えてしまいます。

 「使ってみたい!」という学生の気持ちはワカりますが「仕組み」を概念的に理解していない手法の使用は危険です。

 ボクなら、因子分析については「多くの項目に高く負荷している軸を重要な因子と認めるための前提条件」を、重回帰については「重回帰式から得られる重相関係数Rは何と何の相関なのか」をまず尋ね、正しく答えられないなら使わせません。

 例えば、前者の正答の1つは「項目セットによる対象概念の“均等な網羅”」です。類似した内容を含む項目が多数あればそれが因子になり、かつその因子が多数の項目に高く負荷するのは「当たり前」。逆に、項目セットに含まれていない内容は(それがどれほど重要でも)どんな分析をしようと“出て”はきません。
---------/---------/---------/---------/---------/---------/---------/

3.私は「質問項目」を丁寧に作ること・(先行研究を)読むことにこだわるのですが、学生は時間に追われ「結果」が一番になってしまいます。

 「先行研究」にはそれなりに妥当・適切な測定方法と分析手法が含まれ、また解決すべき問題、達成すべき課題なども言及されているハズです。従って、それらを参考にすることは、自分の研究の意義の確認、並びに妥当な測定と適切な分析の方法の原案確保の最も確実な方法です。それらを踏まえて「可能な進歩・改善」を目指すのが(天才ではない「普通の人」にとっての)研究というものでしょう。2500年前に孔子も言っています:思而不学則殆:思いて学ばざればすなわち殆(あやう)し:Thinking without learning is dangerous.(Confucius)

 他方、概念的理解さえあれば、処理結果は手軽に得られた方がbetter。因子分析を含む多変量解析は以下のfree toolsなどでも可能なようです(‥ただし、ボク自身は使っていません):
 HAD(関西学院大学社会学部の清水准教授作成)
http://norimune.net/had
 College Analysis(福山平成大学経営学部の福井教授作成)
http://www.heisei-u.ac.jp/ba/fukui/analysis.html
 単純な重回帰分析なら Excelでもできますネ:
https://bellcurve.jp/statistics/blog/14015.html
---------/---------/---------/---------/---------/---------/---------/

お気づきの点などあったら書き込みをよろしく!>all

それでは皆さん、お元気で!

※↓頁中ほどの 2011bの右側の青い★が最新版#(←☆のDLはまだ旧版)
http://hp.vector.co.jp/authors/VA021211/
#ver.1.0.8‥recodeやscaleは改訂なし=旧版のまま
 

ブラウザに計算をさせてみよう!:14 その4 算術幾何平均法@標準機能

 投稿者:加藤 厚  投稿日:2017年12月 6日(水)08時00分5秒 em117-55-68-5.emobile.ad.jp
編集済
   円周率を求める方法の1つに算術幾何平均(Arithmetic Geometric Mean)によるアルゴリズムがあります。ガウス(Johann Gauss 1777-1855)の算術幾何平均による完全楕円積分のアルゴリズムとルジャンドル(Adrien-Marie Legendre 1752-1833)の楕円積分に関する関係式に基づくアルゴリズムだそうですが、ボクには何のことやらわかりません。しかし、アルゴリズム(算法=演算手続き)である以上、たとえ理解できなくても使用はできます。しかも「スーパーコンピュータの(略)πの超高精度計算(略)のアルゴリズムがガウス・ルジャンドル法であることが多い」(寒川 2017)とのことで、期待が持てそうです。
http://www.oishi.info.waseda.ac.jp/~samukawa/GaussLeg2.pdf

 以下は、下の画像:上側に示した算術幾何平均によるアルゴリズムの十進BASICプログラム(寒川 2017 p.3)に基づくスクリプトの一例です:
<script>
回数=7;試行=1;a=1;b=1/Math.sqrt(2);s=b*b;t=1
while(試行<回数){
  x=(a+b)/2;y=Math.sqrt(a*b);c=(a-b)/2;t=t*2;s=s+t*c*c
  document.write(試行+": "+((2*a*a)/(1-s))+"<BR※>")
  a=x;b=y;試行++
  }
</script>
※BRは本来は半角(brでもBRでもOK)

 反復計算の回数を 6とし、試行の値が回数=7より小さい間(つまり 6になるまで)の各時点でのπの近似値を試行回数と共に画面表示させます。このスクリプトをコピペして算術幾何平均法@標準機能.htmとして保存し、ブラウザにD&Dすると一瞬で下の画像:下側の結果が示されることでしょう。

 近似値は第4試行(4:)で既に3.14159265‥になり、フーリエ級数で1億項以降まで求めた場合に匹敵する精度が得られています。第5試行で小数第15位まで正しい値が得られる一方、第6試行でも同じ値が示され、11などでも言及した「標準機能で処理できる値」の限界に到達してしまっているようです。ということは、再び任意精度演算ライブラリの出番です。


 

ブラウザに計算をさせてみよう!:13 その3 中学生にも納得できる「面積比と乱数」

 投稿者:加藤 厚  投稿日:2017年11月29日(水)08時00分6秒 em117-55-68-5.emobile.ad.jp
編集済
   フーリエ級数のような「証明には積分が必要で、知らない/忘れた人にはブラックボックス」的な方法ではなく、平均的な中学生にも理解と納得が可能な円周率の求め方に「面積比を乱数で求める方法」※があります。
※「モンテカルロ法」と呼ばれる数値実験(シミュレーション)の一種

 一辺が 1の正方形とそれに内接する円を考えます(下の画像:上参照)。正方形の面積は当然 1、円の面積(=πr2乗)はπ×1/2×1/2=π/4。従って、この円(水色)の面積を何らかの方法で求めれば、面積×4=π。そして、面積は乱数による多数の点で近似が可能です!

 0 以上 1未満の乱数を発生させるMath.random()で正方形中の無作為の位置に一億個の点を打ち、円の中心から点までの距離※が 0.5未満なら円内として数え上げ、途中経過を「試行数:その時点での円周率の近似値」の書式で20回表示するスクリプトの一例は以下の通りです:
※ピタゴラスの定理:直角三角形の斜辺の長さ=√(他の辺1の2乗+他の辺2の2乗)
<script>
回数=100000000;試行=0;円内=0
while(試行<回数){試行++
  x=Math.random()-0.5;y=Math.random()-0.5
  if(Math.sqrt((x*x)+(y*y))<0.5)円内++
  if(試行%(回数/20)==0){
    document.write(試行+": "+4*円内/試行+"<BR※>")
  }
}
</script>
※BRは本来は半角(BRでもbrでもOK)

 コピペして(BRを半角にした上で)円周率@乱数.htmとして保存し、ブラウザにD&Dすると数秒で下の画像:中のような結果※が示されることでしょう。
※乱数に基づくシミュレーションの結果は毎回異なります。

 今回の場合、最初の5百万回までの近似値が3.141508、一億回=最終値も3.14150408で、この程度の試行回数では、点の個数=試行回数が増えてもπの近似値は収束しません。また、途中経過を辿ると、3.141までは一貫して得られているものの、それ以降の桁の値は真の値の3.14159‥を挟んで変動しています。

 今回のスクリプトは、中学生にも理解と納得が可能な手順で「円周率の値は2×2の正方形に内接する円の面積に等しく、その値は3.141~2程度」という内容を示せるものです。しかし、3.141~2という値の精度は、アレクサンドリアのプトレマイオス(下の画像:下参照)が2世紀前半に示した 377/120(≒3.1417)の水準に留まります。無闇に試行回数を増やすことなく、より正確な値を得る方法・工夫は無いものでしょうか?

 

実は「統計」も、ボクたちが最初に学んだ40年前とはずいぶん変わり‥

 投稿者:加藤 厚  投稿日:2017年11月26日(日)22時59分2秒 em1-115-197-234.pool.e-mobile.ne.jp
編集済
  ここ10年ほどは「仮説検定」に加えて「効果量」や「信頼区間」が重視されるようになってきています。

例えば『心理学研究』の「執筆・投稿の手びき」(2015年版)=下のlinkでも「仮説検定はデータ分析の一側面に限られ」るため「研究結果の重要性を評価できるよう効果量とその信頼区間も示す」(p.8) ことが求められたりしています‥その下のlinkは「『教育心理学研究』における p値と効果量による解釈の違い」(2015)。
https://www.psych.or.jp/publication/inst/tebiki20151019_fixed_compress.pdf
https://www.jstage.jst.go.jp/article/jjep/63/2/63_151/_pdf/-char/ja

その意味では SQPも「時代遅れ」なので、現在改訂中の 1.0.8では効果量のいくつか(Cohen's w, η^2, dD)を追加しました(‥年内には 差し替え@Vectorの予定)※。
※1.0.8 はここ↓で本体=sqp.htmのみ開示中(2011bのSQP★:右clickでDLも可)#。
http://hp.vector.co.jp/authors/VA021211/
#I.E.だと桁ズレが発生するので、Chromeなどでの起動がオススメ。

専門的なこと@学会誌はさておき、記述統計と推測統計の基本部分は「自分が実証的に考え、人に客観的≒説得的に伝える『道具』の1つ」として重要なので、「はやりすたり」とは別に自分なりにボチボチやっていきます。

連載中(?)の「ブラウザに計算をさせてみよう!」の最終(15)回でも触れる内容ですが、コンピュータとネットは「個人の能力を拡張することでその活動を支援する道具」、計算の次はクロス表&JavaScript関連、その次はRや Linux関連などであれこれ遊ぶ予定です。;-)
 

卒論ボランティアサポート

 投稿者:唐澤俊英  投稿日:2017年11月26日(日)17時28分0秒 mo220-210-0-45.air.mopera.net
  加藤厚 様

今、加藤君が作成したSQPを使ってデータ分析について学んでいます。

そして、いつもいつも相談にのってもらい感謝しています。

・・・

私は大学時代に統計をしっかり学びませんでした。

「そんな昔のことで困っているの?」「これまで何十年たっているの?」と言われることでしょう。


でも、やはり学ぶ時に学んでおかないと、なかなか難しい。

・・・・

リタイアして、これから何をしたいか! について考えました。

自分の実践活動(子どもや指導者対象の剣道指導、大学生や一般成人対象にした対人関係促進をする指導者養成=実技+調査による評価を含めて、etc.)を行うために調査・データ分析、文章を書く、といったことが必要です。

そこで、大学生の卒論ボランティアサポーターをしながら、いっしょにお勉強することにしました。

身体的にも活動できるのは、あと10年かな。


新宿の風景を見ながら、思いつくまま書いてみました。



 

レンタル掲示板
/3