「統計学関連なんでもあり」の過去ログ--- 042

No.09035　確率の求め方　　【nao】　2009/01/26(Mon) 15:43

下記のようにRで二項数を発生させるときに使う確率の求め方について教えてください。

>rbinom( n, size, prob)

現在，「風速ごとにエサを獲得する確率」をグラフにしたいと考えています(下記データを参照ください)。このときの確率probは，下記の獲得率（現象が起きた数÷標本数）ではないようなのですが，それではどのように求めればいいのかよくわかりません。ご存知の方，教えていただけないでしょうか？すみませんが，よろしくお願いします。

使用するデータ
↓
風速(km/h) 　獲得率(%)　
10         　60%
20         　80%
30         　60%
40         　40%
50         　80%
60         　40%
         
架空の元データ
獲得できた：1，できなかった：0（n=5)
     a b c d e 
10km 1 0 1 0 1
20km 1 1 0 1 1
30km 0 1 1 0 1 
40km 1 0 1 1 0
50km 0 1 1 1 1 
60km 1 1 0 0 0

No.09036　Re: 確率の求め方　　【青木繁伸】　2009/01/26(Mon) 16:39

例えば獲得率60％，n=5 のとき，1,0,1,0,1 のように，1が60％現れるような0/1データを作りたいということですか？
そしてたぶん，そのときには rbinom(5, 1, 0.6) の様にしたのでしょうけど，そのときに出来る5個の0/1データの中に 1 が60％現れないことがあるということですか？

そうだとすれば，「それは，二項*乱*数*なので，正確に 60％になるような乱数列ではないからですよ」というのが答えです。
> set.seed(12345)
> rbinom(5, 1, 0.6)
[1] 0 0 0 0 1　　　1は 20% しかない
> x <- rbinom(1000000, 1, 0.6) # たくさん発生させれば，
> table(x)
x
     0      1 
400405 599595 　　ほぼ 4:6 に近づく
そういうことを聞いているのではない？

No.09037　Re: 確率の求め方　　【nao】　2009/01/26(Mon) 17:55

青木先生ご返信ありがとうございます。

＞例えば獲得率60％，n=5 のとき，1,0,1,0,1 のように，1が60％現れるような0/1データを作りたいということですか？

作りたいのは，「風速が10kmのとき，獲得率が00％」というものです。
虫のエサ獲得率が風速で異なる，さらに足の長さの違うグループ毎でも変わるかどうかを知るために調べたいと考えています。

参考用に画像を添付します（うまく添付できるかわかりませんが・・・）。
このような図を作りたいと考えています。
（後書き：やっぱりエラーが出て添付できませんでした）。

わかりずらくてすみませんでした。

No.09038　Re: 確率の求め方　　【青木繁伸】　2009/01/26(Mon) 18:02

添付画像については，上の「利用上の注意」に書いてあります。
* 添付可能ファイル : PNG
* 最大投稿データ量 : 150 KB
* 画像は横400ピクセル，縦400ピクセルを超えると縮小表示されます。

例えば「風速が10kmのとき，獲得率が60％」って，要するに，n=5 で，1が3個，0が2個ランダムに並んでいるデータということですか？そういうことなら，上に書いたように rbinom ではダメですよ。1 1 1 0 0 を作って，ランダムに並べ替えればよいわけです。

No.09041　Re: 確率の求め方　　【nao】　2009/01/26(Mon) 19:42

青木先生コメントありがとうございます。

私はR初心者で，言葉の使い方なども含めて中途半端な理解しかできていないため，混乱させてしまって申し訳ありません。

青木先生にお教えいただいたやり方は，全くなランダム確率であるのではないかと想像しています。そして私が作りたいのは，自分のデータから推定される確率ということです。
これは異なるように思いますが，どうなのでしょうか？
データをプロットしてみると，風速が小さいときは大きいときに比べてエサ獲得が多いようです。それを示すために，その差について調べたかったのですが，その過程でrbinom()を使うということを知りました。

最初の説明がわかりずらく，お手数をおかけしました。
どうぞよろしくお願いします。

No.09042　Re: 確率の求め方　　【nao】　2009/01/26(Mon) 19:55

参考までに画像の添付です。

No.09043　Re: 確率の求め方　　【青木繁伸】　2009/01/26(Mon) 20:39

> そして私が作りたいのは，自分のデータから推定される確率ということです。

もう少し説明していただかないと，わかりません。

No.09046　Re: 確率の求め方　　【nao】　2009/01/26(Mon) 21:37

お手数をおかけしてすみません。だんだん焦ってきました。
見当違いのことをしているのかもしれません。

以下はこれまで進めてきた行程です。長くなりますがご容赦ください。

説明変数: x
・風速(km/h)

応答変数:nget
・エサを得たか得ていないか（量は無視)

fruit <- glm (cbind(nget, ntotal-nget)~ x, data=data, family=binomial(link=logit))

このfruitモデルから予測される図を描きたいと思っています。
解析の参考にした文献には，

> We analised the probablity of the receiving a fruit using the generalized linear models with binomial error distributions.
とあり，結果を作図したものは，先ほど添付したものです。
この図がわかりやすかったので，自分でもこういう図を描きたいのですが・・・。
この図の軸ラベルに，「probability of・・・」とあったので，これは説明変数（私の場合は風速）ごとにエサ獲得できた確率だと考えました。
しかし，両軸ラベルの範囲は0～1になっているので，単に割り算したわけではないようなので，他にも確率を求める方法があるのだろうと考えました。
調べてみると，ブートストラップというやり方で，乱数を発生させるという方法がありましたので，この方法を使っているのではないかと考え，やり方ネットや文献ではよくわかりませんでしたので，質問させていただきました。

正しいかどうかはわからないのですが，このように考えた次第です。
ですから，一番やりたいことは，glmで得られた結果の予測図ということになるかと思います。

すみません。つたない文章ですが，伝わりましたでしょうか？

No.09048　Re: 確率の求め方　　【青木繁伸】　2009/01/26(Mon) 22:13

glm を使って，ロジットモデルで分析するということなんですか。rbinomなんて，無関係ではないですか。

風速ごとに，餌を取れた個体数 nget，取れなかった個体数 nmiss のデータがあればよいんでしょう。nget+nmiss = ntotal

> d
   x nget ntotal
1 10    2     23
2 20    7     32
3 30   17     45
4 40   12     26
5 50   21     37
6 60   17     19
> ans <- glm(cbind(nget, ntotal-nget)~x, data=d, family=binomial)
> summary(ans)

Call:
glm(formula = cbind(nget, ntotal - nget) ~ x, family = binomial, 
    data = d)

Deviance Residuals: 
      1        2        3        4        5        6  
-0.4559   0.2101   0.6429  -0.3076  -1.0732   1.2585  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.69389    0.47055  -5.725 1.03e-08 ***
x            0.06652    0.01206   5.513 3.52e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 41.1686  on 5  degrees of freedom
Residual deviance:  3.4955  on 4  degrees of freedom
AIC: 28.074

Number of Fisher Scoring iterations: 4
> plot(d$x, fitted.values(ans), ylim=c(0,1))
> points(d$x, d$nget/d$ntotal, pch=19, col="red")

白丸が推定値。赤丸が実測値。

そういうことですか？

No.09049　Re: 確率の求め方　　【nao】　2009/01/26(Mon) 22:38

お付き合いいただき，ありがとうございます。
わかりやすく説明するつもりが意味不明の連発で，申し訳ありませんでした。

こんな風に推定値がプロットできるんですね。
勉強になりました。
「GLMの推定結果をパラメトリック・ブートストラップ尤度比検定にかけた」という説明がねっと上にあったので，そのようなことをしなければダメなのかと思っていました。

最後にもう1つだけお伺いしたいのですが，
rbinom(n, size, prob)の時の，probはどのように算出するのでしょうか？
前に青木先生のHP上で，解説を読んだような記憶があったのですが，見つかりません。
もしも記憶違いであればお許しください。

そして，不勉強により自分でも何をやりたいのかよくわからない段階で，質問などをしてすみませんでした。

No.09050　Re: 確率の求め方　　【青木繁伸】　2009/01/26(Mon) 23:01

> rbinom(n, size, prob)の時の，probはどのように算出するのでしょうか？

算出するものではないですよ。
rbinom 関数は，母比率 prob，試行数 size の二項乱数を n 個作るものです。
例えば，10円玉を5回投げて表の出る回数を記録するような実験をシミュレートするのに使える訳です。10円玉の表の出る確率が prob で，まともな10円玉なら，prob=0.5 でしょう（つまり，裏も表も同じ確率で出るということ）。5回投げるというのが試行回数 size。実験を10回繰り返すということに対応するのが n。結果は 0~5 の数値が10個得られる。
> rbinom(n=10, size=5, p=0.5)
 [1] 5 2 1 3 3 2 2 4 2 2
つまり，母比率は算出するものではなく，仮定する（決める）ものです。

問題：サイコロを10回投げて1または2の目が出る回数を記録する。そのような実験を25回繰り返したとき，どのような結果になるかシミュレーションするにはどうすればよいか。
解答：rbinom(n=25, size=10, prob=2/6)

問題：10円玉を30回投げて，表が出るか裏が出るかを記録するシミュレーションはどのようにするか。
解答：rbinom(n=30, size=1, prob=0.5)

No.09051　Re: 確率の求め方　　【nao】　2009/01/26(Mon) 23:13

青木先生へ

このたびは最後まで丁寧な説明をいただいてありがとうございます。
おかげで，たいへんよくわかりスッキリしました。

先生にはお手数をおかけけしてしまいましたが，とても勉強になりありがたかったです。
これからさらに勉強していきたいと思います。
お世話になりました。