【R言語】統計解析フリーソフトＲ 第6章【GNU R】 [無断転載禁止]©2ch.net
■ このスレッドは過去ログ倉庫に格納されています
# 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)
} >>388
実態調査か何か？
<-と=は挙動が違う場合があるので、使い分けていますが、
代入はどっちかと問われたら、無論 <- または ->
なお、
> 1 -> x
これはエラーにならないけど、
> 1 = x
1 = x でエラー: 代入の左辺が不正 (do_set) です
これはエラー >>389
俺は基本=派。
関数の定義は
z.test <- function(x,n=16,sigma=1){
z=sqrt(n)*mean(x)/sigma
2*pnorm(abs(z),lower=FALSE)
}
と書いている。 ＞1 = x でエラー
当たり前
そんな使い方なんてするかよ
他言語と同じく=一文字の方がすっきりしてイイ こういうのが紛らわしいから、俺は = 推奨。
x <- 1
if(x <- 1) print('YES')
if(x < -1) print('YES') やったことなかったので関数の初期値設定に<-を使うとどうなるかやってみた。
まず、= の場合
> z.test <- function(x,n=16,sigma=1){
+ z=sqrt(n)*mean(x)/sigma
+ 2*pnorm(abs(z),lower=FALSE)
+ }
> z.test(1:3)
[1] 1.244192e-15
<- で初期値設定すると、エラー
> z0.test <- function(x,n<-16,sigma<-1){
Error: unexpected assignment in "z0.test <- function(x,n<-"
> z=sqrt(n)*mean(x)/sigma
Error in mean(x) : object 'x' not found
> 2*pnorm(abs(z),lower=FALSE)
Error in pnorm(abs(z), lower = FALSE) : object 'z' not found
> } 俺は見栄えがいいと思って関数定義には<-を使っているけど
= でも通常に動作する。
> z.test = function(x,n=16,sigma=1){
+ z=sqrt(n)*mean(x)/sigma
+ 2*pnorm(abs(z),lower=FALSE)
+ }
> z.test(1:3)
[1] 1.244192e-15
<- 推奨の人に聞きたいのだけど
<- でないと動作しないってことあるのだろうか？ RStudioつかってりゃ［Alt+-］で簡単入力
それと引数の指定は代入じゃないと思うんだが感覚が違うよかな？ >>396
essだとアンダースコアを入れるを勝手に「<-」になる。
>>394
それは代入ではなく、関数の規定値の設定。
規定値の設定は「=」と決められているので、「<-」はアウト。
ただし、関数を実行するときには「<-」を使うことができる。
> mean(x<-1:10)
[1] 5.5
> x
[1] 1 2 3 4 5 6 7 8 9 10 =に拘り続けて、<-に変える瞬間の君たちがたまらないよ >>398
ありがとうございます。
そんな使い方ができたのですね。
こんなのしか知りませんでした。
> f <- function() {
+ x<<-1:10
+ mean (x)
+ }
> f()
[1] 5.5
> x
[1] 1 2 3 4 5 6 7 8 9 10 一画面しかグラフがないのに
Hit <Return> to see next plot:
と出るの鬱陶しいな、と思っていたら
par(ask=FALSE)
と設定しておけばいいんだな。 同名のdataがあるのでパッケージを::で指定してもうまくいかなかった。
data(netmeta::parkinson)
でなくて
data(parkinson, package = 'netmeta')
とするんだな。 hist(rnorm(100),col=rgb(runif(1),runif(1),runif(1),runif(1)))
だと、1色だけど
hist(rnorm(100),col=sample(colours(),(sample(1:10,1))))
なら1〜10色で表示される。 > hist(rnorm(100),col=rgb(runif(1),runif(1),runif(1),runif(1)))
何色に増やしてもいいがな。
何をやりたいのかな？
hist(rnorm(5000),col=apply(matrix(runif(80),4), 2, function(x){rgb(x[1],x[2],x[3],x[4])})) :::てどういう時に使うんだろ？
ソース読みたくて
library(BayesFactor)
ttestBF_indepSample
では表示されなかったが、
library(BayesFactor)
BayesFactor:::ttestBF_indepSample
だと出てきた。 >>410
グダグダ言っていないで、さっさとヘルプを参照すればいいじゃないか？
?':::' >>411
ありがとう。
stats:::t.test.default
でｔ検定のソースが見られた 関数を値で上書きしてしまった場合に元に戻す方法を教えて頂きたく。
例)
rm <- 1
としてしまった場合にrm関数を元に戻したいです。
よろしくお願いします。 >>413
直後ならRstudioでctrl+zで元に戻す。
再起動してたらこの方法は無理だろな >>413
自分で作ったのではない関数の rm だったら rm(rm) で数値ベクトルオブジェクトが消えて、関数の rm が見えるようになる。 >>413
> rm <- 1
> ls()
[1] "rm"
> rm(rm)
> ls()
character(0)
> rm
function (..., list = character(), pos = -1, envir = as.environment(pos),
[以下略]
深く考えなくてもrm(rm) で元に戻るけど。
rmが数値ではなくて自作関数だったら、少々ややこしいけど、
base::rm()で大丈夫だろう >>413
既に説明されてるけど上書きじゃなくて別オブジェクト扱いになるからrmで消すか明示的にパッケージを::演算子で関数を指定してやればおけ
自作関数だと上書きされちゃうけど rmが自作関数だと戻せないのでは？
> rm = function(x) asin(sqrt(x))
> rm(1)
[1] 1.570796
> rm
function(x) asin(sqrt(x))
> rm = 1
> rm
[1] 1
> rm(rm)
> rm
function (..., list = character(), pos = -1, envir = as.environment(pos),
inherits = FALSE) 同じ名前空間にあるなら上書きされちゃうよ
.Gloalenv（だったかな？）に自作関数があり、そこに上書きして変数にしてるから消したら消えるだけで戻せない ちょっと古いけど、この辺りを読むとなぜそうなるかが分かると思う
（Rの）環境問題について その１。
ttps://qiita.com/kohske/items/325bdf48f4f4885a86f1 >>418
戻せます。
> rm <- function(x) asin(sqrt(x))
> rm
function(x) asin(sqrt(x))
> base::rm(rm)
> rm
function (..., list = character(), pos = -1, envir = as.environment(pos),
[以下略] >>419
>.Gloalenv
case sensitiveなRでそれはないのでは？
> help(".GlobalEnv") >>414-422
みなさまご返信ありがとうございました。
例示したものがbase::rmだったので自作関数のことまでは想像だにしませんでした…(がとても良い勉強になりました。)
自環境でも上書きした関数を元に戻すことが出来ることを確認しました。 >>421
戻したいのが自作関数
rm = function(x) asin(sqrt(x))
だったら上書きされて無理じゃないの？
ゴミ箱から消去したファイルを復活させるソフトもあるから
PCに詳しければ可能かもしれないけれど
Rの機能としては復活させるのは無理と思う。 自作のrm関数はbaseパッケージとは別に定義されるから
base::rm(rm)で自作の方のrmは消えてbase::rmが「復活」します。 ああすまん
自作の方を上書きした場合は、そっちの復活は無理ですね。 >>425
ちょっと意味がわからない。
自作の関数なら定義を再実行するだけなのに、
消える消えないって何が問題なんだ？ >>424
関数だけでなく既定値でも同じ。
> pi
[1] 3.141593
> pi=2
> pi
[1] 2
> rm(pi)
> pi
[1] 3.14159 >>306
知恵袋に漸化式があったので計算スクリプトを書いてみた。
表がでる確率ｐのコインをN回投げてK回以上表が続く確率。
# N: flips
# K: least sequential head
# p: probability of head
seqNp <- function(N=100,K=5,p=0.5){
if(N==K) return(p^K)
q=1-p
a=numeric(N) # a(n)=P0(n)/p^n , P0(n)=a(n)*p^n
for(i in 1:K) a[i]=q/p^i # P0(i)=q for any i
for(i in K:(N-1)){ # recursive formula
a[i+1]=0 # a[i+1]=q/p*(a[i]+a[i-1]+a[i-2]+...+a[i-(K-1)])
for(j in 0:(K-1)){
a[i+1]=(a[i+1]+a[i-j])
}
a[i+1]=q/p*a[i+1]
}
P0=numeric(N) # P0[n] : probability of ending with 0 head when flipped n times
for(i in 1:N) P0[i]=a[i]*p^i # P0(n)=a(n)*p^n
MP=matrix(rep(NA,N*K),ncol=K)
colnames(MP)=paste0('P',0:(K-1))
MP[,'P0']=P0
MP[1,'P1']=p
for(i in (K-2):K) MP[1,i]=0
# MP[k,n] = Pk[n] : probability of ending with k head when flipped n times
for(k in 2:K){ # Pk(n+1)=p*P(k-1)(n)
for(i in 1:(N-1)) MP[i+1,k]=p*MP[i,k-1]
}
ret=1-apply(MP,1,sum)
ret[N]
} >>306
この例で、答えとして何が返ってきて欲しいのか？
2で合ってる？ >>433
1が5回以上出現しているからTRUE（数字なら１）が返り値。 >>434
run.check <- function (x, n=5) {
# x の中に 1 が n 回以上連続していれば TRUE を返す.
chg <- c(TRUE, diff(x) != 0) # 変化があった場所
chgidx <- c(which(chg), length(x)+1) # 変化があった場所の添え字
run.length <- diff(chgidx) # 0や1の連続している個数
true.length <- run.length[x[chg] == 1] # 1の連続している個数
any(true.length >= n) # 連続している個数が n 以上のrunがあるか?
} >>435
whichとdiffを用いての解法のお返事ありがとうございます。
実行時間を>302のスクリプトと比べてみました。
run.check <- function (x, n=5) {
# x の中に 1 が n 回以上連続していれば TRUE を返す.
chg <- c(TRUE, diff(x) != 0) # 変化があった場所
chgidx <- c(which(chg), length(x)+1) # 変化があった場所の添え字
run.length <- diff(chgidx) # 0や1の連続している個数
true.length <- run.length[x[chg] == 1] # 1の連続している個数
any(true.length >= n) # 連続している個数が n 以上のrunがあるか?
}
# 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)
}
> system.time(mean(replicate(10^5,run.check(rbinom(100,1,0.5)))))
user system elapsed
17.56 0.04 17.68
> system.time(mean(replicate(10^5,seqn())))
user system elapsed
9.74 0.07 9.78
後者の方が速いのは1がn個連続したら、TRUEを返して終了するので以後のチェックはしないためだろうと思います。
diff や anyの使い方が勉強になりました。
時間を割いてスクリプトを作成していただいてありがとうございます。<(_ _)> Rstudioで日本語の入力ができないのですが、解決する方法知ってる方いませんか?
具体的には、IME自体をオンにすることができません。
日本語の表示やコピペは問題なくできます。
環境
Linux - Utubntu
RStudio Version 1.0.143 RStudioのQtが古いのでパッチ当てないと日本語入力できないのが現状ですわ
パッチはUbuntu 16.04.LTS用にしかないので他のバージョンならRStudioサーバー使った方が早いかと
なお、パッチはRStudio 日本語とかでグクってみて windows版のRStudioでもたまに日本語入力ができなくなる。
別のソフトでIMEのON/OFFをやってから戻ると直るけど。 Windowsはストアアプリでも同じ症状でるから持病なんじゃないかと思う
Ubuntuはfcitxのパッチ当てとけばWindowsのような症状は出ないので快適 これをやってみた。
http://img-cdn.jg.jugem.jp/d75/1405116/20140609_1105745.jpg
スクリプトは
http://imagizer.imageshack.com/img921/8020/t1VcqA.jpg
小数点以下500桁でみると19個素数があった。
> vip=Vectorize(is.prime)
> (idx=which (vip(N)))
[1] 100 124 150 172 183 202 215 219 255 296 310 314 323 346 400 417 439 441 474
> cat(substring (e,idx[1],idx[1]+9))
7427466391 問題
９９人の囚人がいます。彼らの頭に１〜１００までのナンバーカードが貼りつけられた帽子をランダムにかぶせます。
他人の帽子は見ることができても、自分の帽子は見ることができません。
帽子の数は全部で１００なので、一つ使われずに余ります。
そのナンバーは囚人達にはわからないようにしておきます。
この状況で、囚人たちに一斉に自分のナンバーを宣言させて、全員が正解だったら釈放するという賭けをします。
囚人たちには帽子をかぶせられる前に相談タイムが設けられています。
どういう戦略を取れば、助かる確率を最も高くできるでしょうか？
http://000013.blogspot.com/2010/12/99.html
解答を読んでも数理が理解できないが
確率が0.5になるのはシミュレーションできた。
# http://000013.blogspot.com/2010/12/99.html
inversion <- function(x){
n=length(x)
ret=numeric(n)
for(i in 1:(n-1)){
ret[i] = sum(x[i] > x[(i+1):n])
}
sum(ret) # inversion number
}
is.even= function(x) !inversion(x)%%2 # is inverion number even?
prisoner99 <- function(n=100){
indx=sample(1:n,1) # defective number
X=sample((1:n)[-indx])
Y=numeric(n-1)
for (i in 1:(n-1)){ # select as even permutation
x1=X[-i]
x2=(1:n)[!(1:n) %in% x1] # two numbers unseen for i-th prisoner
tmp=X
tmp[i]=x2[1] ; tmp[n]=x2[2]
Y[i]=ifelse(is.even(tmp), x2[1],x2[2])
}
all(X==Y)
}
mean(replicate(1e3,prisoner99())) >>437
遅レスだが、Ubuntu 18.04LTSならデフォルトのIBus-mozcで日本語入力ができるそうだ。
詳しくはUbuntu Weekly Recipeの第527回を読んでくだされ 覆面算 RED + WHITE = COLOR
どうもＣでやったらうまくいかないのでＲでやってみた。
# RED + WHITE = COLOR
x=c('R','E','D','W','H','I','T','E','C','O','L','O','R')
unique(x)
redwhite <- function(R,E,D,W,H,I,T,C,O,L){
red=10^(2:0)
white=10^(4:0)
color=10^(4:0)
sum(red*c(R,E,D))+sum(white*c(W,H,I,T,E)) - sum(color*c(C,O,L,O,R))
}
x=unique(x) ; x
REDWHITECOLOR <- function(x){
R=x[1]
E=x[2]
D=x[3]
W=x[4]
H=x[5]
I=x[6]
T=x[7]
C=x[8]
O=x[9]
L=x[10]
cat(paste('RED = ',R,E,D,
' WHITE = ',W,H,I,T,E,
' COLOR = ',C,O,L,O,R),'\n')
}
REDWHITECOLOR(unique(x))
library(gtools)
perm=permutations(n=10,r=10,v=0:9)
perm=perm[perm[,1]!=0&perm[,4]!=0&perm[,8]!=0,] # R!=0,W!=0,C!=0
head(perm) ; tail(perm)
n=nrow(perm)
re=numeric(n)
for(i in 1:n){
re[i]=redwhite(perm[i,1],perm[i,2],perm[i,3],perm[i,4],perm[i,5],perm[i,6],perm[i,7],perm[i,8],perm[i,9],perm[i,10])
}
hist(re)
(indx=which(re==0))
REDWHITECOLOR(perm[indx,])
> REDWHITECOLOR(perm[indx,])
RED = 5 8 7 WHITE = 3 9 6 1 8 COLOR = 4 0 2 0 5 # AAB
# × CD
# ------
# CCE
# BFAA
# ------
# BEBDE
f <- function(A,B,C,D,E,F){
(A*100+A*10+B)*D==C*100+C*10+E &(A*100+A*10+B)*C==B*1000+F*100+A*10+A &
(C*100+C*10+E)+(B*1000+F*100+A*10+A)*10==B*10000+E*1000+B*100+D*10+E
}
library(gtools)
perm=permutations(n=10,r=6,v=0:9)
n=nrow(perm)
re=numeric(n)
for(i in 1:n){
re[i]=f(perm[i,1],perm[i,2],perm[i,3],perm[i,4],perm[i,5],perm[i,6])
}
indx=which(re==TRUE)
perm[indx,] 大変初歩的な質問で恐縮ですが、ご教示いただければ幸いです。
PERT分布を学んでいるのですが、ベータ分布のスケーリングの仕方が分かりません。
最小値a、最頻値ｂ、最大値ｃ として、a, b, c から求めたパラメータをa1、a2 とします。
このとき、
dbeta(x, a1, a2)*(c-a)+a
とすると、当然ながら横軸の値の範囲は0から1で、縦軸の値のみ大きくなってしまいます。
横軸の値の範囲がaからcで、縦軸の値は確率密度のままにするには、どうしたらよいのでしょうか？ え？aだけ右に動かしたいならこれでいいのでは？
dbeta(x-a, a1, a2) >>448
447です。ご教示どうもありがとうございました。
a <- 40 # 最小値
b <- 45 # 最頻値
c <- 70 # 最大値
mu <- (a+4*b+c)/6 # PERT分布の平均値
a1 <- b*(mu-a)/(c-a) # パラメータ a1
a2 <- b*(c-mu)/(c-a) # パラメータ a2
pert1 <- function(x) dbeta(x,a1,a2) # ベータ分布布の計算
まではできたのですが、ご教示いただいた方法で
pert1 <- function(x) dbeta(x-a ,a1,a2)
とすると、ｘ の範囲が0から1で、ｙが0の一様分布のようなグラフになってしまいます。 4パラべーたへの変換は下記かな？
https://en.wikipedia.org/wiki/Beta_distribution#Four_parameters_2
a1とa2のパラメータ推定って合ってるんだっけ？僕も正直わかんない
やりたいことはこれだと思う
install.packages("mc2d")してからやってね
----------------------------------------------------
rm(list=ls())
library(mc2d)
a <- 40 # 最小値
b <- 45 # 最頻値
c <- 70 # 最大値
pert4<-function(x) dpert(x,min=a,mode=b,max=c,shape=4)
x1<-40:70
plot(x1,pert4(x1)) >>450
447です。
450様、まさにこれ、これです。どうもありがとうございました。
おかげさまで大変助かりました。
448様にも、ご助言下さいましたこｔ、厚く御礼申し上げます。 lmtestのパッケージをインストールしようとすると不正なマルチバイト文字がありますと出てくるのですがどうすればいいですか？ >>452
Windows 10 （Rの日本語ヘルプはインストールしていない）+ R studioだけど、そのメッセージ出ないな。 職場のSPSS信者たちにRを認めてもらうにはどうしたら良いだろう
「安かろう悪かろう」「俺が知らないものは三流」的な考え方の人が多くて肩身が狭い… >>455
>>452のようなエラーが出たら「無料のものはやっぱりダメね」
有償(SPSS) →信頼できる
ボランティア→いい加減
と言う思考らしい >>454
ユーザー名はアルファベット1文字です
studioインストールしてないんですけど、しないとダメですか？ linuxなら環境変数LANGをCにすれば解決する問題だからwindowsでも似たような対処すれば行けそう。 これかな？
> Sys.setlocale(locale="C") # locale の変更（USの設定にしないと表示が変になる）
[1] "C"
元に戻すには
> Sys.setlocale(locale="Japanese_Japan.932") すみません、教えて下さい。
フィッシャーの正確確率検定で、
適合度（一様性）の検定ってできますか？
できるのは独立性の検定のみ？
Rでいろいろ試したみたんですが、
2群ないとエラーになる。。 >>462
レス感謝です。
そうなんですが、
期待度数が5未満になってしまった場合は、
カイ二乗検定は向かないんですよね？ >>463
こういうので
r1=5;r2=4;n1=10;n2=12
prop.test(c(r1,r2),c(n1,n2),correct=TRUE)
chisq.test(matrix(c(r1,r2,n1-r1,n2-r2),2,byrow=TRUE),correct=TRUE)
警告がでるという話かな？
Warning message:
In chisq.test(matrix(c(r1, r2, n1 - r1, n2 - r2), 2, byrow = TRUE), :
Chi-squared approximation may be incorrect >>463
適合度を見るのに関係ないだろ。
> cnt
[1] 1 0 2 4 3 4 12 6 12 11 21 22 21 37 40 30 59 44 49 47 38 48 36 43 38
[26] 46 33 34 30 21 26 33 17 13 14 12 8 16 7 7 9 6 6 5 3 8 4 4 4 1
[51] 0 2 1 0 0 0 0 2
> prb
[1] 0.0002947255 0.0006464052 0.0012538873 0.0022037831 0.0035720843
[6] 0.0054114364 0.0077413658 0.0105430123 0.0137588816 0.0172972083
[11] 0.0210398698 0.0248524558 0.0285950600 0.0321325300 0.0353432214
[16] 0.0381256526 0.0404028045 0.0421240952 0.0432652773 0.0438266419
[21] 0.0438299769 0.0433147339 0.0423338191 0.0409493611 0.0392287274
[26] 0.0372409836 0.0350539129 0.0327316483 0.0303329196 0.0279098757
[31] 0.0255074215 0.0231629878 0.0209066521 0.0187615256 0.0167443293
[36] 0.0148660914 0.0131329065 0.0115467110 0.0101060386 0.0088067297
[41] 0.0076425765 0.0066058952 0.0056880193 0.0048797146 0.0041715201
[46] 0.0035540195 0.0030180499 0.0025548579 0.0021562071 0.0018144483
[51] 0.0015225575 0.0012741483 0.0010634658 0.0008853653 0.0007352804
[56] 0.0006091856 0.0005035529 0.0004153084
> chisq.test(cnt, p=prb, rescale.p=TRUE, simulate.p.value=TRUE)
Chi-squared test for given probabilities with simulated p-value (based
on 2000 replicates)
data: cnt
X-squared = 65.504, df = NA, p-value = 0.2189
5未満の数があっても関係がない >>466
警告が出るのは独立性検定の場合でしたね。 463です。
カイ二乗検定で、期待度が5未満が多いのは不適切なのかと思い込んでおりましたが、
それは独立性の検定の場合のことであり、
適合度（一様性）の検定の場合は、
気にしなくて良かったのか。。
ならカイ二乗検定でやります。
いろいろ教えて下さり、ありがとうございます。 >>469
いまSPSSを使う理由って何なの？煽りじゃなくマジで >>472
周囲のspssユーザを見ると「指導教員がspssしか使えないので、
自分もspssしか使えない、
新しい（学術系）ソフトは誰かに丁寧に根気よく指導してみらえないと無理」、
これの再生産。 再帰呼び出しに関数名の代わりにRecall使うと処理に要する時間が増えることに気づいた。
フィボナッチ数列の再帰呼び出し
# 関数名での呼び出し
fibo <- function(n){
if(n==1|n==2) return(1)
else fibo(n-1)+fibo(n-2)
}
# Recallを使っての呼び出し
fiboR <- function(n){
if(n==1|n==2) return(1)
else Recall(n-1)+Recall(n-2)
}
> system.time(fibo(30))
user system elapsed
4.27 0.02 4.36
> system.time(fiboR(30))
user system elapsed
6.65 0.00 6.70 cumsumを再帰関数で書いてみた。
何度も試行錯誤
# cumsum with recursive call
cumsumR <- function(v){
res=numeric(length(v))
cumsumR_sub <- function(v,res,i){
res[1]=v[1]
if(i==length(v)) return(res)
else{
res[i+1] = res[i] + v[i+1]
Recall(v,res,i+1)
}
}
cumsumR_sub(v,res,1)
}
cumsumR(1:10)
とりえあず、動作した。
> cumsumR(1:10)
[1] 1 3 6 10 15 21 28 36 45 55 簡略化できた
cumsumR <- function(v,res=NULL,i=1){
res[1]=v[1]
if(i==length(v)) return(res)
else{
res[i+1] = res[i] + v[i+1]
Recall(v,res,i+1)
}
}
> cumsumR(1:10)
[1] 1 3 6 10 15 21 28 36 45 55
> ifelseの動作がどうも納得できないでご教示いただきたいのですが。
これってifelseの仕様でしょうか？バグでしょうか？
> x=2:1
> if(x[1]<=x[2]) x else x[2:1]
[1] 1 2
> ifelse(x[1]<=x[2],x,x[2:1])
[1] 1 >>482
ifelse() はベクトル演算できる関数。
だから、戻り値の要素数は、test部分の要素数に一致する。
testが1つなので、x[2:1]の1つ目の要素が返される。 testの要素を2つにすれば、出力も2つになる。
> ifelse(c(x[1] <= x[2], x[1] <= x[2]), x, x[2:1])
[1] 1 2
さらに３つにすると、helpに書いている通り、yesやno部分は再利用される。
> ifelse(c(x[1] <= x[2], x[1] <= x[2], x[1] <= x[2]), x, x[2:1])
[1] 1 2 1 >>483
ありがとうございました。
バグではなくてそういう仕様だったのですね。
こういう使い方ができるということがわかって勉強になりました。
> x=2:1
> ifelse(c(x[1]!=x[2],x[1]==x[2]),1:2,3:4)
[1] 1 4 ■ このスレッドは過去ログ倉庫に格納されています
