初心者のための R および RjpWiki に関する質問コーナー
新規投稿はできません
青木繁伸 (2005-11-07 (月) 12:40:49)
Macintosh OS 10.4.3 で,R 2.2.0 を使っていますが,R を起動すると vignette という関数(<environment: namespace:utils>)を呼び出して大変時間を食います(今までのバージョンではこのようなことはなかった。少なくともこんなに時間がかかる処理はなかった)。どこかの設定でこれをやらないようにできるのでしょうか。~
現在は,ダミーの vignette 関数を定義してバイパスするようにしております。
どうむ (2005-10-26 (水) 11:08:06)
株式市場の(主に国内)分析等でRをバリバリ使っていますが、「比較的安価で広く出回っている計量経済分析向けの商用ソフト」って、なになのか、気になります。S-PlusとかSASとかSPSSとかSHAZAMとかTSPとか使ったことがありますが、どれも先の要件にあてはまらないような気がします。思い当たるのはS-Plusぐらいでしょうか。なので、とても、気になります。(計量経済分析関連のページ、作っていただいて、意味もなく喜んでいたりします)
[[ ]](2005-10-26 (水) 10:54:37)
Rを使ったICA(独立成分分析)のチュートリアルってどこかにないでしょうか?
山形 (2005-10-26 (水) 10:52:53)
時系列データの解析を行なうにあたり、Box-Jenkins ARIMAを使いたいと考えています。
Rでもこの方法は出来ますか?出来るとしたらどのパッケージをインストールしたら良いのでしょうか?
それとも、普通のARIMAで出来るのでしょうか?
CRAN Task View: Computational Econometrics - [ このページを訳す BETA ] Time series modelling : Classical time series modelling tools are contained in the stats package and include arima() for ARIMA modelling and Box-Jenkins-type analysis. Furthermore stats provides StructTS() for fitting structural time ... cran.r-project.org/src/contrib/Views/Econometrics.html - 12k - キャッシュ - 関連ページみたいのが出てきますが -- 2005-10-26 (水) 11:13:28
ちょろべぇ (2005-10-25 (火) 13:34:32)
現在はじめてUNIXでのR実行にチャレンジしています。
ただ、そこでPDF形式で画像を出力しようとしますと
”バスエラー - コアダンプしました。”
というエラーが出てしまい、異常終了いたします。
またUNIXとLINUXの双方でプログラムをじっこうしてみた所、
UNIXではコアダンプしてしまいますが、なぜかLINUXでは問題なく
動作します。
本件、どのように解決すべき問題なのか、ウェブ等色々調べましたが
まったく見当がつかない状態です。お気づきの点等ありましたら
アドバイスいただければと思います。
====== 以下 詳細情報
・Rのバージョン 2.0.1(UNIXもLINUXも)
・実行したOS ソラリス8(UNIX)、スージーLINUX
・実行したコマンドTestChars <- function(encoding="ISOLatin1", ...) { pdf(encoding=encoding, ...) par(pty="s") plot(c(-1,16), c(-1,16), type="n", xlab="", ylab="", xaxs="i", yaxs="i") title(paste("Centred chars in encoding", encoding)) grid(17, 17, lty=1) for(i in c(32:255)) { x <- i y <- i points(x, y, pch=i) } dev.off() } ## there will be many warnings. TestChars("ISOLatin2")#↑?pdfで出てくるサンプルスクリプト
・読み込んでいるライブラリー Rの標準もの
初心者K (2005-10-23 (日) 18:50:45)
かなり初歩的な質問で申し訳ありません.
毎回同じ大きさグラフをつくりたいのですがデバイスの大きさを指定できません.
どなたかご教示お願いします.
bob3 (2005-10-23 (日) 01:46:36)
連続尺度の変数を名義尺度の変数にコーディングしたいと思っています。
たとえば、以下の例では AREA が1〜5の場合は「P」という区分に、6〜9の場合は「Q」という区分にコーディングしたいのですが、全て「Q」になってしまいます。
x <- data.frame(AREA=1:9) AreaCode <- c('P','Q') ifelse(x$AREA>5, x$A2<-AreaCode[1], x$A2<-AreaCode[2]) # ここでは期待した結果が出力されるのですが、 x # データフレームをみると全て「Q」になっています。
上手なやり方があれば、教えてください。
なお、使っているのは 闇R 2.2.0 です。
よろしくお願いいたします。
> x <- data.frame(AREA=1:9) > AreaCode <- c('P','Q') > ifelse(x$AREA>5, x$A2<-AreaCode[1], x$A2<-AreaCode[2]) [1] "Q" "Q" "Q" "Q" "Q" "P" "P" "P" "P" > x AREA A2 1 1 Q 2 2 Q 3 3 Q 4 4 Q 5 5 Q 6 6 Q 7 7 Q 8 8 Q 9 9 Qifelse の第2,第3オペランドで代入しているから変なことになるようですね。 不等号の件はご愛敬か
ifelse(x$AREA>5, x$A2<-AreaCode[1], x$A2<-AreaCode[2])は冗長なので
x$A2<-ifelse(x$AREA>5, AreaCode[1], AreaCode[2])これで,問題解決。
x$A2<-AreaCode[1+(x$AREA>5)]などとも -- 2005-10-23 (日) 09:36:34
> x <- data.frame(AREA=1:9) > AreaCode <- c('P','Q') > ifelse(x$AREA<=5, x$A2<-AreaCode[1], x$A2<-AreaCode[2]) [1] "P" "P" "P" "P" "P" "Q" "Q" "Q" "Q" > x AREA A2 1 1 Q 2 2 Q 3 3 Q 4 4 Q 5 5 Q 6 6 Q 7 7 Q 8 8 Q 9 9 Qそして、ご回答いただいた方法にするとうまく行きました。
> x <- data.frame(AREA=1:9) > AreaCode <- c('P','Q') > x$A2<-ifelse(x$AREA<=5, AreaCode[1], AreaCode[2]) > x AREA A2 1 1 P 2 2 P 3 3 P 4 4 P 5 5 P 6 6 Q 7 7 Q 8 8 Q 9 9 Q > x <- data.frame(AREA=1:9) > AreaCode <- c('P','Q') > x$A2<-AreaCode[1+(x$AREA>5)] > x AREA A2 1 1 P 2 2 P 3 3 P 4 4 P 5 5 P 6 6 Q 7 7 Q 8 8 Q 9 9 Qところで、名義尺度の水準が三つ以上の場合は、以下のように ifelse を入れ子にする方法のほかに上手いやり方(冗長でない方法)はありますでしょうか。
> x <- data.frame(AREA=1:9) > AreaCode <- c('P','Q','R') > x$A2<-ifelse(x$AREA<=3, AreaCode[1], ifelse(x$AREA>3 & x$AREA<=6, AreaCode[2], AreaCode[3])) > x AREA A2 1 1 P 2 2 P 3 3 P 4 4 Q 5 5 Q 6 6 Q 7 7 R 8 8 R 9 9 Rよろしくお願いいたします。 -- bob3 2005-10-23 (日) 13:22:47
x$A2 <- AreaCode[as.integer((x$AREA-0.00001)/3)+1]のようなのも一つの方法ですね(0.00001 というのは本当は良くない)
x$A2 <- AreaCode[as.integer((x$AREA-1)/3)+1]データに応じた,臨機応変なプログラミングが必要になるでしょうが,for と if-elseif-else で地道にコーディングするのが一番一般的ですね。使い捨てプログラムなら,あれこれ考えて時間をつぶす方がもったいないですから。
x$A2<-ifelse(x$AREA<=3, AreaCode[1], ifelse(x$AREA>3 & x$AREA<=6, AreaCode[2], AreaCode[3]))は冗長で,
x$A2<-ifelse(x$AREA<=3, AreaCode[1], ifelse(x$AREA<=6, AreaCode[2], AreaCode[3]))とするだけでよい -- 2005-10-23 (日) 17:02:02
山形 (2005-10-21 (金) 18:32:18)
Rを用いて時系列解析を行ないたいと考えています。Rは初心者です。
先輩から dynlm というパッケージを用いてやると良いと聞いたのですが、
help(dynlm)を見たときに、Examplesに## multiplicative SARIMA(1,0,0)(1,0,0)_12 model fitted ## to UK seatbelt data uk <- log10(UKDriverDeaths) dfm <- dynlm(uk ~ L(uk, 1) + L(uk, 12))というプログラムがありました。
L(uk,1)やL(uk,12)の1や12は何を表しているのでしょうか?
これと似たようなプログラムを実行したいので、よろしくお願いします。、
花子 (2005-10-20 (木) 14:55:17)
大規模データを用いて数量化理論を実行したいと思っている
Rの超初心者です。
掲示板の記事によると
「Rは実メモリ又は1024Mという上限がある」ようですが、
これは明確な情報としてどこに明記されているのでしょうか?
(疑っているわけではありませんが、私の中でホントに「1024Mの上限があるのか?」で止まっています。)
実際に、試しに1Gから2Gにメモリを増やしてRを起動し、
memory.limit()をみてみましたが、
1Gの時と変わらないのをみると、やはり1024Mが上限なのかなぁという気はしております。
「どこに明記してあるか?」というあいまいな質問ですが、
人からの情報ではなく、なにか説明書のようなものに上限が書かれているのでしょうか?という質問です。
近藤 (2005-10-20 (木) 08:59:58)
例えばdata(sleep) sleep[sleep$g==2,]とすると
>sleep[sleep$g==2,] extra group 11 1.9 2 12 0.8 2 13 1.1 2 14 0.1 2 15 -0.1 2 16 4.4 2 17 5.5 2 18 1.6 2 19 4.6 2 20 3.4 2となります。
この左端の番号 11 - 20 を 1 からふり直すにはどうしたら良いのでしょうか?
リサ (2005-10-14 (金) 08:03:35)
intvlest<-function(x,n,sd,a){ h<-abs(qnorm((1-a)/2)*sd/sqrt(n)) return(c(x-h,x+h)) } conf.intvl<-function(n,k,a=0.95,mu=0,sd=1){ Llim<-rep(0,k) Ulim<-rep(0,k) for(i in 1:k){ x<-rnorm(n,mu,sd) lim<-intvlest(mean(x),n,sd,a) Llim[i]<-lim[1] Ulim[i]<-lim[2] } matplot(cbind(Llim,Ulim)) abline(h=mu) ng<-0 for(i in 1:k){ if(Llim[i]>mu|Ulim[i]<mu)ng<-ng+1 lines(c(i,i),c(Llim[i],Ulim[i])) } return(ng/k) }初めまして。
以上のプログラム、学校の授業で作成したのですが、黒板に書かれたものをただ打ち込んだだけで、さっぱり意味が分かりません・・・・
intvlestっていうプログラムに何かを代入して、グラフを作成しているのは分かるんです。
これが結果的に何を計算しているかも分かるんです。
中心極限定理ですよね?
でも、それぞれの言語の意味が分かりません。
かなり初歩的なことだとは思うのですが、このプログラム1行1行が何の作業しているのか軽くコメントしてくださると嬉しいです。
よろしくお願いします。
conf.intvl <- function(n, k, a = 0.95, mu = 0, sd = 1) { x <- apply(matrix(rnorm(n*k, mu, sd), n, k), 2, mean) h <- qnorm((1-a)/2)*sd/sqrt(n) lim <- cbind(x+h, x-h) matplot(lim) abline(h = mu) sapply(1:k, function(i) lines(c(i, i), c(lim[i,1], lim[i,2]))) return(sum(lim[,1] > mu | lim[,2] < mu)/k) }こんな風にしてみましたが,sapply が,「いやん」な感じ -- 2005-10-14 (金) 15:19:57
たろう (2005-10-13 (木) 17:32:01)
R 2.2.0で,maptoolsのzipファイルをR上でインストールしたのですが,
以下のように表示されてしまいます.改めてzipファイルをXP上でlibraryに
展開しても同じでした.library()と入力すると,mapttolsも表示されるのですが,library(maptools)と入力するとエラーが表示されます.
対策をご教示くださいませ.
utils:::menuInstallLocal()
package 'maptools' successfully unpacked and MD5 sums checked
updating HTML package descriptionslibrary(maptools)
要求されたパッケージ foreign をロード中です
エラー:'%s' が要求したパッケージ '%s' は見つけられませんでした
takahashi (2005-10-11 (火) 20:54:48)
Rとは直接関係ないのでここに書いてみます(説明はRでしてみます).
> x<-c(-3,-2,-1,1,2,3) > y1<-c(0.1,0.2,0.3,0.7,0.8,0.9) > y2<-c(0.2,0.3,0.4,0.6,0.7,0.8) #データは適当ですで,
y1~N1=pnrom(x,0,sd1),y2~N2=pnorm(x,0,sd2) #これはRのコードではありません.
と正規分布に従うことを仮定して,N1,N2に対するsd1,sd2の推定値と信頼区間は得ることが出来ます.ここで,N3~N1+N2としたとき,N3に対するsd3(推定値はsd1+sd2)の信頼区間を得るにはどのような方法があるのでしょうか?
ちなみにあんまり関係ないんですが,↑の過程でpsignifit(http://www.bootstrap-software.com/psignifit/)のR用ラッパー(historyにはto doになっていたものの開発停止でリリースされる見込み無し)を書いたんですが,c側のコードを実行中にたまに落ちます.メモリ使いすぎなのかな?
> y1 <- rnorm(200, mean=5, sd=3) > y2 <- rnorm(200, mean=8, sd=2) > y3 <- y1+y2 > mean(y1)+mean(y2) [1] 12.82907 > mean(y3) [1] 12.82907 > sqrt(var(y1)+var(y2)) [1] 3.726669 > sd(y3) [1] 3.84084となりますが,試行回数を多くするなりしてシミュレーションしてみれば? -- 2005-10-12 (水) 12:04:08
> sim <- function(n, trial=1000) + { + result <- numeric(trial) + for (i in 1:trial) { + y1 <- rnorm(n, mean=5, sd=3) + y2 <- rnorm(n, mean=8, sd=2) + result[i] <- sd(y1+y2) + } + result <- list(sd3 = mean(result), se = sd(result), result=result) + class(result) <- "sim" + return(result) + } > print.sim <- function(result) + { + cat("sd3=", result$sd3, "?tse=", result$se) + } > sim(10) sd3= 3.477083 se= 0.8488288 > sim(200) sd3= 3.601848 se= 0.1748326 > sim(400) sd3= 3.602084 se= 0.1328714とか。。se は sd/sqrt(2n) かなあ。ちゃんと数学で解かないといけないカモね
> result1 <- result2 <- numeric(100) > for (i in 1:100) { + a <- sim(i*10) + result1[i] <- a$sd3 + result2[i] <- a$se} > plot(1:100, result2) > lines(1:100, result1/sqrt(1:100*20), col="red")そんないい加減でいいのか。。 -- 2005-10-12 (水) 12:33:13
-- [[MKR]] &new{2005-10-12 (水) 18:46:59};
るびお (2005-10-05 (水) 23:37:07)
ここのページに書いてあるRuby-RMathlibを使ってみたいのですが、
http://www.okada.jp.org/RWiki/index.php?cmd=read&page=Ruby
そのためには、libRmath.a standalone libraryをCRANからとってこなければならないようなのです。しかし、どうしても見つかりません。
どこでこのライブラリのwindows版もしくはcygwin版等を手に入れることができるのでしょうか?
また、Ruby-RMathlibを使うのに、R自体はインストールする必要はあるのでしょうか?
環境はWindows XP + cygwinです。
hatena (2005-10-03 (月) 16:56:52)
何故か、lme4のパッケージをインストールした後、読み込みができず、利用できません。
library(lme4)と入力すると、
エラー:'%s' が要求したパッケージ '%s' は見つけられませんでした
という表示が出ます。
何か解決法がございましたら、よろしくお願いいたします。
mori (2005-09-28 (水) 16:54:19)
船尾さんのThe R Tipsの練習問題の解答(p.360)にt.test(data, mu=0)の結果として
p値が0.8417なので、母平均が0であるという帰無仮説は棄却される
と書いてあるのですが、p値>0.05の場合は帰無仮説は棄却されないと思うのですがこの記述は間違いでしょうか。アドバイスお願いします。
ひろ? (2005-09-28 (水) 14:17:41)
はじめまして。WINDOWSで、R2.1.1に、Rmapのインストールを試みています。Rmap_1.1.0.zipをダウンロードしてきて、ローカルからのインストールまではうまくいったみたいなのですが、library(Rmap)を実行すると、
「'Rmap' は有効なパッケージではありません。バージョン < 2.0.0 でインストールされたかも知れません」
というメッセージがでてきてしまいます。
ちなみにbinフォルダにprojとshapelibは入れています。
どのように対処すればよいかどなたかご教示いただけないでしょうか。
よろしくお願いします。
ちょろちゃん (2005-09-21 (水) 23:58:08)
R駆け出し中のちょろと申します。
この度はRのソースを自動インデントするエディターを
ご教授いただきたく投稿しました。
現在、勉強もかねてRのスクリプトを作成しております。
そのスクリプトですが、数百行を超えてきたあたりから、
for,ifなどのカッコの対応を追うのが非常に困難になってきました。
そこで、Javaで言うところのEclipseのソース→フォーマットのような
既存のソースをオートインデントするソフト、ないしはエディタが
何か無いかと思い、xyzzyのR-modeの設定を行いましたが、URL http://plaza.umin.ac.jp/~takeshou/xyzzy/rmode.htmlR-modeではカラーリングはばっちりされるのですが、いかんせん
インデンテーション機能が弱いように感じております。
何か、既存のRのソースを整形するよりよいエディターやソフトウェア
などありましたら、ご教授いただきたく考えております。
どうぞヨロシクお願いいたします。m(_ _)m
lost (2005-09-21 (水) 12:10:56)
multiv パッケージがいつのまにかCRAN から消えています。いといろと検索しましたが行方がわかりません。
どなたか、ご教示願います。
R Site Search によれば
multiv is ORPHANED. If you want it you can go the orphaned subdirectory on CRAN
メンテナーがいなくなったか、他に良いものがあるのか、どちらかでしょう。 -- MKR 2005-09-21 (水) 12:18:32
K.Morgen (2005-09-14 (水) 16:41:34)
「Rで自己組織化マップ」を参照し、前段のlibrary(class)から、掲載中のプロセスは、出力まで、確認されました。レーダチャートの数値に、若干の相違がありましたが。 次に、後半の例に従って、library(som)を実行しましたが、”som"のパッケージはない、とのエラーメッセージが出ます。利用可能なRパッケージリストを開いて見ますと、確かに、somの名称のパッケージは見当たりません。
現在、R2.1.1Patchedを使用中ですが、インストールに問題が有ったのでしょうか。
ショーゾー (2005-09-12 (月) 10:09:31)
Mac OS X 10.4.2 + X11 1.1 + R 2.1.1 + Rcmdr 1.1.1の環境です。Rcmdrで日本語が出たのには驚きでしたが,表示されるフォントが汚いです。これをWin環境でのRcmdrのように奇麗なフォントで表示させるにはどうしたらいいでしょう?
Rcmdrの「ツール」->「オプション...」->「デフォルトのフォント」には「1.0」とありますが,ここの修正でどのように例えばOsakaフォントを指定したら良いのでしょうか?
生半可知識ですが,X11のfontsフォルダにKochi-substituteをコピーしてmkfontdirを実行したのですが,fonts.dirに反映されませんでした。よろしくご教示お願いします。
yanagi (2005-09-10 (土) 17:07:37)
g<-expression(2*x^2+3*x+5)
D(g,"x")
y<-D(g,"x")
plot(y,-10,10)
をしましだがerror message "x and y lengths differ"がだます
どうすればいいですか?
curve(2*x^2+3*x+5, xlim=c(-10,10));ではだめですか? -- Akira 2005-09-10 (土) 20:28:22
Fleischmann (2005-09-08 (木) 13:11:16)
昔のAPLのコードをRに変換したいのですが、APLとRの関数の比較表のようなもの、ないでしょうか?
msasaki (2005-09-08 (木) 11:23:54)
6000行×8列の時系列データから、ある列から100個のデータを抽出しようとしています。Samplingを使うとランダムになるため、うまく行きませんでした。関数を探していましたが無いようでして、何か良い方法はございませんでしょうか。皆様におかれましては、大変ご多忙中のこととは存じあげますが、どうぞご指導のほど、よろしくお願い申し上げます。
nori (2005-09-08 (木) 09:48:14)
はじめまして。1000*65のデータの欠損値をKNNで解析を行いたく、今回初めてRにふれている者です。パッケージEMVをインストール後の作業について参考になるページ等ございましたらよろしくお願いいたします。
初心者 ち (2005-09-06 (火) 12:35:22)
c()の引数に16000個の引数を取ることは可能でしょうか?エラーが出ているのですが,要素数の最大値の調べ方が分かりません.google等でしらべたのですがわかりません(^ ^;)
obs1 <- c(4,1,2,1,1,19,5,1,1,1,6,19,1,1,1,10,1,1,4...)
syntax error
Execution halted
eval(parse(text=paste("c(",paste(1:16000,collapse=","),")")))は問題なく通ります
Katsura (2005-09-05 (月) 13:22:09)
連番の変数(m1〜m30)のそれぞれに対数線形モデル(loglin)の結果を代入したのですが、その成分m1$lrt,m2$lrt,...,m30$lrtという値を一括して取り出したいのです。いままでは、いちいちm1$lrt、m2$lrtと30行続けて入力していたのですが、もっと賢い方法がありそうな気がします。ご教唆いただけないでしょうか。
mitsu5 (2005-09-04 (日) 11:34:42)
凡例に mean +/- SD と書きたいのですが、このプラス マイナス を1つの記号としてどうして入力するのでしょうか? 使用しているのは英語版のRです。faqかも知れませんが、うまく探せませんでしたのでよろしくお願いします。
青木繁伸 (2005-09-04 (日) 11:25:25)
横軸の下に,何か描こうとしても,描けません。x <- 0:10 y <- x^2 plot(x, y) text(5, -2, "a")のようにしても,ウインドウサイズにもよりますが,描けません(上の方だけ見えたり,全く見えなかったり)。~
どうしたらいいでしょうか。
TAK (2005-08-29 (月) 20:46:25)
こんばんは。申し訳ありませんが、質問させてください。
Y X1 X2 X3
18.235 13.256 0.632 0.011
30.444 20.668 0.592 0.017
10.529 7.167 0.654 0.018
19.579 13.307 0.668 0.025
27.143 20.413 0.534 0.006
15.000 12.362 0.515 0.010
17.433 13.061 0.448 0.002
30.250 22.299 0.692 0.008
30.000 22.000 0.500 0.007
18.571 13.582 1.311 0.832
31.026 20.697 1.569 0.585
11.028 7.558 1.526 0.702
19.785 14.200 1.138 0.436
27.452 20.470 1.187 0.525
15.598 12.604 0.784 0.525
17.692 13.388 0.684 0.863
30.928 22.752 1.320 0.266
30.086 22.558 0.688 0.434
19.507 14.512 2.239 1.784
31.082 21.684 2.031 0.671
11.577 7.855 1.751 0.888
19.952 14.945 1.900 0.909
27.970 20.639 1.310 0.671
16.577 13.166 1.094 1.473
18.017 13.711 1.623 1.791
上記のデータに対して、目的変数をY、説明変数をX1、X2、X3とした時に、全説明変数で回帰させたところ、X2の偏回帰係数の符号はプラス、Yに対するX2の単相関係数はマイナスとなりました。そこで、マルチコが起こっているものだと考え、http://aoki2.si.gunma-u.ac.jp/R/find_multico.html にあるfind.multico()関数を用いて、X1、X2、X3に対してマルチコのチェックを行ったところ、結果は"no multico"となりました。気になって、VIFも計算したところ、どの変数も皆3より小さい値となっていました。このような場合に、「マルチコという現象は起こっていない」とみなして良いのでしょうか?また、このデータから得られたY=aX1+bX2+cX3+dという重回帰式は、安定した、将来の予測に有効な式となりうるのでしょうか?どなたかお考えを教えてください。よろしくお願い致します。
tolerance VIF X1 0.92483 1.0813 X2 0.44101 2.2675 X3 0.42018 2.3799となり,TAK さんがおっしゃるように X2もX3 もVIF は3以下ですよね。しかし,3 以下なら良いという基準もそんなに確かなことでしょうか?
X1 のみ 0.97442 X1, X3 0.97401 X1, X2 0.97466 X1,X2,X3 0.98074ということです。
阿部 (2005-08-29 (月) 14:34:46)
質問です。代金 緯度 経度 48000 139.77379 35.70892 103000 139.77705 35.69388 68000 139.7803 35.69616 103500 139.7728 35.696のようなデータをRのplotを用いて緯度経度をプロットし、代金によって色分けしたいのですがどのような方法がありますでしょうか。
たとえば代金が、50000未満の点を青、50000以上,80000未満の点を赤、80000以上の点を黒という風に表示したいと思っております。よろしくお願いいたします。
x <- 1:20 color <- ifelse(x < 6, "blue", ifelse(x < 15, "red", "black")) plot(x, x^2, col=color)2行目で color という変数に付値してから使っていますが,直接使ってもよいですよ。
Akira (2005-08-17 (水) 20:06:29)
matrixデータに対して、row毎にt.testを適用し、estimateを取り出そうとしています。
すると、【以下にエラーt.test.default(x) : データは本質的に定数です】とエラーになります。
おそらく、rowの値が全て等しいからだと思うのですが、t.testのヘルプにはそれを回避するオプションがあるか分かりませんでした。ご教授いただけないでしょうか?
今は以下のようにしていますが、上手く行きません。> # このようなデータです。 > x x1 x2 x3 a 7.864811 7.402886 7.633848 b 8.210283 7.838808 8.024545 c 5.624960 5.624960 5.624960 d 3.451917 3.491341 3.471629 e 9.493148 8.908723 9.200936 f 7.759156 8.401847 8.080501 > #関数を定義しました > fun <- function(x, y){ + if(all(x==mean(x, na.rm=TRUE), na.rm=TRUE)){ + x <- mean(as.numeric(x), na.rm=TRUE) + }else{ + x <- t.test(x)[[y]] + } + return(x) + } > > #そしてapplyでrow毎にfunを実行しました > apply(x,1,fun,5) > > #そうすると、以下のエラーが出てしまいます 以下にエラーt.test.default(x) : データは本質的に定数です3行目の値が全て同じなのが、原因と分かりました。
> x.row3 <- x[3,] > x.row3 x1 x2 x3 c 5.62496 5.62496 5.62496 > apply(x.row3,1,fun,5) 以下にエラーt.test.default(x) : データは本質的に定数ですしかし、人為的に作成したデータでは上手く動きます
> test <- matrix(c(1:12, rep(3,3)), ncol=3, byrow=T) > test [,1] [,2] [,3] [1,] 1 2 3 [2,] 4 5 6 [3,] 7 8 9 [4,] 10 11 12 [5,] 3 3 3 > apply(test,1,fun,5) [1] 2 5 8 11 3all(x==mean(x, na.rm=TRUE), na.rm=TRUE)の条件分岐が上手く動いていないと思いますが、Rではx.row3は全て同じ値に見えます。小数点以下の表示していない桁でおかしくなるのでしょうか?
> str(x.row3) `data.frame': 1 obs. of 3 variables: $ x1: num 5.62 $ x2: num 5.62 $ x3: num 5.62 > x.row3[1] x1 c 5.62496 > x.row3[2] x2 c 5.62496 > x.row3[3] x3 c 5.62496よろしくお願いします。
> mean(x, na.rm=TRUE) x1 x2 x3 7.067379 6.944761 7.006070となるのは分かっているのでしょうか。
> mean(x, na.rm=TRUE) x1 x2 x3 7.067379 6.944761 7.006070については、理解しております。x <- mean(as.numeric(x), na.rm=TRUE)としておりますが、それならばx==mean(as.numeric(x), na.rm=TRUE)にするべきでした。
if(all(x==mean(x, na.rm=TRUE), na.rm=TRUE))の部分を
if(identical(all.equal(var(x[2,]), 0), TRUE))にでもかえれば良いはずです。投稿が前後してしまいましたが、本質的に等しい(つまり浮動小数演算の意味で本質的に等しい)データの分散はゼロですから、t 検定など出来るはずがありませんから、エラーになるのは当然です。
if(identical(all.equal(var(x[2,]), 0), TRUE))のような使い方を知りませんでした。
A380 (2005-08-16 (火) 23:40:38)
度分布を計算するRの関数はないでしょうか?
> x <- rbinom(10000, 10, 0.5) > table(x) x 0 1 2 3 4 5 6 7 8 9 10 9 87 439 1142 2140 2483 1993 1162 422 112 11 > hist(x, breaks=0:11-0.5)$counts [1] 9 87 439 1142 2140 2483 1993 1162 422 112 11オンラインヘルプを良く読んで使ってください。 -- 青木繁伸 2005-08-17 (水) 08:10:42
take (2005-08-16 (火) 20:21:32)
お世話になります。
例えば下記のようなリストの表示そのものを保存するコマンドはどなたかご存知ですか?
コピー&ペーストがもっとも手っ取り早いですが、大量にある場合、a b 別々にwrite.table()と保存してDOS窓でcopy/bによって結合するという方法もありますが、R上で行ってしまいたいです。
capture.output()などではベクター形式で表示されてしまいますし、sinkではRで使うには便利ですが、テキストで別のソフトに移すときは不便です。
なお、当方の環境は2.01 WinXPsp1です。どうぞよろしくお願いします
> a <-matrix(1:4,ncol=2); b <-matrix(4:1,ncol=2) > list <-list(a,b) > list [[1]] [,1] [,2] [1,] 1 3 [2,] 2 4 [[2]] [,1] [,2] [1,] 4 2 [2,] 3 1
> a <-matrix(1:4,ncol=2); b <-matrix(4:1,ncol=2) > list <-list(a,b) > class(list) <- "your.class" > print.your.class <- function(x) + { + cat(list[[1]][1,1], list[[1]][1,2], list[[1]][2,1], list[[1]][2,2],"?n") + cat(list[[2]][1,1], list[[2]][1,2], list[[2]][2,1], list[[2]][2,2],"?n") + } > list 1 3 2 4 4 2 3 1cat の出力先をファイルにする,ファイルに追加出力するとか,なんでもかんでも後は自分で色々改造。自分でなんでもかんでも好き勝手にできると言うところが R のよいところ。 -- 青木繁伸 2005-08-16 (火) 21:14:16
econome007 (2005-08-15 (月) 17:19:25)
初心者です。
MASS,nnet,ts,nlsといったパッケージが見あたらないのですが、どこからダウンロードできるのでしょうか。
ラプラス (2005-08-10 (水) 14:29:03)
Rでポアソン方程式を扱える関数はあるのでしょうか?
いろいろ探してみましたが見つからないのですが?
Lost (2005-08-05 (金) 17:17:45)
R のスクリプトで、再帰アルゴリズムを使ったサンプルはないでしょうか?
> fact2 <- function(n) + { + if (n == 0) 1 + else n*fact2(n-1) + } > fact2(8) [1] 40320 > fact2(3) [1] 6 base に本来ある factorial 関数でチェックしてみましょうか・・ > factorial(8) [1] 40320 > factorial(3) [1] 6いかがでしょうか(レポートの手助けだったらやだな)。
else n*fact2(n-1)のところは,
else n*Recall(n-1)のほうがいいらしい。 -- 青木繁伸 2005-08-05 (金) 17:44:25
fact3 <- function(n) { if (n == 0) 1 else n+fact3(n-1) #加算に変更 }R2.1.0 での実行結果は以下のようになりました。
> fact3(998) [1] 498502 > fact3(999) エラー:protect():プロテクションスタックが溢れました結果がオーバーフローしているようではないですが。関数内の"fact2"の替わりに"Recall"を使用すると、引数がより小さな値で、上のエラーが出てしまいます。
> factorial function (x) gamma(x + 1) <environment: namespace:base>
function(a,b=1){ifelse(a==0,b,Recall(a-1,b+a))}とかにすればスタック消費しなくなるはずですけどね 例えばschemeとか -- takahashi 2005-08-06 (土) 03:10:06
issei (2005-08-04 (木) 14:35:12)
標準パッケージのリファランスはどこで得ることができるのでしょうか?
CRAN拡張パッケージについては、ダウンロードするところにあるのはわかるのですが。
spatialを使いたいのですが困っています。
よろしくお願いします。
いしだ (2005-08-01 (月) 16:56:10)
あまりに初歩的なので、「何でも掲示板」に質問しましたが、指摘を受けて、こちらに移動しました。また、しばらく前に、青木先生の掲示板で質問したのですが、解決しませんでした。
MacOSX10.3.9でR2.1.1が起動しません。具体的には、起動させても、ドックの中でバウンドしたまま、予期せず終了してしまいます。何度インストールしてもだめです。
新しいアカウントからログインすると、普通に使えます。また、自分の通常のアカウントでも、ターミナルで"r"と打つと、ターミナルの中で使えます。x11でも同様に使えます。
なぜでしょうか。ライブラリの中のRとついてるものをけしたり、正常に動くパソコンからコピーしてきたりしたのですが、いっこうに、改善しません。何か、解決方法がわかれば、教えてください。
$ /Applications/R.app/Contents/MacOS/Rで何かメッセージが出ないか.
$ otool -L /Applications/R.app/Contents/MacOS/Rついでにライブラリのリンクがどうなるか. にわかMacユーザですんで,あれですが...
$ gdb /Applications/R.app/Contents/MacOS/Rデバッガのめっせーじ
(gdb) run # デバッガで実行なんか出てくる実行可能なら,Rの方が動くが、もし落ちてるようなら
(gdb) bt # とりあえず落ちるならバックトレースなにもかえって来ないようなら
(gdb) [control+C] # シグナルを送って (gdb) bt # とりあえずバックトレースよくわからなかったらすみません.
Trace/BPT trap二文目
/Applications/R.app/Contents/MacOS/R: /System/Library/Frameworks/Cocoa.framework/Versions/A/Cocoa (compatibility version 1.0.0, current version 9.0.0) /Library/Frameworks/R.framework/Versions/2.1.1/Resources/lib/libR.dylib (compatibility version 2.1.0, current version 2.1.1) /System/Library/Frameworks/WebKit.framework/Versions/A/WebKit (compatibility version 1.0.0, current version 1.0.0) /System/Library/Frameworks/Security.framework/Versions/A/Security (compatibility version 1.0.0, current version 177.0.0) /System/Library/Frameworks/ExceptionHandling.framework/Versions/A/ExceptionHandling (compatibility version 1.0.0, current version 4.9.0) /usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 71.1.3)第三文
GNU gdb 5.3-20030128 (Apple version gdb-309) (Thu Dec 4 15:41:30 GMT 2003) Copyright 2003 Free Software Foundation, Inc. GDB is free software, covered by the GNU General Public License, and you are welcome to change it and/or distribute copies of it under certain conditions. Type "show copying" to see the conditions. There is absolutely no warranty for GDB. Type "show warranty" for details. This GDB was configured as "powerpc-apple-darwin". Reading symbols for shared libraries ....... donerunすると
Starting program: /Applications/R.app/Contents/MacOS/R Reading symbols for shared libraries +.................................................................... done Reading symbols for shared libraries . done Reading symbols for shared libraries . done Program received signal SIGTRAP, Trace/breakpoint trap. 0x90a8d318 in _NSRaiseError ()(Rは長時間かけて起動したものの、全く反応無し。 ウィンドウも表示されず、メニューバーも出ない)
#0 0x90a8d318 in _NSRaiseError () #1 0x90a8d1fc in +[NSException raise:format:] () #2 0x90a2d35c in -[NSString stringByAppendingString:] () #3 0x00006ce0 in -[RController doLoadHistory:] () #4 0x00003ce0 in -[RController awakeFromNib] () #5 0x90a5f0e8 in -[NSSet makeObjectsPerformSelector:] () #6 0x92ea2150 in -[NSIBObjectData nibInstantiateWithOwner:topLevelObjects:] () #7 0x92f93c2c in loadNib () #8 0x92eeae24 in +[NSBundle(NSNibLoading) _loadNibFile:nameTable:withZone:ownerBundle:] () #9 0x92f69d28 in +[NSBundle(NSNibLoading) loadNibFile:externalNameTable:withZone:] () #10 0x92f7b51c in +[NSBundle(NSNibLoading) loadNibNamed:owner:] () #11 0x92f69b90 in NSApplicationMain () #12 0x00002dd0 in _start (argc=1, argv=0xbffffe60, envp=0xbffffe68) at /SourceCache/Csu/Csu -47/crt.c:267 #13 0x8fe1a278 in __dyld__dyld_start ()-- いしだ 2005-08-02 (月) 09:40:35
export LANG=C export LC_ALL=C /Applications/R.app/Contents/MacOS/Rでは動かないでしょうか,動くなら,"C"を"ja_JP.UTF-8"に置き換えるといかがでしょうか -- なかま 2005-08-02 (火) 10:50:40
export LANG=ja_JP.UTF-8 export LC_ALL=ja_JP.UTF-8 /Applications/R.app/Contents/MacOS/Rとしないかぎり、コンソールに日本語が表示されません。
生物系院生 (2005-07-27 (水) 13:00:57)
格子データの回帰分析に使用するパラメータ行列を作る際、function機能の中でfor文を使って変数を入れる位置を指定するという作業がうまくいかなかったので質問させて頂きました。
状況を取りうる最小サイズの格子データで説明すると123 456 789といったデータ(数字は変数名)の場合、例えば5のデータはfirst neighborのみを考慮するマルコフ確率場を仮定した場合2,8のデータと4,6のデータによって推測される。このときそれぞれが縦のデータの影響の度合いと横のデータにの影響の度合い関係するパラメータをもっており、それは以下のように表される。(データは周期条件を仮定して、要素内の数字の1が横要素、2が縦要素に関係するパラメータ、行の数字が予測対象の変数名、列の数字が説明変数の変数名に相当する行列)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] 0 1 1 2 0 0 2 0 0 [2,] 1 0 1 0 2 0 0 2 0 [3,] 1 1 0 0 0 2 0 0 2 [4,] 2 0 0 0 1 1 2 0 0 [5,] 0 2 0 1 0 1 0 2 0 [6,] 0 0 2 1 1 0 0 0 2 [7,] 2 0 0 2 0 0 0 1 1 [8,] 0 2 0 0 2 0 1 0 1 [9,] 0 0 2 0 0 2 1 1 0この行列自身は以下のコマンド
el <- 9 #格子データの総要素数 パラメータ行列では行数、列数に相当 col <- 3 #格子データの一行の要素数、すなわち列数 mat <- matrix(0, nrow = el, ncol = el) xiA <- 1 xiB <- 2 #文字列操作では不都合が多いため、パラメータを数値で仮置き for(i in 1:6){ W <- i - 1 if(ceiling(W/col) != ceiling(i/col)){W <- W + col} E <- i + 1 if(ceiling(E/col) != ceiling(i/col)){E <- E - col} N <- i - col if(N < 1){N <- N + el} S <- i + col mat[i,W] <- xiA mat[i,E] <- xiA mat[i,N] <- xiB mat[i,S] <- xiB } for(i in 7:9){ W <- i - 1 if(ceiling(W/col) != ceiling(i/col)){W <- W + col} E <- i + 1 if(ceiling(E/col) != ceiling(i/col)){E <- E - col} N <- i - col if(N < 1){N <- N + el} S <- i + col if(S > col){S <- S - el} mat[i,W] <- xiA mat[i,E] <- xiA mat[i,N] <- xiB mat[i,S] <- xiB }で作ることができるのですが、このコマンドにおいてxiA,xiBの部分を変数のままにした関数を作ってから1,2を代入しようとしてもうまくいかないのです(ちなみに実際のデータを使用する際には400×400の行列を使用する予定になっています)。そのうまくいっていないコマンドは以下のような感じで
el <- 9 #格子データの総要素数 パラメータ行列では行数、列数に相当 col <- 3 #格子データの一行の要素数、すなわち列数 matB <- function(a,xiA,xiB){ for(i in 1:6){ W <- i - 1 #W,Eは横方向に関するパラメータの位置 if(ceiling(W/col) != ceiling(i/col)){W <- W + col} E <- i + 1 if(ceiling(E/col) != ceiling(i/col)){E <- E - col} N <- i - col #N,Sは縦方向に関するパラメータの位置 if(N < 1){N <- N + el} S <- i + col a[i,W] <- xiA #パラメータを代入する位置の指定 a[i,E] <- xiA a[i,N] <- xiB a[i,S] <- xiB } for(i in 7:9){ W <- i - 1 if(ceiling(W/col) != ceiling(i/col)){W <- W + col} E <- i + 1 if(ceiling(E/col) != ceiling(i/col)){E <- E - col} N <- i - col if(N < 1){N <- N + el} S <- i + col if(S > col){S <- S - el} a[i,W] <- xiA a[i,E] <- xiA a[i,N] <- xiB a[i,S] <- xiB } } base <- matrix(0, nrow = el, ncol = el) matB(Base,1,2)for文で要素の位置指定をすることができていないようでした(エラーメッセージは出ず、全角スペースが混じっていないことは検索済)。
手持ちの資料やネット上の解説では、関数定義で行列を扱っている例は少なく、他の言語の知識にも乏しいため行き詰って質問させていただきました。宜しくお願いします。
base<-matB(base,1,2)でよいのでは?forのインデクス云々というのは関係ないと思いますよ -- takahashi 2005-07-27 (水) 13:23:08
el <- 9 #格子データの総要素数 パラメータ行列では行数、列数に相当 col <- 3 #格子データの一行の要素数、すなわち列数 mat <- matrix(0, nrow = el, ncol = el) xiA <- 1 xiB <- 2 #文字列操作では不都合が多いため、パラメータを数値で仮置き for (i in 1:9) { W <- i - 1 if (ceiling(W/col) != ceiling(i/col)) {W <- W + col} E <- i + 1 if (ceiling(E/col) != ceiling(i/col)) {E <- E - col} N <- i - col if (N < 1) {N <- N + el} S <- i + col if (i > el-col && S > col) {S <- S - el} mat[i,W] <- xiA mat[i,E] <- xiA mat[i,N] <- xiB mat[i,S] <- xiB }まだ無駄があるような。 -- 2005-07-27 (水) 16:55:39
> no <- 9 > X <- matrix(1:81, c(no,no)) # テスト用格子データ > N <- W <- E <- S <- vector("list", no*no) # リスト行列の用意 > dim(N) <- dim(W) <- dim(E) <- dim(S) <- c(no, no) > wN <- wW <- wE <- wS <- matrix(0, no, no) # 重み用行列 > > for (i in 1:no) + for (j in 1:no) { + N[[i,j]] <- matrix(c(ifelse(i==1, no, i-1), j), c(1,2)) # 北(上)側格子の座標(周期条件付き) + wN[i,j] <- 2 # その重み + S[[i,j]] <- matrix(c(ifelse(i==no, 1, i+1), j), c(1,2)) # 南(下)側格子の座標(周期条件付き) + wS[i,j] <- 2 # その重み + E[[i,j]] <- matrix(c(i, ifelse(j==no, 1, j+1)), c(1,2)) # 西(右)側格子の座標(周期条件付き) + wE[i,j] <- 1 # その重み + W[[i,j]] <- matrix(c(i, ifelse(j==1, no, j-1)), c(1,2)) # 東(左)側格子の座標(周期条件付き) + wW[i,j] <- 1 + } > > test1 <- function() { # テスト用関数 + for (k in 1:1000) + for (i in 1:no) + for (j in 1:no) { # (i,j) の上下左右の格子の値と重み + c(X[N[[i,j]]], X[E[[i,j]]], X[S[[i,j]]], X[W[[i,j]]]) + c(wN[i,j], wE[i,j], wS[i,j], wW[i,j]) + } + } > gc(); system.time(test1()) # 実行時間 used (Mb) gc trigger (Mb) max used (Mb) Ncells 168482 4.5 350000 9.4 350000 9.4 Vcells 63009 0.5 786432 6.0 280690 2.2 [1] 2.83 0.03 2.89 0.00 0.00 ## 関数形で表現すると > N <- function(i,j) matrix(c(ifelse(i==1, no, i-1), j), c(1,2)) > S <- function(i,j) matrix(c(ifelse(i==no, 1, i+1), j), c(1,2)) > E <- function(i,j) matrix(c(i, ifelse(j==no, 1, j+1)), c(1,2)) > W <- function(i,j) matrix(c(i, ifelse(j==1, no, j-1)), c(1,2)) > wN <- wS <- function(i,j) 2 > wE <- wW <- function(i,j) 1 > test2 <- function() { + for (k in 1:1000) + for (i in 1:no) + for (j in 1:no) { # (i,j) の上下左右の格子の値と重み + c(X[N(i,j)], X[E(i,j)], X[S(i,j)], X[W(i,j)]) + c(wN(i,j), wE(i,j), wS(i,j), wW(i,j)) + } + } > gc(); system.time(test2()) # 但し時間は45倍かかる used (Mb) gc trigger (Mb) max used (Mb) Ncells 167600 4.5 350000 9.4 350000 9.4 Vcells 61585 0.5 786432 6.0 280690 2.2 [1] 125.99 0.07 131.18 0.00 0.00
小山 (2005-07-25 (月) 11:04:15)
uniroot()を使用したく、Rでプログラムを書いています。元データとして約5000行,4列のデータがあります。その4列目のデータを1行目から1つずつ関数に代入させ、uniroot()を用いて関数の最適解を求め、それをまた別の関数に代入し、最終的な解を得たいのですが、繰り返し関数がうまく走りません。
プログラムとしては、まず4列目のデータの個数を調べ(x)、算出した解を格納する配列を定義し、xがデータの個数以下の場合は繰り返し配列に答えを格納していく様にしようと、whileを使用しました。
作ったプログラムは以下の様になっております。no <- length(data[[4]]) ha <- array(data=NA, dim = c(nd, 7)) x <- 0 while (x <= nd){ x <- x+1 a <- data[x, 4] b <- 7 c <- 42 d <- 10 e <- 50 f <- 40 obj <- function(phi,a) 1/cos((phi*pi)/180)-cosh(a*tan((phi*pi)/180)) res <- uniroot(obj,c(0.1,89.9),a=a) phi <- res$root func <- res$f.root depth <- rep(NA,b) i <- 1:b depth[i] <- d+c+L*(((1/tan(phi*pi/180))^2+1)^0.5 -(1/(tan(phi*pi/180))^2+(1-(f+e*(i-1))/L)^2)^0.5) ha[x, ] <- c(depth, a) }しかし配列は綺麗に出てくるのですが、データはすべてNAのままで、計算結果が格納されていません。Rの参考書やネット等で調べましたが、どの部分が間違っているのかわからず、質問させていただくことにしました。つたない説明で申し訳ありませんが、アドバイスよろしくお願いします。
x <- 0 while (x <= nd) { x <- x+1というのはおかしいですね。x が nd のときに,新たなループが生じますが,その中で x は nd より1大きな値になりますよ。そして配列要素のアクセスに失敗するはずですが。
i <- 1:b depth[i] <- d+c+L*(((1/tan(phi*pi/180))^2+1)^0.5 -(1/(tan(phi*pi/180))^2+(1-(f+e*(i-1))/L)^2)^0.5)というところですが,左辺は depth だけでよいのではないでしょうか( [i] はいらないということ。またdepth <- rep(NA, b) も不要になる)
L <- ((en-1)*e+2*f)/2ご指摘いただいた点を参考にして、以下の様に直してみました(長くなるので行間を詰めました)が、やはり出てくるのはNAだけでした。プログラムを実行ごすると 「ERROR: 構文解析エラーです」と出ます。質問なのですが、このプログラム実行後正確に計算されたかを見るには「ha]と入力し実行すれば良いのですよね?
no <- length(data[[4]]) ha <- array(data=NA, dim = c(no, 8)) b <- 7 c <- 42 d <- 10 e <- 50 f <- 40 obj <- function(phi,a) 1/cos((phi*pi)/180)-cosh(a*tan((phi*pi)/180)) x <- 0 while (x < no){ x <- x+1 a <- data[x, 4] res <- uniroot(obj,c(0.1,89.9),a=a) phi <- res$root func <- res$f.root L <- ((en-1)*e+2*f)/2 i <- 1:b depth[i] <- d+c+L*(((1/tan(phi*pi/180))^2+1)^0.5 -(1/(tan(phi*pi/180))^2+(1-(f+e*(i-1))/L)^2)^0.5) ha[x, ] <- c(depth, a) }また、データは以下の様なデータになっております。1列目から3列目までは時間データで、本プログラムでは使用しません。
8 40 10 0.927 8 40 20 0.965 8 40 30 0.962 8 40 40 0.964 8 40 50 0.967
depth[i] <- d+c+L*(((1/tan(phi*pi/180))^2+1)^0.5の前に全角空白が1つ残っていますね。左辺は depth だけ。
no <- length(data[[4]])は
no <- nrow(data)で。 en を適当に定めたら,ha に,一応数値は入ったようですよ。 -- 2005-07-25 (月) 13:51:20
x <- 0 while (x < no) { x <- x+1を,
for (x in 1:no) {とするのです。そうすれば,x を1,2,...,no として順次くりかえすのだなとすぐ分かりますね。no は対象になるのかならないのかとか,余計な心配をしなくてもよい(エラーの種が減る)。
for (j in 1:no){ a <- data[j, 4] if(a > 1.0){ a == 1.0 }else{ a == a } res <- uniroot(obj,c(0.1,89.9),a=a)原因はaの値なのでしょうか。a > 1.0の場合は次のようなデータになっています。
14 15 0 1.015 14 15 10 1.018 14 15 20 1.017 14 15 30 1.014 14 15 40 1.011 14 15 50 1.013
#ref(): File not found: "quad.png" at page "初級Q&A アーカイブ(3)"
# 関数定義 > quad <- function(x) (x-0.1234)*(x-0.5678) # 0,1 の範囲から解を求めようとしたとき > uniroot(quad, c(0, 1)) 以下にエラーuniroot(quad, c(0, 1)) : f() の端点での値が異なった符号を持ちません # つまり,quad(0)= 0.07006652, quad(1)=0.3788665 # 異なった符号を持たない(両方とも正)ということ # 関数値は正,負,正と変化し区間内に二つの解を持つ > uniroot(quad, c(0, 0.1)) 以下にエラーuniroot(quad, c(0, 0.1)) : f() の端点での値が異なった符号を持ちません # つまり,quad(0)= 0.07006652, quad(0.1)=0.01094652 # この場合は,関数値は,いつも正,つまり,この区間内に解はない! > uniroot(quad, c(0, 0.5)) $root [1] 0.1233996 $f.root [1] 1.843689e-07 $iter [1] 7 $estim.prec [1] 6.103516e-05 # つまり,quad(0)= 0.07006652, quad(0.5)=-0.02553348 # 異なった符号を持つということ # 関数値は正,負と変化し区間内に一つの解を持つ #(もしかしたら3個とか5個とか奇数個の解を持つのかもしれないが)よく吟味してみてください。 -- 2005-07-25 (月) 21:50:39
data <- rbind(c(8,40,10,0.927), c(8,40,20,1.20 ), # エラーになるようにした c(8,40,30,0.962), c(8,40,40,0.964), c(8,40,50,1.20 ) ) # エラーになるようにした A <- data[,4]; no <- length(A) b <- 7; c <- 42; d <- 10; e <- 50; f <- 40 obj <- function(phi, a) 1/cos((phi*pi)/180)-cosh(a*tan((phi*pi)/180)) doit <- function(x) { # 一回の試行用の関数 res <- uniroot(obj, c(0.1,89.9), a=A[x]) phi <- res$root # func <- res$f.root # 不要ですね L <- ((e-1)*e+2*f)/2 y <- (1:b)-1 depth <- d+c+L*(sqrt((1/tan(phi*pi/180))^2+1) -sqrt(1/(tan(phi*pi/180))^2+(1-(f+e*y)/L)^2)) return(c(depth, A[x])) } res <- vector("list", no) # 結果を納めるリストを予めサイズ指定で用意 for(i in 1:no) res[[i]] <- try(doit(i), TRUE) # 全試行 > res [[1]] [1] 75.40177 103.94505 131.67010 158.54344 184.53087 209.59754 233.70809 [8] 0.92700 [[2]] [1] "以下にエラーuniroot(obj, c(0.1, 89.9), a = A[x]) : ?n?tf() の端点での値が異なった符号を持ちません?n" attr(,"class") [1] "try-error" [[3]] [1] 69.74911 91.28434 112.08025 132.11997 151.38677 169.86406 187.53553 [8] 0.96200 [[4]] [1] 69.32652 90.34202 110.62861 130.17047 148.95193 166.95752 184.17200 [8] 0.96400 [[5]] [1] "以下にエラーuniroot(obj, c(0.1, 89.9), a = A[x]) : ?n?tf() の端点での値が異なった符号を持ちません?n" attr(,"class") [1] "try-error"
田中 (2005-07-21 (木) 22:41:35)
n×n行列をすでに持っており、各行ごとに最大値を与える添え字を返させ、それをn×2の行列に格納してやりたい(つまりは、ある行について最大値をもつのはここですよというのを、全ての行について記録させたいということです。)のですが、いろいろな方法を試しても上手く出来ませんでした。
解決法を教えていただければと思います。よろしくお願いします。 一応、作ってみたのは以下のようなものですQ <- matrix(0,1,n) for(b in 1:n){ for(c in 1:n){ Q[1,c] <- x[b,c] P[b] <- which(max(Q), arr.ind=TRUE) } }xが調べたいn×n行列で、Pに格納するように作りました。
col=4 MM<-matrix(runif(col^2),ncol=col) SQM<-function(matrix){ ml<-function(n){ncol=length(n);((1:ncol)[n==max(n)])[1]} rbind(apply(matrix,1,ml),apply(matrix,2,ml)) } SQM(MM)
(x <- matrix(rnorm(12), 3, 4)) which(x == apply(x, 1, max), arr.ind=TRUE)本質は二行目。
> (x <- matrix(c(1,2,3,2,1,2,3,3,2,1,2,3),3,4)) [,1] [,2] [,3] [,4] [1,] 1 2 3 1 [2,] 2 1 3 2 [3,] 3 2 2 3 > which(x == apply(x, 1, max), arr.ind=TRUE) row col [1,] 3 1 [2,] 1 3 [3,] 2 3 [4,] 3 4ところで,arr.ind=TRUE のとき,結果の1列目でソートされない結果が返ってくるのは仕様かな。
x <<- mix2_16 #n×n行列 P <- which(x == apply(x, 1, max), arr.ind=TRUE) P <<- P gn <- sprintf("result_mix%g_16.txt",a) write.table(P,file=gn, sep="?t", row.names=FALSE, quote=FALSE) print(P)のようなプログラムををつくり実行したところ、
row col [1,] 4 4 [2,] 5 4 [3,] 8 4 [4,] 9 4 [5,] 13 4 [6,] 15 4 [7,] 7 7 [8,] 19 25 [9,] 1 29 [10,] 2 29のような結果が返ってきました。(n=38のときで11以下省略)行ごとに最大値を与える添え字を出力させたつもりなので、rowが1,2,3・・・,10・・・と並ぶはずなのですがそうではありません。 なぜなのか分かりますでしょうか? よろしくお願いいたします。-- 田中 2005-07-22 (金) 18:19:24
> set.seed(31415); (x <- matrix(runif(20), c(4,5))) [,1] [,2] [,3] [,4] [,5] [1,] 0.9502223 0.0607455 0.9365462 0.77713002 0.1302115 [2,] 0.3357378 0.3579473 0.8376269 0.08262812 0.9671198 [3,] 0.1330718 0.2271160 0.1275758 0.81371308 0.8880396 [4,] 0.4901114 0.4000740 0.7251304 0.30108535 0.5821669 > cbind(1:dim(x)[1], max.col(x)) # 列ごとの最大値位置の添字 [,1] [,2] [1,] 1 1 [2,] 2 5 [3,] 3 5 [4,] 4 3 > rbind(1:dim(x)[2], max.col(t(x))) # 行ごとの最大位置の添字 [,1] [,2] [,3] [,4] [,5] [1,] 1 2 3 4 5 [2,] 1 4 1 3 2 > cbind(1:dim(x)[1], max.col(-x)) # 最小位置が欲しければ -x とするだけ [,1] [,2] [1,] 1 2 [2,] 2 4 [3,] 3 3 [4,] 4 4
sugimoto (2005-07-14 (木) 18:08:02)
初めて質問させていただきます。
R tipsを購入し、昨日から使い始めました。
今、lowess関数を使っています。
この関数を使って変換されたデータと、変換前のデータをつなぐにはどのような方法がありますでしょうか。もしくは、lowess関数に与えるx,yにIDのような値を持たせることは可能なのでしょうか。
大変つたない質問で申し訳ありませんが、よろしくお願いします。
# 変換前のデータをデータフレームにする(マージ用にIDを付加) ( before <- data.frame(ID=1:nrow(cars), cars) ) ID speed dist 1 1 4 2 2 2 4 10 3 3 7 4 ................ # 変換後のデータをデータフレームにする(マージ用にIDを付加) ( after <- data.frame(ID=1:nrow(cars), lowess(cars)) ) ID x y 1 1 4 4.965459 2 2 4 4.965459 3 3 7 13.124495 .................. # ID で紐付けしてデータをマージ result <- merge(before, after, by="ID") # x と dist は重複しているので x を削除 result$x <- NULL # 結果を表示 result ID speed dist y 1 1 4 2 4.965459 2 2 4 10 4.965459 3 3 7 4 13.124495 ..........................
田中 (2005-07-08 (金) 18:50:56)
以前、他の件で質問させて頂きました。
今回は、また別の事柄になるのですが、ファイル読み込みの際のファイル名の指定についての質問です。
具体的には、C:/sp01.txt,C:/sp02.txt,・・・C:/sp30.txtのようなファイルに対してそれぞれ同じプログラムを実行したいのですが、ファイル数が多くなることも考えて、x01 <-read.table("C:/sp01.txt")をfor文で回せないかと考えております。つまり、以下のようにしてfor文で変数?(ただ、分かりやすいように?にしただけです。)を変えてやり、ファイル名の変更を自動化したいわけです。(?の部分を01、02、03・・・と変更するのをfor文で回したい。また、program(a,b)は各ファイルに対し実行したいプログラムのことです。)
kurikaesi <- function(a,b) { for(? in 1:30) { x <- read.table("C:/sp?.txt") y? <- program(a,b) } }自分で調べては見たものの、結論を得ることが出来なかったので何か方法があればおしえてください。 よろしくお願いします。
for (i in 1:2) { fn <- sprintf("test%i.data", i)) x <- read.table(fn) }のような感じ。
> i=1; sprintf("test%02i.data", i) # 「整数二桁、一桁なら最初を0で埋める」の意味 [1] "test01.data" > i=10; sprintf("test%02i.data", i) [1] "test10.data"こんなのもあり(R のループ範囲はベクトル、リストでもよい) -> 青木さんには怒られるかも :-)
> for(file in lapply(1:20, function(i) sprintf("file%02i.data", i))) cat(file, "?n") file01.data file02.data ........... file19.data file20.data
for (file in paste("file", 1:20, ".data", sep=""))のほうが短い。 -- 青木繁伸 2005-07-09 (土) 11:15:09
Riemann (2005-07-04 (月) 20:59:13)
R Bookに載っているExcelとの連携のところを読んで、「これは便利!」と思いやってみたのですが、1つのセルに長いテキストが入っていると途中できれてしまいます。ncharで測るとちょうど255文字になっています。どこかのオプションで制限サイズを変えればよさそうに思うのですが、どこをどういじったよいのかわかりません。ご教示ください。
K (2005-07-04 (月) 00:58:04)
RにはMatlabのSignal Processing Toolboxに相当するような
信号処理に適したパッケージはありますでしょうか?
fftやウェーブレット関係のパッケージはいくつかあるようですが、
ChebychevやButterWorthなど古典的なフィルタ群がないように思います。
「dsp」や「digital filter」、「butterworth」などのキーワードで
検索してみましたが、レスのついていないMLの投稿にヒットするようです。
既存の関数(filterなど)で十分実現可能だから必要ないということなのか
それとも誰も開発していない(需要がない)ということなのか・・。
どなたか非公式パッケージでもご存じでしたら是非教えてください。
Krokodile (2005-06-29 (水) 16:48:38)
R では akima パッケージで、スプライン補間ができますが、
ベジェ曲線による補間はできるのでしょうか?
このサイトを検索しましたが、「ベジェ曲線」はヒットしませんでした。
田中 (2005-06-27 (月) 17:30:29)
以下の質問のご回答よろしくお願いいたします。
k×(m*i)行列Aを作成するときに、1×mの行列Bkm(k、mは添え字で、B行列の要素はk,m,iの関数として与えられる)を合成することで作成したい(つまり、行列Aの要素として行列Bを利用したい)のですが、centroid <- function(i,K,m) { C <- matrix(0,K,m*i) U <- matrix(0,1,m) for (a in 1:K) { for (b in 1:i) { for (w in 1:m) { U[a,b][1,w] <- (a,b,wの3変数による関数) } } } for (a in 1:K) { for (b in 1:i*m) { t <- ceiling(b/m) C[a,b] <- U[a,t][1,(b-(t-1)*m)] } } }としたのですが、エラーがでます。Uの定義の仕方に問題があるのであろうとは思うのですが、調べても(行列に2変数の添え字を付け、さらにその行列の要素を先の2変数+1変数で指定してやる方法について)分からなかったので、助けていただければと思います。
(分かりにくい気がするので、まとめると、C[i,j]=Uij Uij[1,k]=f(i,j,k)という場合の、Uの記述方法について教えていただきたいということです)
よろしくお願いいたします。
C[[i,j]][1,k]のようにアクセスすることでしょうか。-- MKR 2005-06-27 (月) 18:35:22
x <- array(0,c(3,3,3)) for (i in 1:3) for (j in 1:3) for (k in 1:3) x[i,j,k] <- f(i,j,k) x[i,j,k]
PTR【n】 <- matrix(0, n, n) 【】で表される要素はメモリ領域の先頭を指すポインタを表す(^_^;) PTR【2】 は 2×2行列を意味し,PTR【10】は 10×10行列。。こういうのがあったら便利だなと思うこともありますね(あるよ!という答えが返ってきたりして)。 -- 青木繁伸 2005-06-27 (月) 19:13:23
C <- as.list(rep(NA, 3*3)) dim(C) <- c(3,3) # リストを行列化 for (i in 1:3) for (j in 1:3) # リスト行列の各成分に(なんでも)ほうりこむ C[[i,j]] <- runif(i*j) C [,1] [,2] [,3] [1,] 0.6133985 Numeric,2 Numeric,3 [2,] Numeric,2 Numeric,4 Numeric,6 [3,] Numeric,3 Numeric,6 Numeric,9 C[[2,3]][4] # (2,3) 番目のブロック(ベクトル)の4番目の成分 [1] 0.7488104
centroid <- function(i,K,m){ C <- ???? for(a in 1:K){ for(b in 1:i){ for(w in 1:m){ C[[a,b]][w] <- (a,b,wの関数) } } } print(t(C)) }とプログラムを作りました。K×i個のベクトル(要素数m)からK×iのリスト行列(?)を出力したいのですが、Cの宣言の仕方がよく分かりませんでした。(いろいろ調べて試してみたのですが、エラーばかりでてしまいました。)これについて教えていただけませんでしょうか。 たびたびになり申し訳ないのですが、よろしくお願いいたします。-- 田中 2005-06-29 (水) 01:19:36
C[[a,b]][w]という構文は使えません。前者なら使えるはずです。ついでに始めて気づきましたが、リストは行列(風)にアクセスできるだけでなく、dim(C) <- c(3,3,3) 等とすれば配列(風)にアクセスすることもできるようです。なお、リストへのアクセスは純粋な行列・配列へのアクセスに比べれば遅めですから、特に高速性が要求されるときは問題有りかもしれません。-- MKR 2005-06-29 (水) 06:51:01
13m (2005-06-26 (日) 06:26:05)
ある行列もしくはベクトルの要素それぞれに対して、1から任意の整数aまでのx/nの総和を計算するf(x,a)=Σ[n=1:a](x/n)のような関数を適用したいと思い、 以下のようなものを書いてみました。f <- function(x, a) { res <- 0 for(n in 1:a) res <- res + x/n return( res ) }実際に使ってみた様子はこんな感じです
> m=matrix(c(1,2,3,4),ncol=2) > m [,1] [,2] [1,] 1 3 [2,] 2 4 > f(m,10) [,1] [,2] [1,] 2.928968 8.786905 [2,] 5.857937 11.715873これと同じ結果を返す関数を、for文(または要素毎のapply)を使わず、もっと簡潔に、R的に、スマートに書けますでしょうか?すぐにでも書けそうな気がするのですが、私の力不足のせいかなかなかいい方法が思いつきません。なんだかお題みたいな質問になってしまってすみません。
temp <- function(x,a) x*sum(1/(1:a)) temp(m,10) [,1] [,2] [1,] 2.928968 8.786905 [2,] 5.857937 11.715873
g <- function(x, a) { res <- 0 for(n in 1:a) res <- res + 1/(x+n) return( res ) }実行結果
> m=matrix(1:4,ncol=2) > g(m,5) [,1] [,2] [1,] 1.450000 0.8845238 [2,] 1.092857 0.7456349これと同じ結果を返すものを、いかにもfor文を使わずに書けそうな気がするので、もしあれば教えていただきたいと思います。 -- 13m 2005-06-26 (日) 13:13:57
> g <- function(x, n) + { + ncol <- ncol(x) + nrow <- nrow(x) + k <- array(rep(1:n, each=ncol*nrow), dim=c(nrow, ncol, n)) + x <- array(x, dim=c(nrow, ncol, n)) + apply(1/(x+k), c(1,2), sum) + } > m <- matrix(1:4, 2) > g(m, 5) [,1] [,2] [1,] 1.450000 0.8845238 [2,] 1.092857 0.7456349我ながらわかりにくいと思います。for を使ったプログラムの方が遙かに分かりやすい。 -- 青木繁伸 2005-06-26 (日) 20:20:33
> matrix(colSums(1/mapply(seq, m+1, m+5)), dim(m)) [,1] [,2] [1,] 1.450000 0.8845238 [2,] 1.092857 0.7456349
> h <- function(x,n) matrix(colSums(sapply( x, function(x) 1/( x+seq(n) ))), dim(x)) > h(m,5) [,1] [,2] [1,] 1.450000 0.8845238 [2,] 1.092857 0.7456349上のやり方だとかなりアレなので、きちんと別途関数定義して、x,nに関するある程度どんな数列でも扱えるよう工夫してみました。ついでに引数はベクトルにも行列にも対応。
eess <- function(FUN,x,n) { # 引数xのそれぞれの要素に対して一般項FUNに関する1からn項までの総和を求める res <- colSums(sapply(x, FUN, seq(n))) if(is.matrix(x)) res <- matrix( res, dim(x) ) return( res ) } feess <- function(FUN,x,n) { # 上記と全く同様の機能をfor文で書き直したもの FUN <- match.fun(FUN) res <- 0 for(i in seq(n)) res <- res + FUN(x,i) return(res) }使用例とsystem.time()
> gn <- function(x,k) 1/(x+k) # 一般項を適当な関数で定義 > eess(gn,m,5) [,1] [,2] [1,] 1.450000 0.8845238 [2,] 1.092857 0.7456349 > eess(gn,1:5,5) [1] 1.4500000 1.0928571 0.8845238 0.7456349 0.6456349 > eess("^",m,5) # 二項演算子を渡して奇抜な使い方も(Σ[k=1:n]x^k) [,1] [,2] [1,] 5 363 [2,] 62 1364 > system.time(eess(gn,rnorm(10000),100)) [1] 0.51 0.04 0.59 NA NA > system.time(feess(gn,rnorm(10000),100)) [1] 0.06 0.02 0.08 NA NAsapply()の演算は思ったより高速なんですね。驚きました。速度的にはそれほどでなくとも(なんだかんだで最初のfor文が一番速い)、いつか何か別の用途に使えるかもしれません。可読性があるかと言われればあまりないかもしれませんが。 -- 13m 2005-06-27 (月) 02:25:07
坂田 (2005-06-24 (金) 22:05:29)
主成分得点の平面に固体を名前で表示するよりも顔写真を貼り付けたい
と思います(顔に限りませんけど)。このための方法をご教授ください。
巷の本ではプロ野球選手の成績に対して主成分分析を行っているものが
見受けられますが、名前ではなく顔を張りつければインパクトも大きいと思
うのですが。
> library(pixmap) > r <- read.pnm(system.file("pictures/logo.ppm", package = "pixmap")[1]) > x <- runif(10) > y <- runif(10) > plot(x,y) > for ( i in 1:length(x)) addlogo(r,px=c(x[i]-0.05,x[i]+0.05),py=c(y[i]-0.05,y[i]+0.05))
#ref(): File not found: "addlogo.png" at page "初級Q&A アーカイブ(3)"
苺 (2005-06-22 (水) 21:22:30)
はじめまして。
あるデータの平均と、そのなかの特定の1つのデータの図を同じ図にplotするのにpointを使用するときいたのですが、どうすればいいのでしょうか???
> x <- rnorm(100) # 100個の正規乱数を発生させて > plot(1:2, c(min(x), max(x)), type="n") # プロット領域を定義して > points(1.5, mean(x)) # 平均値を描いて > points(1.5, x[3], col="red") # 3番目のデータを赤でプロットしてなんて馬鹿馬鹿しいことをやるわけじゃないんでしょう?
? (2005-06-22 (水) 18:35:56)
重回帰分析をする際にダミー変数を用いました。
しかし、実測値と推定値をplotする際にどうしてもエラーになっていしまいます。どのようにして表現するのでしょうか?
ちなみにダミー変数は数字からアルファベットに変えてやりました。
> x1 <- factor(c("foo", "bar", "bar", "foo", "baz", "baz")) > x2 <- factor(c("hoge", "hoge", "wao", "wao", "hoge", "hoge")) > y <- c( 2, 1, 2, 4, 3, 4) > ans <- lm(y ~ x1+x2) > cbind(y, predict(ans)) y 1 2 2.25 2 1 0.75 3 2 2.25 4 4 3.75 5 3 3.50 6 4 3.50 > plot(y, predict(ans))
#ref(): File not found: "regplot.png" at page "初級Q&A アーカイブ(3)"
というようなことになったが? -- 青木繁伸 2005-06-22 (水) 20:15:21森野 (2005-06-17 (金) 19:40:38)
初歩的な質問で,申し訳ありません。
CRANからRMySQLをダウンロードして,RGuiのパッケージから(ローカルにあるZIPファイルからのパッケージのインストール)を実行しました。
library()で確認したところ,** No title available (pre-2.0.0 install?) **と表示されてしまい,使用する事ができません。
WindowsXPサービスパック1,R 2.1.0を使用しているのですが,どなたかご教授頂けないでしょうか。
山下 (2005-06-15 (水) 17:52:32)
WindowsのVimでRを使用したいのですが、
皆様使用されている方はいませんでしょうか?
エディタ版でVimファイルは見つけたのですが、Syntaxも上手くいかず、
デフォルトの :set Syntax=rを使用しています。
xyzzy並みに使っておられる方がいましたらご教授下さい。
どうむ (2005-06-15 (水) 12:13:12)
LinuxでRを使用しております。
ATLASを使用したいのですが、どのようにしたらいいのでしょうか、、、
以下「うまくいかない」現象です
- おこなったこと(ATLAS)
- http://math-atlas.sourceforge.net/からソースをダウンロード、解凍
- 一般ユーザーで
make- rootユーザーで
make install arch=Linux_HAMMER64SSE2_2
- CPUは、に上からもわかるとおり、AMD64bitです(汗
- 作成された「.a」をパスの通っているディレクトリへ、コピー
- おこなったこと(R)
- バージョン > R-2.0.1
- ソースをダウンロード、解凍
- 一般ユーザーで
./configure --with-blas=f77blas
- この時点でログには
checking for sgemm_ in -lf77blas... no checking for sgemm_... no checking for ATL_xerbla in -latlas... no checking for sgemm_ in -lblas... yes- f77blasを使用することにはなっていません、、、
- 一般ユーザーで
make make check- rootユーザーで
make install
- 案の定lddで実行ファイルを見ても普通のlibblas.soがリンクされている
これを、解決するには
p.s.
Win版ATLASはdllをダウンロードしてきて、既存のものと置き換えるだけでしたので簡単でした(汗
行列演算が、計算の中身にもよりますが、4分の1になったりして、驚いています
Linuxでもこのパフォーマンスを体験したい、、、
make install arch=Linux_UNKNOWNSSE2
./configure --with-blas="-L/usr/lib -lf77blas -latlas" make make check make install
checking for sgemm_ in -L/usr/lib -lf77blas -latlas... yes
test2 <- function ( n=500 ) { A<-array(rnorm(n^2), dim=c(n,n)) B<-array(rnorm(n^2), dim=c(n,n)) C<-array(rnorm(n^2), dim=c(n,n)) D<-array(rnorm(n^2), dim=c(n,n)) BA <- B%*%A return (system.time(A%*%solve(t(BA)%*%BA+C)%*%BA%*%D )) }
> test2() [1] 2.24 0.08 2.33 0.00 0.00 > test2() [1] 2.26 0.06 2.33 0.00 0.00 > test2() [1] 2.27 0.05 2.32 0.00 0.00
> test2() [1] 0.46 0.07 0.54 0.00 0.00 > test2() [1] 0.45 0.06 0.51 0.00 0.00 > test2() [1] 0.43 0.05 0.48 0.00 0.00
make install arch=Linux_UNKNOWNSSE2
# ATLAS make xconfig ./xconfig -F f '-fomit-frame-pointer -O -m64 -fPIC' ? -F c '-fomit-frame-pointer -O -mfpmath=387 -m64 -fPIC' ? -F m '-fomit-frame-pointer -O -mfpmath=387 -m64 -fPIC' make install arch=Linux_HAMMER64SSE2_2 # rootである必要はありません # R (ATLASのスレッドとR-profilingは不仲!) ./configure --disable-R-profiling --with-blas='-L/usr/local/lib -lptf77blas -latlas -lpthread'
> test2() [1] 0.76 0.05 0.57 0.00 0.00 > test2() [1] 0.75 0.07 0.58 0.00 0.00 > test2() [1] 0.75 0.06 0.56 0.00 0.00
学生 (2005-06-14 (火) 18:45:06)
はじめまして、とっても初歩的な質問ですが、
R環境で、入力データ(数値)をソートしたいです、
できれば、出力は順番づけられてほしい
ご存知の方、教えていただければありがたいです。
よろしくお願いいたします
help(sort) # または簡略化された ? sortでもよい今後は自分で調べてもらうことにして,sort の使い方は以下の通り。
> x <- c(3,6,5,2,8) > sort(x) [1] 2 3 5 6 8Net で得られる R の使い方についての資料は,舟尾さんの R-tipsが優れていますが,本を買ってあげてください。
書名(著者・編者等) The R Tips データ解析環境Rの基本技・グラフィック活用集(舟尾暢男著) 発行元 株式会社 九天社 発行年 2005 価格 3500円(税別) ISBN 4-86167-039-xです。 -- 青木繁伸 2005-06-15 (水) 11:27:17
小島 (2005-06-14 (火) 17:40:44)
はじめまして。
heatmap関数を使用して作図をしているのですが、行き詰ってしまいました。
以下の事柄についてご存知の方がいましたら解決方法を教え下さい。
質問1:heatmapの縦横比の変更についてheatmap関数を用いて解析を行うと、その縦横比(X軸とY軸の長さ)が必ず1:1の正方形になるのですが、この比率は変更できないのでしょうか?
扱っているデータがX軸の数の方が圧倒的に多いので、X軸を長くした長方形の図を得たいと思っています。
通常、出力デバイス(png, pdf等)の大きさを変更するとそれに合わせて図も伸縮してくれるのですが、heatmap関数の場合この比率が変更されずに余白が出来てしまいます。
例を挙げると以下のようになります。#サンプルデータ(データをX軸の方が多くなるように加工) data(mtcars) sampledata<-rbind(c(mtcars[,2]), c(mtcars[,3]), c(mtcars[,4]), c(mtcars[,5])) #デバイスの大きさを横長の長方形に指定する pdf("test.pdf", width=12, height=5) heatmap(sampledata) dev.off()marginsをいじれば形は変化するのですが、余白が出来る上、xlabの値が残ってしまいます。
heatmap(sampledata, xlab="XLAB", margins=c(20,5))
質問2:heatmapを使った複数の図の描画について
複数の図を描画する時、
#layoutを用いた画面分割(下図を大きく) mat <- matrix(c(1,2,2), 3, 1, byrow=TRUE) layout(mat) plot(sin) plot(cos)など、plot関数を用いた場合は描画できるのですが、
heatmap関数を用いると、heatmapの図だけ独立して新規に出来てしまいます。mat <- matrix(c(1,2,2), 3, 1, byrow = TRUE) layout(mat) #データは上記と同じものを使用 plot(sampledata) heatmap(sampledata)何か回避する方法はあるのでしょうか?
ヘルプを見たり検索をしてみたのですが、私自身の力不足とheatmapに関する記述が全体的に少なかったため方法が見つかりませんでした。
ご存知の方がいらっしゃいましたらよろしくお願いいたします。
layout(lmat, widths = lwid, heights = lhei, respect = TRUE)という行の,respect を FALSE にします。つまり,
layout(lmat, widths = lwid, heights = lhei, respect = FALSE)とすれば,いいようですが。そして当然使うのは heatmap2(sampledata) のように。以下のようになります(デバイスの大きさを示すため,余分な背景も切り取りました)。
#ref(): File not found: "heatmap.png" at page "初級Q&A アーカイブ(3)"
自己責任で試してみてください。# heatmapでは時間軸の固定が可能 heatmap(sampledata, Colv=NA) # heatmap.2ではColvに欠損値(NA)は入れられないとエラーが出る heatmap.2(sampledata, Colv=NA, trace="none", key=FALSE)以下のように dendrogram を指定すると似たようなグラフは出てくるのですが順序がバラバラになってしまいます。
heatmap.2(sampledata, dendrogram="row", trace="none", key=FALSE)ヘルプの Note: 欄に書いてあった時間軸固定に関する記述がなくなっているので、heatmap.2 ではこのような設定はできないのでしょうか?-- 小島 2005-06-15 (水) 15:35:47
heatmap2 <- heatmap fix(heatmap2) ここで,前記の通りの行を修正後,保存し編集ウインドウを閉じる heatmap2(sampledata)これで,うまく行くと思うのだが。。。 -- 青木繁伸 2005-06-15 (水) 20:35:42
ショーゾー (2005-06-06 (月) 13:20:11)
一つのカテゴリー変数(FとM)と5つの数値変数を説明変数としてyの予測式を重回帰分析(ステップワイズ)で求めようとしています。性別のカテゴリー変数(sex)はfactorとして読み込み,str()でsex : Factor w/ 2 levels "F","M": 1 1 1 1 2 2 ... と確認できました。重回帰分析のsummary()ではこのsexの係数(estimate)が,sexM 6.3007という表示(有意でした)がされました。
ここで,予測値を計算する時はsex=Mのときは6.3007 * 2としてもちいるのでしょうか? sex=Fのときは6.3007 * 1でいいのでしょうか?
チェックのために,Fのときは0,Mのときは1のダミー変数(sexcode)をつくって実行したところ,おなじ結果がでました。このときの予測値の計算の時はsex=Mのときは6.3007 * 1,sex=Fのときは6.3007 * 0となるように思います。すると,factorの内部表現はどうなっているのでしょうか?
> x <- c("A", "B", "C") > as.integer(as.factor(x)) [1] 1 2 3 > x <- c("C", "B", "A") > as.integer(as.factor(x)) [1] 3 2 1という結果を精査すべし。 -- 青木繁伸 2005-06-06 (月) 13:37:25
> x <- as.factor(sample(3, 10, replace=T)) > x [1] 1 3 2 1 1 3 1 3 1 1 Levels: 1 2 3 > y <- rnorm(10) > lm(y ~ x) Call: lm(formula = y ~ x) Coefficients: (Intercept) x2 x3 # x2,x3 は何かといえば,,, -0.04008 0.80732 -0.18625 > x2 <- ifelse(x ==2, 1, 0) # こういう風に解釈してくれている > x3 <- ifelse(x ==3, 1, 0) > lm(y ~ x2+x3) Call: lm(formula = y ~ x2 + x3) Coefficients: (Intercept) x2 x3 -0.04008 0.80732 -0.18625 # で,同じ答えになるこんな風になる。賢いなぁ R は。 -- 青木繁伸 2005-06-07 (火) 00:32:00
> c <- c("Low", "Medium", "High") > f <- factor(c) > f [1] Low Medium High Levels: High Low Medium > as.integer(f) [1] 2 3 1 > f2 <- factor(c, levels=c("Low", "Medium", "High")) > f2 [1] Low Medium High Levels: Low Medium High > as.integer(f2) [1] 1 2 3 >
g (2005-06-06 (月) 10:48:48)
R commanderを使いたくて,パッケージRcmdrを読み込んだら
以下のメッセージが出てしまいました.
同様な状況を解決した人がいましたら解決方法を教えてください.
ちなみにRのバージョンはR2.1.0で
R commandeのバージョンはRcmdr1.0-2です.
また,以下のパッケージも読み込んでいます.
(バージョンもチェック済み)
abind, car (>= 1.0-15), effects (>= 1.0-7), foreign,
grid, lattice, lmtest, MASS, mgcv, multcomp,
mvtnorm, nlme, nnet, relimp, sandwich, strucchange, zoo
*****エラーメッセージの始まり*****local({pkg <- select.list(sort(.packages(all.available = TRUE)))
- if(nchar(pkg)) library(pkg, character.only=TRUE)})
要求されたパッケージ tcltk をロード中です
要求されたパッケージ rgl をロード中です
要求されたパッケージ zoo をロード中です
要求されたパッケージ strucchange をロード中です
要求されたパッケージ sandwich をロード中です
要求されたパッケージ relimp をロード中です
要求されたパッケージ nnet をロード中です
要求されたパッケージ nlme をロード中です
要求されたパッケージ mvtnorm をロード中です
要求されたパッケージ multcomp をロード中です
要求されたパッケージ mgcv をロード中です
This is mgcv 1.3-0
要求されたパッケージ MASS をロード中です
要求されたパッケージ lmtest をロード中です
要求されたパッケージ lattice をロード中です
要求されたパッケージ grid をロード中です
要求されたパッケージ foreign をロード中です
要求されたパッケージ effects をロード中です
要求されたパッケージ car をロード中です
要求されたパッケージ abind をロード中です
以下にエラーparse(file, n, text, prompt) : 構文解析エラーです
エラー:.onAttach は 'attachNamespace' で失敗しました
エラー:'Rcmdr' に対するパッケージもしくは名前空間のロードが失敗しました
ken (2005-06-03 (金) 23:03:02)
はじめまして
plot関数を使って散布図を作成しているのですが、x軸、y軸に対応した項目を直接図の中に表示するにはどうすればよいでしょうか? 例えば、項目 X Y A 3 4 B 7 5 # 投稿法を良く読んで投稿しなさいのような表の散布図を作成した場合、A、Bを直接図の中に表したいと思っています。一通り検索してみましたが、いい方法が見つからなかったので、質問させていただきました。 よろしくお願いします。
> dat <- matrix(c(3,4,7,5), nrow = 2, byrow=TRUE, dimnames = list(c("A", "B"), c("x", "y"))) > plot(dat, pch=rownames(dat))松村俊和 (2005-06-04 (土) 00:33:00)
d <- read.table("test.dat") plot(d, pch=" ") text(d, rownames(d))プロット位置にマークを書いたり,text で書かれる文字列がマークと重ならないようにするにはどうしたらよいかは,オンラインヘルプを参照すればよいでしょう。 -- 青木繁伸 2005-06-04 (土) 21:33:35
Mari (2005-06-02 (木) 13:14:43)
以下のように乱数を発生させ、その平均値の分布のSDがσ/sqrt(n) になることを確かめたかったのですが、なぜかσ/sqrt(n)にはなりません・・。
途中で計算のさせかたを間違っているのかと思いましたが、どうしても分かりません。> dat <- numeric(10000) > for ( i in 1:10000) dat[i] <- mean(rnorm(100,1,3)) > sd(dat) [1] 0.30414 > 3/sqrt(10000) [1] 0.03どなたかご教示いただけませんでしょうか・・・
> x <- matrix(rnorm(1000000,1,3), 10000,100) > sd(apply(x, 1, mean)) [1] 0.29930224896266783ですね。 -- 青木繁伸 2005-06-02 (木) 13:34:23
<ふ> (2005-06-02 (木) 00:04:40)
R2.0.1を使ってます。x <- c(150,160,170,180,190)というデータをつくり、
> hist(x)でグラフをかかせますと、150と160が同じ区間にカウントされます。
help(hist)を読みましたら、
'right = TRUE' (default), では、'(a, b]'で値をカウントするけど、左端の区間は、別、と書いてありました。
それなら、と、> hist(x,xlim=c(140,200))で描画させてみたのですが、左端とは、データの左端なんですね。
次に、> hist(x,right = FALSE)としますと、180と190が同じ区間にカウントされます。
手でグラフを描いたら、こんなふうにはせずに、平らな図にします。このように、最小値なり最大値の扱いを他と変える理由はどこにあるのでしょうか。
よろしくご教示ください。> stem(x) The decimal point is 1 digit(s) to the right of the | 15 | 0 16 | 0 17 | 0 18 | 0 19 | 0 > hist(x,xlim=c(140,200)) > hist(x,xlim=c(140,200),right = FALSE)
hist(x,breaks=15:20*10,right=FALSE)とすれば良いのでは?
Note that xlim is not used to define the histogram (breaks), but only for plotting (when plot = TRUE).しかしまぁ,エラーまでペーストする必要はないでしょう。
RIRISU (2005-05-31 (火) 01:52:02)
クルスカル・ワリスの検定についてお聞きします。1回群 2回群 4回群 8回群 no.1 9 13 19 24 no.2 11 15 19 23 no.3 11 14 21 20 no.4 14 16 22 24 no.5 20 17 19 19 no.6 13 21 20 24このデータから群によって評定の差があるかどうかを検定したいのですが、
> x<-c(9,11,11,14,20,13) > y<-c(13,15,14,16,17,21) > z<-c(19,19,21,22,19,20) > e<-c(24,23,20,24,19,24) > x<-c(x,y,z,e) > g<-factor(rep(1:4,c(6,6,6,6)), + labels=c("1group","2group","4group","8group")) > kruskal.test(x,g)で合っているでしょうか?手計算(表計算ソフトで公式に数値を当てはめるやり方)でやると、カイ二乗値は1494.129となり、自由度は24-1で23になります。
手計算でやると、0.1%の危険率で帰無仮説を棄却するのですが、うえの実行結果ではカイ二乗値=14.4297で自由度=3となり、1%の危険率で帰無仮説を棄却する。となり結果がちがってしまいます。
Sybock (2005-05-26 (木) 18:50:12)
はじめまして、質問させていただきます。
anova を3群のデータに適用し、有意差が出たものについて post-hoc analysis を行いたいと思っております。
中澤港様の「Rによる統計解析の基礎」P.109 を参考に、# 例 firingRatio <- c(52,30,35,41,29,…, 33) types <- c("A", "C", "A", "B", …, "C") # 両ベクトルとも同じ length pairwise.t.test(firingRatio, types, p.adjust.method = "holm")といった形で実行しましたが、この pairwise.t.test という関数の"pairwise"は、t.test(paired = TRUE)とは意味が違うものなのでしょうか。
また、古典的な Fisher の制約つき LSD 法 も行って見たいのですが、RjpWiki 内で Fisher, LSD 等で単語検索を行ったところ、それらしいものを発見できず、yahoo 等で検索しても R における実現方法を探り当てることができませんでした。
知識不足のため、検索キーワードが不適切であった可能性もあり、申し訳ありませんが、上記の件につきましてアドバイスをよろしくお願い致します。
さかな (2005-05-25 (水) 21:35:34)
質問させて頂きます。
例えば、開始日(2002/05/03)と終了日(2004/03/21)を指定した時に、
期間中の一連の日付を生成する良い方法はないでしょうか?
2002/05/02
2002/05/03
....
2004/03/20
2004/03/21
アドバイスよろしくお願いします。
library(chron) seq.dates("05/03/2002", "03/21/2004")。。。-- 青木繁伸 2005-05-25 (水) 23:12:54
x <- seq(as.POSIXct("2002/05/03"),as.POSIXct("2004/03/21"),by = "day") gsub("-","/",as.character(x))
こじろう (2005-05-12 (木) 01:21:32)
Rのカーソルが|ではなく■の形になっていてとても使いにくいです。|の形にしたいのですが、どなたかご教示いただけませんか??
Akira (2005-05-04 (水) 0:01:12)
LaTeXを勉強しろと言われそうですが、グラフと表を一枚のPDFファイルで出力する方法はあるのでしょうか?
LaTeX、Rによるポストスクリプト画像のLaTeXでの利用、mimetexを確認しました。理解不足もありますが、少し違うように思います。
非効率的とはわかりながらも、今は座標を探しながらmtextで表を作成しています。pdf(file="test.pdf") layout(c(1:2)) plot(x=0, y=0, xlim=range(0,4), ylim=range(0,5), pch=".", main="test", xaxs="i", yaxs="i", axes=FALSE, xlab="", ylab="", col="white") for(i in seq(0, 6)){ segments(0.5, i, 14.5, i) } mtext(text="column1", side=3, line = -1, adj=0, at=2, font=2) mtext(text="column2", side=3, line = -1, adj=0, at=3, font=2) mtext(text="row1", side=3, line = -2.2, adj=0, at=1, font=3) mtext(text="row2", side=3, line = -3.4, adj=0, at=1, font=3) mtext(text="row3", side=3, line = -4.5, adj=0, at=1, font=3) mtext(text="row4", side=3, line = -5.7, adj=0, at=1, font=3) mtext(text="data1.1", side=3, line = -2.2, adj=0, at=2, font=1) mtext(text="data2.1", side=3, line = -3.4, adj=0, at=2, font=1) mtext(text="data3.1", side=3, line = -4.5, adj=0, at=2, font=1) mtext(text="data4.1", side=3, line = -5.7, adj=0, at=2, font=1) mtext(text="data1.2", side=3, line = -2.2, adj=0, at=3, font=1) mtext(text="data2.2", side=3, line = -3.4, adj=0, at=3, font=1) mtext(text="data3.2", side=3, line = -4.5, adj=0, at=3, font=1) mtext(text="data4.2", side=3, line = -5.7, adj=0, at=3, font=1) plot(1:10) dev.off()
?documentclass[a4paper]{article} ?SweaveOpts{echo=FALSE} ?usepackage{a4wide} ?usepackage{Sweave} ?begin{document} <<>>= x <- 1:20 y <- x/10+rnorm(20) summary(z <- lm(y ~ x)) @ ?SweaveOpts{echo=true} ?begin{figure}[htbp] ?begin{center} <<fig=TRUE>>= par(mfrow=c(2,2)) plot(z) @ ?caption{linear regression} ?end{center} ?end{figure} ?end{document}
#ref(): File not found: "sweavetest.png" at page "初級Q&A アーカイブ(3)"
?documentclass[a4paper]{article} ?usepackage{graphicx} ?begin{document} ?SweaveOpts{echo=false,pdf=FALSE} <<results=hide>>= library(xtable) hoge<-data.frame(Tokyo=c(100,200), Osaka=c(80,69), Sapporo=c(80,69), Fukuoka=c(80,69), Naha=c(80,69), Kyoto=c(80,69), Hiroshima=c(80,69), Nagoya=c(90,250) ) @ ?section{table and graph} <<results=tex>>= xtable(hoge, caption="table dayo") @ ?setkeys{Gin}{width=1.0?textwidth} <<fig=TRUE,eps=TRUE>>= boxplot(hoge) @ ?end{document}
#ref(): File not found: "test1.png" at page "初級Q&A アーカイブ(3)"
#ref(): File not found: "temp.png" at page "初級Q&A アーカイブ(3)"
原 (2005-05-03 (火) 11:52:12)
東京の地図をmaptoolsで表示させた上に、地価のポイントを表示させたいのですが、やり方がわかりません。
東京の地図は表示できたのですが、地価ポイントの緯度・経度のエクセルファイルのみ(シェープファイルなし)で、東京の上に表示させることはできないのでしょうか。
初心者 (2005-05-02 (月) 22:14:17)
Windows 版 R 2.1.0 ならば,East Asian Languages でインストールした場合に Rconsole と Rdevga を修正する必要がありますが,Mac 版 R 2.1.0 はそのような作業を行う必要は無いのでしょうか?
さとう (2005-05-02 (月) 20:01:21)
debian(woody)においてR-2.0.1をapt-getするためのsources.listファイルの記述はどのようにすればいいのでしょうか。CRANやそのmirrorのサーバーを指定すると,
Ign http://******* Release
のようになりapt-getできません。
# Stable deb http://ftp.de.debian.org/pub/debian stable main contrib non-free deb http://ftp.de.debian.org/pub/debian-non-US stable/non-US main contrib non-free # Testing deb http://ftp.de.debian.org/pub/debian testing main contrib non-free deb http://ftp.de.debian.org/pub/debian-non-US testing/non-US main contrib non-free # Unstable deb http://ftp.de.debian.org/debian unstable main contrib non-free deb http://ftp.de.debian.org/debian-non-US unstable/non-US main contrib non-free
deb http://cran.md.tsukuba.ac.jp/bin/linux/debian stable main
http://cran.au.r-project.org/bin/linux/debian stable/これでも駄目なら deb ファイルをダウンロードし dpkg -i r-base 等と手作業でインストールするか、ソースからコンパイルしたほうがてっとり早いですね。ところで woody 版に固執する理由があるのでしょうか。