初心者のための R および RjpWiki に関する質問コーナー
新規投稿はできません
アール (2006-01-25 (水) 14:17:44)
こんなことをする人はいないので,不具合でも何でもないが,ちょっと焦ったので。
outer 関数で引用する関数が返す値は,引数が絡んでいないとエラーになる。具体的には,定数を返そうとしたらエラーになった。そのようなときには,だましてやるとよいようだ。たまたま0を返す関数を作っていたので,x*0 とかにすればオーケーだった。> x <- y <- c(-1,0,1) > fun <- function(x, y) return(0) > outer(x, y, fun) 以下にエラーouter(x, y, fun) : dim<- : dims [product 9] は object [1] の長さに整合しません > fun2 <- function(x, y) {z <- x+y; return(0)} > outer(x, y, fun2) 以下にエラーouter(x, y, fun2) : dim<- : dims [product 9] は object [1] の長さに整合しません > fun3 <- function(x, y) return(x*0) > outer(x, y, fun3) [,1] [,2] [,3] [1,] 0 0 0 [2,] 0 0 0 [3,] 0 0 0
ying (2006-01-24 (火) 15:50:30)
多変数関数 f(x1,x2,x3,x4)=a*x1+b*x2+c*x3+d*x4 (a,b,c,d:定数)を二次元画面で,x1とx2(任意の二つ変数)を横と縦軸として,グラフを描きたいですけど,どうすればいいですか?Rで描けますか?
#ref(): File not found: "mochi-yaki-ami.png" at page "初級Q&A アーカイブ(4)"
黒く表されているのが結果として描かれる平面。これは,x,y平面に平行で,z 軸と 0 で交わる平面。persp 関数はグリッドを描くが,縮小したので黒くつぶれた。y <- x <- seq(-1,1,length=51) f <- function(x, y) {(x+y)*0} z <- outer(x, y, f) persp(x, y, z, zlim=c(-0.5, 0.5))のようになります。f の定義で単に0を返すように定義すれば良いところを,おかしなことをやっている理由については別件で。 -- 2006-01-25 (水) 13:54:50
ExperimentalDesign (2006-01-24 (火) 14:29:09)
直交表実験の結果は次のとおりです。xが特性値でA〜Hまでが2水準の因子です。Eの因子はありません。交互作用はA:B,A:C,B:C,C:D,D:Fのみ考慮しました。DATA <- structure(list(x = c(16, 38, 41, 25, 38, 58, 43, 45, 19, 27, 22, 38, 5, 41, 32, 28), A = c(-1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1), B = c(-1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1), C = c(-1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1), D = c(-1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1), G = c(-1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, 1, -1), H = c(-1, 1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, 1, -1), F = c(-1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1)), .Names = c("x", "A", "B", "C", "D", "G", "H", "F"), row.names = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16"), class = "data.frame")実験の結果から分散分析を次のようにしたところ
summary(aov(x~.+(A+B+C)^2+C:D+D:F,DATA))G,H,A:C,B:C,D:Fの効果が小さいことがわかったのでこれを無視することにしました。
この結果から、A1B1C1D1F1水準(L1)とA1B2C1D2F2水準(L2)での母平均の推定を行うためにmod <- lm(x~.-G-H+A:B+C:D,DATA) L1<-predict(mod,data.frame(A=-1,B=-1,C=-1,D=-1,F=-1,G=0,H=0), level=.95,interval="confidence",se.fit=T) L2<-predict(mod,data.frame(A=-1,B=1,C=-1,D=1,F=1,G=0,H=0), level=.95,interval="confidence",se.fit=T)として、望む結果を得られました。
この結果を受け、L1とL2の差がどれくらいあるのか推定したいのですがどのようにすれば推定できるのでしょうか?望む結果は fit = -46.0 , lwr = -60.7 , upr = -31.3 です。
> (L1$fit-L2$fit)[1]+c(-1,1)*qnorm(0.995)*L1$residual.scale [1] -60.68449 -31.31551かなぁ? -- 2006-01-24 (火) 14:57:34
> (L1$fit-L2$fit)[1]+c(-1,1)*qnorm(0.995)*L1$residual.scaleで信頼区間を得ることができたので、そういう意味で「できました」と書きました。何かおかしかったでしょうか? -- ExperimentalDesign 2006-01-24 (火) 18:47:14
mitsu5 (2006-01-20 (金) 16:55:33)
きわめて初歩的なことですが、よろしくお願いいたします。
boxplotと散布図を合成したいと思いました。a1=rnorm(30) a2=rnorm(30) a3=rnorm(30) x=c(rep(1,30),rep(2,30),rep(3,30)) boxplot(a1,a2,a3) par(new=T) plot(x,c(a1,a2,a3),xlim=par("usr")[1:2],axes=F)こうすればかなり良い線は行くのですが、X軸の1と3で微妙にずれてしまいます。Y軸はこのままで一致しているようです。しかしこれでは駄目かもしれません。
軸の一致の調整はどうすればよいのでしょうか。よろしくご指導ください。
a1=rnorm(30) a2=rnorm(30) a3=rnorm(30) x=c(rep(1,30),rep(2,30),rep(3,30)) boxplot(a1,a2,a3) #par(new=T) points(x,c(a1,a2,a3))
#ref(): File not found: "boxplot.png" at page "初級Q&A アーカイブ(4)"
今回はこちらで直しておきますが,今後はヘルプを読んで,きれいに投稿しましょう。 -- 2006-01-20 (金) 17:28:28にゃあ (2006-01-19 (木) 14:32:36)
Rでx^2−2x−3=0という方程式を解こうと思います。
Rではpolyroot関数を用いて、polyroot(c(1,-2,-3))で解けると思うのですが、出力が
> polyroot(c(1,-2,-3)) [1] 0.3333333+0i -1.0000000+0iとなって、解がおかしいような気がします。
どなたか検算をお願いします。
> polyroot(c(-3,-2,1)) [1] -1+0i 3-0i
ホームチーム (2006-01-19 (木) 10:15:26)
現在、WindowsとUnixでRのプログラミングを行っております。
バージョンは双方共にR2.0.1です。
Windowsでは、テキストエディターでプログラムを書いて、RのGUIをダブルクリックし、コピペで動作確認をしています。
Unixでは、Windowsで確認したプログラムを、FTPで設置し、R CMD BATCHで実行しています。
プログラムの中にwrite.table(summary_matrix,file="hoge.txt",sep="?t",col.names=NA)という部分があります。 プログラム自身は、自作関数による処理、クラスター解析や主成分分析など色々な処理を行い、途中結果も多数ファイルとして出力し、その最後にサマリーとして該当hoge.txtを出力して終了するような流れになっています。
開発環境であるWindowsXPでは、該当プログラムは全て問題なく動作し、サマリーもテキストファイルとして正しく出力されますが、Unix(ソラリス8)ではなぜか、サマリーとして出したいデータフレームのrow.names情報が、きちんと名前を正しく定義しているにもかかわらず、デフォルトの"1","2","3",,,,という行番号で出力されてしまうバグ?に苛まれております。
UNIXで動作させたときに本当に正しくrow.namesがセットできているか、write.table直前で、cat(row.names(summary.matrix))をしてみても、Win,Unix共に正しく情報はセット出来ていることを確認しています。
本現象について、お気づきの点等ありましたらご教授いただきたく考えております。
やたま (2006-01-19 (木) 10:07:08)
現在Windows環境でR-2.0.1を使っております。
ただ、現状最新のVersionは2.2.1かと思いますので、バージョンのアップを考えています。
そこで、R2.0.1→2.2.1に至るまでは数多くのバージョンアップがあったかと思いますが、その際数々のバグフィックスが等々あったと思いますが、その履歴等を追うにはどこのページを見ればよろしいでしょうか?
Rをインストールしたフォルダ以下にある「CHANGES」テキストファイルを見るしか方法はないでしょうか?
稚拙な質問で恐縮ですが、アドバイスいただければ幸いです。
ying (2006-01-17 (火) 17:57:37)
今svmについて理解している.でも,分からない所があって,教えていただけませんか?
help(svm)からsvmについてたくさん書いてるけど,実際svmのprogramを知りたいです.どうすればいいですか?
実際このようにやってみたら:> edit(svm) function (x, ...) UseMethod("svm") <environment: namespace:e1071>これしかなかった.もっと詳しいプログラムについて知りたいので,教えてください.
> library(e1071) > e1071:::svm.default > e1071:::svm.formula? svm をして,ヘルプおよび,そこからリンクされているウエッブページの内容をよく読むことの方が先決ではないでしょうか。-- 2006-01-17 (火) 18:08:29
しゃけ? (2006-01-16 (月) 15:46:54)
下記のような方法を見つけたのですが、一行目から何をやっているのか全く意味不明です。申し訳ありませんが、誰か詳しく解説していただけないでしょうか?ULONG_MAX = 4294967295以下ばっさり,削除(それにしても,投稿法くらい守ってね)
n<-100000;sum((colSums(matrix(runif(n*2)^2,2))<1))*4/n位でできるしィ -- 2006-01-16 (月) 17:23:28
n<-100000;sum(runif(n)^2+runif(n)^2<1)*4/nもっともっと? -- 2006-01-17 (火) 07:57:38
ying (2006-01-14 (土) 16:42:59)
データirisについてサポートして,グラフを描きたいけど、多次元のデータだから、二次元のグラフを描こうと思う。どうすればいいですか?それと、Rの中にplot.svmについての例:plot(m2, iris, Petal.Width ~ Petal.Length, slice = list(Sepal.Width = 3, Sepal.Length = 4))がある。これはいったいどんな意味ですか?
教えてください!!
## a simple example library(MASS) data(cats) m <- svm(Sex~., data = cats) plot(m, cats) ## more than two variables: fix 2 dimensions data(iris) m2 <- svm(Species~., data = iris) plot(m2, iris, Petal.Width ~ Petal.Length, slice = list(Sepal.Width = 3, Sepal.Length = 4))ですから、plot(m2,...)は「m2についてfomulaで指定する次元について2Dplotする、そのときの色分けはsliceで指定する」といっていると思います(推測)。-- Akira
Akira (2006-01-14 (土) 13:33:59)
Rの引用(Materials&Methodsとして)を何種類か見たのですが統一されていないようにも思います。それとも、お作法などあるのでしょうか?
bioconductorとExcelでもできる簡単な処理ですが、投稿したいと思ってます。
柳 (2006-01-13 (金) 13:02:58)
x<-c(1:3) y<-c(4:6) z<-matrix(0,length(y),length(x)) myfunction<-function(x,y){ for(i in 1:length(x)){ z[,i]<-x[i]+y } }をするとparse errorが出ますがなぜかぜんぜんわかりません
教えて下さい
まいける (2006-01-12 (木) 17:35:58)
はじめまして。
Rにoptimという関数がありますが、そのソースコードを入手することは可能でしょうか?
また、そのソースコードはC等で書かれているものなのでしょうか?
大変初歩的な質問で恐縮ですが、なにとぞご教示よろしくお願い致します。
Moz (2006-01-10 (火) 20:25:23)
x <- factor(c(rep("A",4),rep("B",2),rep("C",3))) x [1] A A A A B B C C C Levels: A B C x <- x[x!="C"] x [1] A A A A B B Levels: A B ClevelからCを消すにはどうしたらよいのでしょうか。table()で集計する
ときにじゃまになります。無理矢理x <- as.factor(as.character(x)) x [1] A A A A B B Levels: A Bとしましたが、もっと正統なやり方があるのではないかと思って質問しました。
x <- x[, drop=T]ですね。 -- びい 2006-01-11 (水) 06:11:20
ひらい (2006-01-09 (月) 00:10:48)
windowsでR2.2.0を使っています。
MySQLの中のデータベースに入っているデータを使って、ユークリッド距離(関数dist)を求めたいのですがどうしたらよいでしょうか?
どなたかご教授お願いします。
aa (2006-01-05 (木) 00:18:37)
Windows XP で,Rを使っています。
長いプログラム(1000行ほど)をRのコンソールに一気に流した際,画面に表示される速度がversionによって全く異なります。
余程複雑な計算でなければ,新しいversion(2.20以降)の方が,画面に表示される時間がかなり遅くなります。
計算の実行速度ではなく,表示時間が遅くなります。
この原因をご存知の方,アドバイスを頂けると幸いです。
本太 (2006-01-01 (日) 23:36:29)
変数が列として表されているCSVファイル(d1.csv)からデータをxに読み込む場合、x<-read.table("C:/R/d1.csv",header=F,sep=",")と記述しています。
変数が行として表されているCSVファイル(d2.csv)を読み込み上記と等しいデータフレームxを作成するためには、どのようにすればよいのでしょうか?
すみませんがよろしくお願いします。----- d1.csv ------ height weight age 172 72.6 25 159 59.8 27 164 78.7 42 179 85.4 35 164 65.7 30 ------------------- ------------d2.csv------------------- height 172 159 164 179 164 weight 72.6 59.8 78.7 85.4 65.7 age 25 27 42 35 30 --------------------------------------
x <- t(read.table("c:/r/d1.csv" ,header=T,sep=","))でいいかと. -- なかま 2006-01-02 (月) 01:12:43
x <- data.frame(t(read.table("c:/r/d1.csv" ,header=T,sep=",")))でしょうか?ただ、要素の属性は変わると思います。 -- Akira 2006-01-02 (月) 13:32:36
変数が列として表されているCSVファイル(d1.csv)からデータをxに読み込む場合、 x<-read.table("C:/R/d1.csv",header=F,sep=",")は
x<-read.table("C:/R/d1.csv",header=T,sep=",")ではないかと思うのはさておき,
height,172,159,164,179,164 weight,72.6,59.8,78.7,85.4,65.7 age,25,27,42,35,30さて,そうであれば,あなたの本当に望むことは(回答例に示されるように)そんなに簡単には求められないので,以下のようになるのではないかと思いますね。
> x <- t(read.table("d2.csv" ,header=F,sep=",")) > name <- x[1,] > x <- x[-1,] > x <- matrix(as.numeric(x), nrow(x)) > colnames(x) <- name > x <- data.frame(x) > x height weight age 1 172 72.6 25 2 159 59.8 27 3 164 78.7 42 4 179 85.4 35 5 164 65.7 30d2.csv を読み込んだときと同じ結果ということは,このようなことを意味するのではないかと,私は考えましたのですが,深読みですか? -- あけましておめでたかったろう 2006-01-02 (月) 23:17:57
にゃあ (2006-01-01 (日) 04:41:36)
Rcmdrでデータのインポート⇒テキストファイルからでCSVファイルを読んだ後、データを表示を選ぶと、データエディタを動かそうとすると、画面がひきつったようになって沢山の残像が出るような画面になってしまいます。
これってバグなんでしょうか?
また、同様に出力した筈のCSVファイルでさえ、データを読み込める場合と読み込めない場合もあります。これもバグなんでしょうか?
じゃんぽけ (2005-12-29 (木) 11:43:38)
Tips等調べて1日格闘してみましたが、わからないポイントがあり投稿しております。
【質問内容】
日付に関する演算を繰り返し行うプログラムを書く予定でして、
1.文字型で格納しておいて取り出す度にDATES()で変換して使う
2.文字型で格納しておいて文字型に変換したものを別の変数に保存して使う
3.最初から日付型で格納しておいて取り出す
で3.の方法で実行する方法がないかなと思い模索している所です。
日付型の変数をデータフレーム内の要素に代入しようとすると、シリアル値が代入されてしまいます。
これを日付型のまま代入する方法はございますでしょうか?
(列をベクトル演算で変換した場合には日付型で代入できるようです。)
また、まだまだ勉強不足のため、下記のような挙動になる要因が何なのかが理解できていない状況です。
何度も質問致しまして恐縮ではございますがご教授頂ければ幸いです。
よろしくお願いします。
【プログラム例】> DATE <- dates("05/12/31",format="y/m/d") > DATE2<- "05/12/31" > > test1 <- data.frame(a=rep(NA,3)) > test1$a[1] <- DATE > test1 a 1 13148 2 NA 3 NA > mode(test1$a) [1] "numeric" > > test2 <- data.frame(a=rep(NA,3)) > test2$a[1] <- DATE2 > test2 a 1 05/12/31 2 <NA> 3 <NA> > mode(test2$a) [1] "character" > > test3 <- data.frame(a=rep(NA,3)) > test3$a[1] <- dates(test1$a[1],format="y/m/d") > test3 a 1 13148 2 NA 3 NA > mode(test3$a) [1] "numeric" > > test4 <- data.frame(a=rep(NA,3)) > test4$a <- dates(test2$a,format="y/m/d") > test4 a 1 05/12/31 2 <NA> 3 <NA> > mode(test4$a) [1] "numeric"
> data.frame(a=rep(dates("12/31/05"),3),b=0) a b 1 13148 0 2 13148 0 3 13148 0 > data.frame(a=dates(rep("12/31/05",3)),b=0) a b 1 12/31/05 0 2 12/31/05 0 3 12/31/05 0
ところで、上記のようにデータフレームにNAなど値を入れる場合であれば列の型を定義できるのですが、要素が全くない場合に事前に型を宣言する方法はあるのでしょうか?
自己レスの修正をしました。 データフレーム作成時に、変数を日付型に宣言しておくことで 空データになっても日付型を保持していますが、 更新時に個別の行を更新しなければ日付型を保持できないようです。
> #データの定義をした後に空データにする > test <- data.frame(a=dates(c(NA,NA),format="y/m/d"),b=c(NA,NA)) > test <- test[FALSE,] > test [1] a b <0 rows> (or 0-length row.names) > > #空のデータフレームに行を追加するとちゃんと日付を保持 > test <- rbind(test,c("2000/01/03",0)) > test a b 1 00/01/03 0 > > #データフレームが1行の場合、全行をまとめて更新すれば日付型で更新 > test$a <- dates("2000/01/22",format="y/m/d") > test a b 1 00/01/22 0 > > #更に1行追加。(実験1用と実験2用のデータを作成) > test1 <- rbind(test,c("2000/01/05",0)) > test2 <- test1 > test1;test2 a b 1 00/01/22 0 2 00/01/05 0 a b 1 00/01/22 0 2 00/01/05 0 > > #実験1:データフレームが2行以上の場合、 #全行をまとめて更新すれば日付型でなく なる > test1$a <- dates("2000/01/22",format="y/m/d") > test1 a b 1 10978 0 2 10978 0 > > #実験1:このデータフレームに日付を文字型で #いれようとすると変数aは文字型に変化 > test1 <- rbind(test1,c("2000/01/03",0)) > test1 a b 1 10978 0 2 10978 0 3 2000/01/03 0 > mode(test1$a) [1] "character" > > #実験2:この場合は、データフレームの行を #それぞれ更新すればOK > > test2$a <- dates(rep("2000/01/22",length(test2$a)),format="y/m/d") > test2 a b 1 00/01/22 0 2 00/01/22 0
じゃんぽけ (2005-12-28 (水) 11:51:41)
一通りTips や Q&A を見てプログラムを書いてみたのですが、もっと簡単にできる方法がないかと思い投稿しております。
以下のように、「FOR文を利用して、条件に合致する行を抽出した上で、新たな変数(列)を追加・代入したい」と考えています。
考えたやり方では、
「1.事前に追加する変数(列)は定義しておく」
「2.抽出条件に当てはまる行がない場合は、IF文でエラーを飛ばす」
ということで解決はできているようなのですが、もっと簡単にできる方法があるのではないかと思っています。
ご教授頂ければ幸いです。宜しくお願いします。
> #その1:新たなdという変数(列)に直接的に値を追加・代入 #テストデータを作成 > test <- data.frame(a=c(1,2,3),b=c(4,5,6),c=c(7,8,9)) > test a b c 1 1 4 7 2 2 5 8 3 3 6 9 > > #合致する行がなくても代入するとエラー > for(i in 4:3){ + test[test$a == i,]$d <- 2 + } 以下にエラー"$<-.data.frame"(`*tmp*`, "d", value = 2) : 置き換えは 1 列ですが、データは 0 列です > > #IF文で合致しない行の場合にスキップしてもエラー > for(i in 4:3){ + if(nrow(test[test$a == i,]) != 0){ test[test$a == i,]$d <- 2 } #こうやれば合致しないものにはエラーはでない + } Warning message: 4 個の変数を与えて、3 個の変数を置き換えようとしています in: "[<-.data.frame"(`*tmp*`, test$a == i, , value = list(a = 3, > > #その2:事前にdという変数(列)を作っておいて値を追加・代入 #テストデータを 作成 > test <- data.frame(a=c(1,2,3),b=c(4,5,6),c=c(7,8,9)) > test$d <- NA > test a b c d 1 1 4 7 NA 2 2 5 8 NA 3 3 6 9 NA > > #合致する行がなくても代入するとエラー > for(i in 4:3){ + test[test$a == i,]$d <- 2 + } 以下にエラー"$<-.data.frame"(`*tmp*`, "d", value = 2) : 置き換えは 1 列ですが、データは 0 列です > > #IF文で合致しない行の場合にスキップすればOK > #ただし、文章としてはやや冗長感がある? > for(i in 4:3){ + if(nrow(test[test$a == i,]) != 0){ test[test$a == i,]$d <- 2 } #こう やれば合致しないものにはエラーはでない + } > > test a b c d 1 1 4 7 NA 2 2 5 8 NA 3 3 6 9 2 >
> ( tmp <- subset(test, a==2) ) # a==2 となる行を抽出 a b c 2 2 5 8 > ( transform(tmp, d=0) ) # データフレーム tmp に変数 d を追加 a b c d 2 2 5 8 0 > transform(subset(test, a==2), d=0) # 上の 2 行をひとまとめに
Shigu (2005-12-27 (火) 18:21:52)
はじめまして、データ分析のためにRの勉強を始めたShiguと申します。
googleと過去ログで調べたのですが、見つからなかったので質問させていただきます。
欠損値がそこらじゅうにあるcsvデータ(400サンプル×130属性)に対して、主成分分析を行いたいのですが、欠損値をna.omitで省いてしまうと、4サンプル×130属性とほとんどデータがなくなってしまいます。
このような欠損値のあるデータに対して、各属性の第一主成分の寄与の大きさを見るために主成分分析を行うにはどのようにすればいいのでしょうか?
どなたかご教授どうぞよろしくお願いします。
鍋 (2005-12-26 (月) 22:24:05)
> mydata [1] "15" "17" "17" "22" "11" "22" "13" "20" "19" "18" "11" "11" "18" "18" "24" [16] "13" "21" "17" "19" "11" "10" "19" "13" "6" "5" "5" "5" "3" "2" "2" [31] "0" "2" > barplot(mydata) Error in -0.01 * height : non-numeric argument to binary operator何でこういうErrorが発生するのでしょうか?どう対応すればいいのでしょうか?わかる人教えてください!!
mydata <- c(15,17,17,22,11,22,13,20,19,18,11,11,18,18,24, 13,21,17,19,11,10,19,13,6,5,5,5,3,2,2,0,2) barplot(mydata)とすれば大丈夫。-- soysoy 2005-12-26 (月) 23:16:45
S. Sugawara (2005-12-25 (日) 14:37:16)
hoge1.sav-hoge50.savと連番になっているspssファイルをまとめてcsvにしたいのですが
library(foreign)
tmp="hoge"
for(i in 1:50)
eval(parse(text=paste(tmp,i,"<-read.spss(", """ , tmp, """ ,i,".sav, use.value.labels=TRUE, to.data.frame=TRUE)",sep="")))
のようにすると"""の部分で文法エラーになります。しかしこれをはずすとreadの中のhoge.savをファイル名と認識してくれないのでエラーとなります
"を含む文字列を作成できれば解決すると思うのですが、その方法を御存じなら教えてください。または他の方法があればよろしくお願いします
> tmp="hoge" > i=12 > text=sprintf("%s%i<-read.spss(?"%s%i.sav?", use.value.labels=TRUE, to.data.frame=TRUE)", tmp, i, tmp, i) > text [1] "hoge12<-read.spss(?"hoge12.sav?", use.value.labels=TRUE, to.data.frame=TRUE)"二つ目の引用符の位置がおかしいと思いましたので直しましたが。 -- 2005-12-25 (日) 15:08:54
Akira (2005-12-24 (土) 00:01:14)
複数のデータに対してヒストグラムを重ねて、データの違いを視覚化しようとしています。ヒストグラムの柱の頂点を結ぶ曲線を描きたいので、plot(density(x))のy軸をFrequencyに変えたいです。#例えばmatというdata.frameがあるとして mat <- data.frame(rnorm(1000), rnorm(1000), rnorm(1000)) for(i in seq(ncol(mat))){ dat <- mat[, i] dat <- density(dat) dat$y <- dat$y * nrow(mat) plot(dat, xlim=c(-10, 10), ylim=c(0, 100), col=i) if(i != ncol(mat)) par(new=TRUE) }などとすればdensity(x)のy軸がFrequencyに変わると思ったのですが、うまくいきません。基本的なことかもしれませんがよろしくお願いします。
dat$y <- dat$y plot(dat, xlim=c(-10, 10), ylim=c(0, 0.5), col=i)ではないでしょうか? -- 2005-12-24 (土) 09:51:25
plot(dat, xlim=c(-10, 10), ylim=c(0, 400), col=i, ylab="Frequency")でしょうけど。 -- 2005-12-24 (土) 22:03:24
velvet (2005-12-23 (金) 13:13:52)
現在以下の環境でRを勉強しているvelvetと申します
OS:windows XP、R:Version2.2.0
現在ある生物の出現が、どのような環境要因によって説明できるについて解析しています。
具体的には以下のようなデータです
100地点において、ある生物が存在したかどうか(A、0(出現しなかった)と1(出現した)のダミー変数です)を従属変数に、その100地点での3つの環境の変数(B、C、D、すべて連続変数です)を独立変数として、ロジスティック回帰モデルで、考えられる全ての独立変数の組み合わせのAICを比較しようと思っています。
このような方法はweb上でいくつか存在しましたが
私の場合、3つの環境変数の間に相関が存在するため、これらの交互作用項を含むモデルも比較しなければなりません。
つまり
glm(A~B, family=binomial) から
glm(A~B + C + D + B*C + B*D + C*D + B*C*D, family=binomial)というモデルまでのAICを計算したいわけです。
どなたか具体的な方法を御存知でしょうか?
宜しく御願い致します
name <- c("B","C","D","B*C","B*D", "C*D","B*C*D") nv <- length(name) str <- character(nv) n <- 2^nv-1 result3 <- character(n) for (i in 1:n) { k <- i m <- 0 for (j in 1:nv) { if (k%%2) { m <- m+1 str[m] <- name[j] } k <- k%/% 2 } result3[i] <- paste(c("glm(A ~", paste(str[1:m], collapse=" + "), ", family=binomial)"), collapse=" ") } print(result3,quote=FALSE)これを,sink でファイルに保存して,しかるべき時に source で読み込んで実行するとか。
[1] glm(A ~ B , family=binomial) [2] glm(A ~ C , family=binomial) : [127] glm(A ~ B + C + D + B*C + B*D + C*D + B*C*D , family=binomial)127通りありますからねぇ。
サイハテ (2005-12-20 (火) 18:00:30)
現在、以下の環境でRの勉強をしているサイハテと申します。
OS:Unix(Sun solaris8)
R:Version2.0.1
Rの使用形態:簡単なスクリプトを”R CMD BATCH”で実行する
非常に稚拙な質問で申し訳ないのですが、Unix上のRでグラフ(plotコマンド)を出力する際、デバイスの制限があるようで、現状PDF形式でしか出力をすることができません。。。
できればR CMD BATCHで、WindowsではメジャーなJPEGやPNG形式で出力したいと考えているのですが、何か方法等ありますでしょうか?
教えていただけると非常に助かります。
よろしくお願いいたしますm(_ _)m
以下、?deviceの引用Devices package:grDevices R Documentation List of Graphical Devices~ Description:~ The following graphics devices are currently available: * 'postscript' Writes PostScript graphics commands to a file * 'pdf' Write PDF graphics commands to a file * 'pictex' Writes LaTeX/PicTeX graphics commands to a file * 'xfig' Device for XFIG graphics file format * 'bitmap' bitmap pseudo-device via 'GhostScript' (if available). The following devices will be available if R was compiled to use them and started with the appropriate '--gui' argument: * 'X11' The graphics driver for the X11 Window system * 'png' PNG bitmap device * 'jpeg' JPEG bitmap device None of these are available under 'R CMD BATCH'.
#ref(): File not found: "hikaku.png" at page "初級Q&A アーカイブ(4)"
実際のプレゼン画面上での pdf はもっときれい。-- 2005-12-20 (火) 19:44:09itoken (2005-12-19 (月) 21:07:52)
ARモデルのあてはめをAICを使って、自分でプログラムし、次数選択したのですが、Rでもともと定義しているAICの次数選択と異なる結果になってしまいました。自分のモデルは
(n-k)×log{Q/(n-k)}+2×k k:次数 n:データ数 Q:残差平方和
というモデルで計算しました。
RでのAICの式はどのようになっているのかを教えていただけると非常に助かります!よろしくお願い致します。
三田 (2005-12-19 (月) 18:18:31)
現在、Rwikiのページを参考に
Meadow+ESSでRのスクリプト作成を行っております。
Version:
Meadow Emacs 21.4.1と表示されています
ESS 2005-09-07版
R 2.1.1
上記のような環境において、Meadow上で変数名として_(アンダースコア)を
入力したいのですが、ESSのショートカットキー
Rwikiのページ:
http://www.okada.jp.org/RWiki/index.php?ESSUsage
において_(アンダースコア)は自動的に”<-”に変更されてしまい、
何とも困っております。
Googleや本Rwikiでも情報収集を試みたのですが、いきづまってしまいました。。ご教授いただければ幸いです(T_T)。
小田 (2005-12-15 (木) 10:04:39)
一通りTips や Q&A を見てみたのですが、力不足で最適なプログラムを思いつかず投稿しています。
どなたかご教授頂けるとうれしく思います。
どうかよろしくお願いします。
【実現したい事】
検索用ベクトルの値を使い、データフレーム内のある列で検索値に該当した全ての行を抽出したいと考えています。>data<-data.frame(LABEL=c("A","B","C","D", "E","F","G","H","I","J"), VALUE=c(1:10)) >data LABEL VALUE 1 A 1 2 B 2 3 C 3 4 D 4 5 E 5 6 F 6 7 G 7 8 H 8 9 I 9 10 J 10 >#検索したい値 >search <- c("A","B","C") >search [1] "A" "B" "C" >#実現したいこと >data[data$LABEL =="A" | data$LABEL =="B" | data$LABEL =="C" ,] LABEL VALUE 1 A 1 2 B 2 3 C 3いくつか方法を考えているのですが、最適なプログラムを思いつけず困っています。
(1) 条件式をpasteで記述後、data[条件式,]を実行> exec <- NA > for(i in 1:length(search)){ + exec <- paste(exec, ifelse(i > 1," | ",""), "data$LABEL =='",search[i],"'",sep="") + } > exec [1] "NAdata$LABEL == 'A' | data$LABEL == 'B' | data$LABEL == 'C'" > data[exec,] LABEL VALUE NA <NA> NA(2) IFELSEを使ってみた
> data[ifelse(search == data$LABEL,TRUE,FALSE),] LABEL VALUE 1 A 1 2 B 2 3 C 3 Warning messages: 1: 長いオブジェクトの長さが短いオブジェクトの 長さの倍数になっていません in: is.na(e1) | is.na(e2) 2: 長いオブジェクトの長さが 短いオブジェクトの長さの倍数になっていません in: "==.default"(search, data$LABEL)どなたかご教授お願いしますm(_ _)m
柴田 (2005-12-14 (水) 17:44:56)
こんにちは。宜しくお願い致します。
OS:AIX5.2
AIXにRをインストールしている段階なのですが、makeコマンド実行時に以下のメッセージが出力されてしまいます。
「make:make 1254-025 既存詳細ファイルが存在するか、またはターゲットを指定する必要があります」
作業内容としましては、CRANにて最新のソース「R-2.2.0.tar.gz」をダウンロードして、FTPにてAIXサーバの/home/operatorにコピーしました。それから、gunzipコマンド、tarコマンドにて上記ディレクトリにR_2.2.0というディレクトリが作成されました。そのディレクトに移動して./configureコマンドを実行した後に、makeコマンドを実行すると上記メッセージが出力されてしまします。
この回避方法をご存知の方がいましたらご教授願いたいのですが、宜しいでしょうか。お願い致します。必要な情報などございましたらお送りいたします。
どうぞ宜しくお願い致します。
R初心者 (2005-12-14 (水) 15:21:25)
重回帰モデルの変数選択を行う上で、
変数増加法を用いたいのですが、
stepAIC()の使い方が良くわかりません。
helpによると、scope = list(upper=・・・
としている辺りが増加法だと思うのですがよく理解できませんでした。
result<-glm(y~???)
result2<-stepAIC(result,scope=list(upper=x1*x2*x3))
とするのでは?とも思いましたが、???の部分に何を入れたら良いのか
わかりませんでした。追加する1つ目の変数はこちらで、
設定しておいたほうが良いのでしょうか?
どなたか、ご指導お願いいたします。
x1<-rnorm(100) x2<-rnorm(100) x3<-rnorm(100) y<-rnorm(100) g<-glm(y~1,gaussian) # 切片から始めたいなら。 library(MASS) stepAIC(g,direction="forward", scope=list(upper=~x1*x2*x3)) # scope=~x1*x2*x3でもよい。
upperの後の~を忘れないように。--KI 2005-12-14 (水) 17:14:48
TOTO (2005-12-13 (火) 20:46:56)
数量化 I 類で総当り法による変数選択(指標はAIC)をしたいのですが、どうすればよいのでしょうか?
as.factorを使って説明変数を変換した後、mle.aic()を使ってみたのですが、アイテム変数単位ではなく、ダミー変数単位で選択されてしまいます。
Akira (2005-12-12 (月) 11:01:01)
bioconductorの質問はここで大丈夫でしょうか?
使用環境はXPpro+闇R2.2.0です。Hmiscを読み込んでからBiobaseを読み込むとBiobaseの読み込みに失敗します。逆は大丈夫でした。# これだとエラーが出てbiobaseが読み込めません library(Hmisc); library(Biobase) # これだとOK。Biobase、Hmiscともに読み込めます。 library(Biobase); library(Hmisc) # エラーメッセージです > library(Hmisc) > library(Biobase) 要求されたパッケージ tools をロード中です 以下にエラーsetMethod("contents", "environment", function(object, all.names) { : no existing definition for function 'contents' エラー:.onLoad は 'Biobase' のための 'loadNamespace' に失敗しました エラー:'Biobase' に対するパッケージもしくは 名前空間のロードが失敗しましたloadNamespaseが変更されているという意味だと思うのですが、読み込みに順番があるのでしょうか?あるとすれば、どこを調べると良いのでしょうか?
よろしくお願いします。
阿部勝延 (2005-12-12 (月) 09:30:20)
R2.2.0(Windows版、Mac版)をインストールしたのですが、「R Book」P50にあるようにプロキシーサーバーの設定をしたのですが、アクセスできません。
状況は
Windows Xp Professional+R2.2.0
プロキシーサーバーがあり、そこからファイアーウォールを抜けていくのでRのプロパティでリンク先"...Rgui.exe”--internet2と設定しています。しかし、Rを立ち上げ、ミラーサーバーの設定でJapan(Aizu)を選択して、パッケージの一覧の取得で落ちます。
Macintosh(MacOSX10.3.9)+R2.2.0
こちらは、ちょっと複雑で、先ほどのWindowsマシンをProxyとして、インターネットに出て行きます。
せっていは、P50にあるとおり、
Sys.putenv("http_proxy"="http://proxyhost:8080")
を実行します。このコマンドは通りますが、パッケージの一覧の取得では、取得できません。
もしいいアドバイスがありましたら、よろしくお願いいたします。
伊藤 (2005-12-11 (日) 11:46:51)
ある条件に従うデータだけにおいて線形回帰をしたいと考えています。x1<-gl(2,10) x2<-rep(gl(2,5),2) x3<-runif(20) y<-rnorm(20) lm(y~x3)このとき、{x1=1,x2=1}の場合だけを回帰したいときはどうすればよいでしょうか。
以下のようなことを試しましたが、うまくいきません。lm(y~x3,subset=list(x1==1,x2==1)) #エラー lm(y~x3,subset=x1==1 && x2==1) # 全部のデータ? lm(y~x3,subset=c(x1==1,x2==1)) # 最初の指定だけ有効?お手数ですが、どうかよろしくおねがいします。
> data <- cbind(x1,x2,x3,y) > data x1 x2 x3 y [1,] 1 1 0.56633985 -0.86012923 [2,] 1 1 0.33230350 0.92239242 [3,] 1 1 0.99880629 0.39366260 [4,] 1 1 0.17977729 -0.36911252 [5,] 1 1 0.14474499 -0.74874814 [6,] 1 2 0.13511632 -0.74440885 [7,] 1 2 0.63733995 -0.28768662 [8,] 1 2 0.20289775 0.61580370 途中省略 [18,] 2 2 0.32580467 1.07307234 [19,] 2 2 0.01826360 1.16538160 [20,] 2 2 0.50006782 1.18701904 > data2 <- subset(data, x1==1 & x2==1) > data2 x1 x2 x3 y [1,] 1 1 0.5663398 -0.8601292 [2,] 1 1 0.3323035 0.9223924 [3,] 1 1 0.9988063 0.3936626 [4,] 1 1 0.1797773 -0.3691125 [5,] 1 1 0.1447450 -0.7487481というところかなぁ。と。-- 2005-12-11 (日) 15:00:17
Omar (2005-12-11 (日) 02:51:25)
中身が多くて、スクロールバーを使っても最初の方の部分が表示されないようなオブジェクトを
すべて見るためにはどうすればよいのでしょうか?
今見ようとしているのは線形モデルオブジェクトです。
Win XP, R 2.2.0を使っています。
g (2005-12-09 (金) 12:12:50)
例えば,a.Rというファイルを作り,中にlibrary(Rcmdr)と記述して,コマンドプロンプト(WinXP SP2を使っています.)から,R CMD BATCH a.Rと打つと,Rcommanderが起動するのですが,すぐにRも終了するので,実際は何もできません.何かオプションを加えて,Rを起動したままにすることができるのでしょうか.バッチ処理の性質上無理なのでしょうか.
実際は,R CMD BATCH a.Rをb.batというファイルに記述し,バッチファイルから起動しようと思ったのですが,うまくいきませんでした.
ご意見いただけると幸いです.
library(grDevices) library(Rcmdr)にすると、動作としては合っているものになるようです。-- 岡田 2005-12-14 (水) 13:06:18
mura (2005-12-07 (水) 11:15:10)
OS: Mac OS X 10.4.3
R: Version 2.2.0 (2005-10-06 r35749)
「環境設定」で文字の色を変更すると、次回以降も有効なままでした。
ところで、メニューバー→フォーマット→フォント→フォントパネルから文字の種類や大きさを変更したら、その場では有効なのですが、次回以降は、また初期設定に戻ってしまいます。
フォントの設定を保存するにはどうしたらいいのでしょうか
にゃあ (2005-12-07 (水) 00:38:30)
パッケージの更新しようとしたら次のエラーメッセージが出てきます。
警告:package 'lattice' is in use and will not be installed
“lattice”ってパッケージは更新できないのでしょうか?
なお、Rは闇2.2.0で、Windows XPで使用しております。
> help(vanilla) No documentation for 'vanilla' in specified packages and libraries: you could try 'help.search("vanilla")' > help.search("vanilla") No help files found with alias or concept or title matching 'vanilla'って出てくるんですが・・・・・・???-- にゃあ 2005-12-08 (木) 05:51:43
サイトウ (2005-12-05 (月) 18:02:07)
WinXPで闇R-2.2.0を使用しています。このバージョンだとRMySQLが対応していないということなのでMySQLにODBCドライバをインストールしてRと連携しようとしています。
そこで、このサイトに載っている手順通りにDBIとRODBC双方のパッケージを読み込んでlibraryの定義をして、次の手順に進むと以下のようなエラーメッセージが出てしまいます。> m <- dbDriver("MySQL") 以下にエラーdo.call(as.character(drvName), list(...)) : 関数 "MySQL" を見つけることが できませんでしたどうしたらよいでしょうか?
どなたかご教授お願いします。
odbcConnect("データソース名","ユーザー名", "パスワード","ケース")という指定になると思います。細かい指定内容はRjpWikiには無いようです。CRAN本家に行ってRODBCの英語のマニュアルを読むか、The R Bookには載ってます。 -- okinawa 2005-12-05 (月) 20:44:59
ミズノ (2005-12-03 (土) 11:48:13)
jaccard係数を導いたのですが、その値を使ってクラスター分析のような樹形図が描けるのでしょうか?
つまり、こちらで距離を指定して樹形図を描けるのでしょうか?
どなたかご指導お願いします。
> x <- matrix(c(0,1,2,1,0,3,2,3,0),3) > x [,1] [,2] [,3] [1,] 0 1 2 [2,] 1 0 3 [3,] 2 3 0 > y <- as.dist(x) > y 1 2 2 1 3 2 3 > ans <- hclust(y) > plot(ans, hang=-1)こんな感じですね。 -- 2005-12-03 (土) 12:41:58
かわぐち (2005-12-02 (金) 10:53:31)
御世話になります。
生存曲線をプロットするには、Rによる統計処理や中澤先生のRによる統計解析の本にlibrary(survival), Surv, survfitあるいはkm.survによる方法がかかれていますが、二つの群をpar(new=T)で重ね合わせようとするとx軸がずれます。
そこでsurvfitを用いplot(res,xlim=c(0,20))とすると軸はあいますが95%信頼区間が邪魔です。
95%信頼区間の点線をのぞく方法はどうしたらよいでしょうか
また、2群ではなく3群、4群となった場合にスマートにプロットさせる方法は
ありますか。
データは”Rによる統計処理”から引用させていただきました。# 富永祐民「治療効果判定のための実用統計学 − 生命表法の解説 −」蟹書房(第3回改訂版)74 ページ,表 4.1 のデータ
# 解析結果は 83-85 ページ
# 1 は A 群,2 は B 群を表す
group <- c(1, 1, 2, 2, 1, 1, 2, 2, 1, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, 2, 2, 1, 1, 1)
event <- c(1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0)
time <- c(2, 20, 5, 1, 3, 17, 2, 3, 15, 14, 12, 13, 11, 11, 10, 8, 8, 3, 7, 3, 6, 2, 5, 4, 2, 3, 1, 3, 2, 1)
res <- survfit(Surv(time,event)~group,conf.type="none") plot(res)ところで、グループ毎の直線を区別するためにイベント発生、打ち切り時のマークを変更するにはどうしたらよいでしょうか。 -- 2005-12-02 (金) 13:19:32
plot(res, lty=c(1,2)) plot(res, col=c(1,2)) plot(res, lwd=c(1,2))などでどうでしょうか?でも、pchはだめみたいですね。僕も教えてほしいです。
蓮井 (2005-11-30 (水) 17:50:27)
ファイル名が、1から40までの数字2つの組合せでできているファイル群があり(0101.txt、0102.txt・・0140.txt、0201.txt・・4040.txt)、
これを順番に開いて計算結果をファイル出力する作業をしたいと考えています。
ただ、このファイル群には、すべての組合せ分1600ファイルはなく、全体では800ファイル程度となっています。ファイル名のリストはありません。
そこで、> for (i in 0:40) { > if (i < 10) I <- paste("0",i,sep="") else I <- as.character(i) > for (j in 1:40) { > if (j < 10) J <- paste("0",j,sep="") else J <- as.character(j) > fname0 <- paste(I,J,sep="") > fname1 <- paste(fname0, ".txt", sep="")として、ファイル名を自動的に作ってファイルを読み込もうとしたのですが、該当するファイルがない場合に「ファイル '0001.txt' を開くことができません 」となって止まってしまいます。
ファイル名がない場合に、その回をスキップするようにするには、どのようにすればよいのでしょうか?
あるいは、別の方法を検討した方がよろしいでしょうか?
よろしくお願いいたします。
filenames <- outer(1:40, 1:40, function(i, j) sprintf("%02i%02i.txt", i, j)) for (str in filenames) { if (file.access(str) == 0) { ファイルが存在するので, 処理を行う } }場合分けで paste というのは,避けたい。 -- 2005-11-30 (水) 19:24:09
Akira (2005-11-30 (水) 13:53:30)
WinXPsp1で闇R-2.2.0を使用しています。Sweaveで画像を貼り付けたいのですが、pdf(..., cidfamily="")のような指定はどのようにすれば良いでしょうか?par(family="")でもうまくいきません。
setHook(packageEvent("grDevices", "onLoad"), function(...) grDevices::ps.options(cidfamily=""))のように設定してください -- なかま 2005-11-30 (水) 15:21:41
ikeike (2005-11-29 (火) 22:56:18)
R初心者@勉強中です。
R 2.2.0にパッケージRcmdr Version 1.1-2を入れていろいろといじっているのですが、
統計量→モデルへの適合、に「非線形回帰」を組み込む(?)
ことは可能でしょうか?
Rでnlsを使えばよいのでしょうが、R Commanderの手軽さが気に入ってしまいました。
Fox教授のHPには、
Fit models
Linear regression
Linear model
Generalized linear model
Multinomial logit model
Proportional-odds logit model
(http://socserv.mcmaster.ca/jfox/Misc/Rcmdr/)
とあり、今のところ対応の予定はなさそうですが。
よろしくお願い致します。
Rcmdr.nls <- function(formula,data=NULL,start, control=NULL, algorithm=NULL, trace=NULL, subset=NULL, weights=NULL, na.action=NULL, model=NULL, lower=NULL, upper=NULL ){ UpdateModelNumber() modelValue <- paste("NonLinearRegressionModel.", getRcmdr("modelNumber"), sep="") if(is.null(data)) data <- ActiveDataSet() else data <- paste(substitute(data)) args <- paste(c(paste("start=", paste(c(substitute(start)))), paste("control=", paste(c(substitute(control)))), paste("algorithm=", paste(c(substitute(algorithm)))), paste("trace=", paste(c(substitute(trace)))), paste("subset=", paste(c(substitute(subset)))), paste("weights=", paste(c(substitute(weights)))), paste("na.action=", paste(c(substitute(na.action)))), paste("model=", paste(c(substitute(model)))), paste("lower=", paste(c(substitute(lower)))), paste("upper=", paste(c(substitute(model))))) [!sapply(list(start,control, algorithm,trace,subset,weights, na.action,model,lower,upper), is.null)],collapse=",") command <- paste("nls(", c(formula), ", data=", data, ",", args, ")", sep="") logger(paste(modelValue, " <- ", command, sep="")) assign(modelValue, justDoIt(command), envir=.GlobalEnv) doItAndPrint(paste("summary(", modelValue, ")", sep="")) activeModel(modelValue) } listNonLinearModels <- function(envir=.GlobalEnv, ...) { objects <- ls(envir=envir, ...) if (length(objects) == 0) NULL else objects[sapply(objects, function(.x) "nls" == (class(eval(parse(text=.x), envir=envir))[1]))] } Rcmdr.nls(density~a*conc,start=list(a=1)) #test code (active data set in #Rcmdr should be `DNase')
oa (2005-11-29 (火) 22:20:26)
仮に以下のようにデータを読み込ませ、condition,sentence,primeを独立変数に、responseを従属変数にして分散分析したい場合には、どのように入力すればいいのでしょうか?condition sentence prime response 1 1 2 1471 1 2 1 812 2 1 2 830 2 2 1 1498 3 1 2 792 3 2 1 1068 4 1 2 753 4 2 1 612
きたむら (2005-11-29 (火) 22:02:46)
ベクトルからある数値を取り出すには
例えば1を取り出すとき、test<-c(1,2,3,4,5,0)
test[-1]
[1] 2 3 4 5 0
となりますが、
0を取り出そうとすると、test[-0]
numeric(0)
となってしまいます。
[-0]は取り出す以外の意味があるのでようか?
また、特定の値が複数ある場合、何個あるかを指定しないで
それら全てを取り除くことは可能でしょうか?
例えば[1,2,0,0,0,3]→[1,2,3]
どなたか、ご指導よろしくおねがいします。
かわぐち (2005-11-27 (日) 21:01:09)
下のようなデータがあり、対数グラフを作成しました。x <- c(3,8,11,15,16) y <- c(5,6,9.6,25.6,29.6) plot(x,y,log="y")初めのデータy=5を基準値として横線を引き、2番目以降のデータ増加量を直線回帰して交点を計算するとx=7.38となりました。
この2直線をグラフに追加するにはどうしたらよいでしょうか。abline(5,0) abline(lm(log(y)~x))ではうまく線が入りません。
よろしくおねがいします。
x <- c(3,8,11,15,16) y <- c(5,6,9.6,25.6,29.6) plot(x,y,log="y") x2 <- x[-1] y2 <- y[-1] ans <- lm(log(y2)~x2) x3 <- c(6, 18) y3 <- exp(ans$coefficients[1]+ans$coefficients[2]*x3) lines(x3, y3, col="red") abline(h=5, col="blue")
#ref(): File not found: "exp.png" at page "初級Q&A アーカイブ(4)"
ってとこでどうでしょ。 -- 2005-11-27 (日) 23:31:41ans2 <- nls(y2~a*b^x2,start=list(a=1.1,b=1.2)) p <- ans2$m$getPars() y4 <- p[1]*p[2]^x3 lines(x3, y4)として,比較してみればよろしいかと。 -- 2005-11-28 (月) 21:37:33
大軽貴典 (2005-11-25 (金) 10:20:00)
はじめまして。
Rでニューラルネットワークなどを使ってみようと思っています。
いま、手元にデータセットがあります。このデータセットは事前にランダムに10個のグループ(各グループは同数)に分けております。グループの情報はカラムを追加させて1,2,..,10というふうに入力しています。
これを用い、90%を学習セットに用い、残りの10%検証用に使おうと思っています。
一度だけなら分けて入力すればよいのですが、グループカラムの情報を使って、1回目はグループ1を学習セットに、2回目はグループ2と繰り返し最終的に10回の作業をしようと思っています。
いきなりすべての工程をコード化できなくてもよいので、うまくあるデータセットからグループカラムの情報を使って学習セットと検証セットの切り出しをする方法はないでしょうか?
データセットは次のような構成になっています。
目的変数 グループ 説明変数1 説明変数2 3 4 5 ...
A 1 1 -3.4
B 3 0 3
C 9 0 5.2
B 10 1 7
どなたかアドバイスお願いできないでしょうか?
よろしくお願いいたします。
for (i in 1:10) { 学習セット <- データフレーム[データフレーム["グループ"] != i, ] 検証セット <- データフレーム[データフレーム["グループ"] == i, ] 何らかの分析 }どうでしょう。
for (i in 1:10) { 学習セット <- subset(データフレーム, データフレーム["グループ"] != i) 検証セット <- subset(データフレーム, データフレーム["グループ"] == i) 何らかの分析 }ですね。 -- 2005-11-25 (金) 10:22:04
model<-svm(subset(training,select=-Object),training$Object)のような命令でNAがあるみたいなエラーが出ています。なお、trainingは、教えていただいた方法で作成した学習セット、Objectとは名義型の目的変数です。
aa (2005-11-23 (水) 17:04:35)
はじめまして。
Rで作成したグラフィックスをPDFで書き出し,TeXでPDFを読み込む際に,エラーが発生します。
対応策をご存知の方,ご示唆を頂けると幸いです。
【作成例 R】setwd("C:/question20051123") #ディレクトリ plot(1:10, main="日本語の入ったPDF") #画像を作成 #画像 をPDFで保存 dev.copy(pdf, "fig1.pdf", cidfamily="Japan1Ryumin", width=12, height=5) dev.off() graphics.off()【作成例 TeX】
?documentclass[12pt,a4paper]{jarticle} ?usepackage[dvips]{graphicx} ?begin{document} %図の挿入 ?begin{figure}[htbp] ?begin{center} ?includegraphics[keepaspectratio=true, height=100mm]{fig1.pdf} ?end{center} ?end{figure} ?end{document}コンパイル時のError Messageは,
! LaTeX Error: Cannot determine size of graphic in (ファイル名) (no BoundingBox)当方の環境は以下の通りです。
OS: Windows XP R: 闇R2.20 TeX: LaTeX2e Ghostscript: Ghostscript 8.52
ebb graph_name.pdfというコマンドだったはず。-- 蓮見 2005-11-24 (木) 01:23:04
bmc -b graph_name.pdfだったような。ほかには mediabb.sty を利用するという方法も。いずれにせよ,R で作成した pdf が pdf として問題のないものになっているのであれば,R ではなく TeX の問題ですから,TeXWiki あたりを参照されては。 -- 2005-11-24 (木) 02:11:14
かわぐち (2005-11-22 (火) 16:09:00)
はじめまして。最近Rを知って感激して統計とともに勉強中です。
さて、生存分析にて型どおりKaplan-Meier,log-rank,Cox回帰分析をしたいと思います。例数は20ー30ほどですが、イベント発生数が群別で0のものや2ー3しかないものがあります。log-rankでは有意差があるのですが、Cox回帰はイベント発生が5個以上必要とのことです。このような場合の群別の比較はどうしたらよいのでしょうか。
よろしくおねがいします。
ねこのて (2005-11-16 (水) 23:37:42)
300行400列のテキスト形式の行列が2つあります。
(細密数値情報と言うデータで、仕様は
http://www.gsi.go.jp/MAP/CD-ROM/saimitu/htmls/format.html にあります。)
この2つのファイルの同じセルについて、数値(土地利用の分類コード)が何から何に変化したかを、場合毎に数えたいと思っています。
(1から2がmケース、1から3がnケース、、、、、)
Rでのデータハンドリングの経験がなく、探せる範囲ではヘルプも探しましたがよくわかりませんでした。
効率的な方法をお願いいたします。
> set.seed(1234) > x <- matrix(sample(1:5, 12, replace=TRUE), 4, 3) > y <- matrix(sample(1:5, 12, replace=TRUE), 4, 3) > # プログラムコンテストは,ここから開始 # 以下で 1000 は,適当な数。必要なら10000でも100000 でも > ans <- table(as.vector(x)*1000+as.vector(y)) > z <-as.integer(names(ans)) > ans <- cbind(floor(z/1000), z%%1000, ans) > colnames(ans) <- c("before", "after", "frequency") > ans before after frequency 1001 1 1 1 1002 1 2 1 2002 2 2 1 3001 3 1 1 3002 3 2 1 4001 4 1 1 4002 4 2 3 4005 4 5 2 5002 5 2 1代替案がいくつか出た段階で,どれが一番計算時間が短いかというコンテストがあるかもしれませんが,そのコンテスト結果が出るまでには,第一案なり何なりが既に答えを出していることは間違いない。(答えが正しいかどうかは別問題。正しい答えを出さないプログラムは問題外)
> x <-read.fortran("S4_2617.TDU", c("i7", "400i2")) > dfx <-data.frame(x) > names(dfx) <-c("ID", paste("var", 1:400, sep="")) > dfx$ID <-NULLとして読込み
> ans <- table(as.vector(x)*1000+as.vector(y))を実施したところ、
> 以下にエラーvector("integer", length) : 指定されたベクトルのサイズが長すぎます、となって実施できませんでした。 どこが間違っているのでしょうか? お手数おかけしますがご指摘いただければ幸いです 。
Aki (2005-11-16 (水) 00:46:20)
自己共分散関数を求めるacfで1次元の配列ではなくn×mの2次元配列を引数にもってくると
グラフにはm×mのグラフが表示されます。
これの非対角要素はccfで求められる相互共分散と同じなのでしょうか?
> ccf function (x, y, lag.max = NULL, type = c("correlation", "covariance"), plot = TRUE, na.action = na.fail, ...) { type <- match.arg(type) if (is.matrix(x) || is.matrix(y)) stop("univariate time series only") X <- na.action(ts.union(as.ts(x), as.ts(y))) colnames(X) <- c(deparse(substitute(x)), deparse(substitute(y))) acf.out <- acf(X, lag.max = lag.max, plot = FALSE, type = type)ということで,引数をまとめて,acf を呼んでいるので acf と同じ結果になるわけでしょう。 -- 2005-11-16 (水) 13:02:30
Akira (2005-11-14 (月) 10:05:32)
ヒストグラムの作図においてy軸に中断線を入れることはできますでしょうか?"中断線"では検索できず、他に良いキーワードを知らないので、上手く見つけられませんでした。
例えば、y軸が1〜100の範囲で1〜20と80〜100で作図したいと思っています。
よろしくお願いします。
muramatsu (2005-11-11 (金) 15:48:14)
あるデータのヒストグラムを描きたいのですが、ばらつきが大きいので、両対数のヒストグラムを描きたいと思っています。
"ヒストグラムと密度の推定"という記事をみましたが、両対数での描き方は載っていませんでした。
plot関数と同様に、hist(data, log="xy")ともやってみましたが、
hist.default(data, log = "xy") : 'x' は数値でなければなりませんというエラーが出てしまいました。しかし、ちゃんと数値データです。
plot(table(data))ではきちんとヒストグラムが描けています。
また、ヒストグラム作成時に区切り幅も対数で決める必要があると思いますが、それはとりあえずおいといて、plot(table(data), log="xy")のようにやってみました。すると、
1: 軸の限界が有限ではありません [GScale(-1.#INF,2.18469,2, .); log=1] 2: Internal(pretty()) の範囲が大きすぎます.. 修正しましたという警告がでて、plotに失敗しました。
どうしたらよいでしょうか?よろしくお願いします。
x <- round(rlnorm(1000),0) plot(table(x))みたいな感じかな。
x <- rlnorm(1000) hist(log(x)) # hist(x) と比較せよみたいなのでよいのでは?横軸の目盛りも,対数目盛りになりますよ。
z <- x <- rlnorm(1000) hist(log(x), xlab="") z <- floor(log10(z)) z2 <- 1:10*10^(log.min <- min(z)-2) if ((n <- max(z)-log.min) > 0) { for (i in 1:n) { z2 <- c(z2, z2*10^i) } } log.z2 <- log10(z2) axis(1, at=log.z2, labels=z2, pos=-30)みたいな感じかな。描かれる横軸二種。上が変換後の値の数値に基づいて目盛った目盛り,下が変換前の数値に基づいて目盛った目盛り。(位置は決めうちしているので,汎用性はない) -- 青木繁伸 2005-11-11 (金) 18:01:49
#ref(): File not found: "pdf.png" at page "初級Q&A アーカイブ(4)"
[1] 83.553719 80.760000 100.750000 ....のようなベクトルデータです。
dist <- density(data) plot(dist, log="xy")下のようにヒストグラムを描こうとすると、
dist <- table(data) plot(dist, log="xy")やはり次のエラーが出ます。
1: 軸の限界が有限ではありません [GScale(-1.#INF,2.18469,2, .); log=1] 2: Internal(pretty()) の範囲が大きすぎます.. 修正しました次のようにしても、
hist(data, log="xy")このようなエラーが出ます。
高水準 plot 関数で,パラメータ "log" を設定できません両対数ヒストグラムを描くにはどうしたらいいのでしょうか? -- muramatsu 2005-11-24 (木) 23:15:28
r <- structure(list(breaks = breaks, counts = log10(counts), intensities = dens, density = dens, mids = mids, xname = xname, equidist = equidist), class = "histogram")に書き換える(counts=counts を counts = log10(counts) にする)。 次に上に示されてことを縦軸にも施す。さっき作った hist2 を呼ぶことに注意。
old <- par(mar=c(5,5,1,1), xpd=TRUE) z <- x <- rlnorm(1000) # hist ではなく hist2 を! frq <- hist2(log(x), xlab="", ylab="") # 横軸を描く z <- floor(log10(z)) z2 <- 1:10*10^(log.min <- min(z)-2) if ((n <- max(z)-log.min) > 0) { for (i in 1:n) { z2 <- c(z2, z2*10^i) } } log.z2 <- log10(z2) axis(1, at=log.z2, labels=z2, pos=-0.27) # 縦軸を描く z <- floor(frq$counts) z2 <- 1:10*10^(log.min <- min(z)-2) if ((n <- max(z)-log.min) > 0) { for (i in 1:n) { z2 <- c(z2, z2*10^i) } } log.z2 <- log10(z2) axis(2, at=log.z2, labels=z2, pos=-4.4) par(old)2つの軸を描く部分が冗長なので,一つの関数にまとめるとよい。このようなヒストグラムは見たことない。 -- 2005-11-25 (金) 00:24:01
#ref(): File not found: "pdf2.png" at page "初級Q&A アーカイブ(4)"
山形 (2005-11-10 (木) 23:16:59)
モデル選択をする場合、AICを基準にモデル選択をするのは知ってますが、
扱っているサンプルが少ないため、AICではモデルを過大評価してしまうことが
あると聞きました。
その場合、補正したAIC(c-AIC)を使うといいということも聞きましたが、この
補正AICをRで行なうことは出来ますか?探してみたのですが、力不足で探せませんでした。よろしくお願いします。
yuta (2005-11-09 (水) 15:41:12)
表題の件で質問です。
以下に例示するような、サンプリングのインターバルが不定で欠損値を含む時系列データを処理する方法はないでしょうか?
なお、行列testmatの1列目は日付のシリアル値(windowsで使用する、1900/1/1=1とする整数値)です。単純にこの列を列ラベルにしたりしてts()で処理することはできるのですが、インターバルについて思うように処理できていないように思います。test<-c({38482,NA,NA,NA,NA,NA, 38531,NA,0.992,NA,NA,NA, 38552,NA,0.988,0.994,0.990,NA, 38562,NA,0.993,0.995,0.992,NA, 38576,NA,0.989,0.990,0.994,NA, 38595,NA,0.994,0.996,0.990,NA, 38608,NA,0.931,0.992,0.990,NA, 38621,0.985,0.956,0.991,0.991,NA, 38637,0.984,0.670,0.980,0.976,NA, 38649,0.989,0.879,0.972,0.983,NA }) testmat <- matrix(test,c(10,6),byrow=T) ts(testmat)使用バージョンはR2.2.0(闇国際化)、OSはWindows2000です。 また、時系列オブジェクトの生成には基本パッケージのts()を使用しています。 宜しくお願いします。