【R言語】統計解析フリーソフトR 第6章【GNU R】 [無断転載禁止]©2ch.net
■ このスレッドは過去ログ倉庫に格納されています
ww.sankei.com/column/ 180307/clm1803070005-a.html ヒトモドキ産経便所ゴキブリは下着ドロの常習犯変態暴論ゴミ便所自殺しなさい今すぐ 欠損値のないデータで解析したモデルをstepAIC()で処理しようとしたところ、TRUE/FALSEが必要なところが欠損値ですとエラーが出たのですが、どう処理すればよいのでしょうか。MuMInのdredge()でも同様になってしまいました。 https://www3.nhk.or.jp/news/html/20190327/k10011863181000.html インフルエンザの新しい治療薬「ゾフルーザ」を投与されたA香港型のインフルエンザ患者30人を調べたところ、70%余りに当たる22人から、この薬が効きにくい耐性ウイルスが検出されたことが国立感染症研究所の調査で分かりました。 調査件数は多くないものの、専門家は現在のような使用を続けると、耐性ウイルスが広がるおそれがあるとして使用基準を見直すべきだと指摘しています。 耐性化率が50%以上である確率は pbeta(0.5,1+22,1+8,lower=F) [1] 0.9946631でいいかな? > (p=1-binom.test(22,30)$p.v) [1] 0.9838752 > binom.test(22,30,conf=p) Exact binomial test data: 22 and 30 number of successes = 22, number of trials = 30, p-value = 0.01612 alternative hypothesis: true probability of success is not equal to 0.5 98.38752 percent confidence interval: 0.5000000 0.8994036 sample estimates: probability of success 0.7333333 yahoo apiで距離取得するコード掲載したサイトないですかね google料金高過ぎて鞍替えなんですが >>680 使ってるQtが古いんや。1.2でだいぶましになった感じだけど入力途中の文字が残る妙な挙動が時々でる。 >>702 書籍は出版される頃には情報が古くなってるから。自分はネットの情報とヘルプだけで十分だな。 RStudioって32ビット版って存在しないの?? ダウンロードしようとしても32ビット版見当たらないし、Windows7,10用の64ビット版セットアッププログラムを起動しようすると32ビットへの選択とかはなくて、単にエラーになっちゃうし 誰か教えてください httpstwitter.com/shotkr16 低脳中卒万引きヒトモドキネトウヨ猿ヒトモドキをさしころせ ここでRをつかった仕事をされてる人いますか? 学生ですか? 仕事あれば教えてほしい >>709 Rを仕事にしてはいないけど、仕事の一環でR使ってる 信頼性、入手しやすさ、扱いやすさのバランス的に自分の環境ではR一択なんだよな >>711 AI目当てでこれから使おうと思って勉強中 今の職場で求められるのは主に統計用途や作図なんでRは手放せないかな twitter.com/tukuhae ゴキブリネトウヨヒトモドキ奇形売春婦肉便器なつこババア滅多刺しにして解体しろ >>712 なるほど。 自分も機械学習やってるけど、R の方がやりやすいよ。ディープラーニングならPythonかもしれないけど データ分析とか統計解析とかだとRかPythonかだね 機械学習よりは統計寄りだとRは仕事でもスタンダード、特にマニアックな統計解析だとRにしかパッケージがない まあ使い慣れてるならRで機械学習やっても全然構わないけど >>715 少し前はRなんてだんだん廃れるみたいな雰囲気あったけど今はむしろ盛り返してデファクトな何時になってるからこの先10年は行きそうだなあ 機械学習特需でpythonの躍進はすごかったけど Rstudio初心者です。 データファイル名を変えてフォルダには変更名で保存されてるのに、読み込むともとの名前のままなんです。どうすれば変更後の名前になるのでしょう? >>719 全く状況がわからん。 主語や目的語、必要な修飾語を省略せずに、分かるように説明しないと誰も助言できない。 いつのまにかどうでもいい国の資格試験の選択科目にまでなったしな なんなんだろう 資格にしがみついたらいつのまにか技術の進展から取り残された国がありましたとさw 0645 ふうL@Fu_L12345654321 学コン1傑いただきました! とても嬉しいです! https://pbs.twimg.com/media/D-IuUuqVUAALnAB.jpg https://twitter.com/Fu_L12345654321/status/1144528199654633477 https://twitter.com/5chan_nel (5ch newer account) RってMac版とWindows版で優劣ないの? 今から始めるならどっちがオススメ? R自体に優劣はないが文字コード絡みの問題が多いWindows環境はおすゝめしない。 >>729 R本体だとグラフィックディバイスに違いはあるけど、実質的に差はない。 でもRは本体だけで使うことは稀で、多種多様なGUIがあるから、 そこでOSの違いによる差が生じる。 面白い問題スレにあったのでシミュレーションしてみた。 # サイコロ # 正6面体のサイコロがある.4面は青色、2面は赤色である. # このサイコロを合計20回振るとき、最も起こりそうな順番はどれか? # 1.赤 青 赤 赤 赤 # 2.青 赤 青 赤 赤 赤 # 3.青 赤 赤 赤 赤 赤 sim <- function(){ a=sample(0:1,20, replace=TRUE, prob=c(4,2)) b=as.character(a) c=paste(b,collapse="") s1=paste(c(1,0,1,1,1),collapse="") s2=paste(c(0,1,0,1,1,1),collapse="") s3=paste(c(0,1,1,1,1,1),collapse="") res=c(grepl(s1,c),grepl(s2,c), grepl(s3,c)) return(res) } k=1e6 re=replicate(k,sim()) mean(re[1,]) mean(re[2,]) mean(re[3,]) 結果は、直感とおり、1が再頻 > mean(re[1,]) [1] 0.124672 > mean(re[2,]) [1] 0.080873 > mean(re[3,]) [1] 0.040564 > grep使わない方法ってあるかな? >>735 ワイの直感的解法 # 1.赤 青 赤 赤 赤 # 2.青 赤 青 赤 赤 赤 # 3.青 赤 赤 赤 赤 赤 だが、以下でも確率同じ。何となく # 1.赤 赤 赤 赤 青 ★ # 2.赤 赤 赤 赤 青 青 # 3.赤 赤 赤 赤 青 赤 ★は赤でも青でもどっちでもOK P(# 1) >P(# 2) >P(# 3) ∵青出やすい R言語等でシミュレーションされ、 自身の確率直感が正しいのを 確認できるとは、素晴らしい。 上記の件、若干の訂正とする # 1.赤 青 赤 赤 赤 # 1'.赤 赤 赤 赤 青 とすると、 # 1と# 1'は、直感で同じ確率と 思ってたが間違えのようだ。 当方のシミュレーションで、 # 1は、0.1248 # 1'は、0.1271 となった。微妙だけど、多分だ。 やっぱり確率計算をコンピュータで モンテカルロシミュレーションのは 素晴らしい。 >>735 【grepは未使用の糞真面目な方法】 # 1.赤 青 赤 赤 赤 についての確率、ほぼ厳密解を得た # 1は、0.124774… だと思う 計算は、モンテカルロ法でない方法 でプログラム、計算した。 で、grepは使用してない。 ちなみに計算誤差は、ほぼ皆無なハズ ソースコードイメージ p01 = (1/3)^4*(2/3) p02 = p01 p03 = p01 * (1 - p01) p04 = p01 * (1 - p01 - p02) p05 = p01 * (1 - p01 - p02 - p03) … p16 = p01 * (1 - p01 - p02 - p03 … - p14) とし、 p01~p16 の合計を算出したところ、 0.124774… となった >>738 レスありがとうございます。 私の直感 # 1.赤 青 赤 赤 赤 # 2.青 赤 青 赤 赤 赤 # 3.青 赤 赤 赤 赤 赤 # 1.を6個に書き換えて #2.と並べると # 1.★ 赤 青 赤 赤 赤 # 2.青 赤 青 赤 赤 赤 ★は赤でも青でもどっちでもOKだから#1.の方が起こりやすい # 2と# 3を比べると # 2.青 赤 青 赤 赤 赤 # 3.青 赤 赤 赤 赤 赤 3個めでは青の方がでやすいので ら#2.の方が起こりやすい よって、P(# 1) >P(# 2) >P(# 3) 分からない問題スレから、 >> 1回3.6%で激レアが出るガチャを10回回した確率って 36%なのでしょうか? それとも0.964*0.964*0.964(略 0.964を10回電卓にかけた数なのでしょうか? 教えてください。 << 百万回シミュレーション p=3.6/100 N=10 sim <- function() any(rbinom(N,1,p)==1) mean(replicate(1e6,sim())) > mean(replicate(1e6,sim())) [1] 0.306628 シミュレーション その2 gacha=c(rep(1,36),rep(0,1000-36)) sim2 <-function() any(sample(gacha,10,replace=TRUE)==1) mean(replicate(1e6,sim2())) > mean(replicate(1e6,sim2())) [1] 0.306904 何故か0.36にならない、どうしてだろ? シミュレーションその3 (処理速度の関係で10万回の平均) gacha=c(rep(1,36),rep(0,1000-36)) sim3 <- function() any(replicate(10,sample(gacha,1))==1) mean(replicate(1e5,sim3())) > mean(replicate(1e5,sim3())) [1] 0.30691 これも0.307弱だな。 何が悪いんだろ? 俺の頭かな?? 3.6%のガチャを10回回して全部外れる確率は (1 - 0.036) ^ 10 ≒ 0.6930592 したがって、 3.6%のガチャを10回まわして1回以上当たる確率は 1 - 0.6930592 = 0.3069408 >>745 ありがとう、 やっぱ、 これだったw 何が悪いんだろ? 俺の頭かな?? 確率統計なんぞ無意味つまりインチキ 事実、結果こそすべて その確率計算の激ナイーブな解法を示す 激レアは、レアだから1個とみなす。 故に、母集団の個数は 1÷0.036 = 27.777… きっと28個だ 手順1) 28枚のカードがある 手順2) 重複しない1~28の番号を振る 手順3) 28枚のカードをシャッフル 手順4) 1~10枚目のどれか1となる確率 と絶対同じハズ、だから、 P(1枚目で当) =27P27 ÷ 28P28 = 1/28 P(2枚目で当) も同様に1/28 P(3枚目で当) も同様に1/28 … P(10枚目で当) も同様に1/28 で此等10個の事象は背反事象だから、 P = 10/28 = 0.357 ∵有効数字3桁と勝手にしちゃう ところてガチャってゲーム何か よくわかんないけど、計算しちゃった 追記というか突然ですが、 そのガチャ3.6%の件、 超幾何分布なのか。二項分布なのか。 確率の小さいとか、母集団が小さい とかだと無視できないと思われる。 確率統計はギャンブル派生数学ぢゃ。 生半可な知識ではカモにされる。 現代の若者たちは、数学特に 確率統計はじめとするギャンブル の能力が特段に欠けており、 R言語等のプログラミング教育で ギャンブルゲームを学習すべきだ。 奇麗事の学問だけの今日の数学ぢゃ カモにされるだけ。 健全な娯楽として賭博系確率統計学 をC R Java Pyson Javascript BASIC の何れかを学校で学習すべきだ。 私のガチャのイメージは、 かつて、昭和の駄菓子屋によくある ガチャガチャすなわちカプセルトイ そのような健全的なギャンブルが 超幾何分布とか二項分布の理解に 役立つのだ。 限られたお小遣の10円玉数枚で、 如何にレアアイテムをGetするかを 子供らは、思考するからだ。 「残り物には福がある」は 確率統計的には正しいのか 子供同士で文学的に議論したものだ。 さて、今のガチャは恐らくは、 デジタルの媒体のスマホゲームだ。 二項分布でよいだろう。 R言語には、 二項分布の密度関数、それの累積関数 はモチロン、それに従う乱数生成を 提供してるようだ。 ゲームの仕組みが複雑化する今、 R言語等の乱数生成プログラミングで これからのデジタル化ギャンブル社会 で、お金より大切な激レアをドンドン 無限にゲットできる人材の増加を 期待できる可能性を秘めている カプセルトイ(ガチャ)1台に1000個カプセルが入っていて36個がアタリ(レアアイテム)とする。 同じカプセルトイが10台ある。カプセル取り出し後は補充されない。 アタリを1個でも手に入れる確率は 1台から10個取り出す場合(G10)と 1台から1個を10台で取り出す場合(G01) ではどちらが高いか? そのシミュレーション rm(list=ls()) N=1000 ; K=36 ;n=10 # アタリ3.6% g=rep(c(1,0),c(K,N-K)) G10 <- function() any(sample(g,n,replace=FALSE)>0) # 非復元(超幾何分布) G01 <- function() any(sample(g,n,replace=TRUE )>0) # 復元(二項分布) mean(replicate(1e6,G10())) mean(replicate(1e6,G01())) シミュレーションと理論値 > mean(replicate(1e6,G10())) ; 1-choose(N-K,n)/choose(N,n) [1] 0.307745 [1] 0.3081121 > mean(replicate(1e6,G01())) ; 1-(1-K/N)^n [1] 0.307295 [1] 0.3069408 1台から10個取り出した方がいいみたい。 >>754 なるほど、 どの台も全部で、1000カプセルで、 どの台も当りが、36カプセルだと G10(1台で10個)取り出す方が 僅かですが、確率良さそうですね。 シミュレーションも理論値も 同様な結論のようであり、 G10(1台で10個)取り出す方が僅かに 有利が分かり何か楽しかったです。 仮に、もしどの台も 250カプセル中当たり9カプセルなら、さらにG10戦略がG01戦略より有利な 感触を掴めました。 >>755 1000個だとシミュレーションの差が微妙で再現性が不安だったけど 250個にすると、差が明らかにつきますね。 > N=250 ; K=9 ;n=10 # アタリ3.6% > g=rep(c(1,0),c(K,N-K)) > > G10 <- function() any(sample(g,n,replace=FALSE)>0) # 非復元(超幾何分布) > G01 <- function() any(sample(g,n,replace=TRUE )>0) # 復元(二項分布) > > mean(replicate(1e6,G10())) ; 1-choose(N-K,n)/choose(N,n) [1] 0.31082 [1] 0.3117069 > mean(replicate(1e6,G01())) ; 1-(1-K/N)^n [1] 0.306782 [1] 0.3069408 rhyper で超幾何分布の乱数を発生させることができたんだな。 これと rbinomを使うと sample関数を使わなくても同じことができた。 N=250 ; K=9 ;n=10 # アタリ3.6% mean(replicate(1e6,rhyper(1,K,N-K,n)>0)) ; 1-choose(N-K,n)/choose(N,n) # 非復元(超幾何分布) mean(replicate(1e6,any(rbinom(n,1,K/N)>0))) ; 1-(1-K/N)^n # 復元(二項分布) 結果 > mean(replicate(1e6,rhyper(1,K,N-K,n)>0)) ; 1-choose(N-K,n)/choose(N,n) [1] 0.311581 [1] 0.3117069 > mean(replicate(1e6,any(rbinom(n,1,K/N)>0))) ; 1-(1-K/N)^n [1] 0.307213 [1] 0.3069408 1台から順次10個取り出すのをイメージすれば残り物には福があるとも言えなくもないな。 > "フランスの数学者パスカル(1623〜1662)が1654年にフェルマーにあてた手紙が、現在の確率 論の始まりだと言われている。当時の有名な賭博師メレがパスカルに以下のような問題を持ち 込み、その問題についてがその手紙のやりとりの中で論じられているそうである。 甲乙二人がおのおの32ピストル(当時のお金の単位)の金を賭けて勝負したとする。 そしてどちらかが先に3点を得たものを勝ちとし、勝った方がかけ金の総額64ピストルをもら えるとする。ところが甲が2点、乙が1点を得たとき、勝負が中止になってしまった。 このとき、二人のかけ金の総額64ピストルを甲と乙にどのように分配すればよいだろうか。 ただし二人の力は互角で、勝つ確率はそれぞれ1/2ずつだとする。" # 先にw(=3)点を得たものを勝ち、甲がA(=2)点、乙がB(=1)点を得たとき、勝負が中止 # 甲の確率を求める gambling <- function(A=2,B=1,w=3,k=1e5){ # k : how many simulate sim <- function(){ while(A < w & B < w){ g = rbinom(1,1,p=0.5) if(g==1){ A=A+1 }else{ B=B+1 } } A > B } mean(replicate(k,sim())) # Pr[A wins] } > gambling(2,1,3)*64 #64ピストル配分 [1] 48.13056 > gambling(2,0,4) # 日本シリーズ [1] 0.80945 >>758 の計算結果は、興味深い結果だ で、 自己流の直感でラフな計算をしてみた。 甲1点、乙1点 ⇒1:1按分 ∵自明 ★ 甲3点、乙1点 ⇒1:0按分 ∵自明 ☆ 安直に★と☆の平均を考えると 甲2点、乙1点 ⇒1:0.5按分 ∵安直 で、総額で64ピストルだから、 甲に、64*1/1.5 = 42.67 ピストル 乙に、64*0.5/1.5 = 21.33 ピストル を配分すれば、思ってしまった。 ラフな計算だと、少々不公平なことが 起こるようだ。 R言語等のシミュレーションは、 超天才のワィの直感より、 さらに厳密な確率計算ができそうだ。 >>758 最後はAが得点して勝負が決定するのを忘れずにAが勝者になるまでのBの点数で場合分け NS <- function(w,A,B,p){ # 先にw点得点した方が勝者、A,B:現在の点数 ,p:甲の勝率 ans=0 for(k in 0:(w-B-1)){ # k: Aが勝者になるまでのBの点数 ans=ans+choose(w-A-1+k,w-A-1)*(1-p)^k*p^(w-A) } return(ans) } > NS(3,2,1,0.5)*64 #64ピストルの分配 [1] 48 > NS(4,1,0,0.5) #日本シリーズで先勝したほうが優勝する確率 [1] 0.65625 NSは日本シリーズの頭文字から命名、No Skinではないww 賭博等の確率計算をイメージすると ワィの脳内は活性化される だから、 ドンドン、妄想全開となるぅぅぅ。 という訳で、妄想内容を記載する 17世紀の地球に行き、胴元になって 賭博ビジネスでガッポリ儲けるゼェ。 17世紀の地球で、胴元になって、 コイントス賭博事業をするゼェ。 で、システムは、2回コイントスし、 表の回数が0回⇒参加者が32万円GET 表の回数が1回⇒胴元に、32万円PAY 表の回数が2回⇒参加者が32万円GET 参加料は、1万円 まんまと、参加者が沢山 1000人いた ある参加者は、以下の様に怪説をした 組み合わせは、3通りだ。 1通り目 表の回数が0回 2通り目 表の回数が1回 3通り目 表の回数が2回 だから、参加者の勝つ確率は2/3だ で、モチロン、胴元ワィの商売は成功 およそ、1000万円儲かっちゃた。 日本シリーズで賭けをする。 日本シリーズは先に4勝したチームが優勝。 勝率はそれまでの(引き分けを除いた)通算勝率に従うとする。 シリーズ開始前の通算成績はA:2勝、B:4勝であった。 今シリーズでAが先勝(第一試合に勝利)した。 この時点でA,Bのどちらに賭ける方が有利か? 分からない問題スレで3進法の小数とかが話題になっていた https://rio2016.5ch.net/test/read.cgi/math/1567866548/845 >> 845 名前:132人目の素数さん[] 投稿日:2019/12/21(土) 05:52:50.74 ID:wum1jR1j ... 10進数で表された5/9を3進数になおせという問題があって、答は0.12だったのですが、 ... << こんなのを作ってみた。 # 小数点付きの数numをN進法で表示する rm(list=ls()) dec2basen <- function(num, N, kmin = 5){ # kmin:最小小数点後桁 int=floor(num) r=int%%N q=int%/%N while(q > 0){ r=append(q%%N,r) q=q%/%N } # rに整数のN進法表示を格納 x=num-floor(num) k=max(nchar(x)-2,kmin) # 同長もしくはkminの長さの小数表示 a=numeric(k) for(i in 1:k){ y=x*N a[i]=floor(y) x=y-a[i] # r . a[1] a[2] a[3] ... a[k] } if(N<=10){ # Nが10以下は数値として表示 cat(r,paste(".",paste(a,collapse = ''),sep=''),'\n',sep='') } else{ # Nが11以上は整数部分と小数部分を数列で表示 print(list(int=r,decimal=a)) } invisible() } > dec2basen(5/9,3) 0.120000000000000 > dec2basen(3.14159265359,7) 3.06636514320 > dec2basen(3.14159265359,8) 3.11037552421 > dec2basen(3.14159265359,9) 3.12418812407 > dec2basen(3.14159265359,14) $int [1] 3 $decimal [1] 1 13 10 7 5 12 13 10 8 1 4 >>764 バグの原因を追求したら、Rの仕様のせいみたいだな。 > (1.2-1)*5==1 [1] FALSE > (1.2-1)*5>1 [1] FALSE > (1.2-1)*5<1 [1] TRUE これがfloor関数が誤値を返してくる原因だった。 > floor(1) [1] 1 > floor((1.2-1)*5) [1] 0 >>765 数値演算についてもっと勉強しな。 たいていの言語でそうなる。 すると、pythonに移植しても同じ結果ということか。 IPython 6.5.0 -- An enhanced Interactive Python. C:\bin\Anaconda3\lib\site-packages\ipykernel\parentpoller.py:116: UserWarning: Parent poll failed. If the frontend dies, the kernel may be left running. Please let us know about your system (bitness, Python, etc.) at ipython-dev@scipy.org ipython-dev@scipy.org""") (1.2-1)*5==1 Out[1]: False (1.2-1)*5>1 Out[2]: False (1.2-1)*5<1 Out[3]: True (1.2-1)*5 Out[4]: 0.9999999999999998 round を使って回避することにした > 0.728*5-3.64 [1] -4.440892e-16 > round(0.728*5-3.64,5) [1] 0 んで、デバッグ版 # 小数点付きの数numをN進法で表示する(62進法まで対応 0-9,a-z,A-Z) rm(list=ls()) dec62 <- function(num, N, kmin = 5){ # kmin:最小小数点後桁 int=floor(num) r=int%%N q=int%/%N while(q > 0){ r=append(q%%N,r) q=q%/%N } # rに整数のN進法表示数列を格納 k=max(nchar(num)-nchar(floor(num))-1,kmin) # 同長もしくはkminの長さの小数表示 a=numeric(k) x=round(num-floor(num),k) # e.g. 7.28-floor(7.28)-0.28 != 0に対応 for(i in 1:k){ y=round(x*N,k) # e.g. 0.728*5-3.64 !=0 に対応 a[i]=floor(y) x=y-a[i] # r . a[1] a[2] a[3] ... a[k] } b=list(integer=r,decimal=a,num=sum(c(int,a)*(1/N)^(0:k))) fig=c(0:9,letters,LETTERS)[1:N] if(N<=62){ # Nが62以下は数値として表示 cat(paste(fig[b$integer+1],sep='',collapse=''), '.',paste(fig[b$decimal+1],sep='',collapse=''),sep='') cat('\n') } else{ # Nが63以上は整数部分と小数部分を数列で表示 print(b[1:2]) } invisible(b) # b$num:検証用 } >>767 Haskellでもやってみた Prelude> (1.2-1)*5 0.9999999999999998 Prelude> 0.72*5-3.6 -4.440892098500626e-16 >>766 レスありがとうございました。 アドバイスに従ってちょっと勉強してみました。 他の言語に移植しても無駄とわかって助かりました。 こうなる理由が理解できました。 (1+1/10-1)*10==1 [1] FALSE > (1-1+1/10)*10==1 [1] TRUE からnまでの順列を列挙するスクリプト。 これをwhile使って高速化するにはどうすればいいだろう? 俺の頭では思いつかない。 perm <- function (n) { v=1:n sub <- function(n, v) { if (n == 1) matrix(v, 1, 1) else { x = NULL for (i in 1:n) x =rbind(x, cbind(v[i], sub(n -1,v[-i]))) x } } sub(n, v) } > perm(4) [,1] [,2] [,3] [,4] [1,] 1 2 3 4 [2,] 1 2 4 3 [3,] 1 3 2 4 [4,] 1 3 4 2 [5,] 1 4 2 3 [6,] 1 4 3 2 [7,] 2 1 3 4 [8,] 2 1 4 3 [9,] 2 3 1 4 [10,] 2 3 4 1 [11,] 2 4 1 3 [12,] 2 4 3 1 [13,] 3 1 2 4 [14,] 3 1 4 2 [15,] 3 2 1 4 [16,] 3 2 4 1 [17,] 3 4 1 2 [18,] 3 4 2 1 [19,] 4 1 2 3 [20,] 4 1 3 2 [21,] 4 2 1 3 [22,] 4 2 3 1 [23,] 4 3 1 2 [24,] 4 3 2 1 > >>771 自己解決 e1071というパッケージに再帰呼び出しなしのpermutationがありました。 library(e1071) permu <- function (n) { z <- matrix(1) for (i in 2:n) { x <- cbind(z, i) a <- c(1:i, 1:(i - 1)) z <- matrix(0, ncol = ncol(x), nrow = i * nrow(x)) z[1:nrow(x), ] <- x for (j in 2:i - 1) { z[j * nrow(x) + 1:nrow(x), ] <- x[, a[1:i + j]] } } return(z) } 可変引数の扱い Cで #include <stdio.h> int main(int argc, char* argv[]) { while(argc--) printf("%s\n", *argv++); return 0; } をRでやるには ... を使って、c(...)もしくはlist(...)で引き出せばいいんだな。 検索してもなかなかわからなかった。 main <- function(...){ argv=c(...) for(i in 1:length(argv)) cat(argv[i],'\n') } >>546 library(Rmpfr) mpfr(factorial(50),5000) one = mpfr(1, 5000) factorial(one*50) > mpfr(factorial(50),5000) 1 'mpfr' number of precision 5000 bits [1] 30414093201713018969967457666435945132957882063457991132016803840 > one = mpfr(1, 5000) > factorial(one*50) 1 'mpfr' number of precision 5000 bits [1] 30414093201713378043612608166064768844377641568960512000000000000 何故か、oneを高精度に設定しておくと、正しい計算をしてくれた。 RのroundのIEEE仕様 > round(2.5) [1] 2 > round(1.5) [1] 2 は丸め誤差が減るという説明だったので実験してみた。 # 四捨五入 vs IEEE with FUN(=mean) for x(=c(-0.9,-0.8,...,0.8,0.9)) comp <- function(n=10,FUN=mean,x=sample((-9:9)/10,n,rep=T),print=T){ X=FUN(x) dif = abs(X-FUN(round(x))) - abs(X-FUN(f45(x))) # round後の実行と四捨五入後の実行の差 if(dif!=0 & print) cat(dif<0,' : ',sort(x),'\n') # 差があれば表示 round後が小さければTRUE return(dif) # 差を返す dif<0:round優位 dif>0:四捨五入優位 } comp() k=1e5 # mean re=replicate(k,comp(print=F)) c(IEEE=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0)) # squared sum comp(FUN=function(x) sum(x^2)) re=replicate(k,comp(FUN=function(x)sqrt(sum(x^2)),print=F)) c(IEEE=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0)) > # 平均 mean > re=replicate(k,comp(print=F)) > c(IEEE=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0)) IEEE 四捨五入 引き分け 0.25235 0.16199 0.58566 > # 平方和 squared sum > re=replicate(k,comp(FUN=function(x)sqrt(sum(x^2)),print=F)) > c(IEEE=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0)) IEEE 四捨五入 引き分け 0.40519 0.01413 0.58068 平均だと小差だが、統計計算の内部処理で使う平方和だと大差がついたので、体感できた。 "隔離中のクルーズ船では船内の換気が共通らしいから運悪く13日後に発症した奴がいるとその近くの部屋のやつは プラスで14日隔離しないといけない。それが今の船内の状況という。 両隣のどちらかが感染したら14日延長、どの部屋も1日で感染する確率pは1% 部屋の配置は長方形でどの部屋にも両隣があるとする。 発症するか、隔離期間が終われば下船できる。 全員定員1の個室として客と乗務員を合わせた人数nは3000人。 クルーズ船から全員下船できる日数の期待値は?" rm(list=ls()) p=0.001 # 1日の感染確率 n=3000 # 収容者人数 d14=14 # 隔離日数 Q0=rep(d14,n) # 必要隔離日数の配列の初期値 # ODA <- function(X){ # One Day After 必要隔離日数の配列 x から1日後の配列を返す Q=X # Xを作業用Qに代入 S=which(Q!=0) # 未感染者susceptibleのindex s=length(S) # 未感染者人数 m=rbinom(1,s,p) # その日の感染数 if(m==0){ # 感染者0なら未感染者の Q[S]=Q[S]-1 # 必要隔離日数を1日減らして返す invisible(Q) }else{ I.idx=sample(S,m) # 未感染者Sからm人が感染、感染者のindex (Infected.idx) id2ext <- function(id){ # 感染者idからの両隣のidを返す plus =ifelse(id==n,1,id+1) # 最大番号なら1を返す,それ以外は1増やす minus=ifelse(id==1,n,id-1) # 1なら最大番号を返す,それ以外は1減らす unique(c(plus,minus)) # 重複を除く } E.idx=as.vector(sapply(I.idx,id2ext)) # 感染者の両隣のidを列挙し E.idx=unique(E.idx) # 重複を除く Extended quarantine.idx Q[E.idx]=d14 # 両隣は必要隔離日数をd14にリセット Q[I.idx]=0 # 感染者は必要隔離日数=0 (リセット者に重複していたら0で上書き) IE.idx=unique(c(I.idx,E.idx)) # 感染もしくは感染者の隣のindex minus1ifnotzero <- function(x) ifelse(x==0,0,x-1) # ゼロでなければ1を減じて返す Q[-IE.idx]=minus1ifnotzero(Q[-IE.idx]) # 感染もなく感染者の隣りでもない人の隔離日数を1日減らす invisible(Q) } } sim <- function(){ AZ=FALSE # All Zero 隔離日数がすべて0か? i=0 # カウンター Q=Q0 # 隔離日数配列初期値 while(!AZ){ # 隔離日数がすべて0になるまで Q=ODA(Q) # One Day After関数を実行 AZ=all(Q==0) # All zero? i=i+1 # カウンター } return(i) # 何日かかったかを返す } k=1e4 # k回のシミュレーションの RE=replicate(k,sim()) mean(RE) # 平均値=期待値 >>778 一日の感染確率。 3000人いて1日3人くらい報告されたから0.1%/日くらいかな 実行したら全員が下船できるまで約1ヶ月になった。 コードにまだバグがあるかも。 Rのcsvの文字化け治らねえ csvそのものをutf-8にしたら今度は不正なマルチバイトがあるってエラー どうすればいいか教えて >>781 どんなデータ? ヘッダーだけじゃなくてセル内も日本語? >>781 エスパーじゃないから、まず確認だけど ・Rを使ってる環境(OS)とRのバージョンは? ・CSVのコード変換に使ってるソフトは? ・CSVを読み込んでる関数と指定している引数は? nCov2019 for studying COVID-19 coronavirus outbreak https://guangchuangyu.github.io/nCov2019/ 後で遊んでみよう >>784 エラーがでた。R 3.6.3にアップデートしたけど > library('remotes') > remotes::install_github("GuangchuangYu/nCov2019", dependencies = TRUE) Downloading GitHub repo GuangchuangYu/nCov2019@master Skipping 1 packages not available: BiocStyle Installing 5 packages: downloader, cowplot, maps, magick, BiocStyle Error: Failed to install 'nCov2019' from GitHub: (converted from warning) package ‘BiocStyle’ is not available (for R version 3.6.3) We also developed a web app (http://www.bcloud.org/e/ ) with interactive plots and simple time-series forecasts. という記述があるんだけど 関数 nls による非線形回帰分析 で作成しているのかなぁ? 面白い問題スレをRでシラミ潰しにカウントしてみた。 "縦n個、横n個のマス目のそれぞれに 1,2,3,...,n の数字を入れていく。 このマス目の横の並びを行といい、縦の並びを列という。 どの行にも、どの列にも、2つの対角線上にも同じ数字が1回しか現れない入れ方は何通りあるか求めよ。(2020京大文系 改) n=4がオリジナル " library(gtools) n=4 (a=permutations(n,n)) # nの順列 r=nrow(a) f<-function(x){ # x=c(1,2,3,4) -> rbind(a[1],a[2],a[3],a[4]) ans=NULL for(i in x){ ans=rbind(ans,a[i,]) } return(ans) } b=permutations(r,n) head(b) ; tail(b) B=apply(b,1,f) diag2 <- function(x){ # 行列xの右上から左下への対角線の配列を返す n=nrow(x) ans=numeric(n) for(i in 1:n) ans[i]=x[i,n+1-i] return(ans) } g <- function(v){ # matri(v,n)で列、対角線の要素がすべて異なればTRUEを返す n=sqrt(length(v)) (x=matrix(v,n)) if(length(unique(diag (x))) < n) return(FALSE) if(length(unique(diag2(x))) < n) return(FALSE) flag=TRUE for(i in 1:n){ if(length(unique(x[,i])) < n){ flag=FALSE break } } return(flag) } counter=NULL for(i in 1:ncol(B)){ if(g(B[,i])) counter=append(counter,i) } length(counter) counter matrix(B[,counter[1]],n) matrix(B[,counter[2]],n) matrix(B[,counter[48]],n) >>789 n=5 にするとメモリー不足でエラーがでた。 サンプリングでやってみる # simulation n=5 sim <- function(n=5){ v=as.vector(t(replicate(n,sample(n)))) # n × n行列要素をベクトル化 if(g(v,dia=TRUE)){ # diagonal latin squareなら print(matrix(v,n)) # 行列化して表示 return(1) } return(0) } k=1e6 sum(replicate(k,sim(5))) 1億回モンテカルロしてようやく、1個みつかった。 [,1] [,2] [,3] [,4] [,5] [1,] 2 5 4 3 1 [2,] 4 3 1 2 5 [3,] 1 2 5 4 3 [4,] 5 4 3 1 2 [5,] 3 1 2 5 4 嫌ね最近は総当たりアルゴリズムで枝刈り使えば実用的な時間で解決できるこ知らない人が多くて versionが3.6,xになって、str2lang とか str2expressionという関数が追加されたんだな。 https://id.fnshr.info/2019/07/20/r3-6-0/ で知った 例 rep("1:6",5) paste(rep("1:6",5),collapse=',') str2lang(paste("expand.grid(",paste(rep("1:6",5),collapse=','),')')) eval(str2lang(paste("expand.grid(",paste(rep("1:6",5),collapse=','),')'))) (s=rep("1:6",5)) (str=paste(s,collapse=',')) (lang=str2lang(paste("expand.grid(",str,')'))) eval(lang) The Clinical and Chest CT Features Associated with Severe and Critical COVID-19 Pneumonia https://journals.lww.com/investigativeradiology/Abstract/publishahead/The_Clinical_and_Chest_CT_Features_Associated_with.98832.aspx の付表にこんなデータがあった # Parameter Total (n=83) Severe/critica l Group (n=25) Ordinary Group (n=58) P Value # Number of lobes Median (interquartile range, IQR) # 5 (4-5) 5 (5-5) 5 (3-5) 0.003 重症群25例 中等症群58例でCTで肺炎像が確認できた肺葉数(1〜5)が 全体で中央値が5 第一四分位数(少ない方から数えて25%の値)が4、第三四分位数(少ない方から数えて75%の値)が5 重症群では中央値5、第一四分位数5、第三四分位数5 中等症群では中央値5、第一四分位数3、第三四分位数5で p値が0.003であるという。 0.0025〜0.0035までとして、 この条件に見合う83例をリストアップすると > re[100,][1:25] [1] 1 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 > re[100,][26:83] [1] 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 [48] 5 5 5 5 5 5 5 5 5 5 5 という類の組み合わせで 19504通りと計算された。 重症群の肺葉数の平均は4.76 中等症群の平均は4.04 と計算された。 タイ値がこれだけ多いデータにWilcox検定を使う意味がわからんな。 ■ このスレッドは過去ログ倉庫に格納されています
read.cgi ver 07.5.5 2024/06/08 Walang Kapalit ★ | Donguri System Team 5ちゃんねる