学問・資格

2017年3月28日 (火)

サンプルサイズの決定(1つの母平均の検定)をカシオの高精度計算サイトkeisan.casio.jpの自作式として作った。

統計の計算式のUPシリーズ第六弾。このサンプルサイズを決めることを一つの目標としていろいろな自作式をUPしてきたのでした。まずは一つの母平均の検定の場合。参考文献は「サンプルサイズの決め方」(永田靖さん)で、この本では近似を使っておられますが、カシオのサイトは高精度に非心t分布が計算できるので厳密に計算しなおしました。
こちらがリンク。
サンプルサイズの決定(1つの母平均の検定)
こんな感じで計算できます。

1つの母平均の検定時に、効果量(Δ=(μ-μ0)/σ 平均の差が標準偏差の何倍か?)と有意水準を与えたとき、必要なサンプルサイズを計算します。 帰無仮説:μ=μ0で、対立仮説としてはμ≠μ0、μ>μ0、μ<μ0の3種類が選べます。

検出力
有意水準
効果量
Δ=(μ-μ0)/σ
対立仮説の種類
μはμ0と等しくない μはμ0より大きい μはμ0より小さい

nが
13
のとき 検出力は
91.07
nが
12
のとき 検出力は
88.29
求めるサンプルサイズ
13
参考文献:
統計勉強シリーズ:

「2群の平均の差のt検定」をカシオの高精度計算サイトkeisan.casio.jpに自作式として作ってみた。

「母比率の信頼区間の推定」計算式をカシオの高精度計算サイトにUP!出口調査で何人投票してたら当選確実になるか?のような。

「1つの平均のt検定」を計算する自作式をカシオの高精度計算サイトkeisan.casio.jpにUP!

「検出力曲線(1つの平均の場合)」を描画する自作式をカシオの高精度計算サイトkeisan.casio.jpにUP!

2017年3月27日 (月)

「検出力曲線(1つの平均の場合)」を描画する自作式をカシオの高精度計算サイトkeisan.casio.jpにUP!

統計勉強シリーズ第五弾。
今回は検出力を計算してみよう。ついでにグラフにする。まずは1つの平均の場合。
こういうのはG*POWERというソフトで計算することが多いそうですが、今回は「サンプルサイズの決め方」に載っている方法を導入(ただし、カシオのサイトが非心t分布を精度よく計算できるので近似なしで)。
リンクはこちら。
検出力曲線(1つの平均の場合)
説明は、

検出力曲線(1つの平均の場合の1-β)を計算します。帰無仮説はμとμ0が等しく、対立仮説はμがμ0と違う場合、μがμ0より大きい場合、小さい場合の3種類が選べます。サンプルサイズと有意水準α、効果量Δ(=(μ-μ0)/σ)、つまり平均の差が標準偏差の何倍か?)の最大・最小値・分割数を入力すると検出力曲線がグラフにできます。

サンプルサイズ
有意水準
対立仮説の種類
μはμ0と異なる μがμ0より大きい μがμ0より小さい
効果量最小値
効果量最大値
分割数
のような感じで、実行した後のグラフはこんな感じで。
μ>μ0のとき
Ken1
μ≠μ0のとき
Ken2
参考文献:
統計勉強シリーズ:

「2群の平均の差のt検定」をカシオの高精度計算サイトkeisan.casio.jpに自作式として作ってみた。

「母比率の信頼区間の推定」計算式をカシオの高精度計算サイトにUP!出口調査で何人投票してたら当選確実になるか?のような。

「1つの平均のt検定」を計算する自作式をカシオの高精度計算サイトkeisan.casio.jpにUP!

「1つの平均のt検定」を計算する自作式をカシオの高精度計算サイトkeisan.casio.jpにUP!

統計を勉強するシリーズ第四弾。
2つの平均の差のt検定をもう作った後なので順番が前後しますが、1つの平均のt検定も作ってみました。
リンクはこちら。
1つの平均のt検定
こんな感じで計算します。データは奥村先生のサイトで紹介されていたものを使っています。
標本平均
標本不偏分散
サンプルサイズ
帰無仮説用の平均μ0
有意水準
%


t値
3.67992
自由度
9
p値
0.0050761
帰無仮説(μ=μ0)
棄却される
信頼区間下限[
0.8977
,
3.7623
]信頼区間上限
参考文献:
統計勉強シリーズ:

「2群の平均の差のt検定」をカシオの高精度計算サイトkeisan.casio.jpに自作式として作ってみた。

「母比率の信頼区間の推定」計算式をカシオの高精度計算サイトにUP!出口調査で何人投票してたら当選確実になるか?のような。

2017年3月22日 (水)

「母比率の信頼区間の推定」計算式をカシオの高精度計算サイトにUP!出口調査で何人投票してたら当選確実になるか?のような。

カシオの高精度計算サイトに統計関係の自作式を勉強のため作ってみるシリーズ第三弾。
今回は母比率の信頼区間の推定です。
リンクはこちら。

母比率の信頼区間の推定

参考にしたのは文献の入門 統計学と青木先生のページ、
で計算方法は、
のWaldとAgresti-Coullを選択できるようにしてます。
画面イメージはこんな感じ:

選挙の出口調査の得票率から合格確実かどうか?を調べるような母比率の信頼区間を推定します。
標本数(出口調査の人数)、条件を満たした数(ある候補者に投票した人数)、信頼区間のパーセント(95%など)と推定方法(WaldとAgresti-Coullより選択)を入力すると信頼区間が計算されます。

標本数
人 or 個
条件を満たした数
人 or 個
信頼区間のパーセント
%
推定方法
Agresti-Coull Wald

母比率
70
%
信頼区間の下限 [
57.6001
% ~
82.3999
% ] 信頼区間の上限
参考文献:
統計勉強シリーズ:

「2群の平均の差のt検定」をカシオの高精度計算サイトkeisan.casio.jpに自作式として作ってみた。

2017年3月20日 (月)

21世紀に初めて見つかった三体問題の解が13個ある。いくつかアニメにしてみた。

久しぶりに数値計算シリーズ。

昔、このScienceの記事と

Physicists Discover a Whopping 13 New Solutions to Three-Body Problem

arxivの元論文のリンク、

Three Classes of Newtonian Three-Body Planar Periodic Orbits

を読んだ。たまたま今日、その話を知り合いとしていたのでアニメでどんな軌跡か見てみよう。

*計算はルンゲクッタ8次のDormand & Prince (DOP853)のExcel VBA移植版を使っています。

http://sci.tea-nifty.com/blog/2011/04/dormand-prince8.html

http://sci.tea-nifty.com/blog/2011/04/dormand-princ-1.html

Twitterに動画を連続投稿してます。

2017年3月17日 (金)

「2群の平均の差のt検定」をカシオの高精度計算サイトkeisan.casio.jpに自作式として作ってみた。

統計、特に推定と検定を再勉強している(理由はまた書きます、と前にも言っていて書いてない、、、)。でその勉強の足跡をカシオの高精度計算サイトに自作式をUPすることで残しておこうと思った。
少し前に、

ABテストに必要なサンプル数の計算式をカシオの高精度計算サイトkeisan.casio.jpに自作式として作ってみた。

というのを作ったが、今回は平均の差のt検定。
リンクはこちら。

2群の平均の差のt検定

こんな感じで入力する。
標本1の個数
標本2の個数
標本平均1
標本平均2
標本不偏分散1
標本不偏分散2
有意水準
%
等分散仮定
等分散を仮定できる 仮定できない(Welch)
と、こんな感じで出てくる。
tの値
-3.11581
自由度
23
p値
0.004862
平均が等しいという帰無仮説
棄却される
信頼区間[
-2.27979
-0.46047
]
Rの出力に合わせてみた。奥村先生のサイトのデータを使わせてもらいました。
Rと数値があってるのでたぶん間違ってないはず、、、Welchの検定もつけてます。
---
参考文献:

2017年2月16日 (木)

FPU(Fermi-Pasta-Ulam)問題、非線形振動子でどうやってエネルギーがモードに分配されるか?をシンプレクティック8次で計算[解説編]

先ほど 「FPU(Fermi-Pasta-Ulam)問題、非線形振動子でどうやってエネルギーがモードに分配されるか?をシンプレクティック8次で計算。エネルギーではなくて、実際の変位をGIFアニメにしてみた。あまり見たことないのでは?詳細はのちほどブログにて。」 というのをTweetした。

ではその詳細を。この手の計算は昔いっぱいしてこのブログでも取り上げていたのですが、さっき

というPainleve先生のTweetを見て再掲しようと思った。Nasaの女性”コンピュータ”を描いた話題の映画Hidden Figuresも作られたことだし。 では過去の記事の再掲:

非線形の再帰現象のさきがけとなったフェルミ・パスタ・ウラム問題。

http://oldweb.ct.infn.it/~rapis/corso-fsc2/FPU-recent-review-Ruffoetal05.pdf

一次元調和振動子を弱く非線形結合(2次や3次)して、最初に1つのモードにエネルギーを与えたら、全部のモードに分配されるかな?と思ってやってみたら最初のモードに戻ってしまったというもの。

例えば2次で結合する以下の方程式で、、、

d2xi/ dt2 = (xi+1 + xi-1 - 2xi) + α[(xi+1 - xi)2- (xi - xi-1)2]

α=1/4、N=64としてシンプレクティック8次でモード毎のエネルギーを計算(もちろんExcel VBAを使用)した結果がこれ。

最初のモードからいろんなモードにはエネルギーは行くものの、いつか戻ってくるという。

Fpunew01

*FermiたちはMANIACというコンピュータでleap frog法で計算したそうです。

で、これのエネルギーでなくて変位をGIFアニメにしたのが最初に書いたTweetでした。

シンプレクティック8次のお話はこちら。

http://sci.tea-nifty.com/blog/2009/04/8excel-vba-7680.html

2017年1月 3日 (火)

誕生日のパラドクスはよく知られてますが、ではN人中K人の誕生日の一致する確率の計算は? #大人のピタゴラスイッチ #ピタゴラスイッチ

昨日見た、大人のピタゴラスイッチで、クラスに40人いたらどれか2人は誕生日が一致する確率が89.1%というのが出ていた。

これを計算するには例えばカシオの高精度計算サイトを参照。

誕生日が一致する確率

案外高いので誕生日のパラドクスとか言われている。

---

では、これを拡張して、、、

N人中、k人の誕生日が一致する確率は?というのは結構難しい。

しかし、非常によく合う近似式が知られている。

”Methods for Studying Coincidences” 、PERSl DlACONlS and FREDERICK MOSTELLER,

Journal of the American Statistical Association
December 1989, Vol. 84, No. 408, Applications & Case Studies

の中の、、、

Birthday_problem2013


という式。

R言語のpbirthday, qbirthdayという関数に使われていたりする。

これ:http://stat.ethz.ch/R-manual/R-patched/library/stats/html/birthday.html

では図示すると、、、

50%を超えるのは

Birthday_problem20132

2人一致:23人、 3人一致:88人、4人一致:187人、5人一致:313人、6人一致:459人とかになる。

本当にこの近似式いいの?という疑問があるので、メルセンヌツイスタでN=1000000でモンテカルロシミュレーションをしたのがこれ。

Birthday_problem20133

なかなかの精度。
--

でこの式をお手軽に計算したいときは、カシオの高精度計算サイトkeisan.casio.jpにアップしている式をどうぞ。

誕生日の一致する確率(N人中k人以上)

誕生日が一致する確率(N人中2~5人、グラフ表示)

6個のサイコロを投げて、どれかの組み合わせで10になる、という確率が98%になるという話、実際にプログラム組んで計算。列挙リスト付。 #大人のピタゴラスイッチ #ピタゴラスイッチ

昨日、2017/1/2に大人のピタゴラスイッチを見ていた。

その中で面白い話が出ていた。

それは、6個のサイコロを投げて、その中で10になる組み合わせができる確率は98%ある、というもの。例えば番組中では2個で10になる例(5,5と4,6)が出ていた。
これ、自分でも調べてみようとExcelのVBAで計算してみた。
あまりにひどいプログラムなのでソースは自粛、、、でも結果だけ示す。
・2個で10になる(残りの4個は問わず(なので4個内に10になる組み合わせがあることもある)
→ 28363通り。
・3個で初めて10になる。
→ 14928通り。
・4個で初めて10になる。
→ 1917通り。
・5個で初めて10になる。
→ 340通り。
・6個で初めて10になる。
→ 126通り。
なので、10になる組み合わせは全部で45674通り。
総組み合わせは6^6=46656通りなので、
確率は、、、97.895%、、、なるほど98%だ!
列挙したリストはこちら。

2017年1月 1日 (日)

#2017年の数学 ペル方程式 x^2 - N*y^2 =1は N=2017にするとすごく桁数の多い解になるが、2015,2016, 2018・・・は桁数が少ない。(PARI/GPで計算)

Pell's equationというのがあって、それはx^2 - N*y^2 =±1の自然数解(x,y)で、その最小解を求めるPari/GPプログラムを以前にかいて、N=10000まで計算していた。

今年2017年はどうなの?ということで調べた。ついでの周りの年代も。
これ:
N,    x,   y

2015, 404, 9 

2016, 449, 10
2017, 22691017898615873418283839489716246568157231499338273 ,505243842362839347335084683756885179819279000763128 
2018, 56280003, 1252834 
2019, 674, 15
2020, 809, 18
今年だけすごく桁数が多くなってる!2020までは少ないな。
次にすごくおおくなるのは?2053年でした。
2053,250607337495961345281030961012407293763134459221296801448080649,	5530944591957664391168044355584108756032547977899443803262020	

より以前の記事一覧

無料ブログはココログ
2017年3月
      1 2 3 4
5 6 7 8 9 10 11
12 13 14 15 16 17 18
19 20 21 22 23 24 25
26 27 28 29 30 31  

最近の記事

最近のコメント

フォト