Ｒコミュの初歩的な質問

Tomoさん
コメントありがとうございます。
もう少し勉強してみます。

並び替え関連でもう一つ質問させて下さい。
例えば、
> dataa<-c(1,0,0,0,0,0,1,0,1,1,0,0,0)
> table(dataa)
dataa
0 1
9 4
上記のようにdataaのテーブルが作成されますが、
この出力を簡単に
1 0
4 9
と並べ変えるようなコマンドはありますか？
t()コマンドが良いかと思ったんですが、駄目でした。
> t(table(dataa))
宜しくお願い致します。

> ののさんさん
質問の意図を履き違えてたらすみません。
行列の列の順番を入れ替える操作のことだとすると、
dataa[,2:1]
でどうでしょうか。

すみません、質問させて下さい。
ある行列Xを共分散行列にし、それを正規化するには
どう書けばよろしいのですか？

続けてで申し訳ありません。
QR分解について教えて頂けないでしょうか。

ある行列X
X =

2 3 1 7 2
6 8 0 1 6
4 0 7 0 9
1 3 2 0 0


をQR分解して、

X = Q*R となるような X と同じ次元の上三角行列 R とユニタリ行列 Q を出力するにはどうしたら出来ますか？

私はmatlabをRに移植する作業をしているのですが、
これをmatlabで書くと、[Q,R] = qr(X) で、結果は以下のようになりました。

Q =

-0.2649 -0.2000 0.2067 -0.9204
-0.7947 -0.4000 -0.3963 0.2267
-0.5298 0.8000 0.2786 0.0412
-0.1325 -0.4000 0.8501 0.3160

R =

-7.5498 -7.5498 -4.2385 -2.6491 -10.0664
0 -5.0000 4.6000 -1.8000 4.4000
0 0 3.8568 1.0506 0.5428
0 0 0 -6.2160 -0.1099

宜しくお願い致します。
matlabでは
http://infoshako.sk.tsukuba.ac.jp/ShakoDoc/MATLAB5/jhelp/techdoc
/ref/qr.html

http://www.mathworks.co.jp/access/helpdesk_archive_ja_JP/r13/help/toolbox/matlab/ref/index.html?/access/helpdesk_archive_ja_JP/r13/help/toolbox/matlab/ref/func_b12.html&http://www.google.co.jp/search?q=qr%E5%88%86%E8%A7%A3%E3%80%80matlab&lr=lang_ja&ie=utf-8&oe=utf-8&aq=t&rls=org.mozilla:ja-JP-mac:official&client=firefox-a
を参考にしました。

Zenwさん
ありがとうございます。
dataaは行列ではないので、
tableの後に[2:1]を挿入することで目的がかないました。

> table(dataa)[2:1]
なぜこんなことをしたかったかというと、
binom.testの引数にtableでも受け付けてくれることに気が付いたんですが
number of successが反対になってしまうのをひっくり返す方法を模索していました。
とりあえず、以下の方法で目的が達成できました。
> binom.test(table(dataa)[2:1])
ありがとうございました。

>たりぃさん
QR分解については
Q = qr.Q(qr(X))
R = qr.R(qr(X))
で求まるようです。ちょっと冗長ですが。

＞Zenw さま
QR分解、出来ました　
ありがとうございました

93で質問させて頂いた　たりいです。
補足で書かせてください。　
まず始めに

行列X
[,1] [,2] [,3] [,4] [,5]
[1,] 2 3 1 7 2
[2,] 6 8 0 1 6
[3,] 4 0 7 0 9
[4,] 1 3 2 0 0
を共分散行列にすると、
> cov(X)
[,1] [,2] [,3] [,4] [,5]
[1,] 4.9166667 3.8333333 -0.1666667 -2.0000000 6.916667
[2,] 3.8333333 11.0000000 -8.6666667 0.3333333 -1.833333
[3,] -0.1666667 -8.6666667 9.6666667 -4.3333333 7.500000
[4,] -2.0000000 0.3333333 -4.3333333 11.3333333 -4.666667
[5,] 6.9166667 -1.8333333 7.5000000 -4.6666667 16.250000
になりました。
これは、

N>1 ここで N は、観測値の数の場合、N-1 によって正規化したものです。これによって cov(X) は、観測値が正規分布で得られる場合、共分散行列の最良の不偏推定となります。N=1 の場合、cov は N によって正規化します。
（http://www.mathworks.com/access/helpdesk_ja_JP/help/techdoc/index.html?/access/helpdesk_ja_JP/help/techdoc/ref/cov.html&http://www.google.co.jp/search?hl=ja&client=firefox-a&rls=org.mozilla%3Aja-JP-mac%3Aofficial&q=cov++matlab+%E4%BB%8A%E6%97%A5%E5%88%86%E6%95%A3&btnG=%E6%A4%9C%E7%B4%A2&lr=lang_ja&aq=f&oq=　からの引用です。）

私が計算したいのは、
　matlabではcov(X,1)で表現されるものなのですが、
これはN によって正規化し、観測値の平均に関する 2 次モーメント行列を作成たものだそうです。

実際に　matlab　で計算すると、

>> cov(X,1)

ans =

3.6875 2.8750 -0.1250 -1.5000 5.1875
2.8750 8.2500 -6.5000 0.2500 -1.3750
-0.1250 -6.5000 7.2500 -3.2500 5.6250
-1.5000 0.2500 -3.2500 8.5000 -3.5000
5.1875 -1.3750 5.6250 -3.5000 12.1875

でした。
これを　R　でかくにはどうしたら宜しいのですか。
長文申し訳ありません。　宜しくお願い致します。

>たりぃ さん
僕は統計の知識がまったく無いのでなんともいえないのですが、
たりぃさんの内容と
[R-tips -基本統計量の算出]http://cse.naro.affrc.go.jp/takezawa/r-tips/r/59.html
を参照したところ、どうやら
cov(X) は var(X)
cov(X,1) は var(X)*(N-1)/N
で表現できるようです。Nは観測値の数です。

＞Zenw　さま
度々　ありがとうございます。
行列Xの行数を　 ncol(X) で求め、
N<-ncol(X) で入れてみると出来ました。
URLもありがとうございます。　頑張ります。

申し訳ないのですが　また行列の並び替えに関して質問させてください。

行列X
> X
[,1] [,2] [,3] [,4] [,5]
[1,] 2 3 1 7 2
[2,] 6 8 0 1 6
[3,] 4 0 7 0 9
[4,] 1 3 2 0 0 を

[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 0
[2,] 2 3 1 0 2
[3,] 4 3 2 1 6
[4,] 6 8 7 7 9
のように各列について昇順に並べ替えるにはどうしたらできるでしょうか。
X[order(X[,1]),] #一列目について並び替えた行列をつくり、
(X[order(X[,1]),])[,1] #その1列目だけ取り出す、

この作業を各列について繰り返し

matrix(c((X[order(X[,1]),])[,1],(X[order(X[,2]),])[,2],(X[order(X[,3]),])[,3],(X[order(X[,4]),])[,4],(X[order(X[,5]),])[,5]),ncol=ncol(X))
とすれば何とかできたのですが、これでは長くなってしまううえ、一般の場合に対応しないと思っています。

私なりに、
> for (i in 1:ncol(X)){
+ a=(X[order(X[,i]),])[,i]
+ print(a)
+ }
と書いてみたのですが
[1] 1 2 4 6
[1] 0 3 3 8
[1] 0 1 2 7
[1] 0 0 1 7
[1] 0 2 6 9　と逐一表示されてしまいました。

>たりぃ　さん
Rにはapplyという関数が用意されてます。
今回のことでしたら
、
apply(X,2,sort)

でいけるのではないでしょうか。

ちなみに
applyは行列の一部に関数を作用させるもので、
apply(X, MARGIN, FUN)
の形で使います。
X は　操作したい行列です。
MARGIN は　作用させたい要素のことで、1のとき行、2のとき列、c(a,b)のときX[1,b]になります。
FUN　は作用させる関数です。今回は昇順に並び替えるとのことで[sort]がいいと思います。

＞しょうさま
　お礼＆報告が遅くなりすみません。
書いて下さったようにやると出来ました。cbindで列ベクトルの結合を行うのですね
ありがとうございました

＞Zenw　さま

本当に何度もお世話になっております。　
apply(X,2,sort) で出来ました。詳しい解説もありがとうございます
関数applyをこれから色々応用させていきたいと思います

度々申し訳ありません　質問させて下さい。

関数
SELF <- function(X,Y,beta,r,metric,kNN){
略
return(list(T=T,Z=Z))
}
で結果として出力された T,　Zを、他の関数の中に呼び出して使うにはどうしたら出来ますか　関数呼び出しの方法を教えて頂けないでしょうか。

106: たりぃ さん

retv <- SELF(X,Y,beta,r,metric,kNN)
another.function(T=retv$T, Z=retv$Z)

というようなことをしたい、ということでしょうか?

107: とおる　さま
ありがとうございます。　
とおるさまが書いて下さったのは、関数SELFをretvと置き換え、another.function　の入力引数にする、というものでしょうか。

すみません、少し補足をさせて下さい。
具体的には、 demo_SELFと言う関数内に、SELFでの結果T,Zを使っています。
このdemo_SELFはplotをするためのもので、

plot(cbind( -T[1] , T[1])*100,cbind(-T[2] , T[2])*100,col='green',type='l',lty=2,ann=F
)
・・・・（略）
のように線を引きたいと考えております。このdemo_SELF内でのSELFの関数＆T、Zの呼び出しは可能でしょうか。

108: たりぃ さん

> 関数SELFをretvと置き換え、another.function　の入力引数にする、というものでしょうか。

関数の戻り値にretvという名前をつけ、another.functionの引数にしました。

> 出力された T,　Zを、他の関数の中に呼び出して使うにはどうしたら出来ますか

ということでしたので。


以下のような感じのことをしたいのでしょうか?

myfun1 <- function(x) {
T <- 1:x
Z <- rev(T)
list(T=T, Z=Z)
}

myfun2 <- function(x) {
retv <- myfun1(x)
T <- retv$T
Z <- retv$Z
plot(T, Z)
}

myfun2(10)

109 :とおる　さま
ありがとうございます。
関数の戻り値がretvだったのですね、失礼致しました。
はい、まさに109で書いて下さったような事がしたかった事です。
応用頑張ります。

いつもお世話になっております。
立て続けで本当に恐縮なのですが、どなたかご教授して頂けないでしょうか。

データの欠損値の際の処理に悩んでいます。
今、例えばtest1.txt　という、

158 51
162 55
111 72
173
64　　この様なデータがあったとします。

これを欠損値をなくして、
158 51
162 55
111 72

とし、新たなtxtとするにはどうしたら良いでしょうか。
実際に扱うデータは１万行以上有り、手作業で出来ません。
一般にどのように書けば宜しいのか教えて頂けると嬉しいです。宜しくお願い致します。

＞みっひーさん
78への遅いコメントです。

力技でなんとかStataとほぼ同じ出力を返す信頼区間を求める関数が出来ました。
もちろんStataのoptionコマンドまでは対応していません。

ci<-function(x){
dig <- as.numeric(options("digits"))
options(digits=7)
mcc <- match.call()
mct <- as.character(mcc[2])
tts <- t.test(x)
s.e <- sd(x)/sqrt(length(x))
cat(" Variable | Obs Mean Std. Err. [95% Conf. Interval]\n")
cat("-------------+---------------------------------------------------------------\n")
cat(ifelse(nchar(mct)>12, abbreviate(mct, minlength=12, method=c("both")),
formatC( mct, width=12)), " |",
formatC( length(x), digits=7, width=11) ,
formatC( mean(x), digits=7, width=12) ,
formatC( s.e, digits=7, width=12) ,
rep(" ", 8), sep="")
cat(tts$conf.int, "\n", sep=" ")
options(digits=dig)
}

多分、素人くさい関数だと思いますが、Rの良い勉強になりました。
改善点があれば是非教えて下さい。

>たりぃ さん
えっと、そのtxtのデータ量によるかと思います。
もしRで扱える程度の量、つまりread.tableなどでRの変数にデータを格納できる程度の量であるなら、以下のようにして、欠損値をはじけます。
※いま"test1.txt"のあるフォルダは"c:/test/"にしてあります。

x<-read.table("c:/test/test1.txt",sep=" ",header=F)

x<-x[(1:length(x[,1]))[is.na(x[,2])==F],] #欠損値(NA)の除去

write.table(x, "C:/test/output.txt", quote=F, col.names=F, append=T)

こんな感じでうまくいくかと。
もし、Rに格納できないほど巨大なデータだとしたら、たぶんＣ言語なりで処理したほうが早そうですね。

追記、
ファイル分割してデータ量を抑えてやればいいだけの話でした……

111: たりぃ さん

他にも

x <- na.omit(x)
x <- na.exclude(x)
x <- x[complete.cases(x), ]

とか。

>とおる さん
ごり押しでやってしまった・・・。そんな便利なものがあるとは。勉強になります。

116: Zenw さん。

第2列に NA がある列を取り除くなら Zenw さんの方法でも良いですが、
x[!is.na(x[,2]), ]
と書いた方がシンプルかな。

たくさんの列のうち、どこかに NA のある行を取り除く場合、
is.na() を使うのであれば
x[apply(is.na(x), 1, sum) == 0, ]
という方法も。

＞Zenwさま　とおるさま　

111 で質問させて頂いたたりいです。
遅くなり、申し訳有りません。おかげさまで解決しました。
ご丁寧に、また何通りも教えて頂きありがとうございました。
とても勉強になりました。

すみません　メモリ増設について質問です。
今、一万行の行列を読み込んでデータ解析をしているのですが、


エラー： サイズ 827.9 Mb のベクトルを割り当てることができません
R(218,0xa088a720) malloc: *** mmap(size=868114432) failed (error code=12)
*** error: can't allocate region
*** set a breakpoint in malloc_error_break to debug
R(218,0xa088a720) malloc: *** mmap(size=868114432) failed (error code=12)
*** error: can't allocate region
*** set a breakpoint in malloc_error_break to debug

とでて、グラフなどの結果が表示されません。
（小さなデータだと上手くできるのですが・・・）
Rでのメモリ増設方法をご存知の方、ご教授していただけないでしょうか。
　read.tableで、データの読み込みまでは出来るようです。

対策として、プログラムの途中でデータの大きい物は
rm(X)などして軽くしていっているのですが、それでも上手くいきません。

このコミュニティのトップで紹介されている参考サイトの1番にあげられているRjpWikiで、「メモリ」で検索してみると、以下のページがヒットしました。

Windows版RのFAQ
http://www.okada.jp.org/RWiki/?Windows%C8%C7R%A4%CEFAQ

関係する部分を引用します。


2.7 R が使うメモリーのサイズに制限があるように見えるのですが. (There seems to be a limit on the memory it uses!)

実際に制限があります.

その上限値は, コマンド・ラインのフラグ --max-mem-size で設定されます (「Windows 版 R をインストールするには? (How do I install R for Windows?)」参照). デフォールト値は, 物理的な RAM のサイズと 1Gb の小さいほうの値です. 10M 以上の任意の値に設定することができます (10M 未満では R は動作しません). ですが, 注意してください: Windows には (ほとんどのバージョンで) ユーザーの仮想メモリーに 2Gb の最大値があり, その一部はプロセスによって予備に取っておくことができますが, 使用することはできません.

メモリー・マネジャーの働き方のため, 空きメモリーがあっても R がそれを利用できない可能性があります.

メモリーの用法についての情報は, ?Memory および ?memory.size をご利用ください. 動作中の R のセッションの中で memory.limit を呼ぶことにより, メモリーの上限値を大きくすることができます.

フラグ --max-mem-size の値があまりに大きすぎると R の起動に失敗することが知られています: その上限値は Windows 2000 Professional では約 1.7Gb のようです.

別のメモリー・マネジャーを使うように R をコンパイルすることができます. そのメモリー・マネジャーは, 大きいサイズのメモリーを使うのに優れていますがかなり遅くなります (ある種の作業では R が数倍遅くなります).



参考になれば幸いです。

>とおるさま
ありがとうございます。参考にさせて頂きました。

現在Windows VistaでR2.6.0を使っています。
新しいパッケージをダウンロードしたところ、
「R2.7.2で作られた」とあってRをアップデートしなければなりません。

パッケージはツールバーから更新ができるようですが、
Rそのもののアップデートはどのようにしたらよいのでしょうか？
非常に初歩的な質問で恐縮ですがご存じの方がいらっしゃいましたらよろしくお願いします。

＞きんちゃん　さん

ありがとうございます。
新しいバージョンをインストールするのですね。

パッケージのダウンロードが大変だ…。

Stataで論理値を要素に持つ変数(例として変数logistic)に対して
>tab(logistic) を実行すると次のようなきれいな表ができます。
logistic | Freq. Percent Cum.
------------+-----------------------------------
0 | 6,260 60.37 60.37
1 | 4,110 39.63 100.00
------------+-----------------------------------
Total | 10,370 100.00
Rでも実現したかったので関数を書いてみました。
>statatab(logistic)
logistic | Freq. Percent Cum.
------------+-----------------------------------
0 | 6260 60.37 60.37
1 | 4110 39.63 100.00
------------+-----------------------------------
Total | 10370 100.00
不満は2x2にしか対応していないことです。
３桁毎のカンマも出ませんがそれは良いと思っています。
上手く動いているようですが、初心者なので何か間違い・アドバイスなどあればお教えください。
#####
#
# 論理型の変数に対して（のみ）、Stata風の2x2 Tableを作成する。Stataでは >tabulate(変数名) or tab(変数名）
#
#####

statatab <- function(x){ # xは論理型の変数
mcc <- match.call() # as.character(mcc[1])にstatatab
mct <- as.character(mcc[2]) # as.character(mcc[2])に変数名xが入る
tab <- table(x) # 度数による2x2 table
tabr <- 100*prop.table(table(x)) # ％による2x2 table
cat(
if(nchar(mct) > 11){ # 変数名xが11文字を超える場合には省略名を使用
abbreviate(mct, minlength=11, method=c("both"))
} else {
formatC(mct,width=11) # 変数名xが11文字以下の場合にも11文字分のスペース確保
},
"| Freq. Percent Cum.\n")
cat("------------+-----------------------------------\n")
if(tabr[1]==100){ # 変数xの要素が100%同じの場合
cat(ifelse(x[1]==0, " 0 |", " 1 |"), # 変数xの要素が全て0か1かを判断
sprintf("%10.0f %10.2f %10.2f\n", tab[1], tabr[1], tabr[1]))
cat("------------+-----------------------------------\n")
cat(" Total |", sprintf("%10.0f %10.2f", sum(tab),100.00))
} else { # 変数xの要素が100%同じでない場合
cat(" 0 |", sprintf("%10.0f %10.2f %10.2f\n", tab[1], tabr[1], tabr[1]))
cat(" 1 |", sprintf("%10.0f %10.2f %10.2f\n", tab[2], tabr[2], 100.00))
cat("------------+-----------------------------------\n")
cat(" Total |", sprintf("%10.0f %10.2f", sum(tab),100.00))
}
}
整形した文が崩れないように書けないかな。

すみません、
Rでロジスティック回帰モデルの予測確率の信頼区間を
Bootstrap法で求めるプログラムというのは存在いたしますでしょうか。

もし御存知のかたがいらっしゃいましたら、ご教示願えませんでしょうか。

ロジスティック回帰分析で、
glm.fit:数値的に０か１である確率が生じました
と警告メッセージが出てしまいます
何が悪いのか分からないので、教えといただきたいです。