【R言語】統計解析フリーソフトＲ 第6章【GNU R】 [無断転載禁止]©2ch.net
■ このスレッドは過去ログ倉庫に格納されています
>>278
>argsは見れる関数と見れない関数があるのであんまり役に立たない
だったら、例が悪いよ。見られない関数の例をあげないと
> ヘルプで見ようと思ったら書いてなかった
こっちだってわざわざ下請け関数のオプションまでは書かないよ。
ヘルプを見て未解決ならソースを見てくれ。 rm(list=ls())
で R のメモリーをクリア 先週長々と質問させていただいたものです。その後、試行錯誤して行った結果十解決しました。
以前は「summary」としか入れていなかったのですが、「dd$summary」と入力したところ
全体のデータが出力されました。それ以後「dd$summary$dir.area」等で出せるようになりました。
大変初期的な場面で躓いておりましたが、現在何とか進んでおります。
色々とお答えくださった方々、とても参考になりました。本当にありがとうございました。 N(=100)回コインをなげてｎ(=5回)以上続けて表がでる確率。
seqn<-function(n=5,N=100,p=0.5){
rn=rbinom(N,1,p)
count=0
for(i in 1:N){
if(rn[i] & count<n){
count=count+1
}
else{
if(count==n) {return(TRUE)}
else{
count=0
}
}
}
return(count==n)
}
mean(replicate(10^5,seqn()))
> mean(replicate(10^5,seqn()))
[1] 0.81085
案外、高い確率になった。 pooledVariance <- function(...) {
args = list(...)
n.args=length(args)
ss2=0
df=0
for(i in 1:n.args){
ss2 = ss2 + var(args[[i]])*(length(args[[i]])-1)
df = df + (length(args[[i]])-1)
}
ss2/df
}
effectsize <- function(y1,y2){
diff=mean(y1)-mean(y2)
var=(var(x1)*(length(x1)-1)+ var(x2)*(length(x2)-1))/(length(c(y1,y2))-2)
sd=sqrt(var)
diff/sd
}
library(effsize)
cohen.d() >>303
単なる二項分布。
コインが５回続けて表がでたら、0.5^5 <0.05なのでイカサマコインといわれちゃいそうなんだが、
１００回やってみると案外、５回表が続くので確率を計算しようと思ったが、解析的にできる頭がないので
シミュレーションしてみた。 >302は
> rbinom(100,1,0.5)
[1] 0 1 0 0 1 0 0 0 1 1 1 0 1 0 1 1 0 1 1 0 1 1 0 0 0 1 0 1 1 1 1 0
[33] 0 1 0 0 0 0 0 0 1 0 1 1 0 0 0 1 0 0 1 0 1 1 1 1 1 1 1 1 1 1 1 0
[65] 0 0 1 1 0 0 1 0 1 0 0 0 0 0 0 1 0 0 1 1 1 1 1 0 1 1 0 0 1 1 0 0
[97] 0 0 1 1
で５回以上１が連続するときTRUEを返す関数なのだが
もっと簡単にやれないかなぁとは思っている。
rep(1,5) %in% rbinom(100,1,0.5)は １個ずつ評価されてTRUEが５個返ってくるだけ。
文字列にしてgrepを使うとなんとかなりそうな気がしないでもないのだけど.... すまん、教えてほしいだけど
分析するために初めてＲをインストールしようと思って、このスレのあるように公式サイト行ったら、esetが「JS/Redirector.NAV トロイの木馬」を検知したんだが…；
どうしたらいいだ… ホームページダメとかわけわかんねぇ・・・とんでもねぇな
ありがとう >>309
CRANはここ
https://cran.r-project.org
プロジェクトの方はESETの誤検知っぽいんだよな >>310
してねぇよｗｗ
totalvirusで調べると一件引っかかるし、Redirector検知だから、しょうがないね col=rgb(runif(1),runif(1),runif(1),runif(1))で色指定すると
走らせるたびに色がちがっておもしろい。
hist(rnorm(100),col=rgb(runif(1),runif(1),runif(1),runif(1))) stanやJAGSのコードでgamma関数を使おうとして
y = gamma(x)
と、書いたらエラーになった。
stanだと y=tgamma(x)、JAGSだとy=exp(loggam(x))で動作した。 ある大学の入学者男女の比率は１であるという帰無仮説を検定する課題が花子と太郎に課された。
花子は５０人を調査できたら終了として入学者を５０人をみつけて１８人が女子であるという結果を得た。
帰無仮説のもとで
５０人中１８人が女子である確率は 0.01603475
これ以下になるのは５０人中０〜１８人と３２〜５０人が女子の場合なので
両側検定して
> sum(dbinom(c(0:18,32:50),50,0.5))
[1] 0.06490865
> binom.test(18,50,0.5)$p.value
[1] 0.06490865
で帰無仮説は棄却できないと結論した。
http://i.imgur.com/XDIp9rM.png
一方、本番と十八番が好きな太郎は一人ずつ調べて１８人めの女子がみつかったところで調査を終えることにした。
１８人めがみつかったのは花子と同じく５０人めであった。
帰無仮説のもとで
１８人がみつかるのが５０人めである確率は0.005772512
これ以下になるのは２３人以下５０人以上番めで女子１８人めがみつかった場合なので
両側検定して
pnb=dnbinom(0:999,18,0.5)
> 1 - sum(pnb[-which(pnb<=dnbinom(50-18,18,0.5))]) # < 0.05
[1] 0.02750309
http://i.imgur.com/K3T7utr.png
で帰無仮説は棄却される。
どちらの検定が正しいか、どちらも正しくないか？
検定する意図によってｐ値が変わるのは頻度主義統計の欠陥といえるか？
花子の横軸は裏口入学者数、太郎の横軸はサンプル数なので
サンプルでの裏口入学率を横軸にして９５％信頼区間を示す。
花子の検定での信頼区間は0.36〜0.72で18/50を含む、p=0.06491
http://i.imgur.com/SeTLk8K.jpg
太郎の検定での信頼区間は0.375〜0.72で18/50を含まない、p= 0.0275
http://i.imgur.com/tNzlfxe.jpg
主観である、検定の中止の基準の差でｐ値や信頼区間が変化するのは変だという批判である。 （訂正）
ある大学の入学者男女の比率は１であるという帰無仮説を検定する課題が花子と太郎に課された。
花子は５０人を調査できたら終了として入学者を５０人をみつけて１８人が女子であるという結果を得た。
帰無仮説のもとで
５０人中１８人が女子である確率は 0.01603475
これ以下になるのは５０人中０〜１８人と３２〜５０人が女子の場合なので
両側検定して
> sum(dbinom(c(0:18,32:50),50,0.5))
[1] 0.06490865
> binom.test(18,50,0.5)$p.value
[1] 0.06490865
で帰無仮説は棄却できないと結論した。
http://i.imgur.com/XDIp9rM.png
一方、十八という数字が好きな太郎は一人ずつ調べて１８人めの女子がみつかったところで調査を終えることにした。
１８人めがみつかったのは花子と同じく５０人めであった。
帰無仮説のもとで
１８人がみつかるのが５０人めである確率は0.005772512
これ以下になるのは２３人以下５０人以上番めで女子１８人めがみつかった場合なので
両側検定して
pnb=dnbinom(0:999,18,0.5)
> 1 - sum(pnb[-which(pnb<=dnbinom(50-18,18,0.5))]) # < 0.05
[1] 0.02750309
http://i.imgur.com/K3T7utr.png
で帰無仮説は棄却される。
どちらの検定が正しいか、どちらも正しくないか？
検定する意図によってｐ値が変わるのは頻度主義統計の欠陥といえるか？
花子の横軸は女子数、太郎の横軸はサンプル数なので
サンプルでの女子の割合を横軸にして９５％信頼区間を示す。
花子の検定での信頼区間は0.36〜0.72で18/50を含む、p=0.06491
http://i.imgur.com/SeTLk8K.jpg
太郎の検定での信頼区間は0.375〜0.72で18/50を含まない、p= 0.0275
http://i.imgur.com/tNzlfxe.jpg
主観である、検定の中止の基準の差でｐ値や信頼区間が変化するのは変だという批判である。 kainokousiki<-function(a,b,c){return (-b+sqrt(b^2-4*a*c))/(2*a)} #解の公式
kainokousiki(1,-5,6)
でrunすると3じゃなくて6を返すんだけど、どこが間違ってる？ 機械学習をきっかけにPythonに逆転された感じだね 『Rを使った〜』だとPythonじゃないのかよって思うよね ずっとRだけでPython触ったこと無いけど、覚え直す価値ある？
環境構築からもう面倒なイメージ Python自体は:と直後のインデントさえ気を付ければかなり簡単
3系は数が全て小数扱いなので楽
Anacondaというパッケージでインストールすれば、今流行りのJupyter Notebookという開発環境で対話的にコーディングできる（Rも使える） アナコンダてのがRぽくできるのね、ありがとう
dplyrやggplot2みたいに素人でも簡単便利だといいんだけど
Pythonはオブジェクト志向ぽいしすぐ諦めそう >>335
誤爆？
何をブラウズする環境が無くなったの？ RStudio Severあるで。Dockerで使うのがいいんじゃないかな？！ 統計学とウェブ解析を交えて実践的な勉強と練習を
したいのですが、おすすめな書籍やサイトはありますか。
実際に解析ツールや分析ツールを用いて
自分で分析解析してから
解答を見て適切な手順や方法、考察を
解説してくれるものが良いです。
統計学は統計検定2級の知識はありますが
ウェブ解析はテキスト読んだだけです。 R studioは日本語コメント書く度にIMEが無効になったりカーソルがずれたり黒文字の予測変換が黒背景と重なって見えなくなったりと散々だわ 会社規定なんでしかたない
UNIX環境使えるのが羨ましい Docker for WindowsでRStudioサーバー動かせば？ >>342
IME無効になるのは俺だけじゃなかったと知ってほっとした。 ストアアプリも同じ症状でるからRStudio固有の問題でなくWindows環境の不治の病だと思ってる Run any R code you like. There are over three thousand R packages preloaded.
https://rdrr.io/snippets/ Fisher test検定時に
p<2.2e-16
と表示されるんですが、これより小さい値の指数桁数を正確に表記する方法教えて下さい。
例えば5.8e-35となるようにです。 fisher.test関数の返り値はlist型で、その中にp.valueという名前でp値が格納されているから$演算子を使って直接参照するか、broom::tidy関数に返り値を渡してdata.frame形式で出力してやれば見れる >>352
>>354が言う直接的な参照
> fisher.test(matrix(c(1,120,130,2),2))$p.value
[1] 1.691912e-69 352です。よく分かりました！
ありがとうございます！ >>356
技術的な助言をしたけど、学術的に言えば、
p < 0.01 は全て p < 0.01 として、具体的なp値を考える意味はないと思うよ。
一部の例外的な研究分野を除いて(e.g., 遺伝統計学)。 はい、まさにその例外的な分野で使おうとしてます。ありがとうございます。 Rのガンマ関数はいくつでオーバーフローするかやってみた。
> i=1
> while(gamma(i)!=Inf){
+ i=i+1
+ }
Warning message:
In gamma(i) : value out of range in 'gammafn'
> i
[1] 172
> gamma(172)
[1] Inf
Warning message:
value out of range in 'gammafn'
> gamma(171)
[1] 7.257416e+306 matplot()で折れ線グラフ描いたときに、X軸をカテゴリで示したいのですが、
可能でしょうか？
例えばtemp <- c("0時間","8時間","24時間","48時間")として、
matplot()の引数にtempをとるやり方です。
他にもやり方あれば教えてください。 >>371
matplot(..., xaxt="n")
axis(1, at=seq(along=temp), lab=temp) >>372
遅くなりましたがありがとうございました。
できました！ 特定の長方形の中に複数の長方形を最小面積で敷き詰める平面充填に関するパッケージってありませんかね # jonckheereテストを書いてみた
jonckheere <- function(L,
alternative = c("two.sided", "increasing", "decreasing"),
cat=TRUE){
# L : list of vectors A1,A2,...,Ak, with assumed tendency
How.Many.Greater.Pairs <- function(A,B){ # How many pairs of A[i] > B[j], count as 0.5 when equal,
n.a = length(A)
n.b = length(B)
how.many.greater.pairs = 0
for(i in 1:n.a){
for(j in 1:n.b){
how.many.greater.pairs = how.many.greater.pairs+ifelse(A[i]==B[j],0.5,A[i]>B[j])
}
}
return(how.many.greater.pairs)
}
Sum.of.Greater.Pairs <- function(L){ #L=list(A1,,,,Ak),A1 < A2 < A3,..,< Ak : vector
k = length(L)
comb = combn(1:k,2) # possible combinaition of pairs to compare
n.comb = ncol(comb) # how many combinations
J = 0 # sum of greater pairs
for(i in 1:n.comb){
J = J + How.Many.Greater.Pairs(L[[comb[1,i]]],L[[comb[2,i]]])
}
return(J)
}
J = Sum.of.Greater.Pairs(L)
n = sapply(L,length)
N = sum(n)
EJ = (N^2-sum(n^2))/4
VJ = (N^2*(2*N+3)-sum(n^2*(2*n+3)))/72
Z = (J-EJ)/sqrt(VJ)
alternative = match.arg(alternative)
p.value = switch(alternative, 'two.sided' = 2 * min(pnorm(Z), pnorm(-Z), 0.5),
'increasing' = pnorm(Z),
'decreasing' = pnorm(-Z))
if(cat){
cat( 'p.value = ', p.value,'\n')
cat('alternative hypothesis: ' ,alternative,'\n')
}
invisible(p.value)
} ■ このスレッドは過去ログ倉庫に格納されています
