初心者のための R および RjpWiki に関する質問コーナー

新規投稿はできません。




encodingをグローバルに設定したい

akira? (2007-12-31 (月) 22:26:19)

WinXPとubuntuを並用していています。ubuntuでSHIFT_JISを指定せずに読み込んで、grepするとエラーになってしまいます。
仮にhogeの中身を

test\n
テスト\n
test\n

とします。

事前にencodingがわかっていれば、

> con <- file("hoge.tsv", open="r", encoding="SHIF_JIS")
> x <- readLines(con)
> close(con)
> grep("hoge", x) <- SHIFT_JISを指定しないと In grep("hoge", x) :  入力文字列 2 はこのロケールでは不適切です となります
1 2

としています。

これを、UTF-8とSHIFT_JISを区別せずに使いたいのですが、グローバルに設定を変える方法、もしくは検索キーワードを教えてほしいです。

localeのLC_CTYPEを変更すれば良いと思うのですが...

> sessionInfo()
R version 2.6.1 (2007-11-26)
i486-pc-linux-gnu

locale:
LC_CTYPE=ja_JP.UTF-8;LC_NUMERIC=C;LC_TIME=ja_JP.UTF-8;LC_COLLATE=ja_JP.UTF-8;
LC_MONETARY=ja_JP.UTF-8;LC_MESSAGES=ja_JP.UTF-8;LC_PAPER=ja_JP.UTF-8;
LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=ja_JP.UTF-8;
LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

loaded via a namespace (and not attached):
[1] rcompgen_0.1-17
  • options(encoding="cp932") とかは, いかがでしょか. 尚, SHIFT_JISを使ってはいけません. 時々, cp932のデバッグに `localedef -i ja_JP -c -f SHIFT_JIS -A /usr/share/locale/locale.alias ja_JP.SJIS' 等と言うような副作用満点な行為を行うことがありますが, 思い留まった方が健康には良いでしょう.(Winあげるのめんどいんで) -- なかま 2008-01-07 (月) 22:25:30
  • ありがとうございます。これですね。この記事を知ってはいたのですが、特に目立つエラーもないので、いまひとつ意味もよくわからないまま、使ってました。最近、Rまわり(?)の作業が多く、本筋以外で悩むことが多いです… -- akira? 2008-01-10 (木) 07:12:18

計数データのPoissonGLMによるモデルからカウントされる0の度数を求めるコマンドについて

idaten? (2007-12-27 (木) 04:25:41)

0が多い計数データ、solder.balanceというデータをGLMでPoisson回帰しました。

> P=glm(formula = skips ~ ., family = poisson, data = solder.balance2)

そのモデルからカウントされる0の数を調べたいときに、色々調べていましたら、

"the obserbved zero counts are compared to expected number of zero counts for the likelihood model:"
> round( sum( dpois( 0,fitted(P) ) ) )

として得られると書いてありました。
fitted(P)は推定されたPoisson分布の平均パラメータλですから、このコマンドの意味は各λで0をとる確率の和ということになると思うのですが、これでなぜ0をカウントすることができるのでしょうか?対数尤度関数として考えるのでしょうか?
そもそも訳が間違っていて、このコマンドは0を得るのでしょうか?
もし情報不足、不適な記述がありましたら、申し訳ございません。指摘していただけると助かります。

  • これは当てはめで得られたポアソン分布パラメータに対するポアソン分布の0確率を求める命令に見えますが。その(理論)確率は一般に実数になりますから丸め関数で整数に落とすのでしょう。わからないのは sum 関数をとっていることですが、もしかして fitted(P) は複数の推定値からなるベクトルになっているのでしょうか。 -- 2007-12-27 (木) 07:12:52
  • 説明不足で申し訳ありません。確かにfitted(P)は720の推定値からなるベクトルになっています。データの詳しい説明もなしに申し訳ありませんでした。全てのゼロ確率を足してから丸めるのと、各ゼロ確率を丸めるたのを足すのでは意味も値もぜんぜん異なると思うのですが、この記述は前者ですが、後者がだめな理由がわかりません -- idaten? 2007-12-27 (木) 09:35:30

マッピングのためのgridの設定

R 初心者? (2007-12-23 (日) 22:13:19)

こんにちは。
空間統計学のためにRを使っている者ですが、次の質問に対して皆さまのヘルプを頂けると嬉しいです。

アメリカ(127.5W-67.7W, 24N-52N)領域を4x5のgridに区切り、それぞれの値をgridごとに計算し、その値の空間時間分布を解析しようとしています。
すなわち、gridの中央点は070W42N,070W46N,...となります。
しかし、この地理情報をshape fileとして与えなければ絵を描くことができません。
そのようなgrid情報はArcIGSで書かなければならないのでしょうか。
もしそうであれば、ArcGISを一度も使ったことがないものに対して、簡単なアドヴァイスを頂けると嬉しいです。
どうかよろしくお願いします。

  • RjpWiki カテゴリー別見取り図からたどれるGIS関係記事は既に見たのですか。 -- 2007-12-24 (月) 09:50:57
  • shape fileとして与えなくてはいけないという思い込みはどこから来たのでしょうか?また、shape fileの扱いにArcGISが必要という思い込みもどこから来たのでしょうか?gridは最も単純なので、spパッケージを使うまでのない気もしますが、幾何変換やいろいろ処理が必要ならspを使えばよいと思います。 -- 谷村 2007-12-29 (土) 00:56:19
  • ご回答頂きありがとうございます。確認が遅れ申し訳ございません。 krigはできるものの中央点と境界を与えるだけではgridを一つずつ色を塗ることができないので、ArcGISで情報を与えなければならないと考えました。どのようにしてgridに色塗りができるか教えて頂けると幸いです。 -- R 初心者? 2008-01-05 (土) 09:52:00
  • http://r-spatial.sourceforge.net/gallery/ ここに例題がありますので、勉強してみてください。ここにあるソースを見るとわかりますが、X,YからSpatalPolygons?クラスを作成して、その四角いポリゴンに色を塗っているようです。R内部では、ShapeFileを直接扱っているのではなく、SpatalPolygons?クラス等のオブジェクトに変換して扱います。もちろんArcGISなどで四角いポリゴンがタイル状に並んでいるものを作成して、それをRに読み込んで色を塗ることはできるかもしれませんが、Gridの大きさを変更するたびにShapeFileを作り直さなければなりませんので現実的ではありません。 -- okinawa 2008-01-07 (月) 09:01:43
  • okinawaさま、ご回答頂きありがとうございます。今から勉強します。 -- R 初心者? 2008-01-10 (木) 20:39:42
  • 書き込みに気がつきませんでした。2kmグリッド(32x16)を作成して色を塗る単純な例を示します。 -- 谷村 2008-01-12 (土) 16:06:44
    library(sp)
    grd <- GridTopology(cellcentre.offset=c(932668,1350376),cellsize=c(2000,2000),cells.dim=c(16,32))
    grd.sp <- as.SpatialPolygons.GridTopology(grd)
    grd.d <- data.frame(runif(512), row.names=getSpPPolygonsIDSlots(grd.sp))
    grd.spd <- SpatialPolygonsDataFrame(grd.sp,grd.d)
    spplot(grd.spd)
sp-grid-20080112.png

正規分布とロジスティック分布の密度関数の重ね描きは?

idaten? (2007-12-21 (金) 04:26:30)

 正規分布N(0,1)とロジスティック分布の密度関数の違いを図で確認するために、この2つを一つの図で重ねがきしたいのですが、どのようにしたらいいのでしょうか?それぞれの密度関数はかけても、x軸の範囲が違ったりするので、一つの図できれいに表現することができません。どなたかコマンドを教えていただけると助かります。よろしくお願いします。

  • xlim=c(-3,3) 等でx軸作図範囲を共通にしたら。 -- 2007-12-21 (金) 07:04:44
  • これで
    plot(x <- seq(-4, 4, .1), dlogis(x), ylim=c(0, .5), type="l", ylab="Density")
    points(x, dnorm(x), type="l", col="red")
    どうでしょう -- 2007-12-21 (金) 08:25:34
  • 上の2行目は,points ではなく lines を使うと type="l" などとする必要はなくなります。 -- 2007-12-21 (金) 08:52:57
  • 丁寧なお返事ありがとうございます。うまく図が描けました。しかし図を描いてあまり二つの図が似ていないのにはびっくりしました。 -- idaten? 2007-12-22 (土) 01:31:11
  • リンク関数をLogitとProbitにしたときのGLMで推定結果(ResidualDeviance?)に差がなかったことを、リンク関数の密度関数が酷似しているからだと思っていたのですが、もう少し考えてみる必要がありそうです。どうもありがとうございました。 -- idaten? 2007-12-22 (土) 01:37:24

あるベクトルから条件に合致した要素を取り出す際に、条件が複数の場合はどうすればいい?

AobaAoba? (2007-12-20 (木) 18:24:29)

どうぞよろしくお願いします(WinXPで、バージョン2.5.1使用です)。

下記のxのベクトルのうち、yのベクトルの要素に含まれるものだけを取り出して、zというベクトルを作りたいのです。

x <- c(2,2,3,6,2,0,6,7,9)
y <- c(2,7)

そこで、単純には以下のようにすれば良いことは初心者ながらにわかります。

z <- x[x == 2 | x == 7]

しかし、yがもっと長いベクトルの場合はどうすれば良いのでしょうか?
私が思いつけるのは「forループを使って、yの要素を1つずつ取り出してxにあるかどうかチェックし、ベクトルを結合していく」という以下のような方法です。

z <- c()
for(i in 1:length(y)){
z <- c(z,x[x == y[i]])
}

もっと簡単な方法があればよろしくお願いいたします。

  • x[x%in%y] -- 2007-12-20 (木) 18:45:07
  • えっ?それだけで良いんですか?あれ、確かにあってる。どうもありがとうございました。 -- AobaAoba? 2007-12-20 (木) 18:55:25

データフレーム名の自動生成or変換

うち? (2007-12-20 (木) 17:13:04)

irisのデータを例にして説明します。
種名ごとに抽出したデータフレームを自動で生成し,そのデータフレーム名を種名と対応させたいと考えております。
下記のような方法を試してみましたが,当然ながら全くダメでした。データフレーム名の自動生成方法,またはデータフレーム名変換方法を御教示いただければ幸いです。

種名 <- levels(iris$Species) #種名を抜き出す
種数 <- length(種名) #種数を取得
for(i in (1:種数) ) {
  paste(種名[i]) <- subset(iris, Species==種名[i]) 
}
# 「関数 "paste<-" を見つけることができませんでした」とエラーが出る。
for(i in (1:種数) ) {
  df <- subset(iris, Species==種名[i])
  names(df) <- paste(種名[i]) 
}
# これだとデータフレーム名はそのままで,列名が変わってしまう。
  • 実際に何がやりたいかによるんだろうけど,別に新しくデータフレームなんか作らない方が良いのでは?名前の違うたくさんのデータフレームができたら,今度はそれらを統一的に使うのが面倒になるだけ。
    以下のような使い方があるのはご存じ? -- 2007-12-20 (木) 19:00:43
    a <- split(iris, iris[,5])
    lapply(a, summary)
    summary(a[["versicolor"]])
  • 上と同意見ですが、どうしても別にしたいなら、list()で代用してはだめですか?
    dat <- lapply(unique(iris$Species), function(x) subset(iris, Species==x))
    names(dat) <- unique(iris$Species)
    として、
    dat[["setosa"]]
    とかする -- akira? 2007-12-20 (木) 19:08:16
  • こんなことがやりたい? -- 2007-12-20 (木) 20:31:27
    > for(i in (1:種数) ) assign(種名[i],  subset(iris, Species==種名[i]) )
    > ls()
    [1] "setosa"     "versicolor" "virginica"  "種数"       "種名"
  • 皆様,ご回答有り難うございます。
    また,私の投稿を読みやすく修正していただき有り難うございました。
    assign を教えて下さった方,どうもありがとうございます。assign だと変数名と値を同時に代入できるんですね!
    しかし,上のお二方のご指摘のようにデータ処理方針を考え直す必要があるようですね。私のやりたいことをもう少し詳しく説明してみます。
    例はirisですが,底びき網に入った魚の話です。
    底びき網には魚が溜まる場所が数箇所あり,場所によって種組成・サイズ組成が違います。
    大量の魚が獲れるため,魚の溜まる場所(サンプル区分)ごとに抽出サンプリングをしており,サンプル区分ごとに抽出率は異なります。データは下記のような感じになります。
    iris2 <- iris
    iris2$区分 <- rep(c("A","B"))
    iris2$抽出率 <- ifelse(iris2$区分=="A", 0.8, 0.4 )
    そのまま度数分布表を作成したのでは網全体のサイズ組成がゆがんでしまうので,抽出率を元にして種ごとに網全体のサイズ組成を復元したいというのが本来の目的です。
    種ごとサンプル区分ごとにデータフレームを抜き出して度数分布表を作って,抽出率の平均値の逆数を度数分布表に乗じてやれば元の組成を復元できると考えたのですが,種数×サンプル区分はかなりの数になりますので,データ処理の自動化を狙って先の質問をさせていただきました。
    上のお二方のご提案を参考にして問題解決ができるのか,考えてみます。
    R初心者なので,すぐにお返事できないと思いますが,どうかご容赦下さい。 -- うち? 2007-12-20 (木) 21:15:21
  • paste() 関数は文字列及び文字列化した数値を分離文字をはさんで結合した文字列を作る関数ですから、文字列が一つだけなら使う必要はありません。参考までに paste() 関数を使った連番変数名の作成とそれへの値の付値の例を紹介します。 -- 2007-12-21 (金) 08:48:57
    > for (i in 1:10) assign( paste("mydata",i,i+5,sep="."), i:(i+5) )
    > ls()
     [1] "mydata.1.6"   "mydata.10.15" "mydata.2.7"   "mydata.3.8" "mydata.4.9" 
     [6] "mydata.5.10"  "mydata.6.11"  "mydata.7.12"  "mydata.8.13" "mydata.9.14" 
    > mydata.4.9
    [1] 4 5 6 7 8 9
  • 解決できました。 -- うち? 2007-12-21 (金) 20:04:11
    皆様に頂いたアドバイスを参考にして,下記のようなコードで自分がやりたかった処理を実現できました。
    やはり,ご指摘のように大量のデータフレームを作ってから処理を行うのは非効率だと分かりました。
    今回は大変勉強になりました。どうも有り難うございました。
    # データの生成
    iris2 <- iris
    iris2$区分 <- as.factor(rep(c("A","B")))
    iris2$抽出率 <- ifelse(iris2$区分=="A", 0.8, 0.1 ) 
    
    種名 <- levels(iris2$Species)
    種数 <- length(種名)
    区分名 <- levels(iris2$区分)
    区分数 <- length(区分名)
    
    階級下限 <- seq(0,10,by=0.5)[-length(seq(0,10,by=0.5))]
    
    # 種別・区分別の度数分布を作製して抽出率補正した結果を1つのデータフレームにまとめる。
    for( i in (1:種数)){
     for( n in (1:区分数)){
       temp <- subset(iris2, Species==種名[i] & 区分==区分名[n])
       ans <- hist(temp$Sepal.Length, br=seq(0,10,by=0.5), right=FALSE, plot=FALSE)
       補正度数 <- ans$counts/mean(temp$抽出率) # ここで抽出率補正
       階級下限 <- cbind(階級下限, 補正度数) 
       colnames(階級下限)[ncol(階級下限)] <- paste(種名[i],区分名[n],sep=".")
       補正サイズ組成 <- data.frame(階級下限) # データフレーム化。重複した列をまとめてくれる。
    }}
    
    # 補正後の度数分布を種ごとに合体して,全体?のサイズ組成を復元して作図。
    pdf(file="補正サイズ組成.pdf", width=7, height=3) 
    for( i in (2: (length(補正サイズ組成))) ){
      ifelse( i%%2==0, {
      title <- substring(colnames(補正サイズ組成[i]), 1, nchar(colnames(補正サイズ組成[i]))-2 ) #タイトルを種名にする
      barplot(補正サイズ組成[[i]]+補正サイズ組成[[i+1]], names.arg=補正サイズ組成[[1]], main=paste(title, "補正サイズ組成",  sep="") )
      }, next) 
    }
    dev.off()
  • こんなふうな図を添付したかったのでしょう。やっていることに比べてプログラムは長すぎる気がしますね。もっと単純に記述できそうに思いますが,やってみようという気が起きません(^_^) -- 2007-12-21 (金) 22:14:48
    iris78975839789.png
  • 私の汚いソースを直す気にならないのは当然ですよね。すっきりしたソースを書けるように勉強します。画像添付は、最初jpgで失敗してpngでやり直したのですが上手くいきませんでした。上のグラフをアップしていただいた方と同様のソースを書いたつもりだったのですが。。 -- うち? 2007-12-21 (金) 22:42:56
  • こんな風になるかな(短くはならないけど,少しは見通しが良くなる?) -- 2007-12-22 (土) 00:16:36
    iris2 <- iris
    iris2$区分 <- as.factor(rep(c("A","B")))
    種名 <- levels(iris2$Species)
    種数 <- length(種名)
    brks <- seq(0,10,by=0.5)
    split.df <- split(iris2,
    	list(iris2$Species, iris2$区分)) # データファイルを2要因で分割
    ans <- sapply(split.df, function(temp) # 集計
        hist(temp$Sepal.Length, br=brks, right=FALSE, plot=FALSE)$count)
    dim(ans) <- c(20, 3, 2) # 配列に割り付け
    ans2 <- ans[,,1]/0.8+ans[,,2]/0.1 # 重み付けの合計
    pdf(file="補正サイズ組成.pdf", width=7, height=3) # 出力ファイルの用意
    for (i in 1:種数) { # for を使わないとかえって面倒
    	barplot(ans2[,i],
    	  names.arg=brks[-length(brks)], # x軸の名前
    	  main=paste(種名[i], "補正サイズ組成",  sep="")) # 表のタイトル
    }
    dev.off()
  • こんなにスッキリするとは!ありがとうございます。split(x,f)のfにリストを与えると複数の要因で分割できるのですね。その他の技についてもじっくり解読させていただきます。 -- うち? 2007-12-22 (土) 14:10:13

ARIMAモデルの次数選択

kuri? (2007-12-20 (木) 02:14:02)

ARIMAモデルの次数選択をAICを選択基準として求めるプログラム(すべての組合せを行い、最小のAICのときの次数を返す)を作ってみたのですが、途中でエラーが出て停止してしまいます。
エラーが出ても無視して、全通りを実行することはできますか?

  • ?try -- 2007-12-20 (木) 05:17:27

微分できない関数のoptim()での最適化

Katoh? (2007-12-18 (火) 12:37:52)

はじめまして。
optim()関数について質問です。

こちらの「汎用最適化関数 optim() の使用法」ページによると、
>微分できないような関数では"Nelder-Mead" 法 や"SANN" 法 しか選択肢は無い。
とありますが、
実際に準ニュートン法や共役勾配法でも微分できない関数の最適(と思わしき)値を見つけることができます。
どのようにして微分値を求めているのでしょうか?それとも求めずに何か他の手を使って最適化を行っているのでしょうか?

よろしくお願いします。

  • 関数値だけを使うアルゴリズムですから、微分値は仮にそれが存在しても使いません。具体的な詳細は自分で調べましょう。 -- 2007-12-18 (火) 12:50:08
  • >関数値だけを使うアルゴリズム  これはNelder-Meadや焼きなまし法のことですよね。それとも準ニュートン法や共役勾配法がoptim()内では微分値を使わないアルゴリズムとして実装されているのでしょうか? (質問の仕方が悪くて誤解されたのかな・・・) -- Katoh? 2007-12-18 (火) 13:46:54
  • (誤解されました ...)。 optim の準ニュートン法や共役勾配法オプションは、与えられた関数が微分可能かどうか判断できません。導関数が明示的に与えられれば、それで微分値を計算します(与えられた導関数が間違っていても)。いっぽう、導関数を明示的に与えなければ、内部的に数値微分(イメージとしては微分商)で微分値を推定します。この場合も、本当にそこで微分できるかどうかは(数値微分が失敗しない限り)どうでもよいことです。おそらくあなたが試した例は解がたまたま微分できる範囲内にあったような例では無いですか。なお、このことは help(optim) を読めば皆書いてあります。 -- 2007-12-18 (火) 16:06:21
  • ありがとうございます。そしてごめんなさい、ちゃんとhelp読みます。「微分値は仮にそれが存在しても使いません。」というのは、(関数の)微分値が、という意味なんですね。導関数を与えても使わないのかと誤解しました。 これまた誤解かもしれませんが、ということは「微分できないような関数では"Nelder-Mead" 法 や"SANN" 法 しか選択肢は無い。」とあるのは<良い最適値を見つけるためには>という前提のもと書かれているということでよろしいでしょうか。 -- Katoh? 2007-12-18 (火) 16:23:09
  • 「良い最適値」とは厳密に言えば形容矛盾ですね。一方で数値誤差を除いても、いかなる最適化手法も局所的最適値しか見つけられないのも事実で、最後の吟味はユーザの責任になります。あなたの得たらしい解がそうなのかどうかは ご自分で判断してください。 -- 2007-12-18 (火) 22:19:04

barplotで、特定のbarだけをかえたいのですが

[[<ふ>]] (2007-12-15 (土) 22:10:24)

たとえば、
x <- matrix(seq(10,20,2),6,2)
のようなデータを用意して、

barplot(t(x))
とすると、各barは、同じ二色(defaultでは、glay scaleですが)で塗り分けられます。

これを、(例えば)右から二本目のbarの色を、"lightblue","green"に変更するには、どのようにしたらいいでしょうか。

(bar内が分割されてない場合は、そのbarだけ色を変えることはできています。)

よろしくお願いいたします。

  • やはり、そういう引数ないですよね。なら、ちからずくで、と思ってはいたものの、どうやらればいいのかわかりませんでした。ありがとうございます。助かりました。 -- [[<ふ>]] 2007-12-15 (土) 23:49:44
  • xx<-t(x); barplot(xx,col=c("blue","red")); par(new=T); xx[,c(1:4,6)]<-NA; barplot(xx,col=c("green", "black"),ylim=par()$usr[3:4]) みたいなのも思い出してあげてください. -- なかま 2007-12-17 (月) 10:22:54

ftableで作成したクロス集計で合計がゼロの行を消したい

Transalp? (2007-12-14 (金) 18:32:49)

お世話になります。
要素が4つ以上あるデータフレームの値をtableで集計し、ftableしてクロス集計表を作りましたが、行合計がゼロのデータが集計表に表示され、困っております。

hoge.txtの中身
所属	名前	月度	出席
営業部	営業太郎	01月	はい
総務部	総務花子	02月	はい
営業部	営業次郎	03月	はい
購買部	購買好男	03月	はい
経理部	経理京子	04月	はい
技術部	技術冴子	04月	はい
営業部	営業太郎	04月	はい

Rコンソールから

df <- read.delim("hoge.txt", header=TRUE)
tbl <- table(df[,1], df[,2], df[,3])
ftbl <- ftable(tbl)

ftable実施後

                01月 02月 03月 04月
                                   
営業部 営業次郎     0    0    1    0
       営業太郎     1    0    0    1
       技術冴子     0    0    0    0
       経理京子     0    0    0    0
       購買好男     0    0    0    0
       総務花子     0    0    0    0
技術部 営業次郎     0    0    0    0
       営業太郎     0    0    0    0
       技術冴子     0    0    0    1
       経理京子     0    0    0    0
       購買好男     0    0    0    0
       総務花子     0    0    0    0
経理部 営業次郎     ・・・以下略

営業部には営業太郎と営業次郎しかいないはずなので、

tbl <- tbl[rowSums(tbl) > 0,, drop=FALSE]

としましたが「次元が違う」とエラーになってしまいます。そもそも、こういうやり方ができるのかどうかをここ2〜3日調査しているのですが、ftableオブジェクト(?)に関する情報があまり見つけられなくてはまっております。

この「行合計がゼロ」のftableを削除する方法をご教授下さい。
実行環境はWindowsXP SP2 R-2.6.1 です。

よろしくお願い申し上げます。

  • プログラムを書いてしまえばよいのだが。集計次元数が決め打ちなんだけど,以下のようにすれば一応はそれらしいものが。 -- 2007-12-14 (金) 21:35:37
    > df2 <- data.frame( table(df[,1], df[,2], df[,3]))
    > df2 <- df2[df2[,ncol(df2)] > 0, ]
    > a <- reshape(df2, idvar=c("Var1", "Var2"), timevar="Var3", direction="wide")
    > a[is.na(a)] <- 0
    > colnames(a) <-c(colnames(df)[1:2], dimnames(tbl)[[3]])
    > rownames(a) <- 1:nrow(a)
    > a
        所属     名前 01月 02月 03月 04月
    1 営業部 営業太郎    1    0    0    1
    2 総務部 総務花子    0    1    0    0
    3 営業部 営業次郎    0    0    1    0
    4 購買部 購買好男    0    0    1    0
    5 技術部 技術冴子    0    0    0    1
    6 経理部 経理京子    0    0    0    1
  • おなじく、決め打ちなんですが
    x <- apply(ftbl,1,sum)>0
    y <- cbind(rep(attr(ftbl, "row.vars")[[1]], each=length(attr(ftbl, "row.vars")[[2]])), 
               rep(attr(ftbl, "row.vars")[[2]], times=length(attr(ftbl, "row.vars")[[1]])))
    z <- data.frame(y[x,],ftbl[x,])
    colnames(z) <- c("所属","名前",unlist(attr(ftbl,"col.vars")))
    z
        所属     名前 01月 02月 03月 04月
    1 営業部 営業次郎    0    0    1    0
    2 営業部 営業太郎    1    0    0    1
    3 技術部 技術冴子    0    0    0    1
    4 経理部 経理京子    0    0    0    1
    5 購買部 購買好男    0    0    1    0
    6 総務部 総務花子    0    1    0    0
    とか、 集計だけ欲しいなら、
    df$"name" <- apply(df[,1:2],1,paste,collapse="-")
    x <- matrix(0, nrow=length(unique(df$"name")), ncol=length(unique(df$"月度")),
                dimnames=list(unique(df$"name"),unique(df[,3])))
    for(i in rownames(x)) for(j in colnames(x)) x[i,colnames(x)%in%df$"月度"[df$"name"%in%i]] <- 1
    x
                    01月 02月 03月 04月
    営業部-営業太郎    1    0    0    1
    総務部-総務花子    0    1    0    0
    営業部-営業次郎    0    0    1    0
    購買部-購買好男    0    0    1    0
    経理部-経理京子    0    0    0    1
    技術部-技術冴子    0    0    0    1
    とか -- akira? 2007-12-15 (土) 00:40:09
  • ftable は該当因子の全ての組合せに対するクロス集計表を作りますから、ご希望の形の出力は無理でしょう。存在する因子の組み合わせ自体を単一の因子にすれば問題はほぼ解決します。例えば次のようにする。 -- 2007-12-15 (土) 09:50:05
    # df に所属変数と名前変数を(全角空白をはさんで)結合した新しい変数を加える
    # 新しい変数のラベルは「所属及び名前」
    > df2 <- transform(df, 所属及び名前=paste(所属,名前, sep=" ")) 
    > df2
        所属     名前 月度 出席    所属及び名前
    1 営業部 営業太郎 01月 はい 営業部 営業太郎
    2 総務部 総務花子 02月 はい 総務部 総務花子
    3 営業部 営業次郎 03月 はい 営業部 営業次郎
    4 購買部 購買好男 03月 はい 購買部 購買好男
    5 経理部 経理京子 04月 はい 経理部 経理京子
    6 技術部 技術冴子 04月 はい 技術部 技術冴子
    7 営業部 営業太郎 04月 はい 営業部 営業太郎
    > ftable(df2, row.vars=5, col.vars=3) # 第5列を行、第3列を列とする集計表
                    月度 01月 02月 03月 04月
    所属及び名前                             
    営業部 営業次郎         0    0    1    0
    営業部 営業太郎         1    0    0    1
    技術部 技術冴子         0    0    0    1
    経理部 経理京子         0    0    0    1
    購買部 購買好男         0    0    1    0
    総務部 総務花子         0    1    0    0
  • それを許すなら,ftable より xtabs が好きだなあ。ftable は表頭の表示が気にくわない。 -- 2007-12-15 (土) 11:58:15
    > df2 <- transform(df, 所属及び名前=paste(所属,名前, sep=" ")) 
    > xtabs(~所属及び名前+月度, df2)
                     月度
    所属及び名前 01月 02月 03月 04月
      営業部 営業太郎    1    0    0    1
      営業部 営業次郎    0    0    1    0
      技術部 技術冴子    0    0    0    1
      経理部 経理京子    0    0    0    1
      総務部 総務花子    0    1    0    0
      購買部 購買好男    0    0    1    0
  • なるほど、こういった形でプログラムすればいいんですね。Rの振る舞いがいまひとつ判らずどういう風にコードを書けばいいかわからなくて模索しておりました。みなさんお休み中にもかかわらずありがとうございました。 -- Transalp? 2007-12-17 (月) 09:30:52

ケースの選択方法

mikasa? (2007-12-13 (木) 20:10:39)

SPSSからの転向に向けて、Rコマンダーを使って勉強中です。
(1) データセットから特定のデータだけを分析対象にする方法(SPSSで言うところのケースの選択)ですが、「データセットの部分集合を抽出」で新たなデータセットを作る以外に方法はないでしょうか?
(2) データセットを分割して一回の指示で分割したデータセットそれぞれの分析をすることは可能でしょうか?(SPSSで言うところの、ケースの分割)
よろしくお願いします。

  • ローマ数字や丸付き数字は使わない方がよいですよ。
    (1) の方は,質問の意図が分かりませんね。新たなデータセットを作る必要なんかないというか,普通は R はオンメモリで処理しますからね。
    irisデータセットのSpeciesがsetosaのものの統計量を求めたいなら,以下の一行でよし。
    summary(iris[iris[,5]=="setosa",])
    iris[,5] の3種ごとに統計量を求めたいなら,以下の一行。
    by(iris, iris[,5], summary)
    (2) の方は,iris データセットのSpeciesで3分割し,おのおのについて統計量を求めるってことが,split, lapply でできます。以下の一行でよし。
    lapply(split(iris, iris[,5]), summary)
    まあ,同じことをやるのでも,何通りもあるということで。
    Rコマンダーなど使わない方がよいかも。 -- 2007-12-13 (木) 20:42:16
  • (1) 「データセットの部分集合を抽出」しかないのでは?
    (2) 難しいかと・・・.ただ,たとえばグラフの場合は「層別」にグラフを描けるものもありますので,「層別」で代用という手もあります. -- 2007-12-14 (金) 01:20:13
  • SPSSとRではケース選択や分割の考え方(操作)が全くちがうので、同じやり方で考えると「なんだ、こんなこともできないの?」となるかもしれません。「元SPSS使い」の私は「ケースの選択」の操作自体が大嫌いだったので、この際「Rのやり方」を勉強したほうがいいと思います。 -- okinawa 2007-12-14 (金) 09:17:25

値を文字列として加えるには?

AobaAoba? (2007-12-11 (火) 17:03:50)

エクセルを卒業してRの世界に足を踏み入れたばかりのものです。
さきほどからこのサイトやらなにやらを3時間ばかり彷徨っているのですが、それでもわからず、たまらず質問させて頂きます。

やりたいことはもう少し複雑なのですが、以下のことができれば後は自分でできるはずなのです。


以下のように変数に数字を入れて、それぞれの値を代入したいのです。
x_1 <- 1
x_2 <- 2
x_3 <- 3
. . . x_10 <- 10


そこで次のようなことを書けば良いと思ったんですが、うまくいきません。

for (i in 1:10){
x_i  <- i
}


x_iのところのiの前後に何か入れるんじゃないかと思うのですが(エクセルのVBAではそうだったので)、どうすれば良いのでしょうか?
お手数おかけしてしまって本当にすみませんが、どなたかどうぞお願いいたします。

  • 何をしたいのかよくわからないのですが、 x[i]<-i じゃだめなのでしょうか。Rは変数の中にベクトル(連続した複数の値)を入れることができます。「どうしても変数名に番号を貼り付けて、その変数に値を代入したい」のであれば、eval(parse(text)=***)を使わないと動きませんね。(前にも同じような話がでたような・・・)-- okinawa 2007-12-11 (火) 17:17:53
    for (i in 1:10){
    Command<-paste("x_",i," <- ",i,sep="")
    eval(parse(text=Command))
    }
  • コメントありがとうございます。本当は、xがデータフレームになっていて、
    x[[i]]
    つまり「i番目の列のベクトル」に新たに名前をつけたいのです。なので、新しい名前に列の番号が付けられないかなと思った次第です。いただいたコメントで解決しそうです。ありがとうございました。 -- AobaAoba? 2007-12-11 (火) 17:28:54
  • どういたしまして。次回からは、実行環境(OS,Rのバージョン)を明示して投稿してください。RはWin,Mac,Linux版があり、OSやバージョンによって回答が異なることがあります。 -- okinawa 2007-12-11 (火) 17:40:18
  • 聞きたいことを素直に,詳しく説明して聞けばよいのです。「こんなことだろう」と勝手に思って全然別のことを聞いているんじゃないですか??
    データフレームの列の名前を付けたいのですよね。だったら,データフレームの名前をxとすれば,以下のようにすればいいのでは? -- 2007-12-11 (火) 17:41:34
    colnames(x) <- paste("x", 1:ncol(x), sep="_")
  • colnamesという関数があるのですね!言い訳のつもりではないのですが、そんなものの存在を知らないので、どうやったら良いかといろいろ回り道をしているところです。でも、これは良いことを聞きました。ありがとうございました。素直に質問する、実行環境を明示する件、両方とも心したいと思います。お二方ともどうもありがとうございました。 -- AobaAoba? 2007-12-11 (火) 18:10:19
  • まあ、forをやってるあたりでR魔人の先生方から「何か来そう」な感じはありましたが、これからもめげずに質問してください。 -- okinawa 2007-12-11 (火) 18:30:03
  • コメントありがとうございます。足掻いた上でどうしてものときは、また質問させて頂きます。 -- AobaAoba? 2007-12-11 (火) 18:50:54

Bray-Curtis Distanceの関数について

R.,M.? (2007-11-22 (木) 10:56:38)

初めて質問させていただきます。
昨日、クラスター解析の準備としてBray-Curtis Distanceの計算を試みたのですが以下のエラーメッセージが出て計算できません。

> dist.BC(x)
 以下にエラー dimnames(x) <- dimnames : 
   'dimnames' はリストでなくてはなりません 

そこでdimnames(x)を実行すると”行の番号”と”変数名”が正常に出てきます。ちなみにdist.BC関数を含むpackage"clusterSim"中にある他の関数(例えばdist.SM)では問題なく計算できています。また、packageのヘルプを実行したところ以下の警告が出ています。

> help(package="clusterSim")
Warning message:
In readLines(f) :
   'C:/PROGRA~1/R/R-26~1.1RC/library/clusterSim/INDEX' で不完全な最終行が見つかりました 

データは観測点における各生物分類群の現存量(対数変換済み)でデータの規模は173測点、52分類群のデータとなっています。他の関数(ユークリッド距離など)では問題なく計算できていますので、取り込んだデータの形式に問題は無いと思うのですが・・・

原因と修正方法をお教えていただけませんか。

  • dist.BC のソースコードでそういったエラーが出そうな箇所といえば wynik<-matrix(nrow=nr,ncol=nr,dimnames=names(x)) ぐらいです。names(x) の実行結果はどうなりますか. 最後の質問はヘルプドキュメントが壊れているということでしょう. もともとこわれているのか、インストール時にうっかり壊れたかは知りません. -- 2007-11-22 (木) 12:28:28
  • ありがとうございます。> names(x)を実行したところ52個の変数が正常に書き出されました([1] "変数1” "変数2" "変数3"・・・)。加えて先ほどpackageをアップデートしたのですが改善されていませんでした。 -- R.,M.? 2007-11-22 (木) 14:14:48
  • dimnames として names(x) を使うなら、後者は list(c("row1","row2",...),c("col1","col2",...)) の形であるべきですね. -- 2007-11-22 (木) 14:51:42
  • 14:51:42 さんのアドバイスにより,dist.BC の2行目を
       wynik <- matrix(nrow = nr, ncol = nr, dimnames = list(rownames(x),rownames(x)))
    とすればよいのでは? -- 2007-11-22 (木) 17:55:24
  • 毎度のことですが、この手の質問は問題点が再現できる最小限の大きさのデータを添えれば、回答者も無駄な憶測をせずにすむんですがね。 -- 2007-11-22 (木) 22:31:51

T,FとTRUE,FALSE

初心者? (2007-11-16 (金) 20:50:15)

TRUE, FALSEは上書きは不可能なんですが,それらのエイリアスであるT,Fは上書き可能で,いろいろな値が代入できてしまいます(こないだこれでドツボにはまりました).

プログラミングの際にはT,Fは使わずにTRUE,FALSEを使うのが好ましいということなんでしょうか

  • そういうことです。Rを使う上での注意事項で,必ず取り上げられる事項です。その背景がよく分かっていればT,F を使っても,一向差し支えはない。 -- 2007-11-16 (金) 21:08:26
  • 「背景がよくわかっていれば」、TRUEとしてゼロ以外の任意実数(Inf, -Inf も可)、FALSE としてゼロ(0, 0.0, 0+0i 等)を使ってもよい、のだけれど、中級以上でないとお薦めできません. -- 2007-11-16 (金) 21:45:25
  • TRUEもFALSEも補完されるので、spell outする癖をつけた方がよいかも。ちなみに、ヘルプには、'TRUE' and 'FALSE' are part of the R language, where 'T' and 'F' are global variables set to these.って書いてあるね -- 2007-11-20 (火) 16:41:53

R 言語でベクトル解析

Euler? (2007-11-16 (金) 20:31:13)

 てな、Webぺーじないかな?

  • ぐぐりなさい -- 2007-11-16 (金) 21:10:08
  • Rのシステム関数はほとんどベクトル化されていますから、それと知らずにベクトル解析していることになるんですよ(笑)。 -- 2007-11-16 (金) 21:48:39

boxとプロットの間の隙間

初心者? (2007-11-16 (金) 16:16:59)

いつもお世話になっています.

>theta <- seq(0,2*pi,0.01)
>plot(theta,sin(theta),xlim=c(0,2*pi));grid()

このようにプロットした際,boxの大きさはxlimで指定した範囲+アルファに設定され,プロットしたデータとboxの間に隙間ができてしまいます(つまり,boxの両端がx=0,x=2*piになるようにしたい).この隙間を調整するにはどのようにすれば良いのでしょうか.

  • なぜそう言うことをしたいのか。
    以下のようにすれば取りあえずの目的は達成できますか?
    bty="n" で枠を描くのを止めて,rect で自分で枠を描きます。 -- 2007-11-16 (金) 16:25:16
    theta <- seq(0,2*pi,0.01)
    plot(theta,sin(theta),type="l", bty="n")
    rect(0, -1, 2*pi, 1)
  • ある行列データをimageで,行列データの各列の積分値のデータをplotで書き,それらのプロットをx軸を共通にして重ねたかったのですが,plot側のx軸のスケールが変わってしまい困っていました.plot(...,xaxis="i")とすることで解決できました.ありがとうございました -- 初心者? 2007-11-16 (金) 16:45:33
  • 少し不正確では?
    > theta <- seq(0, 2*pi, 0.01)
    > oldpar <- par(xaxs="i")
    > plot(theta, sin(theta), xlim=c(0, 2*pi), type="l")
    > par(oldpar)
    xaxsはpar()の引数です. -- 2007-12-12 (水) 10:46:53
  • > xaxsはpar()の引数です.
    その通りではあるが,par() の引数の多くは plot の引数として与えることができますよ。
    plot のオンラインヘルプを見ればわかりますが,plot が ... 引数を持ち,その説明で
    ... Arguments to be passed to methods, such as graphical parameters(see par).
    実際に,
    > plot(theta, sin(theta), xlim=c(0, 2*pi), type="l", xaxs="i")
    > plot(theta, sin(theta), xlim=c(0, 2*pi), type="l")
    を試してみれば,違いがわかるでしょう。 -- 2007-12-12 (水) 12:59:34
  • なるほど.こんな便利な方法があったのですね.ありがとうございます. -- 2007-12-13 (木) 09:05:54

インストーラなしのWin版Rのインストール

内部統制? (2007-11-14 (水) 10:30:43)

 Win版Rをインストーラを使わないでインストール方法はあるでしょうか?

  • ソースから自分でコンパイルすれば?(なんでそいういうやり方をしないといけないの?) -- 2007-11-14 (水) 10:47:03
  • どこぞのWinの[Program Files]の中の[R]フォルダを丸ごとコピーする。(レジストリに依存しないパッケージのみ) -- okinawa 2007-11-14 (水) 13:39:02
  • Rjpwikiの過去記事にUSBでWin版Rを持ち歩く方法があります。これはinstall不要だったと思います。あとは、内部統制さんがご自分で探して下さい。それぐらいはできるでしょ? -- 2007-11-14 (水) 13:41:43

2次輸送問題

OR mania? (2007-11-13 (火) 11:14:42)

 Rには2次輸送問題(QTP:Quadratic Transportation Problem)を解く関数はあるのでしょうか?

 こちらの「数理計画」のページを見ましたが見当たりませんでした。

barplotのlegendについて

K? (2007-11-12 (月) 14:05:32)

barplotのlegendについて教えて下さい。
グラフをプレゼン用に使おうと、par(ps=22)として文字サイズを大きくしました。

par(ps=22)
par(mar=c(5,5,3,3))
barplot(sum, names.arg=breaks2,
        col = c("darkred", "darkgreen", "lightblue", "pink"),
        xlab="胸高直径階(cm)",  ylab="立木本数(本/ha)", main="")
legend("topright", c("アカマツ", "ヒノキ", "スギ", "その他" ),
       fill=c("darkred", "darkgreen", "lightblue", "pink"))

そして、以下のようなヒストグラムを得ました。

histgram_k2.png

この凡例中の文字列の前の四角のサイズを大きくしたいと思い、legendのマニュアルを見たのですが方法が分かりませんでした。どなたか方法をご教授いただければ幸いです。よろしくお願いします。

  • 制御する引数がないので,legend 関数の中の,rect2 の呼び出しの所をいじるしかないかなあ -- 2007-11-12 (月) 14:54:40
  • legend 関数の lty, lwd パラメータを適当に指定してみたらそれらしくなりませんか。 -- 2007-11-12 (月) 16:38:45
  • なるほど,fill の代わりに col を使い(fill で指定しているものを col で指定) と lwdを大きめにして,y.intersp を少しあければなんとか見られますね。 -- 2007-11-12 (月) 17:00:29
  • 早速のコメント、ありがとうございます。 頂いたアドバイスを参考に、
    legend("topright", c("アカマツ", "ヒノキ", "スギ", "その他" ),
            col=c("darkred", "darkgreen", "lightblue", "pink"), 
            lwd=15, y.intersp=2)
    とすることで、以下のような凡例を得る事が出来ました。
    legend_K1.png
    四角から角の丸い太線になったのは、fillの代わりにcolを用いたためですよね。出来れば角の丸い太線ではなく、四角く見えるようにしたいのですが、無理でしょうか。 ともあれ、お陰様でPowerPoint?でのプレゼンなどには十分使えるグラフを得る事が出来ました。どうもありがとうございました。 -- K? 2007-11-13 (火) 09:45:21
  • おおっ!アドバイス通りにやってみたら、見事に欲しかったグラフが得られました。Rではこんな技も使えるのですね。勉強になります。助かりました。どうもありがとうございました。 -- K? 2007-11-14 (水) 13:49:49

Mac OS X 10.5にRのインストールができません

suzuk? (2007-11-11 (日) 11:37:40)

(Rがインストールされていない)Mac OS X 10.5 (Leopard)にRのインストールを試みましたが、以下のメッセージが表示されてしまいます。 「このソフトウェアのより新しいバージョンがシステムにインストールされています。この古いバージョンをインストールする必要はありません。」 R-2.6.0、R-2.5.0、R-2.4.1の何れのバージョンにおいても同様の現象が認められます。アドバイスをよろしくお願いします。

  • CRANのR for Mac OS Xというページhttp://cran.r-project.org/bin/macosx/から.dmg形式のインストーラーをとってきて、インストールしましたが、とくに問題もなくインストールできました。ちなみにダウンロードしたのはこれ。R-2.6-branch-43408.dmg (latest version) -- fujiwara? 2007-11-11 (日) 16:40:09
  • R-2.6-branch-43408.dmg (latest version)はLeopardに正常にインストールできました。ありがとうございました。 -- suzuk? 2007-11-11 (日) 18:39:34

キュー、スタックの実現

kbi? (2007-11-11 (日) 09:38:12)

google等の検索エンジンとサイト内検索を試したのですが
有効な情報が見つからなかったので質問させてください。
Rにおいて、基本的なデータ型としてキューやスタックは
用意されていないのでしょうか。
待ち行列などのモデルをシミュレーションしてみたいと思っています。
調べるのに有効なキーワードだけでも教えていただけると幸いです。

  • R本体にはキューやスタックという機能はありません。貢献パッケージにもなさそうです.困ったときは google よりもトップページから辿れるR専門の検索エンジン Rsitesearch を使いましょう。例えばキーワード "stack pop push" で検索すると次のようなスタック関数がMLに投稿されていました.要するにリスト構造を使えばキューもスタックも単純なリスト操作です。 -- 2007-11-11 (日) 12:42:17
    > Stack <- function(){  # Rには既に stack という名前の関数があるので注意
    +  it <- list()
    +  res <- list(
    +              push=function(x){
    +                it[[length(it)+1]] <<- x
    +              },
    +              pop=function(){
    +                val <- it[[length(it)]]
    +                it <<- it[-length(it)]
    +                return(val)
    +              },
    +              value=function(){
    +                return(it)
    +              }
    +              )
    +  class(res) <- "stack"
    +  res
    +}
    > print.stack <- function(x,...) print(x$value())
    > push <- function(stack,obj) stack$push(obj)
    > pop <- function(stack) stack$pop()
    > bz <- Stack() 
    > bz
    list(0)
    > push(bz, 1:10)
    > bz
    [[1]]
     [1]  1  2  3  4  5  6  7  8  9 10
    > push(bz, rnorm(10))
    > bz
    [[1]]
     [1]  1  2  3  4  5  6  7  8  9 10
    [[2]]
     [1]  1.1519118  0.9921604 -0.4295131  1.2383041 -0.2793463  1.7579031
     [7]  0.5607461 -0.4527840 -0.8320433 -1.1665705
    > pop(bz)
     [1]  1.1519118  0.9921604 -0.4295131  1.2383041 -0.2793463  1.7579031
     [7]  0.5607461 -0.4527840 -0.8320433 -1.1665705
    > bz
    [[1]]
     [1]  1  2  3  4  5  6  7  8  9 10
    > pop(bz)
     [1]  1  2  3  4  5  6  7  8  9 10
    > bz
    list()
  • なるほど、これは使いやすそうです。自力でリストから構造を作ればいいのですね。検索方法も分かりました。ありがとうございました。 -- kbi? 2007-11-11 (日) 21:53:55

パッケージのインストールができない

tomo? (2007-11-09 (金) 20:27:05)

Rの初心者です。 R2.4.1を使用しています。
パッケージのインストールができません。ダウンロードはしていますが,解凍ができていないようです。

> utils:::menuInstallPkgs()
URL 'http://prs.ism.ac.jp/bioc/1.9/bioc/bin/windows/contrib/2.4/Heatplus_1.4.0.zip' を試しています
Content type 'application/zip' length 436881 bytes
開かれた URL
downloaded 426Kb
~
以下にエラーzip.unpack(pkg, tmpDir) : 
ファイル 'C:/Program  Files/R/R-2.4.1/library/file57d3669a/Heatplus/chtml/Heatplus.chm' を開くことができません

違うパソコンからだと,同じ操作でインストールできました。パソコンの設定が悪いのでしょうか。 よろしくお願いいたします。

  • まずは,最新バージョンにバージョンアップしましょう。うまくいく行かないは,それからのお話です。バージョンのふるい状態でなんだかんだいっても,普遍性がないです。それと,「違うパソコン」との違いは(ああ,おもしろい)なんなんでしょうか? -- 2007-11-09 (金) 22:31:31
  • Vista の場合は,R を起動する際に「右クリック」→「管理者として実行」で起動してください. -- 2007-11-10 (土) 01:07:48
  • 皆さま ありがとうございました。違うパソコン=XP できないパソコン=Vista と言うのをお伝え忘れていましたが,適切なご助言,ありがとうございました。無事,インストールできました。 -- tomo? 2007-11-10 (土) 12:22:42

ウィンドウサイズの拡大

クラスタ? (2007-11-07 (水) 20:30:56)

hclust()を使い樹形図を作成しました。
しかし、データ数が600個あり文字が重なって読めません。
ウィンドウサイズを拡大しようにも画面サイズより広がりませんでした。
なんとか読めるように文字が重ならない大きさまで拡大したいです。
ご教示お願いします。

  • 初級Q&A アーカイブ(6)の『クラスター分析におけるフォントサイズ変更方法について』を見た上で書いてる? -- 2007-11-07 (水) 21:51:46
  • 読んでませんでした。同じ問題を抱えた人の質問でしたが、フォントサイズではなくウィンドウサイズの拡張限界の変更方法は、ないでしょうか? -- クラスタ? 2007-11-07 (水) 23:16:38
  • 600個のデータを一度に見よう(一部分だけ見てもそれはそれで意味不明になる可能性大)というのが無茶なのでは。example(dendrogram) で見られるような文字グラフィックスで眺めてみたらどうでしょう。 -- 2007-11-08 (木) 00:12:34
  • こんな表示方法もあるんですね、ウィンドウの拡張は諦めて、この方法で試行してみます。ありがとうございます。 -- クラスタ? 2007-11-08 (木) 00:33:27
  • これと同じ質問にここ数週間内に答えたんだけど,元記事が消去されたのか,なくなっていますね。
    それはさておき,画像を画面一杯に限って重ならないようにフォントを小さくしても,今度は,フォントが読めなくなるだけ。アラース。
    画面に描かれたものを見るだけでよいという局面は少ないわけで,何らかの形でプリントするのが最終形?とすると,なにも,グラフのサイズはコンピュータの画面のドット数には限定しなくても良いということ。
    いったん画像ファイルに保存すれば,それをそのまま印刷するには巨大すぎても,そのファイルを画像処理プログラム(フォトショップとかその他なんでも)で読み込めば,部分を拡大してみることもできるわけ。普通の大きさのフォントを使っても,重ならないように大きな描画領域をとればよいだけ。
    つまり,横幅(縦幅も)を十分大きく取って描画し,それをファイルに落とすようにしておけばよいだけ。
    以下を実行して,できるファイルをご覧下さいませ。 -- 2007-11-08 (木) 00:54:04
    x <- matrix(rnorm(1200),600)
    pdf("cluster.pdf", height=10, width=150)
    plot(hclust(dist(x)), hang=-1)
    dev.off()
  • おお、できました。これを求めていたんです。お二方ありがとうございました。 -- クラスタ? 2007-11-08 (木) 01:40:35

関数 cut() の引数 right について

実験万? (2007-11-01 (木) 14:40:22)

ある数値データベクトルの度数分布を求めるときに、階級の端点に一致する数値がある場合、その数値を、下側と上側、どちらの階級に含めるべきか、という問題を検討しています。そこで、次のような実験をしてみました。

まず、0.1 から 9.9 まで、0.1 刻み、99 個の数値からなるデータベクトルを用意します。このデータベクトルについて、最初の階級の下限 0、最後の階級の上限 10、階級幅 1 の条件で度数分布を求めます。

実験 (1)
階級の端点に一致する数値がある場合、その数値を「下側」の階級に含めることにすれば、次のような結果を得ます。

命令:
x <- seq( 0.1, 9.9, 0.1 )
table( cut( x, breaks=seq(0,10,1), right=TRUE ) )
結果:
(0,1]  (1,2]  (2,3]  (3,4]  (4,5]  (5,6]  (6,7]  (7,8]  (8,9] (9,10]
   10     10      9     11     10     10     10     10     10      9


実験 (2)
階級の端点に一致する数値がある場合、その数値を「上側」の階級に含めることにすれば、次のような結果を得ます。

命令:
x <- seq( 0.1, 9.9, 0.1 )
table( cut( x, breaks=seq(0,10,1), right=FALSE ) )
結果:
[0,1)  [1,2)  [2,3)  [3,4)  [4,5)  [5,6)  [6,7)  [7,8)  [8,9) [9,10)
    9     10     10     10     10     10     10     10     10     10


それぞれ、最初と最後の階級の度数をみると、期待したとおりの結果になっています。

しかし、ひとつ問題が見つかりました。それは、実験 (1) の結果で、(2,3] の度数が 9、(3,4] の度数が 11 であることです。本当は、それぞれ 10 であるはずだと思うのです。このような結果になる原因について、ご教示いただけませんでしょうか。

  • 計算機で扱う実数は,「近似値に過ぎない」からです。近似値でないのは,整数と,小数部が1/2のべき乗の和で表される数値(たとえば 0.125 とか0.07670455など)の場合のみです。0.1 などは,2 進数で表すと循環小数になりますので,正確な値を保持できるわけがないということです。
    あなたの例でも,全体を10倍して小数を含まないデータにしてみると以下のようになります。
    要するに,正確な値が欲しければ整数を扱うようにするのがよいと言うことです。
    ちなみに,日本では「○○以上□□未満」ということがあるので,左側の級限界を含むのが常識になっています。 -- 2007-11-01 (木) 14:48:42
    > x <- seq( 1, 99, 1 )
    > table( cut( x, breaks=seq(0,100,10), right=TRUE ) )
    
      (0,10]  (10,20]  (20,30]  (30,40]  (40,50]  (50,60]  (60,70]  (70,80] 
          10       10       10       10       10       10       10       10 
     (80,90] (90,100] 
          10        9
  • ありがとうございます。勉強になります。これからは、必要に応じて、10^n を掛け算して、整数にしてから、処理することにします。 -- 実験万? 2007-11-01 (木) 15:08:42

0.999999999

ビットマップ? (2007-11-01 (木) 12:18:08)

R version 2.5.1です。

set.seed(100)
vo <- round( rnorm(40,100,5), 1)


table(diff(sort(vo)))

とすると

                0 0.0999999999999943  0.100000000000009  0.200000000000003 
                4                  8                  5                  1 
0.299999999999997  0.300000000000011  0.399999999999991  0.400000000000006 
                5                  1                  1                  1 
              0.5  0.599999999999994  0.700000000000003                  1 
                7                  1                  1                  1 
 1.20000000000000              2.800                3.5 
                1                  1                  1 


という結果が得られました。
0.0999999999999943 と 0.100000000000009 はどちらも 0.1 と扱ってほしいのですが、そうするにはどうしたらよいですか?

  • 必要な段階でroundすること。良いか悪いかはともかく,table(round(diff(sort(vo)),1)) とした結果をみてみればいかが?
    この質問も,この上の質問と解答に関連しています。 -- 2007-11-01 (木) 14:42:36
  • そうですね.round でゴリおししてみます.ただ,何桁必要かがデータによってマチマチなので,それをユーザーが指定するオプションにしないといけないことになりそうです. -- ピットマッブ? 2007-11-02 (金) 03:09:53

グラフの軸を太くしたい

surg? (2007-10-28 (日) 07:40:31)

Rで作成したグラフを、Microsoft Wordの書類に挿入しようと思っております.WindowsとMacintoshの双方で利用したいため,グラフを大きめのbitmapで作成し,Word上で縮小表示しようとしております.

png("test.png", width=324*4, height=324*4)
par(cex=4, lwd=4)
hist(rnorm(1000))
dev.off()

上記のようにしてWord上で15%程度に縮小して印刷すると,グラフの軸だけがかすれてしまいます.軸を太くしたいのですが,何か方法はありますでしょうか? help(par) は読んだのですが,それらしい項目は見当たりません.

  • lwd -- okinawa 2007-10-28 (日) 10:08:41
  • par じゃなくて hist の引数でね -- 2007-10-28 (日) 10:34:03
  • ありがとうございます.ご指摘のようにlwdをparではなくhistの引数で渡したところ,希望どおりの動作となりました. -- surg? 2007-10-28 (日) 13:58:01
  • histはこれで大丈夫だったのですが,plotではだめでした.
    png("test.png", width=324*4, height=324*4)
    par(cex=4, lwd=4)
    plot(rnorm(1000), rnorm(1000), lwd=4)
    dev.off()
    こうするとグラフのほとんどの構成要素は太くなるのですが,「軸の目盛り」だけが太くなりません. help(plot), help(plot.default)にはそれらしい記述は見当たりませんでした.-- surg? 2007-10-28 (日) 15:15:32
  • 実際に呼び出される plot.default に問題があるように思います
    取りあえずは,plot.default にパッチを入れて使うとか -- 2007-10-29 (月) 11:03:54
    plot.default <- function (x, y = NULL, type = "p", xlim = NULL, ylim = NULL, 
        log = "", main = NULL, sub = NULL, xlab = NULL, ylab = NULL, 
        ann = par("ann"), axes = TRUE, frame.plot = axes, panel.first = NULL, 
        panel.last = NULL, asp = NA, ...) 
    {
    # 実際は以下のようになっているのだが,これでは lwd は Axis に渡らないのでは?
    #    localAxis <- function(..., col, bg, pch, cex, lty, lwd) Axis(...)
    # 次に示すように lwd を取ると,期待通りの動きをする
    # col とか bg, pch, cex ,lty も同じになると思うが,
    # この関数を作った人は,plot の軸を描くときに col, bg, pch, cex, lty は
    # 弄るべきではないと思ったのかもね
         localAxis <- function(...) Axis(...)
  • ありがとうございます.plot.defaultに手を入れることで,無事目的を果たせました.
    仰るとおり,作者は他のパラメータと同じく,Axis()のlwdもpar()で制御するつもりだったのでしょうね. -- surg? 2007-10-29 (月) 12:35:03

PLSとCCAは同等か

yishii? (2007-10-24 (水) 23:04:38)

初心者なので教えていただきたいことがあり、投稿させていただきます。あるケモメトリックスの入門書に、PLS(Partial least squres)法と正準相関解析(Canonical correlation alalysis,CCA)が同等だと書いてあったのですが、それは正しいのでしょうか。もしそうだとすれば、この両者の関係を論じた文献があれば、ご教示いただければ幸いです。よろしくお願いします。

  • 統計プロパーの質問はこのWikiでは歓迎されません。 -- 2007-10-25 (木) 00:07:54

R commanderでの生存時間解析(cmprsk)における多変量解析(crr)について教えてください

takeiteasy? (2007-10-24 (水) 00:22:52)

先日、Rcommanderでのcmprskを教えていただいたものです。おかげさまで、cumincに関しては、Rcommanderで解析可能になりました。ありがとうございます。たびたび申し訳ありませんが、cmprskにおいての多変量解析であるcrr(competing risks regression analysis)をRcommanderに機能実装しようと試みましたが、うまくいきません。
Rcommanderハンドブックの6章のCox回帰を参考にやってみましたが、根本的な理解が不足していて、できません。
申し訳ありません。

  • Rでも動かないのでしょうか?何を試して、何がダメだったか、説明して下さい。
    takeiteasy?さんはコマンド入力でRを使用されているようですので、そのスクリプトをRcmdrに組み込めば済むような気がします。私はRcmdrもcmprskも使わないので、トンチンカンならすみません。実装できたら教えてくださいね。 -- akira? 2007-10-24 (水) 07:21:23
  • crr のヘルプを見ますと、crr のモデルの指定の仕方は「Cox回帰」のようなモデル指定ではありませんねぇ。。。
    ↓で定義しました関数 Mycuminc() を改造してみてはいかがでしょうか?? -- 舟尾? 2007-10-24 (水) 22:52:31
  • ありがとうございます。Mycuminc()を改造してやってみます。 -- takeiteasy? 2007-10-25 (木) 19:29:06

e1071パッケージのsvmメソッドの損失関数について

codedir? (2007-10-19 (金) 12:51:21)

e1071パッケージのsvmメソッドの回帰における損失関数はε-インセンシティブ損失関数で線形だと思うのですが、この損失関数を線形だけでなく2次関数に変更したりできるように、プログラムを触りたいんです
損失関数を定義しているプログラムを探したんですが、見当がつきません。プログラム内のepsilonが出てくる部分を見てみても、epsilonを使って計算しているようには見えません。
よろしくお願いします

e1071:::svm.default
e1071:::svm.formula
  • svmtrain って,Cで書かれたプログラムを呼んでるんでは?
    そのソースプログラムを読めばよいかと。 -- 2007-10-19 (金) 17:36:08
    cret <- .C("svmtrain", as.double(if (sparse) x@ra else t(x)), 中略, as.double(epsilon), 後略
  • ありがとうございます!google code searchで検索したら、e1071/src/Rsvm.c がヒットしました。おそらくこの中のsvmtrainだと思います
    通常このソースはe1071のどこに保存されているのでしょうか?e1071のフォルダを一通り調べたのですがわかりません。どうやって調べたらわかるのでしょうか?よろしくお願いします -- codedir? 2007-10-19 (金) 22:27:15
  • なんだかなあ。? svm してでてくる中で,参考文献として出ているウェッブページなんかを探せば,ソースファイルがあるんじゃないかと(底にはないかも知れないけど,どこかには公開されているはず)。私はそんな熱意がある訳じゃないから苦労して探そう何滴はないけど,あなたほどの熱意があるなら,簡単にみつけられるんじゃ? -- 2007-10-20 (土) 00:18:26
  • 近くにある、「Rに関して包括的な情報をまとめたサイト」 = CRANへ行って、探してみましょう。windows用のバイナリでけではなく、ソースコードが置いてあります。 -- ジェットラグ? 2007-10-20 (土) 06:56:53
  • お返事が大変おそくなって申し訳ありません。ありがとうございました! -- codedir? 2007-11-3 (土) 17:32:49
  • すみません、追加の質問になるんですがお願いします。e1071のソースをダウンロードしたら、srcフォルダにRsvm.c、svm.cpp、svm.hがありました。e1071に実装されているのはLibsvmだそうで、Libsvmのページに行き、損失関数をいじるためにソースを修正するところもわかりました。また「他言語からの利用」を参考にして、CのソースをRで読めるようにするところまでたどりつきました。e1071のsrcの中にあったsvm.cppとsvm.hから、損失関数を書き換えたsvm2.dllをつくり、svmメソッドを実際に走らせると「C シンボル名 "svmtrain" はパッケージ "svm2" の DLL の中にありません」となってしまいます。(svm2は自作のパッケージです)Rsvm.cの中には"svmtrain"があるんですが、計算のソースであるようには見えませんし、Rsvm.cではDLLファイルが作れません。svm.cppの中には計算部分があるように思えますが、"svmtrain"というルーチンがみあたりません。どこを改良したら、C部分の損失関数をいじったRのプログラムをはしらせれるようになるのでしょうか?すみませんが、お願いします-- codedir? 2007-11-17 (土) 17:36:39
  • 自己解決しました!svm.dllをビルドする時にsvm.hとsvm.cppとRsvm.cをビルドすればいいだけの話でした。 -- codedir? 2007-11-20 (火) 14:39:34

Rcommanderでcmprskはできますか?

takeiteasy? (2007-10-17 (水) 23:46:49)

医学データの解析に使用したいのですが、
いままで、JMPなどで統計解析していましたが、
cmprskはJMPなどでは解析できません。
コマンド入力などは素人でなかなか難しいです。
Rcommanderなどを用いてcmprskは解析できるでしょうか?

  • R には cmprsk というパッケージがあるものの、Rcommander 経由で使えるかどうかは知りません(おそらく使えないでしょう)。 -- 2007-10-17 (水) 23:54:37
  • R commander では cmprsk を含め,生存時間解析の機能はサポートされておりません
    (R Commander に cmprsk の機能を実装すれば解析できます^^). -- 舟尾? 2007-10-18 (木) 01:05:46
  • ありがとうございます。コマンド入力でがんばってみます。 -- takeiteasy? 2007-10-18 (木) 01:28:56
  • 言いだしっぺの法則….パッケージ Rcmdr の etc フォルダに「MyProgram?.R」を作成してください.
 ### MyProgram.R の中身
 Mycuminc <- function(){
  initializeDialog(title="cumulative incidence function")
  timeBox   <- variableListBox(top, Numeric(), title="時間変数")
  eventBox  <- variableListBox(top, Numeric(), title="イベント/打ち切り変数")
  groupBox  <- variableListBox(top, Factors(), title="群変数")
  Var1      <- tclVar("0")
  Var1Entry <- tkentry(top, width="6", textvariable=Var1)
  onOK <- function(){
    XXX <- getSelection(timeBox)
    YYY <- getSelection(eventBox)
    ZZZ <- as.factor(getSelection(groupBox))
    AAA <- as.numeric(tclvalue(Var1))
    if (length(XXX) != 1 || length(YYY) != 1 || length(ZZZ) != 1) {
      errorCondition(recall=Mycuminc, message="変数を1つずつ選択してください。")
      return()
    }
    closeDialog()
    library(cmprsk)
    if (is.na(AAA)) {
      command <- paste("cuminc(", 
        ActiveDataSet(), "$", XXX, ", ",
        ActiveDataSet(), "$", YYY, ", ",
        ActiveDataSet(), "$", ZZZ, ", ",
        "cencode=0)", sep="")
    }
    else {
      command <- paste("cuminc(", 
        ActiveDataSet(), "$", XXX, ", ",
        ActiveDataSet(), "$", YYY, ", ",
        ActiveDataSet(), "$", ZZZ, ", ",
        "cencode=",           AAA, ")", sep="")
    }
    doItAndPrint(command)
    tkfocus(CommanderWindow())
  }
  OKCancelHelp(helpSubject="cuminc")
  tkgrid(getFrame(timeBox), getFrame(eventBox), getFrame(groupBox), sticky="w")
  tkgrid(tklabel(top, text=gettextRcmdr("打ち切り")), Var1Entry, sticky="e")
  tkgrid(buttonsFrame, columnspan=2, sticky="w")
  tkgrid.configure(Var1Entry, sticky="w")
  dialogSuffix(rows=2, columns=1)
}

Mycifplot <- function(){
  initializeDialog(title="cumulative incidence function")
  timeBox   <- variableListBox(top, Numeric(), title="時間変数")
  eventBox  <- variableListBox(top, Numeric(), title="イベント/打ち切り変数")
  groupBox  <- variableListBox(top, Factors(), title="群変数")
  Var1      <- tclVar("0")
  Var1Entry <- tkentry(top, width="6", textvariable=Var1)
  onOK <- function(){
    XXX <- getSelection(timeBox)
    YYY <- getSelection(eventBox)
    ZZZ <- as.factor(getSelection(groupBox))
    AAA <- as.numeric(tclvalue(Var1))
    if (length(XXX) != 1 || length(YYY) != 1 || length(ZZZ) != 1) {
      errorCondition(recall=Mycuminc, message="変数を1つずつ選択してください。")
      return()
    }
    closeDialog()
    library(cmprsk)
    if (is.na(AAA)) {
      command <- paste("tmp <- cuminc(", 
        ActiveDataSet(), "$", XXX, ", ",
        ActiveDataSet(), "$", YYY, ", ",
        ActiveDataSet(), "$", ZZZ, ", ",
        "cencode=0); plot(tmp)", sep="")
    }
    else {
      command <- paste("tmp <- cuminc(", 
        ActiveDataSet(), "$", XXX, ", ",
        ActiveDataSet(), "$", YYY, ", ",
        ActiveDataSet(), "$", ZZZ, ", ",
        "cencode=",           AAA, "); plot(tmp)", sep="")
    }
    doItAndPrint(command)
    tkfocus(CommanderWindow())
  }
  OKCancelHelp(helpSubject="cuminc")
  tkgrid(getFrame(timeBox), getFrame(eventBox), getFrame(groupBox), sticky="w")
  tkgrid(tklabel(top, text=gettextRcmdr("打ち切り")), Var1Entry, sticky="e")
  tkgrid(buttonsFrame, columnspan=2, sticky="w")
  tkgrid.configure(Var1Entry, sticky="w")
  dialogSuffix(rows=2, columns=1)
}

MySetup <- function(){
  initializeDialog(title="cmprsk のセットアップ")
  onOK <- function(){
    closeDialog()
    install.packages('cmprsk')
  }
  OKCancelHelp(helpSubject="install.packages")
  tkgrid(buttonsFrame, columnspan=1, sticky="w")
  dialogSuffix(rows=1, columns=1)
}

次に、同フォルダにある Rcmdr-menus.txt の一番最後に以下を追記して下さい.

### Rcmdr-menus.txt の一番最後に追記
### メニュー追加・ここから
   menu  MyMenu   topMenu  ""                       ""         ""  ""
   item  topMenu  cascade  "メニュー"               MyMenu     ""  ""
   item  MyMenu   command  "CIF(計算)"            Mycuminc   ""  ""
   item  MyMenu   command  "CIF(プロット)"        Mycifplot  ""  ""
   item  MyMenu   command  "cmprsk のセットアップ"  MySetup    ""  ""
### メニュー追加・ここまで

まずは R commander を起動し,「メニュー」→「cmprsk のセットアップ」を行った後,
サンプルデータを R Commander でアクティブにしてから,「メニュー」→「CIF(計算)」などを実行してください.

### サンプルデータ
x <- data.frame(time=c(1:10,1:10), censor=floor(3*runif(20)), group=c(rep("A",10),rep("B",10)))

R Commander への機能追加については 説明スライドの最後の方R Commander本の6章を参照ください.ちなみに,コマンドでやられる場合は・・・

x <- data.frame(time=c(1:10,1:10), censor=floor(3*runif(20)), group=c(rep("A",10),rep("B",10)))
library(cmprsk)
( tmp <- cuminc(x$time, x$censor, x$group) )
plot(tmp)
cif.gif

です.-- 舟尾? 2007-10-18 (木) 02:02:03

  • 舟尾先生、詳細な説明ありがとうございます。やってみたいと思います。 -- takeiteasy? 2007-10-18 (木) 11:48:47

音声データをスペクトル解析したい

trantila? (2007-10-16 (火) 18:28:30)

音声データ(WAV形式です。)をスペクトル解析したいのですが、どのようにしたら良いのでしょうか?
ここのサイトのスペクトル解析のところは見たのですがいまひとつよく分かりませんでした。
Rでは音声データを解析対象として扱えないのでしょうか?
もし、そうであれば音声データを時系列データに変換する方法が必要なのでしょうか?

よろしくお願いします。

RESTFul な Web Service で Rを利用

tantan? (2007-10-15 (月) 11:48:37)

RWeb なような対話型の web service ではなくて、CGI の Get や Put などで、入力パラメータを指定して、サーバサイトで R を実行し、戻り値を XML 形式にして出力するような Web アプリってないでしょうか?

  • おれもほしい!(言いだしっぺの法則は無しにして) -- okinawa 2007-10-15 (月) 13:21:39
  • RApache+XMLでなんとかなるんでは? -- なかま 2007-10-15 (月) 17:56:50
  • RApache は Windows 用はないようですね。 -- 2007-10-15 (月) 19:21:24
  • これは? http://data.vanderbilt.edu/~koyamat/brew/magicsquare.brew -- コインランドリー? 2007-10-16 (火) 07:56:23

RからWSDLを読み込みたい

akira? (2007-10-14 (日) 17:17:22)

なかまさんのRからのBBDJアクセスメモを試したところ、エラーがでてしまいました。proxyは変更していません(<-これか?)。実のところ、SSOAPがよくわかりません。

> library("SSOAP")
 要求されたパッケージ XML をロード中です 
 要求されたパッケージ RCurl をロード中です 
> ddbj <- processWSDL ("http://xml.nig.ac.jp/wsdl/GetEntry.wsdl")
> iface <- genSOAPClientInterface(def=ddbj)
 以下にエラー parse(text = txt) :   予想外の 入力  です  以下中: 
 "       server = .defaultServer, .convert = .operation@returnValue, .opts = list(), ..., 
         nameSpaces =  &#65533;" 

とエラーになってしまいました。
関数を読み込むと、

> processWSDL
function (fileName = "KEGG.wsdl", handlers = WSDLParseHandlers(fileName), nameSpaces = character()) {
    library(XML)
(中略)
    types = processWSDLTypes(doc[["types"]], doc)
(以下、省略)

と、"processWSDLTypes"を見つけられず、トレースできませんでした。 で、exampleを確認すると、

>  tmp = processWSDL(system.file("examples", "KEGG.wsdl", package = "SSOAP"))
Warning messages:
1: In function (node)  :
  skipping import node with no schemaLocation attribute
2: In function (node)  :
  skipping import node with no schemaLocation attribute

となり、KEGG.wsdlの読み込みでエラーになったようで、

> tmp = processWSDL("http://soap.genome.jp/KEGG.wsdl")
Warning messages:
1: In function (node)  :
  skipping import node with no schemaLocation attribute
2: In function (node)  :
  skipping import node with no schemaLocation attribute

と、大元にアクセスしてもだめでした。
一つのスレで2つの疑問が生じていますが、

ddbj <- processWSDL ("http://xml.nig.ac.jp/wsdl/GetEntry.wsdl")
tmp = processWSDL("http://soap.genome.jp/KEGG.wsdl")

のエラーを回避したいです。

> sessionInfo()
R version 2.6.0 (2007-10-03) 
i486-pc-linux-gnu 
locale: LC_CTYPE=ja_JP.UTF-8;LC_NUMERIC=C;LC_TIME=ja_JP.UTF-8;
LC_COLLATE=ja_JP.UTF-8;LC_MONETARY=ja_JP.UTF-8;LC_MESSAGES=ja_JP.UTF-8;
LC_PAPER=ja_JP.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;
LC_MEASUREMENT=ja_JP.UTF-8;LC_IDENTIFICATION=C
attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     
other attached packages:
[1] SSOAP_0.4-6 RCurl_0.8-3 XML_1.93-2 
loaded via a namespace (and not attached):
[1] rcompgen_0.1-15

です。

  • SSOAPの0.4-3をお使いください. 直しかけましたが, 挫折しました... -- なかま 2007-10-14 (日) 20:52:56
  • ありがとうございます。DDBJにアクセスできました!
    ただ、SSOAPの0.4-3でも、example(processWSDL)はこけますね。これはexampleにあるkegg.wsdlに"schemaLocation attribute"がないっていうエラーなんでしょうか? SAXとか、DOMとか、よくわからないんです...スミマセン
    RSiteSearch?("processWSDL")すると、"http://www.ncbi.nlm.nih.gov/entrez/eutils/soap/eutils.wsdl"でもこけてるみたいですね。
    > ncbi <- processWSDL("http://www.ncbi.nlm.nih.gov/entrez/eutils/soap/eutils.wsdl") 
     以下にエラー ans[[ns]] <<- o :  一つ未満の要素を選択しようとしました 
     追加情報:   50 件以上の警告がありました (警告を見るには warnings() を使って下さい) 
    > warnings()
     警告メッセージ: 
    1: In FUN(X[[2L]], ...) :
      Failed to handle node HasLinkOut of type attribute&simpleType in processWSDLType
    2: In FUN(X[[2L]], ...) :
      Failed to handle node HasNeighbor of type attribute&simpleType in processWSDLType
    3: In FUN(X[[1L]], ...) :
    (以下、略) 
    うーん、DDBJだけじゃなく、KEGG APIもeUtilsも使いたかったのですが... -- akira? 2007-10-14 (日) 21:18:15
  • 最新のSSOAP_0.4-6でも,options(useFancyQuotes?=FALSE)で`ddbj'なら動きます(dQuoteがエンコードにより異なる仕様になったのが原因でした). ちょっと, これに取り組む余力はありませんです. -- なかま 2007-10-15 (月) 11:18:18
  • ありがとうございます。今はWinXPで使っています。SSOAP_0.4-4ならoption指定無しでddbjは動くようです。 -- akira? 2007-10-15 (月) 12:30:17

グラフのフォントサイズを設定したい

田中? (2007-10-14 (日) 15:41:02)

R2.6.0,Windows XPを使っております。
barplot()を使って棒グラフを描き,PowerPoint?のスライドに貼り付けてプレゼンをしようと考えています。
PowerPoint?に貼り付けてみると,文字が少し小さく,フロアからは見えにくいようです。
x軸,y軸のラベルや図のタイトルのフォントサイズを変更することはできるでしょうか。
どなたか対処法をご存知の方がいらっしゃいましたら,お返事をどうぞよろしくお願い致します。

  • この手の質問は過去の質問履歴から検索するか、グラフィックパラメータが吉。きっとお望みの結果が得られると思いますよ。下の質問と同じですが、調べたこと、でも分からなかったこと、を書いてくれると答え様もあるのでは? -- akira? 2007-10-14 (日) 16:06:46
  • akira様,早速のコメントありがとうございました。"フォントサイズ""設定"などで検索していたら全く出てこなかったのですが,アドヴァイスの通りグラフィックパラメータのページを拝見しましたら,cex.label等で基準のサイズから拡大・縮小できることが分かりました。"ラベル""拡大"などのキーワードで検索したら良かったようですね。どうもありがとうございました。 -- 田中? 2007-10-14 (日) 17:59:55

線分の描き方について

真矢? (2007-10-11 (木) 18:27:55)

座標(X,Y)を始点とした傾きm、長さLの線分の書き方を教えてください。
R2.4.1を使っています。

  • ?segments -- 2007-10-11 (木) 18:50:55
  • この手の質問は検索すれば一発で出てきますよ。あと、リンク集の初心者コースも良いと思いますよ。質問する前に、頑張って調べた人に天使は微笑みます。調べたこと、でも分からなかったこと、を書いてくれると答え様もあるのでは? -- akria? 2007-10-14 (日) 16:11:34

ESSからRcmdrを日本語で呼び出したい

<ふ>? (2007-10-02 (火) 21:57:19)

(連続で失礼します。)

essをMeadowで使ってます。
(WindowsXP、R2.5.1、Rcmdr1.3-5、ess-5.3.4、Meadow3)
このessのR processからlibrary(Rcmdr)を実行すると、langが渡らないのか、起動したRcmdrのメニューが日本語になりません。

そこで、Rのインストールディレクトリの /etc/Rconsole を開いて、

## Language for messages
language = ja
のように、強制的にjaにしてみました。

これで、起動するRcmdrのメニューなどは、めでたく日本語になったのですが、「データセットの表示」「データセットの編集」では、日本語が表示されません。
やはり、Tcl/Tkにちゃんとlangを伝えないとだめなようです。

Meadowの中のなにかを設定すればいいのではないかと思おうのですが、よくわかりません。

どなたかアドバイスをお願いいたします。

  • essd-r.el を見てみて下さい. たぶん, default-process-coding-system が windows の時に, undecided-dos なってる筈です. essのチームに伝わってる(R教授曰く、必ず伝える)と思うのですが, Rtermまわりで欧州近辺(win9xとかで)で問題があるのかもしれませんね. -- なかま 2007-10-02 (火) 23:40:28
  • なかまさん、ありがとうございます。essの設定では、http://www.okada.jp.org/RWiki/?ESSConfig を参考にさせていただいていてます。euc-japan-unix の部分は、sjis-unixにしてますので、ご指摘の部分は、undecided-dos にはなってません。ESS本家のMLをのぞいてみます。うまく現象を表現できるかどうか不安ですが....。 -- <ふ>? 2007-10-03 (水) 10:56:22
  • Rconsoleのlanguageはデフォルトのまま変更してませんが、環境変数LANGをja_JP.SJISに設定していると、最初からRコマンダーのメニューは日本語になってますし、「データセットの表示」「データセットの編集」で、Excelで編集した(日本語の混じった)csvファイルの中身が表示されましたよ。動作環境はESS以外は<ふ>さんと同じです(ESSだけ5.3.6)。5.3.6からJGRのように関数の入力中に引数のリファレンスがミニバッファに表示されるようになったので便利です。 -- 2007-10-03 (水) 13:17:49
  • たぶん、Meadowが環境変数LANGを設定している(私のは"JPN"になってました)ので,JPNなぞ知らんと言われている為のようです. Rを起動したら、`Sys.unsetenv("LANG")'としてから,library(Rcmdr)としてみて下さい. WindowsのRではアクティブコードページの932を元に日本語処理をしますが,TCLはLANGを優先(しかし"JPN"は知らない)してしまうからです. -- なかま 2007-10-03 (水) 13:50:13
  • 引数にLANGを与えればいいので、こんな感じでどうでしょう -- なかま 2007-10-03 (水) 14:19:19
    (set-default-coding-systems 'cp932)
    (set-terminal-coding-system 'cp932)
    (set-keyboard-coding-system 'cp932)
    (set-buffer-file-coding-system 'cp932)
    (require 'ess-site)
    (define-key ess-mode-map "\177"   'delete-char)
    (setq ess-ask-for-ess-directory nil)
    (setq ess-pre-run-hook
          '((lambda () (setq S-directory default-directory)
    	  (setq default-process-coding-system '(cp932 .   cp932))
    	  )))
    (setq inferior-R-args "LANG=")
  • なかまさん、ありがとうございます。LANGをja_JP.SJISにしただけでうまくいってしまいました。Sys.ensetenvはやってません。上にいただいたlispは後ほど試してみます。 -- <ふ>? 2007-10-03 (水) 22:16:48

datasetにhelpをつける方法を知りたいのですが

<ふ>? (2007-10-02 (火) 14:58:10)

RcmdrでR2.5.1を使っています。

データの集計は、Excelでおこなって、それをRのDatasetとして読み込ませたあと、数値を因子に変換したり、あれこれ修正をして、.Rdaで保存し、それをもとに分析を行っています。

この.Rdaでsaveされたデータセットに、helpをつける方法がわかりません。

Packageで提供されているDatasetには、helpがついていて、便利なものですあから(というよりも不可欠..)、一度整形したデータを使いまわすためにデータの定義など、きちんとしたドキュメントをデータにつけておきたいと思っております。

ご教示ください。

  • パッケージを作るを参照のこと。今回の場合,プログラムを含まないデータだけのパッケージを作ればよい。でも,結構面倒でしょう?日本語で作れるかどうかもよく分からないし,(たぶん大丈夫だろうけど)データセットの説明書をちゃんと作ればよいだけでは?(Word ででも何ででも) -- 2007-10-02 (火) 17:21:32
  • help(comment) -- 2007-10-02 (火) 17:36:41
  • help(comment)でオブジェクトにコメントする方法は、やってみます。ありがとうごじざいます。それと、やはり、packageなんですね。http://cran.r-project.org/doc/manuals/R-exts.pdf をちゃんと読んでやってみます。データセットの説明書はつくってあるのですが、やはり、一体化したいので。Tryしてまた報告します。 -- <ふ>? 2007-10-02 (火) 21:26:21

折れ線プロットは?(plot)

kd? (2007-09-25 (火) 03:46:53)

折れ線プロット(すなわちポイント点は無印)を描こうと思っております.
(Mac版 R-2.5.1.dmg を使っております.)
しかし下記を試しても折れ線プロットを描いてくれませんでした.何か方法はあるでしょうか?

plot(x$V2~x$V1,pch=-1,lty=1) ---> 丸印でポイント点をプロットするだけ

plot(x$V2~x$V1,pch="",lty=1) ---> 無印,無線

  • あのー、検索すれば一発ですが..."plot" "折れ線"とかしてみました?調べてすぐわかることを質問する人は怒られますよ。
    答えを教えてもkdさんのためにならないので、検索の仕方を教えますと、googleで"CRAN+plot+折れ線"とか、グラフの引数の疑問なら"グラフィックパラメータ"が吉。
    他、リンク集の初心者部分も熟読しましょうね。
    わからないことがわからないと、何から手を付けて良いのかすらわからないものですよね。学校の勉強と同じですね。 -- 2007-09-25 (火) 07:38:37
  • まず何より先に,ヘルプを読むべし。
    ? plot しましょ。掲示板に質問して答えを待つなんて,まだるっこしい。
    ... Arguments to be passed to methods, such as graphical parameters (see par). Many methods will accept the following arguments:
    type
    what type of plot should be drawn. Possible types are
    "p" for points,
    "l" for lines,
    "b" for both, -- 2007-09-25 (火) 10:51:07
  • 久しぶりに?plotしたんですが、type="h"とか"s"って、昔からありましたっけ?いや、まあどうでもいいんですが。 -- 2007-09-25 (火) 11:32:43
  • どれくらい前のことを「昔」というかによりますが,動いているシステムの一番古いものがR1.8.1ですが,そのときからありましたよ。 -- 2007-09-25 (火) 11:46:24
  • あ、そうでしたか。ありがとうございます。 -- 2007-09-25 (火) 12:25:49
  • 古い R をインストールして調べてみました。R 1.1.0 では,p,l,b,o,h,s,S,n で(少なくとも R 0.61.3 から),R 1.2.0 から p,l,b,c,o,h,s,S,n と c が増え現在に至るということみたいですね。なお,R 1.0.0 は古いperlがなくてインストールできないようでした。 -- 2007-09-25 (火) 12:58:56

igraphについて

あひる? (2007-09-20 (木) 23:42:45)

現在、大学でigraphを用いた授業を受けております。家でもじっくり課題ができるようにと早速Rをダウンロードしたのですが、どうしたらigraphが使えるようになるのかがわかりません。Googleで検索したページでigraphのダウンロードはできたようなのですが、Rでlibrary(igraph)と打ってもエラーになってしまいます。どうしたらigraphが使えるようになるのでしょうか?どなたかご教授頂けたら幸いです。

  • どのバージョンとか、OSは何だとか、どんなエラーが出るか、どこの大学か、書かないと怒られますよ。 -- 2007-09-21 (金) 01:18:33
  • どこの大学かは書かなくてもいいので(笑)、インストールしたOS名とRのバージョンを教えてください。質問の「感じ」では、パッケージインストールのやり方をまちがえているような気がします。 -- okinawa 2007-09-21 (金) 09:18:51
  • > Googleで検索したページでigraphのダウンロードはできた  と書いているところを見ると,ちゃんと「インストールできていない」とみた。 -- 2007-09-21 (金) 10:34:18
  • 大学の授業なら、担当教官に聞くのが一番早いですよ。 -- 2007-09-21 (金) 10:58:02
  • 今夏休みで,かつ向学心に燃える学生なんでしょう -- 2007-09-21 (金) 11:12:29
  • すみません、あまりこうゆうのには慣れていないもので。。OSはWidouwsXPでRはR2.5.1です。 -- あひる? 2007-09-22 (土) 03:25:06
  • まず、Rを起動してください。次にメニューから[パッケージ][パッケージのインストール...]を選びます。するとダウンロードするサイトが出てきますので、Japan(Tukuba)を選んでください。その後パッケージのリストが出てきますので、igraphを選んで[OK]を押すと自動的にigraphが読み込まれます。特殊なパッケージ以外はすべてこの方法でインストールを行ってください。 -- okinawa 2007-09-22 (土) 07:54:02
  • 上の方法で無事解決しました。どうもありがとうございますm(_ _)m -- あひる? 2007-09-23 (日) 23:44:21

R Consoleの表示容量について

カイロ? (2007-09-20 (木) 22:25:29)

Rを利用して3日目です.医療データを用いたROC曲線を用いて予後推定のcut-off値を求める作業をしていますが,データの取得で躓いています.RjpWiki内,Google,等を検索しても分かりませんでしたので,こちらでお伺い申し上げます.大変初歩的な質問と思いますが,ご教示いただけましたら幸いです.ROC曲線作成には青木先生のライブラリからsourceを使用させて頂いております"http://aoki2.si.gunma-u.ac.jp/R/ROC.html".

図のように,R Console上で数千行ある計算結果が得られるのですが,スクロールバーを一番上まで持っていっても全ての結果を表示できません.GUIプリファレンスでbuffer bytes等を増やしてもうまくゆきません.R Console上で表示,あるいは別ファイルに書き出し,などして全ての計算結果を取得する方法はどうすれば宜しいでしょうか.

宜しくお願いいたします.

当方の環境:
R version 2.5.1 (2007-06-27)
i386-pc-mingw32
locale:
LC_COLLATE=Japanese_Japan.932;LC_CTYPE=Japanese_Japan.932;
LC_MONETARY=Japanese_Japan.932;LC_NUMERIC=C;LC_TIME=Japanese_Japan.932
attached base packages:
[1] "stats" "graphics" "grDevices" "utils" "datasets"
[6] "methods" "base"

consolegamen.png


  • sink 関数を調べましょう。sink("foo.bar"); 作業; sink() でコンソールに表示されるであろう内容が "foo.bar" というファイルができます。 -- 2007-09-20 (木) 23:01:30
  • できました! ".bar"を拡張子とするとうまく読めませんでしたが,適当なファイル名(R.logとしました)で指定すると目的の値を取得することができず,ずっとここで作業が止まっておりました. Prompt replyありがとうございました. -- カイロ? 2007-09-20 (木) 23:46:08
  • ところで,sink関数で吸い出された計算結果が必ずR本体のあるbinフォルダに書き出されるのですが,これを任意のフォルダに設定することができますでしょうか.指定した作業ディレクトリに書き出されると便利なのですが. 重ねての質問すみません.-- カイロ? 2007-09-20 (木) 23:51:46
  • setwd("作業ディレクトリ") -- aa? 2007-09-21 (金) 07:42:58
  • sink("c:/ファイル名") で直接ディレクトリ指定して落とせます。日本語windowsの¥は使えませんのでご注意を。なおファイル名の拡張子は.txtとかにするといいですよ。 -- okinawa 2007-09-21 (金) 09:12:31
  • ご教示頂いた方法でうまくゆきました.aaさん,okinawaさん,ありがとうございました. -- カイロ? 2007-09-21 (金) 14:47:28

nlme ライブラリ,lme のデフォルト

キュロシ? (2007-09-20 (木) 08:32:11)

Mixed Effects Models in S and S-PLUS を勉強中のキュロシです.
作者(Pinheiro & Bates)に訊きなさい,という種の質問かもしれません.

> library(nlme)

として nlme パッケージを使います.いくつかのデータセット (groupedData) がありますが,ここでは Orange を使って説明します.

> f1 <- lme(Orange)

とすると,以下の結果が得られます.

> f1
Linear mixed-effects model fit by REML
Data: Orange 
Log-restricted-likelihood: -139.9061
Fixed: circumference ~ age 
(Intercept)         age 
17.3996502   0.1067703 
Random effects:
Formula: ~age | Tree
Structure: General positive-definite
           StdDev      Corr  
(Intercept)  1.89732541 (Intr)
age          0.02402639 -1    
Residual    10.00622121       
Number of Observations: 35
Number of Groups: 5

つまり,lme の argument である,fixed, random を data を指定することによって,はしょっています. とすると,

> f2<-lme(fixed=circumference ~ age, data=Orange, random=~age|Tree)

とすることで同じ結果が得られそうなものですが,得られません. f1 のモデルが 収束しない f2 と何が違うのかがわかりません.何かヒントになりそうなことがありましたら,是非とも教えてください.よろしくお願いいたします.

  • lme(Orange) と同じ結果になるのは lme(fixed=circumference ~ age, data=Orange) ですよ。
    lme(fixed=circumference ~ age, data=Orange, random=~age|Tree) と同じではないですね。 -- 2007-09-20 (木) 09:26:02
  • その場合の random effects は何なんでしょうか? -- キュロシ? 2007-09-20 (木) 13:18:42
  • 私には分かりません。help に書いてあるかも。そもそも Mixed Effects Models in S and S-PLUS (って,本ですかね)を読めば分かるのではないかと思いますが,私はその勉強をするつもり(必要)はないので。 -- 2007-09-20 (木) 13:50:40
  • Pinheiro に訊いたところ, lme(Orange) とすると,まず lmList で random effects の分散共分散行列の推定を行うそうです.それで,このやり方だと収束して,他のやり方だと収束しないことがある,ということでした.ただ,Orange データには非線形モデルを使わなくっちゃいけない,と怒られちゃいました. -- キュロシ? 2007-09-26 (水) 00:18:00

回帰分析についての質問

netbean? (2007-09-18 (火) 16:46:17)

Rの初心者です。他ツールで解析経験がありますが、Rはまだ数時間の実践しかやっていません。
step関数を用いてステップワイズ変数選択を実施する際、下記エラーが出ました、

> step(RegModel.2)
Start:  AIC= 189.44
。。。
object$residuals + predict(object) 中で警告がありました:
         長いオブジェクトの長さが
         短いオブジェクトの長さの倍数になっていません
以下にエラーlm.fit(x[, jj, drop = FALSE], y, offset = offset) : 
        矛盾した次元です

回帰のモデル精度が非常に低く、stepwiseを期待したいところが、ここで固まって進めない状況です。どなたか対処法をご存知でしたら、教えてください。ありがとうございます。

  • エラーメッセージの通りです。object$residualsとpredict(object)の長さが違うためのエラーです。x[,jj,drop=FALSE]という次元はないというエラーです。おそらく文法間違いです。対処法は「正しい文法を使う」ということでしょうか? -- 2007-09-18 (火) 17:11:01
  • 早速の回答、ありがとうございます。step関数の文法関しては、step(RegModel?.5)で単純コマンドを入れました。他のデータセットで同じコマンドでは、動いたため、なぜですか(Step関数の中身?それともデータの問題?)。初歩的な質問、申し訳ありませんが、教えてください。 -- netbean? 2007-09-18 (火) 17:26:50
  • あなたが何をどうやったのかがわかりませんので,回答不能でしょう。実際に使ったデータが大きすぎるとか,データを公開できないという場合には,エラーを再現できる最小限の架空データでよいからデータを添えて,あなたのやった操作が分かるプログラムを添えて質問すべきでしょう。そうでないと「エラーメッセージの通りです」以外のコメントは付きません。 -- 2007-09-18 (火) 18:00:47
  • RegModel? が何なのかわからないのに答えようがある? -- 2007-09-18 (火) 18:01:45
  • それにしても「使いはじめてから数時間」で質問するという度胸に感(寒)心 -- 2007-09-18 (火) 18:12:23
  • RegModel?は、線形回帰の結果です。RegModel? <- lm(X1_G~B+BP+CH+CHP+CHPR, data=dataset) -- netbean? 2007-09-18 (火) 18:33:30
  • Data公開準備中です。Rの三日間セミナーで集中勉強したが、殆んど、コマンダーを使いました。また、自分のデータの解析は、まだ数時間の実践しかないとの意味でした。初歩的ばっかりで申し訳ありませんでした。 -- netbean? 2007-09-18 (火) 18:45:27
  • 他のデータセットで動いて手持ちデータで動かないということは、step関数ではなく「他のデータセットと手持ちのデータセットの違い」が問題じゃないですか? -- 2007-09-18 (火) 20:46:10
  • RegModel?? <- lm(X1_G~B+BP+CH+CHP+CHPR, data=dataset) と書かれても,何の情報もないということが分からないのでしょうか?あなたが質問を受ける側で,ここに書かれている情報だけで何が分かるかと,反省してみたらいかが?(というか,そう言うことができないのでしょうね) -- 2007-09-18 (火) 21:07:03
  • > 「使いはじめてから数時間」で質問する というのは,なにも卑下する必要はありません。質問するのに経験は不要。むしろ,経験がないからこそ質問はあるもの。しかし,質問の仕方というのは注意すべし。上の方に色々注意事項が書いてあるでしょう?そう言うことに無頓着で,指摘されても分かっているのか分かっていないのか無視するという態度はほめられたものじゃない。 -- 2007-09-18 (火) 21:10:35
  • 質問者本人じゃないですが、そんなに寄ってたかってつつかないでもいいと思うけど・・・。自分が初心者の頃を思い出してみましょうよ。 -- 2007-09-18 (火) 21:11:14
  • つついているわけではなくて。ちゃんとした回答を得ようと思うときにはこういう情報が必要ですよといっているだけでしょう?それが不満なら,回答を放棄すればよいだけ。あなたが初心者だったとき,誰かにアドバイスを求めたかった(求めたとき),この質問者のようにしましたか?この質問者のように,足りないところを指摘されて,なおかつそれに答えないという対処をしましたか?あなたはそんな質問をすることも必要なくRの達人になったのかも知れないけど。この質問者に対して潜在的な回答者はだれも「ついついたりしてない」し「軽んじてもいない」し,あなたが思っているような態度で接しているわけではない。一刻も早く問題解決をしたいならどうしたらよいか,少々辛口かも知れないが,その方策を呈示しているに過ぎないと思うが(注意事項は書かれているわけだし)??あなたのような見方こそ,質問者を特殊視しているのでは?あなたが質問者に同情するなら,今の情報だけで質問者に回答(の候補)を差し上げればいかが?
    ちょうど,この下の質問「判別分析の変数選択」を見てみたらいかが?指摘に応じて情報を提供したら(まあ,回答も辛口だけど)直ぐに答えは出たじゃないですか。今回の質問者は質疑応答の2巡目を終わっているんですよ。次が3巡目。 -- 2007-09-18 (火) 21:28:12
  • この掲示板って役にたつけどかなり雰囲気悪いよね.こわーい.こわーい.恋をしたら優しくなれるよ. -- 2007-09-19 (水) 00:00:25
  • 恐いですよ(笑)。常連回答者はいささか無神経な質問に内心うんざりしていますからね。しっかり勉強した上で、それでもわからなければ必要情報をしっかり与えた上で質問しましょう。それぐらいなら自分で調べてやる、となればなおさら結構。 -- 2007-09-19 (水) 00:42:37
  • 他の掲示板で、ここの場所の批判(袋だたきに会う、雰囲気悪い)を目にしたことがありますが、回答者は真剣に、そして正確に答えを見つけてあげたいと思っています。(その解析が人の生き死にを左右することだってあるからです)そのため「教育的な意味を含めた辛口な回答」もあります。私自身何度も「ムッ」ときたことがありましたが、言葉の裏に「教育的指導」を読み取り今日に至っています。 -- okinawa 2007-09-19 (水) 08:39:42
  • みなさん、いろんな意見とコメント、ありがとうございます!初めての質問、必要な情報さえ提供できず、質問方が良くなかったことを反省して、これからはまず良い質問者になりたいと考えております。この問題について、自分でもう一度突っ込んでみます。(データの問題と思うがデータの公開方法がうまく見つけず、また、面倒くさいがとりあえず回避案をやってみた)。Rはひとりで勉強しており、壁か多く質問さえもうまくできないくらいに思っちゃったりしますが、みなさんは、一番最初の時、ように勉強しましたか。たぶん、どこかに答えがあると思いますが、なかなか辿り着けなくて、、、 -- netbean? 2007-09-19 (水) 11:06:24
  • うまくいくほうといかないほうのデータを見比べて、少しずつうまくいかないほうを単純化してみるとか・・・、再現できなくなったらその変更点が問題。単純化させても再現できたて、その(単純な)データを見て調べても原因がわからなければ、データを添付して質問する。周りに使える人がいないと時間もかかって大変だけど、しょうがない。まずは本のチュートリアルを一通りやってデータの形式や構造になれるのがよいかと。 -- 2007-09-20 (木) 03:00:05
  • 最初は他の人のスクリプトを書き換えながらこつを覚えるのが良いと思います。私の経験では思いついた方法でも、どこかに似たスクリプトがあったりしました。 -- 2007-09-20 (木) 07:43:44

判別分析の変数選択

mayu? (2007-09-15 (土) 22:49:29)

discの線形判別関数のソースを使い、判別分析をしています。
自分が選んだ変数だけで解析したい場合どうすればいいのでしょうか?
変数選択が出来るソースを探していたのですが見当たらず質問させて頂きました。
自分でプログラムを組めばいいのかもしれませんが、プログラムが分からず、
勉強不足ですみませんが、よろしければ回答お願いします。

  • disc とは何ぞや。「ソースを使い、判別分析」とは何ぞや.赤の他人にもわかるように質問しましょう. -- 2007-09-15 (土) 23:21:01
  • > 自分が選んだ変数だけで解析したい場合どうすれば
    そもそも,自分で選んだ説明変数を引数に指定するのでは?
    自分の望まない変数まで使われてしまうというのは被害妄想では?
    二群判別なら,step(lm(ホゲホゲ))を使う。 2007-09-16 (日) 00:59:53
  • 回答ありがとうございます。 説明不足ですみません。  discは、青木先生のホームページから頂いてきた線形判別関数で、それを使い判別分析の勉強をしています。
       自分で選んだ説明変数を指定して解析しようとするとエラーが出てしまいました。
    > result <- disc(iris[,2,4], iris[5])
    以下にエラーdata[OK, ] :  次元数が正しくありません
    指定の方法が間違っているのでしょうか。  step(lm())ですか、勉強してみます。ありがとうございます。 -- mayu? 2007-09-16 (日) 23:29:35
  • disc( iris[c(2,4)], iris[5] ) . それと「変数選択」という言葉は統計学では「なんらかの基準で最適な変数の組合せを選択する」という意味で使われることが多く、単に「自分で選択した変数の組合せ」という意味で使うと誤解を招きます.disc の解説にある「変数選択は行えません」とは前者の意味です。 -- 2007-09-17 (月) 07:45:33
  • この質問に関してデータフレームの添字操作に関する不可解な挙動に気づきました。 -- 2007-09-17 (月) 08:06:23
    > x <- data.frame(a=1:3, b=rnorm(3), c=runif(3), d=3:1)
    > x[,2]
    [1] -0.3922563  1.3033630  0.1997175
    > x[,2,3]                             # これがエラーにならない                
    [1] -0.3922563  1.3033630  0.1997175
    > x[,2,3,4]
     以下にエラー`[.data.frame`(x, , 2, 3, 4) :  使われていない引数 (4)
    ああなるほど、x[,2,drop=3] つまり x[,2,drop=TRUE] と解釈されたんですね.
  • result <- disc(iris[,c(2,4)], iris[,5]) または,result2 <- disc(iris[c(2,4)], iris[5]) ではないですか?添え字指定を復習すべし。
    たいていの場合は,自分の使い方のミスを反省すべし。 -- 2007-09-17 (月) 21:11:48
  • お礼が遅くなりすみません。皆様、ご回答ありがとうございました。勉強不足だと痛感しましたが、前に進むことが出来ました。ありがとうございます。 -- mayu? 2007-09-18 (火) 21:46:34

数値データ、非数値データが混在するデータのカウント方法について

terada? (2007-09-14 (金) 21:20:20)

こんにちは。数値データ、非数値データが混在するデータのカウント方法について質問させてください。

> dat
      記号1    記号2     記号3       記号4
90         B1         A1          A2           A2
91         B1         A1          A1           A1
92         B1         A1          A1             
93         B1         A1          A1           A1
94         B1         B1          B3            K
95                                              K                               
96         B0         B0          B1           A1
97         B0                              
98         B1         A1                    
99          3                              
100                    3                              
101         2                              
102         2                              
103         1      

以上のようなデータ(14行4列のデータで、空いているところはnullです。)がある場合に、このデータフレームのすべての要素、B0、B1、B3、A1、......... 、1、2、3がそれぞれ何個ずつあるかカウントし、かつこのデータフレームにはない要素、B2、A3は0であるとカウントさせ、その数値データをベクトルの形で得たいのですが、何か方法があれば教えてください。table関数を使って色々やってみたのですがデータフレーム単位の結果がどうしても出せません。
よろしくお願いします。

すみません。質問の仕方が不適切でした。訂正します。 以下のような結果が得たいという意図で質問しました。

A1  A2  A3  B0  B1  B2  B3   K  1  2  3
11   2   0   3   7   0   1   2  1  2  2

再度アドバイスをいただければ幸いです。よろしくお願いします。

  • 例えば(簡単な例で、質問者も意味がわかる最小限の例を提示すると解答しやすい) -- 2007-09-14 (金) 21:48:17
    > ( d <- data.frame(x=c("a","",3,1), y=c(1,2,"","a")) )
      x y
    1 a 1
    2   2
    3 3  
    4 1 a
    > ( ud <- unlist(d) )             # リスト、データフレームをベクトル化する定番の方法
    x1 x2 x3 x4 y1 y2 y3 y4 
     a     3  1  1  2     a 
    Levels:  1 3 a 2
    > levels(ud) <- c(levels(ud),"d") # データフレームに登場しない要素を加える
    > ud
    x1 x2 x3 x4 y1 y2 y3 y4 
     a     3  1  1  2     a 
    Levels:  1 3 a 2 d
    > table(ud)
    ud
      1 3 a 2 d 
    2 2 1 2 1 0   # 先頭の2は空文字列が2回登場した、という意味
  • ついでながら。数値と文字が混在するデータフレームでは、数値は実際は文字列化(より詳しくいうと因子化)されていますから、正確には数値ではありません。更に注意すると空文字列 "" は NULL ではありません。NULL は文字通り何もない,空文字列は「空な文字列」がある。 -- 2007-09-14 (金) 21:51:58
  • うーん,このやりかたではうまくいかない? それとも空項目は無視したいというのなら。-- 2007-09-14 (金) 23:11:37
    > table(ud)[-1]  # 但し空文字水準がいつも先頭にくるかどうかは知りませんから,要確認
    ud
    1 3 a 2 d 
    2 1 2 1 0
  • 最初の質問を書き換えませんでした?一度コメントが付いた記事を書き換えるのはルール違反です.書き換えずに、新たに追加してください. -- 2007-09-15 (土) 00:19:13
  • 重ね重ねすみません。大幅な変更ではないので問題ないと思い、書き換えてしまいました。最初の質問に戻したほうがよろしいでしょうか? -- terada? 2007-09-15 (土) 00:30:36
  • ルール違反というのは,それに対するコメントが頓珍漢になる可能性がある(実際そうなっていますね)から。今後は気をつけましょう。で、肝心なことにまだ答えていませんが、問題は解決したのですか. -- 2007-09-15 (土) 06:27:38
  • 今日やっと解決しました。アドバイスありがとうございました。 -- terada? 2007-09-17 (月) 17:22:52
  • ついでにどうしたら解決したのか報告していただくと、他の方にも参考になりますね.お願いします. -- 2007-09-17 (月) 19:20:44
  • 2度目にアドバイスしていただいた方法で解決することができました。コメントが遅くなってすみません。 -- terada? 2007-09-20 (木) 15:44:04

ベクトルから特定の要素を削除するには?

青葉? (2007-09-13 (木) 13:40:42)

昨日から始めた初心者です。壁にごつんごつん当たりながら進んでます。

xから、yの中にある要素を数も考慮して取り除きたいのですが、どうすればよいでしょうか?
x <- c(1,1,1,2,2,2,3,3,3)
y <- c(1,2,2)
この2つから、(1,1,2,3,3,3)というベクトルを作りたいのです。
また、xから要素の種類だけ特定して取り除くことが可能でしょうか?
例えばx <- c(1,1,1,2,2,2,3,3,3)から2という数字だけを指定して、すべての2を取り除く方法です。
検索の仕方が悪いのかもしれませんが、苦戦しております。どなたかどうぞよろしくお願いします。

  • この程度のことは自分でなんとかして初めてRへの理解が進むと思いますので、がんばってください。 -- okinawa 2007-09-13 (木) 13:44:37
  • コメントありがとうございます。ヒントだけいただいてよろしいでしょうか?forループを使うしなないのでしょうか?そうするとなんとか出来そうです。 -- 青葉? 2007-09-13 (木) 14:05:22
  • (okinawa さんの意見に内心同感ながら) x[-which(x == 2)], x[x != 2]。 -- 2007-09-13 (木) 14:12:10
  • x[!x%in%2]とか。
    教えて君にならないように頑張ってください。上手な検索方法を身に着けることも大切です。 -- 2007-09-13 (木) 14:44:43
  • 最初の例は for ループを素直に使うべき例? -- 2007-09-13 (木) 14:50:17
    > xx <- x
    ## match は最初のマッチ位置を返すので好都合
    ## xx を順次縮小変更していることに注意
    > for (i in y) xx <- xx[-match(i,xx)] 
    > xx
    [1] 1 1 2 3 3 3
  • みなさまコメントありがとうございます。小さな壁は自分で乗り越えてるつもりなんですが、もうちょっと努力するよう頑張ってみます。forループを使ったコメントありがとうございます。私がやったら、30行くらいかかってしまいました。 -- 青葉? 2007-09-13 (木) 15:10:06
  • 奇策を弄す -- 2007-09-13 (木) 15:26:04
    > x <- c(1,1,1,2,2,2,3,3,3)
    > y <- c(1,2,2)
    > result <- merge(data.frame(a=table(x)),
    +                 data.frame(a=table(y)),
    +                 by.x="a.x", by.y="a.y", all=TRUE)
    > rep(as.numeric(levels(result[,1]))[result[,1]],
    +     ifelse(is.na(result[,3]), result[,2], result[,2]-result[,3]))
    [1] 1 1 2 3 3 3
  • rep(unique(x),c(table(x) - table(factor(y,levels=unique(x))))) -- 常夏? 2007-09-17 (月) 09:02:43

ベクトルの要素の種類数を知るには?

青葉? (2007-09-13 (木) 00:13:47)

今日使い始めたばかりの初心者です。
ベクトルの要素の数を知るにはlengthで良いことは分かったのですが、種類数を知る方法が見つからず1時間ほど苦戦しています。
例えば、x <- c(1, 1, 2, 2, 2, 3)だったら、要素の数は6ですが、要素の種類数は、1と2と3の3つです。この3という値を知りたいのです。
どなたかどうぞよろしくお願いします。

  • length(unique(x)) 他にも一発の関数はありそうだけど。 -- takahashi? 2007-09-13 (木) 00:33:58
  • ありがとうございます!これでやっと先に進めます。感謝です! -- 青葉? 2007-09-13 (木) 00:40:39
  • R 本体には一発関数はないと思います。takahashi さんの例が一番正統派ですが、要するに一意化すれば良いわけですから、副作用を利用した次のようなやりかたも。 -- 2007-09-13 (木) 10:30:40
    > factor(x)
    [1] 1 1 2 2 2 3
    Levels: 1 2 3
    > nlevels(factor(x))
    [1] 3
    > union(x,NULL)
    [1] 1 2 3
    > length(union(x,NULL))
    [1] 3
    > duplicated(x)
    [1] FALSE  TRUE FALSE  TRUE  TRUE FALSE
    > sum(!duplicated(x))
    [1] 3
    > table(x)
    x
    1 2 3 
    2 3 1 
    > names(table(x))
    [1] "1" "2" "3"
    > length(names(table(x)))
    [1] 3
    > sum(table(x) == 2)  # 丁度二回登場した要素数
    [1] 1
  • 副作用ならsummary(factor(x))とか -- 2007-09-13 (木) 12:08:09
  • 答えは3にならないといけないので,summary(factor(x)) ではなく length(summary(factor(x))) でしょうね。他に,length(table(x)); length(levels(factor(x))) -- 2007-09-13 (木) 12:25:39
  • 他の方々もコメントありがとうございます。予想外に多くの関数を知ることができて得した気分です。 -- 青葉? 2007-09-13 (木) 13:11:39

tcltk のロードができません

yosito? (2007-09-11 (火) 18:21:39)

R 初心者です。環境は MacOS X (10.4) にて、R.app を使っています。
tcltk は標準で package に含まれているのですが、いざロードしようとしますと、
以下のようなメッセージが出てロードできませんでした。

要求されたパッケージ tcltk をロード中です
Tcl/Tk インターフェース ロード中... Error in fun(...) : couldn't connect to display ":0"
Error :  .onLoad は 'tcltk' のための 'loadNamespace' で失敗しました
エラー: パッケージ 'tcltk' をロードできませんでした

どなたか対処法をご存知でしたら、教えてください

  • Tcltkが必要なパッケージを起動させる場合は、X11を起動した後にRを起動し、R.appの画面からSys.putenv("DISPLAY"=":0")を実行してください。 -- okinawa 2007-09-11 (火) 18:25:28
  • Rを起動した後でも,X11が必要になったら,ツールバーの「X11を起動」をクリックすればおーけー。ツールバーにないときは「ツールバーをカスタマイズ」 -- 2007-09-11 (火) 18:30:07
  • どうもありがとうございました。X11 を起動させておく必要があったのですね。無事 tcltk がロードできました。とても迅速な対応に感謝、感謝です。 -- yosito? 2007-09-11 (火) 21:53:52

nlme ライブラリ内の plot のマージンを変更する方法

キュロシ? (2007-09-11 (火) 04:42:46)

Pinheiro & Bates の "Mixed-Effect Models in S and S-PLUS" という本で勉強しています.

library(nlme)

とした後

plot(Dialyzer,outer=~QB)

を実行すると,下図のようなきれいな図が描けます.

groupedData1.jpg


ひとつだけ大きな外れ値があったとします.すなわち

Dialyzer$rate[1] <- 500

です.

ここで,先ほどと同じ図を描こうとすると

plot(Dialyzer,outer=~QB)

マージンが妙に大きい奇妙なプロットになってしまいます.

groupedData2.jpg


このマージンを自由に変更したいのですが,いまひとつやり方がわかりませんでした.
ヘルプファイルは

?plot.nffGroupedData(nlme)

だと思うのですが,該当するヘルプがわかりませんでした.
このマージンはどのように変更できるのでしょうか?

  • マージンって何ですか?ylimのことですか?plot(Dialyzer,outer=~QB, ylim=c(0,100))じゃだめですか? -- 2007-09-11 (火) 07:14:31
  • そんなに極端な外れ値がある場合,全部を含めた図を描こうとするよりも先に,外れ値を含めて描くことの意味,外れ値の生じた原因を考える方が先だと思います。それでもなおかつ描きたいならば,軸を対数目盛りにする機能があればそれを使う,なければ,y軸に取る変数を前もって対数を取っておくとか。まあ,対数グラフというものに対する賛否もいろいろありますが。 -- 2007-09-11 (火) 09:31:59
  • ヘルプがあるとすれば見方は ?plot.nffGroupedData?  でしょう。 -- 2007-09-11 (火) 09:49:19
  • ylim の値によって,プロットの見た目が上記の例のように変わってしまうを防ぐ手立てが知りたいです.最初の例では「でっぷり」しているのに,二番目の例では,極端にスリムになってしまっています.outer margin でしたっけ? -- きゅろし? 2007-09-11 (火) 10:31:23
  • コメントに従って help を見ていますか?
    plot(Dialyzer,outer=~QB, log="y", asp=1)
    のように,log パラメータと asp パラメータを使うのが吉かと。
    それにしても,極端な外れ値ですねぇ。
    それと,グラフファイルはjpegは不適切ですよ。pngがよいでしょう(Windowsだとどっちもどっちかな?) -- 2007-09-11 (火) 11:23:32
    nlme-plot.png
  • ありがとうございました.別に外れ値でなくてもよかったのですが,ylim の値によって, asp が変わる,というのはどういう仕様なのでしょうか?理由がいまひとつわかりませんが・・・ -- キュロシ? 2007-09-11 (火) 14:10:57
  • > ylim の値によって, asp が変わる
    普通はその方が好ましい(少なくとも作者はそう考えた)と言うことでしょう。
    しかし,場合によっては困るのでasp引数を使うと言うことでしょう。
    理由が知りたい場合は,パッケージの作者にお聞き下さい(^_^;) -- 2007-09-11 (火) 14:26:25
  • パッケージ作者はこんなおかしなデータのことまで考慮していないのでしょう。どうせ例示用にせよ、もう少しましな例をお願いします。 -- 2007-09-11 (火) 15:59:35

特定の行と列を削除した行列の作り方

KK? (2007-09-10 (月) 11:57:13)

ある行列Aの特定の行および列を削除した行列Bを作りたい場合、どのようにすればよいでしょうか。
よろしくお願いいたします。

  • 添え字の指定の仕方を調べましょう行列Tips大全などをまずは参照すべきでは? -- 2007-09-10 (月) 12:43:15
  • 検索不足でした。ありがとうございます。 -- KK? 2007-09-10 (月) 12:53:15

日付+時刻の時系列のプロットでの扱い

たけ? (2007-09-01 (土) 15:54:42)

たとえば2006/01/01 00:00:00から2007/09/01 15:48:32までの時系列で欠損値-999、行の欠落はないデータを扱う際、解析自体は適当に行番号でできるのですが、plotする際に、x軸に日付をどころどころ表示したプロットを出力したいです。ts()でも年+月までであればそれっぽいものができますが、この日付+時刻の状態を扱ってプロットにもっていけるパッケージやtsオブジェクトの作成方法をご存知の方がいましたらお教いただけますでしょうか。環境は2.51 w32です。
どうぞよろしくお願いします。

  • 例えば -- 2007-09-01 (土) 19:42:36
    gnp <- ts(cumsum(1 + round(rnorm(365), 2)))          # 例示用時系列
    Dt <- seq(as.Date("2001-1-1"), len=365, by="1 day")  # 2001-1-1 から一年間の日付
    DDt <- sub("2001-","",Dt)                            # "2001-" を取る
    plot(gnp, xaxt="n")                                  # x-軸注釈無しでまずプロット
    Dy <- 1+50*(0:7)                                     # 注釈を付ける日付. 1-1 から50日おき
    Dc <- DDt[Dy]                                        # 注釈文字列ベクトル
    axis(1,at=Dy, labels=Dc)                             # 最後に注釈を付ける
    あとは Dy (必要に応じて DDt も) を好きなように変えればよい。
  • ありがとうございます。テキスト処理に不慣れでしたが、何日かおきに手作業で注釈をつけるのもこんなに短くできるのですね。参考にしていじります。 -- たけ? 2007-09-01 (土) 22:48:02

2つの因子の重複

KEN? (2007-08-30 (木) 17:53:43)

2つのデータ(data1とdata2)の各因子(1〜6)の一致をみたいと思っています。一致したデータ数(今回の場合は3行目の1つ)とその割合(今回の場合は1/5)を計算するにはどうすればよろしいでしょうか。Excelに限界を感じ先週よりRの勉強を始めました。超初心者ですが、よろしくお願いいたします。

	data1	data2
1	A	B
2	B	A
3	B	B
4	B	A
5	A	B

使用環境:R version 2.4.1 (2006-12-18) i386-pc-mingw32

  • たとえば -- 2007-08-30 (木) 18:01:59
    > x
      data1 data2
    1     B     B
    2     A     B
    3     B     A
    4     A     A
    5     A     A
    > sum(x$data1 == x$data2)
    [1] 3
    > mean(x$data1 == x$data2)
    [1] 0.6
    注:アーカイブに投稿されていましたので、移動しました.
  • どうもありがとうございました。早速試してみます。 -- KEN? 2007-08-30 (木) 18:20:45
  • 試してみましたところエラーがでました「以下にエラーOps.factor(x$data1, x$data2) : 因子の水準セットが異なっています」。データに間違いがないか何度かためしてみたのですが同じ結果でした。データの読み込みミスなどで何か形式が誤っているのでしょうか。 -- KEN? 2007-08-30 (木) 19:22:53
  • 老婆心ながら。x$data1 == x$data2 は各行毎に一致するかどうかを示す論理値ベクトルを返します.論理値 TRUE, FALSE は数値が要求される局面ではそれぞれ整数値 1,0 とみなされますから,上のような操作ができるわけです. -- 2007-08-30 (木) 19:29:07
  • コマンド str(x$data1) と str(x$data2) の結果を投稿してください。恐らく data1 と data2 の因子水準が異なるのでは. -- 2007-08-30 (木) 19:31:08
  • > エラーがでました データフレームはちゃんと作られているのかな? -- 2007-08-30 (木) 19:34:46
  • factor の比較は注意が必要。 -- 2007-08-30 (木) 19:38:38 比較するfactor ベクトルの水準数が等しくないとエラーになる
    > x <- data.frame(data1 = factor(sample(letters[1:6], 10, replace=TRUE)), 
    + data2  = factor(sample(letters[1:6], 10, replace=TRUE)))
    > x
       data1 data2
    1      d     f
    2      c     d
    3      f     f
    4      a     b
    5      c     a
    6      d     d
    7      d     e
    8      f     c
    9      c     a
    10     e     e
    > x$data1 == x$data2
     以下にエラーOps.factor(x$data1, x$data2) :  因子の水準セットが異なっています
    > levels(x$data1)
    [1] "a" "c" "d" "e" "f"    ← "b" がない
    > levels(x$data2)
    [1] "a" "b" "c" "d" "e" "f"
    文字型のまま比較した方がよい
  • 本来の値に戻して sum( levels(x$data1)[x$data1] == levels(x$data2)[x$data2] ) ではどうですか。 -- 2007-08-30 (木) 20:01:53
  • ↑ character にして比較しているので○。もし,データファイルから読んでいるなら,read.table の as.is 引数を調べれば,文字型で読み込めるので,それもためせば◎? -- 2007-08-30 (木) 21:24:25
  • str(x$data1)は「Factor w/ 4 levels」、str(x$data2)は「Factor w/ 3 levels 」でした。水準数が異なっていたようでした。一番初めに実データを添付しておけばよかったのですがデータが大きすぎたために一部のみのとさせていただきました。お手数おかけしました。〜 -- KEN? 2007-08-31 (金) 10:13:05
  • sum( levels(x$data1)[x$data1] == levels(x$data2)[x$data2] ) でうまくいきました。read.table の as.is 引数についても試してみます。 -- KEN? 2007-08-31 (金) 11:47:12
  • 念のため。 levels(x$data1)[x$data1] は仮に x$data1 が数値でも文字列にしてしまいます。もし本来数値であるなら as.numeric(levels(x$data1))[x$data1] を使います。単に一致性のチェックだけなら文字列のままでも良いわけですが。 -- 2007-08-31 (金) 12:37:13
  • 因子Tips大全 参照。 -- 2007-09-02 (日) 07:52:26

heatmapのlegend

みゅー? (2007-08-29 (水) 11:02:21)

かなり調べたつもりなのですが、heatmapのlegendのつけ方がわかりません。heatmapのカラー表示がどのレベルに対応するのかを明示したいので、方法が書いてあるページなどをお教えいただけないでしょうか。

  • 例えば次の例など参考になりますか.これは Zoonekynd 氏の Statistics with R 中にある例です。元マニュアルをみてください.他にも参考になる例があるはず。 -- 2007-08-29 (水) 11:58:35
    # Polar plot to spot seasonal patterns
    x <- as.vector(UKgas)
    n <- length(x)
    theta <- seq(0, by=2*pi/4, length=n)
    plot(x * cos(theta), x * sin(theta),
         type = "l",
         xlab = "", ylab = "",
         main = "UK gas consumption")
    abline(h=0, v=0, col="grey")
    abline(0, 1, col="grey")
    abline(0, -1, col="grey")
    circle <- function (x, y, r, N=100, ...) {
      theta <- seq(0, 2*pi, length=N+1)
      lines(x + r * cos(theta), y + r * sin(theta), ...)
    }
    circle(0,0, 250, col="grey")
    circle(0,0, 500, col="grey")
    circle(0,0, 750, col="grey")
    circle(0,0, 1000, col="grey")
    circle(0,0, 1250, col="grey")
    segments( x[-n] * cos(theta[-n]),
              x[-n] * sin(theta[-n]),
              x[-1] * cos(theta[-1]),
              x[-1] * sin(theta[-1]),
              col = terrain.colors(length(x)),
              lwd = 3)
    text(par("usr")[2], 0, "Winter", adj=c(1,0))
    text(0, par("usr")[4], "Spring", adj=c(0,1))
    text(par("usr")[1], 0, "Summer", adj=c(0,0))
    text(0, par("usr")[3], "Autumn", adj=c(0,0))
    legend("topright", legend = c(1960, 1973, 1986),
           fill = terrain.colors(3))  # このあたり
  • 助かります。ありがとうございます。早速試してみます。 -- みゅー? 2007-08-29 (水) 13:11:48
  • R Graphics Manual はRのパッケージ中の Example 画像を一挙に集めて見せるサイトです.多すぎて探すのは大変だけれど、探せば報われる。これなんかいかにも -- 2007-08-29 (水) 23:00:36

Rのコンパイラー

あつ? (2007-08-28 (火) 18:31:23)

http://www.cs.uiowa.edu/~luke/R/compiler/
RのコンパイラーをWindowsのマシンに入れようとするのですがうまく行きません。
tar.gz形式なので、一度展開してzip形式にしてGUI画面から読ませると、昔のRしか動かないよ。。と言った感じのメッセージが出ます。
仕方なく

R CMD INSTALL  compiler_0.1.-3.tar.gz

としても、テンポラリーの位置が合わないのか展開されせん。tar,gz,perlは、入っています。
どうやったら、これがR-2.5.1で動くか、どなたかお教え願えませんでしょうか?

  • エラーが良く分からないのですが、古いRが残っていて、新しいRのpathが通ってないとか -- 2007-08-28 (火) 19:59:32
  • R byte compiler は久しく開発が停止しています.そもそも Windows 版 R で動くかどうか怪しいですね.どうしても動かしたければ、Vmware + Ubuntu 等の仮想 Linux 環境を試してみればどうでしょう.問題なくインストールできました. -- 2007-08-28 (火) 21:33:43
  • 古いRは、残っていません。> > library(compiler)以下にエラーlibrary(compiler): ’compiler’は有効なパッケージではありません。バージョン<2.0.0でインストールされたかも知れません  と出ます -- あつ? 2007-08-28 (火) 21:46:31
  • Linuxはないので、sparcのsolaris10で試してみたら、そのままOKでした。 やはりwindowsだと駄目なのでしょうか? -- あつ? 2007-08-28 (火) 21:54:00
  • 解凍したパッケージのR/cmp.Rを、そのままwindowsの画面に貼り付けたら動作しました。要するにパッケージ作り方そのものがwindows用では無いようです。パッケージの作り方を知らないので、とりあえず作業スペースを保存して使っています。パッケージの作り方のページを見ましたが、よくわかりませんでした。 -- あつ? 2007-08-29 (水) 12:57:19

変数の交換

とも? (2007-08-27 (月) 17:29:39)

中間変数を使わずに内容を交換することは可能ですか。例としてa1とa2の内容を交換する方法をお示しします。

a1 <- c(1,2,3)
a2 <- c(2,3,4)
z <- a1
a1 <- a2
a2 <- z
  • できないでしょう。 -- 2007-08-27 (月) 18:38:24
  • 結局明示的にせよ暗黙にせよ中間変数を使わない限りスワップはできないでしょうね.次のような関数を定義すれば明示的に中間変数を使わないですみますが。-- 2007-08-27 (月) 18:44:59
    > swap <- function(a,b) {tmp <- b; b <<- a; a <<- tmp}
    > a=1:5; b=rnorm(5)
    > swap(a,b)
    > a 
    [1]  0.25953799 -0.13130845  0.51087979 -0.67557507  0.63444793 
    > b 
    [1]  1  2  3  4  5 
    なお、例えば行列の列や行の交換は次のようにすれば(明示的に中間変数を使わずに)できます。
    > X[c(m,n),] <- X[c(n,m),] # m 行と n 行の交換
  • "中間変数"を使わない例(笑) -- 2007-08-27 (月) 19:10:48
    > a <- 1:3; b <- rnorm(3)
    > a <- list(a,b); b <- a[[1]]; a <- a[[2]]
    > a
    [1]  1.23461645 -0.08733515  0.62929149
    > b
    [1] 1 2 3
  • listはデータ構造が違うので,ベクトルの中間変数ではないが,「"中間変数"を使わない例」というのはねえ(^_^;) ~そのようなものが許されるなら,a,bをファイル出力して,b,aの順で読むとか,いろいろトリックありそう。 -- 2007-08-27 (月) 20:52:45
    a <- 1:3
    b <- rnorm(3)
    write(a,"temp")
    write(b,"temp2")
    b <- scan("temp")
    a <- scan("temp2")
  • list の例の類例ですが -- 2007-08-27 (月) 21:07:54
    a <- 1:3; b <- rnorm(3)
    a <- rbind(a, b)
    b <- a[1,]
    a <- a[2,]
    まあ,いずれも,トリビアルな例
  • RubyではできるのでRではどうかなと思っていましたが、できないということが分かって有益でした。ありがとうございました。 -- とも? 2007-08-27 (月) 22:33:13
  • スワップを行うシステム関数があるかどうかというご質問だったのでしょうが,確かに R にあってもよさそうで無い関数ですね.あと Inc <- function(n) n <<- n+1 といった類の動作をするシステム関数もあってもよさそうなのになぜか無いですね. -- 2007-08-27 (月) 23:50:24
  • 「Rubyではできるので」 裏で中間変数を使ってるんじゃないの? -- 2007-08-28 (火) 10:16:41
  • Rにはスワップのための関数や演算子の使い方があるのかなと思って質問したものです。Rには無いということを調べるのは難しかったもので。 -- とも? 2007-08-28 (火) 12:57:11
  • a1 <- a1+a2 a2<- a1-a2 a1 <- a1-a2 でできませんかね. -- 2007-08-28 (火) 13:03:02
  • この問題が「クイズとして」出されるときに裏技として解説される技ですね。数値オブジェクトのときには有効です。+,- 演算子が使えない非数値オブジェクトのときにはできませんね。また,以下のような場合にも注意が必要です。 -- 2007-08-28 (火) 13:23:47
    > a1 <- 1e10; a2 <- 1e-10
    > a1
    [1] 1e+10
    > a2
    [1] 1e-10
    > a1 <- a1+a2; a2<- a1-a2; a1 <- a1-a2
    > a1
    [1] 0     ← これは困る
    > a2
    [1] 1e+10

関数内の変数が展開されない

kutakuta? (2007-08-24 (金) 20:46:58)

よろしくお願いします。 R version 2.5.1 (2007-06-27) powerpc-apple-darwin8.9.1 を使用しています。 ggplotというライブラリをつかって箱ヒゲ図+jitterplotを作ろうとしています。 行列"data"の中に"class"というカラムと"colx"(x=1-N)というカラムがあって、classごとの箱ヒゲ図を作るのが最終目的です。
例: http://bg9.imslab.co.jp/Rhelp/R-2.4.0/src/library/ggplot/man/images/ggjitter_004.png

p <- ggplot(data,aesthetic=list(y = class, x = col1))
ggjitter(ggplot(p))

とするとうまく目的の図が得られるのですが、

i <- 1
p <- ggplot(data,aesthetic=list(x = class, y = paste("col",i,sep=""))
ggjitter(ggplot(p))

とするとうまくいきません。 summary(p)で変数の中身をみてみると、xの値が望みの"col1"ではなく"paste("col",i,sep="")"というように変数展開されていない事が分かりました。 いろいろ試してみたのですが、どうしてもこの問題が解決できません。 どなたかヒントをいただけないでしょうか。

  • paste("col",i,sep="")で変数名が作れるとは思えません。eval(parse(text=***))あたりを使わないといけないのでは? -- okinawa 2007-08-24 (金) 21:34:09
  • 必要なデータの定義などもちゃんと付けて質問しましょう。文字列のcol1 ではなく,オブジェクトの col1 は,eval(parse(text=paste("col",i,sep=""))) で作ることができると思いますよ。 -- 2007-08-24 (金) 21:34:59
  • つまりあなたは x = col1 ではなく x = "col1" と指定したことになっていますね。単純に以下のような指定では駄目なのですか. -- 2007-08-24 (金) 21:47:56
    x = data[[i+1]]  # データフレームなら
    x = data[,i+1]   # 行列なら
  • ご返答ありがとうございます。質問で十分な情報を書かなかったために意図が伝わっておりませんでした。
    dataにはクラス名(0/1の2クラス)が"class"という添字付きで入っているカラムと、正の連続値が"col1"という添字付きで入っているカラムが含まれます。クラス0/1ごとにデータの分布を見ようとしています。
    ggplot関数でx=...の...の部分はカラムの添字を代入するようです(だからx = col1とするとうまく行きました)。なのでxにはオブジェクトではなくカラム名の文字列を代入するものと考えました(文字列もオブジェクトなのでここに矛盾があるのかなと思います)。そこでx=paste(...)とやってみましたがうまくいきませんでした。pasteを使わず、単純にvar <- "col1"; p <- ggplot(... x = var ...)としてもxの値が"col1"にならずに"var"となってしまいます。カラムの添字を指定するのに文字列を代入しているあたりが基本的な事がよくわかっていない部分だと思うのですが何か気付かれる事はあるでしょうか? -- kutakuta? 2007-08-24 (金) 22:07:48
  • かなり意味不明(失礼)。まずカラムなどといういい方はどこかの世界で普通なんでしょうか。列ラベルもしくは変数ラベル(名)でしょうか。またそもそものオブジェクトは行列ではなくデータフレームなんでしょうね.変数名は文字列(names()関数で操作できる)ですが,変数自体は文字列ではありません。ですから x=col1 は適正な指定ですが,x="col1" は x の値を文字列 "col1" にするという意味になってしまいます。一方でデータフレームの列を変数名で取り出す x$col1 という方法以外に、変数名文字列で取り出す x$"col1" (もしくは, x[["col1"]、但し xcol1? は駄目 ) という方法も用意されていますからややこしい(R はかなり御都合主義変態言語です)。 結局あなたの混乱の源は変数名と変数名文字列の混同からきているようにみえるのですが。-- 2007-08-24 (金) 23:32:57
  • 小さいセットで良いので手持ちのデータとあなたがしたことを正しく再現して見せてください。お力になりたいのですが、意味不明です。十分な情報と冗長な文章は別物です。 -- 2007-08-24 (金) 23:46:08
  • あなたのやったことを再現できないんですよ
    > p <- ggplot(data,aesthetic=list(y = class, x = col1))
    > ggjitter(ggplot(p))
    > とするとうまく目的の図が得られるのですが、
    と書いているんだが,そもそもこの前の方のdata がどういう風に作られているかわからない。せいぜい
    data <- data.frame(class=factor(rep(1:4, each=10)),
                       col1=rnorm(40), col2=rnorm(40))
    ぐらいかなと思うが,これじゃエラーメッセージが出てしまい動かんのよね -- 2007-08-25 (土) 07:38:27
  • kutakutaさん、まず str(data) の結果を見せていただけますか? -- okinawa 2007-08-25 (土) 10:02:10
  • 質問の仕方があまりにも拙く、余計なお手数をおかけしてしまっているようで大変恐縮しております。もう一度私の質問を説明するためにサンプルコードを書いてみましたので、これで説明し直させてください。
    ## load library
    library(ggplot)
    
    ## create data
    data <- data.frame(matrix(c(rep(0,10),rep(1,10),rnorm(60)),ncol=4))
    data[,1] <- factor(data[,1])
    colnames(data) <- c("class","col1","col2","col3")
    
    ## do each process (it works).
    pdf("col1_doeach.pdf")
    p <- ggplot(data,aesthetic=list(y = col1,x = class))
    ggjitter(ggboxplot(p))
    dev.off()
    
    pdf("col2_doeach.pdf")
    p <- ggplot(data,aesthetic=list(y = col2,x = class))
    ggjitter(ggboxplot(p))
    dev.off() 
    
    pdf("col3_doeach.pdf")
    p <- ggplot(data,aesthetic=list(y = col3,x = class))
    ggjitter(ggboxplot(p))
    dev.off()
    
    ## use variable (doesn't work)
    column_name <- "col1"
    pdf("col1_var.pdf")
    p <- ggplot(data,aesthetic=list(y = column_name,x = class))
    ggjitter(ggboxplot(p))
    dev.off()
    
    ## use 'for' (doesn't work).
    for (i in 1:3){
     column_label = paste("col",i,sep="")
     p <- ggplot(data,aesthetic=list(y = column_label,x = class))
     pdf(paste(column_label,".pdf",sep="")
     gjitter(ggboxplot(p))
     dev.off()
    }
    以上のコードのうち、"do each"の部分は意図したように動きます。"use variable"の部分は動きません。
    ご回答者の方が指摘されている通り、列のラベルを文字列で指定しようとしているところがまちがっているように思います。
    なぜ変数に列のラベルを入れる事にこだわっているかと言うと、サンプルコードの最後にある"use for"の部分のように複数の行に関するboxplotを作りたいからです。
    質問内容、質問の仕方ともに未熟な部分があり不愉快な思いをされている方もいらっしゃるかもしれませんが、どうかご回答よろしくお願い致します。-- kutakuta? 2007-08-25 (土) 13:56:29
  • 本質的な解決ではないかもしれませんが。ggplotが吐き出すオブジェクトを修正すると動きます。kutakutaさんが例示なさっている,動いた例と動かない例のggplotが吐き出すオブジェクト違いを確認してみて下さい。-- aa? 2007-08-25 (土) 16:14:36
    for (i in 1:3){
     column_label = paste("col",i,sep="")
     p <- ggplot(data,aesthetic=list(y = column_label,x = class))
     p$ylabel <- column_label
     p$defaults[[1]] <- uneval(as.formula(paste(column_label, "~class",collapse="")))[[1]]
     ggjitter(ggboxplot(p))
    }
  • ご返答どうもありがとうございます。
    ご案内いただいた通りにやってみてsummary(p)で確認するとpの中身が期待していたものになっていました。ありがとうございます。
    ただ、私の環境では最後のggjitter(ggplot(p))の部分で
    以下にエラーas.data.frame.default(x[[i]], optional = TRUE) :
          クラス "function" はデータフレームに強制変換できません
    というエラーが出てしまいました
    私自身いただいたソースの意味をまだきちんと理解できていませんので(特に下から3行目)、もう少し勉強したいと思います。
    週末にもかかわらずご返答ありがとうございました。 -- kutakuta? 2007-08-25 (土) 16:54:45
  • p <- ggplot(data,aesthetic=list(y = get(column_name),x = class)) -- surg? 2007-08-25 (土) 16:52:25
  • ご回答ありがとうございます。ご指摘いただいた方法でやってみますと、pの中身が
    > summary(p)
    Title:     
    Labels:    x=class, y=get(column_label)
    Data:      class, col1, col2, col3 [20x4] 
    Defaults:  y=get(column_label), x=class
    Scales:    
    Grobs:      
    Facetting: ". ~ ." 
    Margins:   FALSE 
    という風にget(column_label)の部分が展開されないままとなってしまいます。-- kutakuta? 2007-08-25 (土) 16:54:45
  • p <- ggplot(data,aesthetic=list(y = get(column_name),x = class)); p$ylabel <- column_name -- surg? 2007-08-25 (土) 17:33:17
  • 最初の質問で aesthetic=list(y = class, x = col1) になってたよね(変だと思ったのだけど)
    で,要するにaesthetic=list(y = col1, x=class) なんであって,col1 の1の部分を i=1 などのように変数に入っている数値を使って指定したいと言うことなんでしょう?最初にokinawaさんと私が言ったように,eval と parse を使うということが解決法でしょう。やってみなかったのですか?つまり以下のようになるということ。
    library(ggplot)
    ## create data
    data <- data.frame(matrix(c(rep(0,10),rep(1,10),rnorm(60)),ncol=4))
    data[,1] <- factor(data[,1])
    colnames(data) <- c("class","col1","col2","col3")
    i=3
    p <- ggplot(data,aesthetic=list(y = eval(parse(text=paste("col", i,sep=""))),x = class))
    ggjitter(ggboxplot(p))
    最初から,data を付けて質問すれば,直ぐ解決したことだと思います。。。未だに解決していないのですか? -- 2007-08-25 (土) 18:33:27
  • ご返答ありがとうございます。ご指摘の通りはじめに示した式の中にまちがいがありました(訂正させていただきました)。先日教えていただいた方法はすぐに試させていただきましたがうまくいきませんでした。結果のpの中身を示します。
    > summary(p)
    Title:     
    Labels:    x=class, y=eval(parse(text = paste("col", i, sep = "")))
    Data:      class, col1, col2, col3 [20x4] 
    Defaults:  y=eval(parse(text = paste("col", i, sep = ""))), x=class
    Scales:    
    Grobs:      
    Facetting: ". ~ ." 
    Margins:   FALSE 
    なにぶん初心者なものでプログラミングのスジ自体が良くないのかもしれません(私が初心者の頃に書いたCやJavaのコードは今見ると基本的なグランドデザインがダメなために余計なことばかりしていますから)。
    私自身スジの悪さは認識しておりますので、ご回答者様が不愉快な思いをされていることも理解できます。質問の仕方にも問題がありましたし。これ以上ご迷惑をかけるのは心苦しいので、今回の質問に関しましてはここで取り下げさせていただきたく存じます(ここまでのやり取りは他のご回答者様からいただいたものもありますので私自身はこのままにさせていただきます)。申し訳ありませんでした。 -- kutakuta? 2007-08-25 (土) 23:49:59
  • kutakutaさん,お役に立てず申し訳ございませんでした。なお,私がコメントさせて頂いた部分の「ggplotが吐き出すオブジェクト違いを確認」とは,str(p)とすることを意味します。kutakutaさんが,理解できなかったとおっしゃっていた,下から3行目の部分は,ggplot.defaultから,参考に作成した物です。他の方の方法でも,同一かと思いますが,ggplotの作成するオブジェクトのリストの成分の成分名を,目的のものにすることがポイントかと思います。kutakutaさんの,ご質問のおかげで,library(ggplot) の存在を知り,大変勉強になりました。お役に立てず恐縮でしたが,誠にありがとうございました。 -- aa? 2007-08-26 (日) 01:25:06
  • おかしいですね。なぜ動かないのでしょう。
    別に回答者は不快に思っているわけではないと思いますよ。こういう風に質問してくれたら回答しやすいのにと,もどかしくは思っておりますが。それと,提案された方法を試したら,その結果を直ぐにコメントするようにすれば,やりとりがスムーズに行くと思いますが。先ほど示したプログラムをそのままコピーして,コンソールにペーストしてもうまく動かなかったと言うことなんでしょうか。
    > library(ggplot)
     要求されたパッケージ grid をロード中です
     要求されたパッケージ reshape をロード中です
     要求されたパッケージ RColorBrewer をロード中です
    
     次のパッケージを付け加えます: 'ggplot'
    
    
    	The following object(s) are masked _by_ .GlobalEnv :
    
    	 alpha 
    
    > ## create data
    > data <- data.frame(matrix(c(rep(0,10),rep(1,10),rnorm(60)),ncol=4))
    > data[,1] <- factor(data[,1])
    > colnames(data) <- c("class","col1","col2","col3")
    > i=3
    > p <- ggplot(data,aesthetic=list(y = eval(parse(text=paste("col", i,sep=""))),x = class))
    > ggjitter(ggboxplot(p))
    > summary(p)
     Title:     
     Labels:    x=class, y=eval(parse(text = paste("col", i, sep = "")))
     Data:      class, col1, col2, col3 [20x4] 
     Defaults:  y=eval(parse(text = paste("col", i, sep = ""))), x=class
     Scales:    
     Grobs:      
     Facetting: ". ~ ." 
     Margins:   FALSE 
    > str(p)
    List of 11
     $ data       :'data.frame':	20 obs. of  4 variables:
      ..$ class: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
      ..$ col1 : num [1:20]  0.205  1.523  0.946 -0.496  0.890 ...
      ..$ col2 : num [1:20]  0.0696 -1.5451  1.1064 -1.4103 -0.3079 ...
      ..$ col3 : num [1:20] 0.527 0.828 0.972 0.352 1.200 ...
     $ grobs      : list()
     $ scales     : list()
      ..- attr(*, "class")= chr "scales"
     $ defaults   :List of 2
      ..$ y: language eval(parse(text = paste("col", i, sep = "")))
      ..$ x: symbol class
     $ title      : chr ""
     $ fixedaspect: logi FALSE
     $ xlabel     : chr "class"
     $ ylabel     : chr "eval(parse(text = paste(?"col?", i, sep = ?"?")))"
     $ formula    : chr ". ~ ."
     $ margins    : logi FALSE
     $ facet      :List of 1
     $ : num 0
     - attr(*, "dim")= int [1:2] 1 1
     - attr(*, "dimnames")=List of 2
      ..$ : chr "value"
      ..$ : chr "value"
     - attr(*, "rdimnames")=List of 2
      ..$ :'data.frame':	1 obs. of  1 variable:
      .. ..$ value: Factor w/ 1 level "value": 1
      ..$ :'data.frame':	1 obs. of  1 variable:
      .. ..$ value: Factor w/ 1 level "value": 1
     - attr(*, "idvars")= chr "value"
     - attr(*, "class")= chr "ggplot"
    こんな風になって,以下のような図が描けましたけど。 -- 2007-08-26 (日) 19:56:25
    ggplot99.png
  • たぶん望む物が2つ重なって質問されているのではないでしょうか。(1)y軸ラベルがcol1,col2,col3になっててほしい(2)作画をループで回したい。で、(1)ですが、
    i<-3 # = より <- を使いましょう
    text<-paste("p <- ggplot(data,aesthetic=list(y = ","col",i,", x=class ))",sep="")
    eval(parse(text = text))
    ggjitter(ggboxplot(p))
    のようにするとy軸ラベルが望む物になるとおもいます。また、(2)ですが、試したところggplot(),ggboxplot(),ggjitter()のオブジェクト構造が思いのほか複雑で、単純なループやベクトルの代入ではうまく複数処理ができないようです。(plot()関数と同じ仕様だと私自身勘違いしていたようです。ggplot()なんて紛らわしい名前なので)-- okinawa 2007-08-26 (日) 21:40:46
  • まあ,確かに ggplot は変な仕様ですね。ついでに言えばデータポイントをgitter で描くというのはあまりお勧めではないです。ランダムなデータポイントの配置が意味ありげなものに見えることがあるのと,逆に意味もなく散らばりすぎることがあること。上の作図例で,データポイントがx軸方向に無意味にぶれて描画されているのがやなかんじ。普通のボックスプロットにpoints 関数でデータポイントを重ねる方がましなような気がする。付値に=を使ったのは意図的ではないが,まあ,かんべんして! -- 2007-08-26 (日) 22:17:18
  • boxplot と points を使って描くと以下のようなもの。ラベルとか,外れ値の記号・色を変えるとかは更に工夫が必要。 -- 2007-08-26 (日) 22:40:30
    set.seed(1234)
    data <- data.frame(class=factor(rep(0:1,each=10)), col1=rnorm(20), col2=rnorm(20), col3=rnorm(20))
    layout(matrix(1:3, ncol=3))
    par(cex=1.2)
    for (i in 1:3) {
    	boxplot(eval(parse(text=paste("col", i, sep="")))~class, data, outline=FALSE)
    	points(jitter(as.numeric(data$class)), eval(parse(text=paste("data$col", i, sep=""))), col="blue", pch=20, cex=2)
    }
    ggplot98.png
    なお,上のような場合だと,
    boxplot(data$class, data[,i+1], outline=FALSE)
    points(jitter(as.numeric(data$class)), data[,i+1], col="blue", pch=20, cex=2)
    とする方が,簡明・明快
  • good job! 質問と乖離しましたが、ggplotを使わないほうがいいみたいですね。教訓「packageを使う前に通常の方法で実現できないか試せ」ですか・・・。kutakutaさんはggplotから入られたので、かえって難しくなってしまったのかもしれませんね。 -- okinawa 2007-08-27 (月) 08:01:14
  • うまく行きました。皆さまありがとうございました。 -- kutakuta? 2007-08-27 (月) 08:22:23

Rコマンダーの日本語表示について

J? (2007-08-23 (木) 20:55:38)

Mac OS 10.4.10 上で R ver.2.5.1 とRcmdr 1.3-5 を使っています。
Rコマンダーの日本語メニュー表示がやや読みにくいのですが,これはfontの設定等で改善できるものでしょうか?

  • どういう風に読みにくいのでしょうか?また,もし「fontの設定等で改善できるものでしょうか」ということなら,fontの設定などを変更してみたのでしょうか(何をどのように変更してみたのでしょうか)。回答に値する質問をしましょうよ。 -- 2007-08-23 (木) 22:40:19
  • 言葉足らずで失礼しました。メニュー等の日本語表示で文字の輪郭が粗くなっており,細かい部分が判読し辛いという意味でした。
    オプションでdefault fontを変えてみようと試みましたが失敗しておりました。また自分で調べてみます。ご迷惑おかけしました。 -- 2007-08-24 (金) 09:39:56
  • 我流ですが...
    1. ttfを持ってくる
       http://sourceforge.jp/projects/efont/ から sazanamiあたりがいいかと.
       tar jxvf sazanami-20040629.tar.bz2
    
    2. 自前用のフォントの格納ディレクトリを作成
       mkdir -p /usr/X11R6/lib/X11/fonts/japanese
    
    3. フォントのコピー
       cp sazanami*.ttf /usr/X11R6/lib/X11/fonts/japanese/
    
    4. /usr/X11R6/lib/X11/fonts/japanese/fonts.dirの作成
    ==================== ここから ==============================
    32
    vl=y:sazanami-mincho.ttf               -sazanami-mincho-medium-r-normal--0-0-0-0-c-0-iso8859-1
    vl=y:sazanami-mincho.ttf               -sazanami-mincho-medium-r-normal--0-0-0-0-c-0-jisx0201.1976-0
    vl=y:sazanami-mincho.ttf               -sazanami-mincho-medium-r-normal--0-0-0-0-c-0-jisx0208.1983-0
    vl=y:sazanami-mincho.ttf               -sazanami-mincho-medium-r-normal--0-0-0-0-c-0-jisx0212.1990-0
    vl=y:ai=0.2:sazanami-mincho.ttf        -sazanami-mincho-medium-i-normal--0-0-0-0-c-0-iso8859-1
    vl=y:ai=0.2:sazanami-mincho.ttf        -sazanami-mincho-medium-i-normal--0-0-0-0-c-0-jisx0201.1976-0
    vl=y:ai=0.2:sazanami-mincho.ttf        -sazanami-mincho-medium-i-normal--0-0-0-0-c-0-jisx0208.1983-0
    vl=y:ai=0.2:sazanami-mincho.ttf        -sazanami-mincho-medium-i-normal--0-0-0-0-c-0-jisx0212.1990-0
    vl=y:ds=y:sazanami-mincho.ttf          -sazanami-mincho-bold-r-normal--0-0-0-0-c-0-iso8859-1
    vl=y:ds=y:sazanami-mincho.ttf          -sazanami-mincho-bold-r-normal--0-0-0-0-c-0-jisx0201.1976-0
    vl=y:ds=y:sazanami-mincho.ttf          -sazanami-mincho-bold-r-normal--0-0-0-0-c-0-jisx0208.1983-0
    vl=y:ds=y:sazanami-mincho.ttf          -sazanami-mincho-bold-r-normal--0-0-0-0-c-0-jisx0212.1990-0
    vl=y:ds=y:ai=0.2:sazanami-mincho.ttf   -sazanami-mincho-bold-i-normal--0-0-0-0-c-0-iso8859-1
    vl=y:ds=y:ai=0.2:sazanami-mincho.ttf   -sazanami-mincho-bold-i-normal--0-0-0-0-c-0-jisx0201.1976-0
    vl=y:ds=y:ai=0.2:sazanami-mincho.ttf   -sazanami-mincho-bold-i-normal--0-0-0-0-c-0-jisx0208.1983-0
    vl=y:ds=y:ai=0.2:sazanami-mincho.ttf   -sazanami-mincho-bold-i-normal--0-0-0-0-c-0-jisx0212.1990-0
    vl=y:sazanami-gothic.ttf               -sazanami-gothic-medium-r-normal--0-0-0-0-c-0-iso8859-1
    vl=y:sazanami-gothic.ttf               -sazanami-gothic-medium-r-normal--0-0-0-0-c-0-jisx0201.1976-0
    vl=y:sazanami-gothic.ttf               -sazanami-gothic-medium-r-normal--0-0-0-0-c-0-jisx0208.1983-0
    vl=y:sazanami-gothic.ttf               -sazanami-gothic-medium-r-normal--0-0-0-0-c-0-jisx0212.1990-0
    vl=y:ai=0.2:sazanami-gothic.ttf        -sazanami-gothic-medium-i-normal--0-0-0-0-c-0-iso8859-1
    vl=y:ai=0.2:sazanami-gothic.ttf        -sazanami-gothic-medium-i-normal--0-0-0-0-c-0-jisx0201.1976-0
    vl=y:ai=0.2:sazanami-gothic.ttf        -sazanami-gothic-medium-i-normal--0-0-0-0-c-0-jisx0208.1983-0
    vl=y:ai=0.2:sazanami-gothic.ttf        -sazanami-gothic-medium-i-normal--0-0-0-0-c-0-jisx0212.1990-0
    vl=y:ds=y:sazanami-gothic.ttf          -sazanami-gothic-bold-r-normal--0-0-0-0-c-0-iso8859-1
    vl=y:ds=y:sazanami-gothic.ttf          -sazanami-gothic-bold-r-normal--0-0-0-0-c-0-jisx0201.1976-0
    vl=y:ds=y:sazanami-gothic.ttf          -sazanami-gothic-bold-r-normal--0-0-0-0-c-0-jisx0208.1983-0
    vl=y:ds=y:sazanami-gothic.ttf          -sazanami-gothic-bold-r-normal--0-0-0-0-c-0-jisx0212.1990-0
    vl=y:ds=y:ai=0.2:sazanami-gothic.ttf   -sazanami-gothic-bold-i-normal--0-0-0-0-c-0-iso8859-1
    vl=y:ds=y:ai=0.2:sazanami-gothic.ttf   -sazanami-gothic-bold-i-normal--0-0-0-0-c-0-jisx0201.1976-0
    vl=y:ds=y:ai=0.2:sazanami-gothic.ttf   -sazanami-gothic-bold-i-normal--0-0-0-0-c-0-jisx0208.1983-0
    vl=y:ds=y:ai=0.2:sazanami-gothic.ttf   -sazanami-gothic-bold-i-normal--0-0-0-0-c-0-jisx0212.1990-0
    ==================== ここまで ==============================
    
    5. フォントパスの追加
    xset +fp reash の追加
    --- /etc/X11/xinit/xinitrc.orig
    +++ /etc/X11/xinit/xinitrc
    @@ -24,6 +24,9 @@
         xmodmap "$usermodmap"
     fi
     
    +xset +fp /usr/X11R6/lib/X11/fonts/japanese
    +xset reash
    +
     # start some nice programs
     
     xterm &
    
    otfだとTTCapのイタリックが使えないので, ttfがR的にはおすすめかも. -- なかま 2007-08-24 (金) 14:02:19
  • RにおけるMacOSXでのtcltkはX11を使います. X11のフォントの設定はデフォルトでは古典的なフォントしか無いので変更したい人が多いかと思います. fontの設定をしているような方なら,部分的に(フォントパスを追加したいとか)質問もされるでしょうが, この問題を回答に値する質問の仕方というのは難しいように思います. むしろ,「OSXでRcmdrのフォントの変更方法は?」としか言えないように思います. 質問者の方は出来れば「一般の方にも理解可能な日本語!」に編集していただければ迷惑にはなりません. -- なかま 2007-08-24 (金) 14:22:55
  • 仰る通り,私が知りたかったのは「OSX上でRcmdrのフォントを変更する方法」でした。 ご教示頂いたやり方は,私には少し敷居が高いのですが勉強してみます。ありがとうございました。 -- J? 2007-08-27 (月) 20:17:00

ファイル出力加工

Bun? (2007-08-23 (木) 20:21:21)

R初心者です。分かりづらい標記でしたら、ご指摘の程よろしくお願いいたします。

以下のプログラムで、クロス集計の結果を外部出力しました。

x<-c(1,2,1,2,1,2)
y<-c(1,1,1,2,2,2)
z<-data.frame(x,y)
z[,1] <- factor(z[,1], levels=1:2, labels=c("男", "女"))
XTBL<-table(z[,1],z[,2])
write.table(XTBL, file="C:/XTBL.txt", append=F, quote=F, col.names=T)

下記は、出力結果

   1 2
男 2 1
女 1 2

上記の結果ファイルをエクセルで読み込み、作業したいのですが、Rからの出力のデータファイル時に、男女の上に変数名を付けて、出力することは可能でしょうか?

  • どのような結果を期待しているのかよくわからない。以下のようなことなのか? -- 2007-08-23 (木) 22:55:28
    > x<-c(1,2,1,2,1,2)
    > y<-c(1,1,1,2,2,2)
    > x <- factor(x, levels=1:2, labels=c("男", "女"))
    > y <- factor(y, levels=1:2, labels=c("あり", "なし"))
    > XTBL<-table(x,y)
    > XTBL
        y
    x    あり なし
      男    2    1
      女    1    2
    > 上記の結果ファイルをエクセルで読み込み、作業したいのですが
    なんで,Excel なんかで読んで作業するんだろう
  • fix(XTBL)でデータエディタで見ると、row.namesという変数名が見て取れますが、Rから出力した結果を、エクセルで読み込むと、1行目が列づれを起こしてしまいます。エクセルで読み込む際に、列づれを起こさないように加工することはできますでしょうか。 -- Bun? 2007-08-23 (木) 23:54:18
  • 次のようにすればそれらしくなりますが。 -- 2007-08-24 (金) 00:19:41
    > z <- data.frame(x,y)
    > z[,1] <- factor(z[,1], levels=1:2, labels=c("男", "女"))
    > XTBL <- table(z, dnn=c("sex",""))
    > XTBL
        
    sex  1 2
      男 2 1
      女 1 2
  • ご回答ありがとうございます。テキストファイルの生成時に、依然ラベルが吐き出されずに列づれを起こしている状態でした。 -- Bun? 2007-08-24 (金) 11:49:47
  • ?write.tableでcol.namesのoptionを見ましょう。Bunさんの要望はおそらく
    write.table(z, firle="hoge.csv", col.names=NA, sep=",")
    です。ExcelでできるならExcelですれば良いと思います。私は面倒なので使いませんが。
    それにしてもhelpを読まない人が多すぎます... -- 2007-08-24 (金) 12:23:51
  • ↑(内心)同感ですが、一方で計算機ソフトのヘルプ文章は、ある程度知っている人でないとそもそも意味がわからないことが多いのも真実。 -- 2007-08-24 (金) 14:13:24
  • 解決できました。どうもありがとうございました。 -- Bun? 2007-08-24 (金) 14:39:06

Rコマンドファイル内からディレクトリ作成は出来ますか?

みゅー? (2007-08-23 (木) 16:21:27)

R初心者につき言葉足らずな質問になっていたらごめんなさい。ご指摘ください。

大量のファイルが生成されるため種類ごとに出力を分けたいと思いました。
シェルコマンドのmkdirに対応するようなRコマンドはありますか?
みつけられませんでした。
(ls()やrmは使えるようですね。)

単に、./以下にhogeディレクトリを作成していない状態のまま"./hoge/filename.pdf"とやってももちろん./hogeは自動生成しませんでした。
先に./hogeを用意しておかないとだめ、
ディレクトリ生成はできないということでしょうか。

  • shell("mkdir hoge")でたいがいのシステムではいけるんじゃないかな。ちなみにlsやrmはシェルコマンドのそれとは違うのでご注意を。 -- takahashi? 2007-08-23 (木) 16:42:38
  • ?dir.create -- なかま 2007-08-23 (木) 16:48:10
  • dir.createでいけました ありがとうございます  -- みゅー? 2007-08-23 (木) 17:07:16
  • shell は通りませんでした。 関数"shell"を見つけることが出来ませんでした といわれてしまいます。 どうすればこの方法がうまく通るでしょうか。 また、他にもこのようにUnixコマンドを通すことができる関数はありますか? Rコマンドファイル内からUnixコマンドを通せるようになるととても便利になるのですが・・・。 -- みゅー? 2007-08-23 (木) 17:22:32
  • 追加自己解決です。windowsならshell()でいけますね。unixならsystemですね。 -- みゅー? 2007-08-23 (木) 17:25:27
  • Windowsでもsystem()は使えます。DOS系コマンドはshell()でそれ以外はsystem()のような気がします。<未確認 -- 2007-08-23 (木) 19:36:30

変数名の修正

Bun? (2007-08-23 (木) 15:22:23)

既に変数名を含んだデータセットについて、一部の変数名を新しい変数名にしたい場合、どのような関数、プログラムを利用すればよいのでしょうか。

  • 簡単な具体例を挙げてください。回答者に余分な推測をさせないこと。 -- 2007-08-23 (木) 15:32:11
  • colnames(df)[4] <- "new.name" -- 2007-08-23 (木) 15:39:36
  • 元のデータファイルの変数名を書き換える -- 2007-08-23 (木) 15:46:41
  • ありがとうございます。2番目の方の回答で解決できました。 -- Bun? 2007-08-23 (木) 19:07:02

関数実行時の入力文字列の取得方法

r_user? (2007-08-18 (土) 10:14:36)

たとえば、f(x)と入力して関数fを実行すると、関数にはオブジェクトxの中身が渡されますが、
中身だけでなく、文字列xも取得して利用できるような関数を定義したいと考えています
引用符を付けてf("x")のようにするのではなく、普通にf(x)と入力して、
xの中身に加えて、"x"という文字列も取得して利用できないものでしょうか。

  • f<-function(...)print(paste(match.call())[-1]) -- takahashi? 2007-08-18 (土) 11:39:35
  • 少し趣旨が違いますが,Rで仮引数の(未評価)実引数を取り出すには substitute 関数を使います.但し結果は文字列ではなく language オブジェクトです。 -- 2007-08-18 (土) 11:40:16
    > foo <- function(x) {substitute(x)}
    > y <- foo(1:10)
    > str(y)
     language 1:10
    > eval(y)
     [1]  1  2  3  4  5  6  7  8  9 10
  • ご教授ありがとうございます。match.call()、substitute()ともに、これからお世話になりそうです。色々と応用可能性を試してみます。 -- r_user? 2007-08-18 (土) 12:03:43
  • もう少し簡略に(htest 属性を持つ stats の関数で常用されています) -- 2007-08-19 (日) 07:08:05
    bar <- function(x) return(deparse(substitute(x)))

条件によるIDの振り分けについて

suzuki? (2007-08-17 (金) 17:17:46)

はじめまして。こんにちは。
データによるIDの振り分けについて質問させてください。

DATA1	 DATA2	 DATA3
1	 30	 1
1	 30	 1
1	 30	 1
1	 30	 2
1	 30	 2
1	 31	 1
1	 31	 1
1	 31	 1
1	 31	 1
2	 30	 1
2	 30	 2
2	 30	 2
2	 30	 2
2	 30	 3
2	 30	 3


以上のようなデータがあって、以下のようにDATA1、DATA2、DATA3をキーとして全ての値が同じであれば、同じIDとするように
IDを振り分けていきたいのですが、各IDの行数がランダムであるためうまく処理できる方法が見つかりません。
ある程度、手作業でできなくもないのですが、20万行ほどのデータであるため、何とか自動化できないものかと考えています。
アドバイスよろしくお願いします。

DATA1	 DATA2	 DATA3	 ID
1	 30	 1	 1
1	 30	 1	 1
1	 30	 1	 1
1	 30	 2	 2
1	 30	 2	 2
1	 31	 1	 3
1	 31	 1	 3
1	 31	 1	 3
1	 31	 1	 3
2	 30	 1	 4
2	 30	 2	 5
2	 30	 2	 5
2	 30	 2	 5
2	 30	 3	 6
2	 30	 3	 6
  • 美しくないしコストも考えてないけど、dはデータフレームで id<-1;d$ID<-1:dim(d)[1];aggregate(d$ID,list(d$DATA1,d$DATA2,d$DATA3),function(i){d[i,]$ID<<-id;id<<-id+1}) -- takahashi? 2007-08-17 (金) 17:50:32
  • 以下のようにするのはいかが?このやり方だと,フィールドが幾つあっても簡単。区切り文字は適切に選ぶ必要がある。 -- 2007-08-17 (金) 21:08:59
    > y <- apply(x, 1, paste, collapse="/")
    > z <- names(table(y))
    > x["ID"] <- sapply(y, function(i) which (i == z)) 
    > x
       DATA1 DATA2 DATA3 ID
    1      1    30     1  1
    2      1    30     1  1
    3      1    30     1  1
    4      1    30     2  2
    5      1    30     2  2
    6      1    31     1  3
    7      1    31     1  3
    8      1    31     1  3
    9      1    31     1  3
    10     2    30     1  4
    11     2    30     2  5
    12     2    30     2  5
    13     2    30     2  5
    14     2    30     3  6
    15     2    30     3  6
  • もし DATA1,DATA2,DATA3 が全て非負整数値であまり大きくなければ次の方法が使えます. -- 2007-08-17 (金) 22:16:25
    > a <- max(x)+1
    > xx <- transform(x, ID=as.factor(x[,1]*a^2+x[,2]*a+x[,3]))
    > xx
       DATA1 DATA2 DATA3   ID
    1      1    30     1 1985 # 数値に拘らなければこれだけでもOK
    2      1    30     1 1985 
    3      1    30     1 1985
    4      1    30     2 1986
    5      1    30     2 1986
    6      1    31     1 2017
    7      1    31     1 2017
    8      1    31     1 2017
    9      1    31     1 2017
    10     2    30     1 3009
    11     2    30     2 3010
    12     2    30     2 3010
    13     2    30     2 3010
    14     2    30     3 3011
    15     2    30     3 3011
    > levels(xx[,4]) <- 1:nlevels(xx[,4]) # もし1,2,3,...というID番号の方がよければ
    > xx      
       DATA1 DATA2 DATA3 ID  # ただしIDは因子です
    1      1    30     1  1
    2      1    30     1  1
    3      1    30     1  1
    4      1    30     2  2
    5      1    30     2  2
    6      1    31     1  3
    7      1    31     1  3
    8      1    31     1  3
    9      1    31     1  3
    10     2    30     1  4
    11     2    30     2  5
    12     2    30     2  5
    13     2    30     2  5
    14     2    30     3  6
    15     2    30     3  6
    > xx[,4] <- as.integer(xx[,4]) # もし因子でなく整数値の方が良ければ
    時間比較例:
    > set.seed(1)
    > x <- data.frame(DATA1=sample(1:10, 2e5, rep=TRUE),
                      DATA2=sample(20:40, 2e5, rep=TRUE),  
                      DATA3=sample(10:20, 2e5, rep=TRUE))
    > a <- max(x)+1
    > system.time({xx <- transform(x, ID=as.factor(x[,1]*a^2+x[,2]*a+x[,3])); 
                   levels(xx[,4]) <- 1:nlevels(xx[,4])})
       ユーザ  システム      経過 
        0.196     0.124     0.929 
    > system.time({y <- apply(x, 1, paste, collapse="/"); z <- names(table(y));
                  x["ID"] <- sapply(y, function(i) which (i == z)) })
       ユーザ  システム      経過 
      147.873     0.656   154.912
  • 本当にありがとうございます。自分が組んでいたプログラムより圧倒的に早く処理できました。 -- suzuki? 2007-08-18 (土) 14:24:38

Rcmdr QQプロット

KT? (2007-08-15 (水) 18:25:40)

Rcmdrで2つのデータをQQプロットに同時にのせたいと思います。可能でしょうか?

  • RcmdrはRにGUI環境をつけただけと思うので、Rでできるならできるのでは?情報が少なく、こんな一般的な返事しかできません。 -- 2007-08-15 (水) 21:00:23
  • Rcmdr でできるかどうかは問題じゃなくて,R でできるかどうかが重要だが,要求がはっきりしないので回答できない。どんなデータがあって,どんなQQプロットが望みなのか。(まあ,QQプロットは分かるんだけど)(Rcmdrでということから,Rプログラムを書くのは難しいということなのかなあとか) -- 2007-08-15 (水) 23:17:54
  • 色をかえて、2つのデータを同時に1つの図の上に載せたいということでしたら、Rcmdr上では無理です。1つずつ描かせたときにスクリプトウィンドウに出ているqq.plot(...) 文をR Consoleにコピペして実行し、1つ目のプロットと2つ目のプロットの間に > par(new=T) という一文を実行してあげればかけることはかけますが、軸も色も重なってしまうので望むものではないでしょう。 -- 2007-08-16 (木) 11:00:08
  • それらしいものを書くには、スクリプトウィンドウに出ているqq.plot() ではなく、qqnorm()関数を直接コンソールから打ち込むほうがいいでしょうね。データが変数a, b に入っているとして、
    > qqnorm(a,xlim=c(-3,3),ylim=c(min(c(a,b)),max(c(a,b))),col="red")
    > par(new=T)
    > qqnorm(b,xlim=c(-3,3),ylim=c(min(c(a,b)),max(c(a,b))),col="blue")
    とかにすると、軸を一致させたプロットになります。 でも、直線にだいたいのっているかどうかを比べてみたいだけなら、
    > qqnorm(a,xlim=c(-3,3),yaxt="n",col="red")
    > par(new=T)
    > qqnorm(b,xlim=c(-3,3),yaxt="n",col="blue")
    というのもありかもしれませんね。-- 2007-08-16 (木) 11:08:08
  • どうもありがとうございました。Rcmdrでは無理とのことですので、 -- KT? 2007-08-17 (金) 14:48:18
  • どうもありがとうございました。早速試してみます。 -- KT? 2007-08-17 (金) 14:54:31

RのグラフをPower point上で変更

sh? (2007-08-08 (水) 01:40:14)

Rのグラフ(軸ラベルのフォント、プロット記号)をPower point上で変更することは可能でしょうか?

  • Windowsならば、emf(メタファイル)でCOPY&PASTEすれば可能でしょう。ppでグループ解除すればオブジェクトがバラケます。(あなたの実行環境を明示してください。RはWindows,MacOSX,Linuxで動作します。回答者はあなたの環境がわからないので、的確な回答ができません) -- okinawa 2007-08-08 (水) 09:06:09
  • emfファイルで保存するにはどうすればよろしいでしょうか?実行環境はMac OS X 10.4.9です。 -- sh? 2007-08-11 (土) 10:38:02
  • macosxでemfは使えません。 -- okinawa 2007-08-11 (土) 10:47:38
  • macosxにおいて、RのグラフをPower point上で変更する方法はないでしょうか? -- sh? 2007-08-12 (日) 02:35:04
  • Power Point なんか使うのやめたら? -- 2007-08-12 (日) 11:15:05
  • 前コメント者のような極端な意見はともかく,フォントやプロット記号なら,labels()とかaxes()とかplot()とかのパラメータ(詳しくはヘルプ参照)でいくらでも変更が効くので,事前にR上で変更したものを使ったほうが幸せになれると思うのですが?? -- 2007-08-13 (月) 00:24:26
  • 確かに,Rで図を作成する際は,可能な限りラベルやプロット記号をRで処理してから,Power Pointにコピーした方が楽ですよ。Power Pointを全く使わないというのは,多くの方には無駄な時間を要することになりますので,あまり気にしない方が良いです。 -- 2007-08-13 (月) 00:37:19
  • Macに限らずLinuxでもそうですが、ベクター型のグラフィック出力の基本はps(eps),PDFだと思います。それを修正するのであれば、それなりのツール(CANVAS,イラストレータ,Acrobat,GhostScript?など)を使わないとできませんので、画像変換の知識とお金が必要になります。そういう意味で、R内で処理をしてPNGやJPEGに出力し、pptに貼付けした方が安上がりということだと思いますよ。 -- okinawa 2007-08-13 (月) 08:44:58
  • ちなみに,Macなら,PDFで出力して貼り付けるほうがきれいです(JPEGというのはグラフには向かないし,PNGはジャギーです)。MacのPowerPoint?だとpdf画像ファイルを貼り込むと,背景の色は(白であっても)透明にならない(白として貼り込まれる)が,PowerPoint?代替品のKeynoteだと,ちゃんと透明になるので好都合。 -- 2007-08-13 (月) 15:22:53
  • pptにPDFインポートできるとは知らなかった・・・。keynoteではできることは確認してましたが。手元にmac版pptがないと答えるのがつらいですね。 -- okinawa 2007-08-13 (月) 18:43:06
  • とどめというか,pdf 関数の引数の family, fonts をよく調べるとよいでしょう。 -- 2007-08-13 (月) 22:42:45
  • epsをtextediterで直接編集という手もありますね。玄人向けか・・・。 -- okinawa 2007-08-14 (火) 08:18:04
  • ベクター型のグラフィック出力というと、SVG出力で Inkspace はどうかな? -- 2007-08-14 (火) 17:43:03
  • inkscapeですよね。RSvgDevice?パッケージでdevSVG("***")でSVG作成して・・・。これはすごい!ちゃんと日本語も使えますね。PDFにも出力できるし、これいいですよ。使えますよ。mac_ppt上で編集という本来の目的とは離れてるけど。 -- okinawa 2007-08-15 (水) 09:02:26
  • 結局,質問者は,何をどうしたかったんでしょうねえ(何の返答もないし,わからんなあ) -- 2007-08-16 (木) 00:06:50
  • RのグラフをPower point上で変更したかったのですが、事前にRで処理してからPower Pointにコピーすることにします(何の返答もなく、大変失礼しました) -- sh? 2007-08-16 (木) 03:23:17

同じ要素を持つ行を選択するには

ななし? (2007-08-08 (水) 01:14:02)

次のようなデータがあったとき

> x <- matrix(c("a","b","c","c","e","d","a","c","c","c","c","c"), c(4,3), byrow=T)
> x <- data.frame(x)
> x
  X1 X2 X3
1  a  b  c
2  c  e  d
3  a  c  c
4  c  c  c

X1が等しくかつX3が等しい行を選択したい(この場合は1行と3行)のですがどのようにすればよいのでしょうか?

>  which( x[,1]==x[,1] & x[,3]==x[,3])
[1] 1 2 3 4


>  x[,4] <- 1:nrow(x) #自分と同じ行を除けばいけるかも…
>  which( x[,4]!=x[,4] & x[,1]==x[,1] & x[,3]==x[,3])
integer(0)


> which(duplicated(x[,1]) & duplicated(x[,3]))
[1] 3 4

では(当然ながら)意図する結果を得られませんでした。

単純にforで1行ずつ回すしか無いのでしょうか?

  • x[is.element(x[,1], x[,3]),] でどうでしょう -- aa? 2007-08-08 (水) 01:45:06
  • ありがとうございます。ですがこれは偶然一致しただけですね。 -- ななし? 2007-08-08 (水) 02:03:28
    > x <- matrix(c("a","b","f","c","e","d","a","c","f","c","c","f"), c(4,3), byrow=T)
    > x <- data.frame(x)
    > x
      X1 X2 X3
    1  a  b  f
    2  c  e  d
    3  a  c  f
    4  c  c  f
    でも同じなのですがx[is.element(x[,1], x[,3]),]ではうまくいきません
  • この場合 is.element(x[,1],x[,3]) は "a","c","a","c" の各々が集合 c("c","d","c","c") 中にあるかどうかを示す論理値をあたえますから、たまたまこの例では正しい結果になったというだけですね。結局総当たりで調べるしか仕方がない問題ですが、工夫で for ループを無くすなら:-- 2007-08-08 (水) 02:18:32
    > y <- outer(x[,1], x[,1], FUN="==") * outer(x[,3], x[,3], FUN="==")
    > diag(y) <- 0
    > y
         [,1] [,2] [,3] [,4]
    [1,]    0    0    1    0
    [2,]    0    0    0    0
    [3,]    1    0    0    0
    [4,]    0    0    0    0
    > which(y == 1, arr.ind=TRUE)
         row col
    [1,]   3   1     # 1,3 が求める行という意味
    [2,]   1   3     # 重複するのが目障りだが
  • 重複が目障りというなら diag(y) <- 0 のかわりに,y[upper.tri(y, diag=TRUE)] <- 0 とするだけ -- 2007-08-08 (水) 08:10:05
  • なるほど外積を使うのですね。条件が三つになったときも外積を追加すれば問題なくいけました。ありがとうございました。 -- ななし? 2007-08-08 (水) 10:15:18
  • 外積を使うプログラムの計算所要時間はデータフレームの行数の二乗に比例しします。以下に示すプログラムは行数に比例します。 -- 2007-08-08 (水) 12:26:31
    y <- paste(x[,1], x[,3], sep="@") # sep の設定には注意
    z <- y %in% names(table(y))[table(y)>1]
    x[z,]
    一回の処理に要する秒数
       行数 外積法  別解
    1     4  0.002 0.001
    2     8  0.001 0.001
    3    16  0.001 0.000
    4    32  0.001 0.001
    5    64  0.004 0.001
    6   128  0.009 0.001
    7   256  0.030 0.001
    8   512  0.346 0.001
    9  1024  1.075 0.002
    10 2048  3.931 0.003
    11 4096 13.471 0.005
  • いったん行を一つの文字列にして比較するわけですね。ありがとうございました。 -- ななし? 2007-08-08 (水) 14:05:11
  • 既に賞味期限切れだが。複数変数を連結して一文字にするというアイデアは素晴らしいが、示唆された方法では重複ケースが一括して取り出され、具体的にどことどこが同じかわかりにくい困難がある。アイデアを活かして次のような方法を考えてみた。因子化すると自動的に同じ値がグループ化されるという機能(8月17日の記事参照)を使っている。 -- 2007-08-21 (火) 22:18:59
    > set.seed(1);x <- paste(sample(letters[1:5],3e4,1), sample(letters[1:5],3e4,1), sep="")
    > df <- data.frame(a=x[1:1e4], b=x[(1e4+1):2e4], c=x[(2e4+1):3e4]) # テスト用データフレーム
    > y <- paste(df[,1],df[,2],df[,3],sep="`") # 3変数を連結
    > y[y %in% setdiff(y, y[duplicated(y)])] <- NA # 一回しか登場しない文字列をNA値で置き換え
    > y <- as.factor(y) # 因子化
    > levels(y) <- 1:nlevels(y) # 水準を変更(NA値は保存されることを注意)
    > y[1:20]
     [1] <NA> <NA> <NA> <NA> <NA> <NA> 1855 1694 <NA> <NA> 525  399  1292 671  <NA>
    [16] <NA> 1599 <NA> 838  <NA>
    2148 Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 ... 2148
    > which(y==3)        # ID番号が3である行番号
    [1]  910 1105 7367
    > df[which(y==3),]   # チェック
          a  b  c
    910  aa aa eb
    1105 aa aa eb
    7367 aa aa eb
    > system.time({y <- paste(df[,1],df[,2],df[,3],sep="`");
                   y[y %in% setdiff(y, y[duplicated(y)])] <- NA;
                   y <- as.factor(y); levels(y) <- 1:nlevels(y)})
       ユーザ  システム      経過 
        0.020     0.004     0.026 

データハンドリング 縦横変換

蔵正? (2007-08-07 (火) 09:56:23)

Aのデータを読み込み、Bのようにデータ変換したいのですが、どうすれば良いでしょうか、ご教示いただければ幸いです。

A----------------------------------------------------------------
PATNO VISIT SUBSCALE QUESTION RESP
1103 Baseline Total Score 99 0.571428571
1103 Month 1 Total Score 99 1.214285714
1103 Month 6 Total Score 99 1.4
1103 Month 9 Total Score 99 0.547619048
1104 Baseline Total Score 99 2.048780488
1104 Month 1 Total Score 99 1.095238095
1104 Month 3 Total Score 99 1
1104 Month 6 Total Score 99 0.904761905
1104 Month 9 Total Score 99 0.880952381
1105 Baseline Total Score 99 1.571428571
1105 Month 1 Total Score 99 1.666666667
1105 Month 3 Total Score 99 0.976190476
1105 Month 6 Total Score 99 0.928571429
1105 Month 9 Total Score 99 1.023809524
  :

B----------------------------------------------------------------
PATNO Baseline Month 1 Month 3 Month 6 Month 9
1103 0.571428571 1.214285714 1.142857143 1.4 0.547619048
1104 2.048780488 1.095238095 1 0.904761905 0.880952381
1105 1.571428571 1.666666667 0.976190476 0.928571429 1.023809524
  :

  • reshapeでピボットテーブルを見よ。その他reshape関数などいろいろある。 -- okinawa 2007-08-07 (火) 10:24:35
  • その前にこのデータはかなり不規則ですから、整形しないといろいろ困難がありそうです。(1) Baseline はひとつの文字列なのに、"Month 1" は二文字。(2) Month 3 がある場合と無い場合がある(これはまちがいでは無く、そのままですか)。 -- 2007-08-07 (火) 11:29:05
  • 整形するくらいなら,readLines で読んで splitchar や grep 使って,行列を埋めていくという戦略を取る方がよいかも。
    ただし,データの規則を十分に把握していないとプログラムできない。
    Baseline と Month 1,Month 3,Month 6,Month 9 で欠測値もあるとか,ひょっとして,Month 12やMonth 15やそれ以上の月のものがあったり,Month 5 なんてものもあり得るのかどうかとか。同一患者のデータは一箇所にかたまっているのかとか。 -- 2007-08-07 (火) 11:36:49
  • xtabsではだめでしょうか?たしかrjpwikiのどこかでピポットテーブルのような使い方をみた覚えがあります。例題や欠測の場合の処理などはヘルプをご覧ください。 -- ぼう? 2007-08-07 (火) 17:59:24
  • なるほどね。テキストエディタなんかを使ってデータを全文置換して以下のようにまとめる -- 2007-08-07 (火) 18:23:03
    PATNO VISIT RESP
    1103 Baseline 0.571428571
    1103 Month1 1.214285714
    1103 Month6 1.4
    1103 Month9 0.547619048
    1104 Baseline 2.048780488
    1104 Month1 1.095238095
    1104 Month3 1
    1104 Month6 0.904761905
    1104 Month9 0.880952381
    1105 Baseline 1.571428571
    1105 Month1 1.666666667
    1105 Month3 0.976190476
    1105 Month6 0.928571429
    1105 Month9 1.023809524
    そのあとデータフレーム(df とでも)に読んで,
    > ( ans <- xtabs(RESP~PATNO+VISIT, df) )
          VISIT
    PATNO   Baseline    Month1    Month3    Month6    Month9
      1103 0.5714286 1.2142857 0.0000000 1.4000000 0.5476190
      1104 2.0487805 1.0952381 1.0000000 0.9047619 0.8809524
      1105 1.5714286 1.6666667 0.9761905 0.9285714 1.0238095
    PATNO=1103 の Month3 は実際は欠測値なんだけど 0 が入っている(0というのが有効な測定値でないなら,0をNAに置き換えるだけでよいね。
    > ans[ans==0] <- NA
    > class(ans) <- "matrix"
    > ans
          VISIT
    PATNO   Baseline   Month1    Month3    Month6    Month9
      1103 0.5714286 1.214286        NA 1.4000000 0.5476190
      1104 2.0487805 1.095238 1.0000000 0.9047619 0.8809524
      1105 1.5714286 1.666667 0.9761905 0.9285714 1.0238095
    attr(,"call")
    xtabs(formula = RESP ~ PATNO + VISIT)
  • xtabでももちろんできます。私も以前はxtabを使ってましたが、面倒くさくなって今はreshapeパッケージのcast()melt()を使ってます。 -- okinawa 2007-08-07 (火) 18:44:10

バックスラッシュを置換できますか?

ぼう? (2007-08-06 (月) 16:51:00)

winのファイルパス名が入ったテキストをいじる際、バックスラッシュを置き換えたいと思います。以下のようにすると「??」を置換することはできますが、エスケープシーケンスを減らしても、「?」1文字だけを置換することができません。私の正規表現の理解が不十分なのかもしれませんが、どなたかご教示いただければ幸いです。環境はWinXp? R2.51です。

> x <- "aa?bb??cc"
> gsub("????", "/", x)
[1] "aa?bb/cc"    #本当は"aa/bb??cc"としたい
  • ¥b は二文字ではなく一文字なんだから,¥とはマッチしないでしょうね。 -- 2007-08-06 (月) 17:31:10
  • Rでは¥t(Tab)¥n(Return)のような使い方をしますので、上記の文字列の場合は [a][a][¥b][b][¥¥][c][c] と認識されます。ですから、[¥bb]=>[/bb]に変換する工夫をすればいいのではないでしょうか。 -- okinawa 2007-08-06 (月) 17:58:16
  • なるほど、次の文字に対してかかってしまい?bなどとなってしまうところを見落としていました。うまくフォルダ名のフォーマットにあわせて変換できるよう工夫してみます。ありがとうございます。 -- ぼう? 2007-08-06 (月) 19:21:46

中央値でソートして箱ヒゲ図を表示するには

sh? (2007-08-04 (土) 13:43:20)

グループの中央値でソートしての箱ヒゲ図を表示するにはどうすればよろしいでしょうか?

boxplot(count ~ spray, data = InsectSprays, col = "lightgray") # 1
boxplot(decrease ~ treatment, data = OrchardSprays, log = "y", col = "bisque") # 2
# 1を2のように表示したい
  • count ~ spray というモデル式は因子 spray の水準の順に箱型図を書きます.help(reorder) を見てください.因子 spary の水準を count の中央値の順に並べ替えればよいのでは。-- 2007-08-04 (土) 15:44:02
    > attach(InsectSprays)  
    > spray2 <- reorder(spray, count, median) 
    > spray
     [1] A A A A A A A A A A A A B B B B B B B B B B B B C C C C C C C C C C C C D D
    [39] D D D D D D D D D D E E E E E E E E E E E E F F F F F F F F F F F F
    Levels: A B C D E F       # 最初の水準の並び
    > spray2
     [1] A A A A A A A A A A A A B B B B B B B B B B B B C C C C C C C C C C C C D D
    [39] D D D D D D D D D D E E E E E E E E E E E E F F F F F F F F F F F F
    attr(,"scores")
       A    B    C    D    E    F 
    14.0 16.5  1.5  5.0  3.0 15.0 
    Levels: C E D A F B             # count の中央値によって水準を並べ替えた
    > boxplot(count ~ spray2, col = "lightgray")
    # 一行で書けば
    > boxplot(count ~ reorder(spray, count, median), data = InsectSprays, col = "lightgray") 
    # 逆順にソートするなら、count を -count にするトリックで
    > boxplot(count ~ reorder(spray, -count, median), data = InsectSprays, col = "lightgray")
  • ありがとうございました。大変勉強になりました -- sh? 2007-08-05 (日) 00:57:03
  • 新たなオブジェクトを作るのは,ときとして考え物。at 引数を活用すべきではないでしょうか。 -- 2007-08-06 (月) 11:58:19
    attach(InsectSprays) # 当然と思って省略したが。。。
    boxplot(count ~ spray, data = InsectSprays, col = "lightgray", at=rank(by(count, spray, median)))
  • 正しくは(あらかじめ attach(InsectSprays?) を実行していない限り) -- 2007-08-06 (月) 13:35:03
    boxplot(count ~ spray, data = InsectSprays, col = "lightgray",  
            at=rank(by(InsectSprays$count, InsectSprays$spray, median)))
  • attach(InsectSprays?) したなら data 引数は冗長。 -- 2007-08-06 (月) 13:45:48

降順のソート関数

*? (2007-07-31 (火) 09:10:01)

降順のソートの関数はあるのでしょうか?
探しているんですが昇順しか見当たりません。
探し不足かもしれませんが、よければ回答お願いします。

  • help(sort)と入力しましたか? decreasingでは不十分なケース? -- なかま 2007-07-31 (火) 09:42:01
  • rev(sort(x)) でも可。 -- 2007-07-31 (火) 13:01:27
  • 回答ありがとうございます。実行することが出来、とても助かりました。 helpも確認不足でした。これから気をつけます。ありがとうございました。 -- *? 2007-07-31 (火) 18:24:58

関数の出力を表示させない

RYUI? (2007-07-30 (月) 16:37:13)

Rのバージョンは2.5.1です.
非計量多次元尺度構成法を行っているのですが,そこで次元数PとストレスSの折れ線グラフを描く関数を作ってみました.

library(MASS) # isoMDS関数を使用するため
#類似度データ
ruijiData <- data.frame(
  c(0,1,5,3),
  c(1,0,2,4),
  c(5,2,0,5),
  c(3,4,5,0)
)
#次元数とストレスのプロット関数
plot.mds.stress <- function(data) {
  stress <- NULL
  for(i in 1:c(nrow(data)-1)) {
    mds <- invisible(isoMDS(dist(data),k=i))
    ps <- c(i,mds$stress)
    stress <- rbind(stress, ps)
  }
  plot(stress,main="ストレスのプロット",xlab="次元数 P",ylab="ストレス S",type="o")
}

plot.mds.stress(ruijiData) を実行すると,図と以下の(isoMDS関数による)出力が得られます.

initial  value 13.343941 
iter   5 value 0.313496
final  value 0.006694 
converged
initial  value 0.000000 
final  value 0.000000 
converged
initial  value 0.000000 
final  value 0.000000 
converged

ここで,上記の出力が表示されないようにしたいのですが,invisible関数ではできませんでした.
どのようにしたら,出力を消すことができるのでしょうか? よろしくお願いいたします.

  • 出力させないと言うだけなら,for の前に sink("適当なファイル名") ,for の後に sink() とすれば? -- 2007-07-30 (月) 17:27:14
  • 上の出力はおそらく optim() 関数の経過出力で、本当に収束したかどうかチェックするためのものです。消してしまって本当に大丈夫ですか。どうしても消したければ例えば(少なくとも以下の例ではうまくいく)。最後に sink() を実行しないとコンソールへの全ての出力が消えてしまいますからご注意。 -- 2007-07-30 (月) 17:24:17
    > res <- optim(sq, distance, genseq, method="SANN",
    +    control = list(maxit=6000, temp=2000, trace=TRUE))
    sann objective function values
    initial       value 29625.000000
    iter     1000 value 14815.000000
    iter     2000 value 14815.000000
    iter     3000 value 14815.000000
    iter     4000 value 13541.000000
    iter     5000 value 13482.000000
    iter     5999 value 13323.000000
    final         value 13323.000000
    sann stopped after 5999 iterations
    >  sink("all.Rout")                 # ファイル名は適当
                          # 以下全てのコンソール出力はファイルにリダイレクトされる     
    > res <- optim(sq, distance, genseq, method="SANN",
    +    control = list(maxit=6000, temp=2000, trace=TRUE))
                       # 何も出力されない   
    > sink()               # sink 関数によるリダイレクト終了 
    > res
     $par
     [1]  1 19 16  8 13 15  2  9 12 14  5 18  4  3 11  7 20 10  6 17 21  1
    
    $value
    [1] 12919
    
    $counts
    function gradient 
        6000       NA 
    
    $convergence
    [1] 0
    
    $message
    NULL
  • 別解:capture.output を使う。簡単にするためにシュガーコートする。 -- 2007-07-30 (月) 17:36:57
    plot.mds.stress2 <- function(data) junk <- capture.output(plot.mds.stress(data)) # シュガーコート関数
    plot.mds.stress2(ruijiData)
  • やはり収束の善し悪しを見ないで結果(プロット)だけ得る、というやりかたは危険だと思う。 -- 2007-07-30 (月) 18:05:48
  • 収束の善し悪しをプログラムで判定して(上のプログラムで言えば junk に入っている文字列を解析して),収束しなかったときにだけ(短い)warning を出すようにすればよい -- 2007-07-30 (月) 18:11:49
  • help を良く読みましょう!trace という引数があるでしょう? -- 2007-07-30 (月) 18:21:30
    mds <- isoMDS(dist(data),k=i, trace=FALSE)
    ヘルプに曰く,
    Side Effects
    If trace is true, the initial stress and the current stress are printed out
    every 5 iterations.
  • 皆様,お返事が遅れ大変申し訳ございません.皆様の回答により,満足のいく結果が得られました.また,収束を見ないで分析する危険性など,ご指摘抱きまして,誠にありがとうございました! -- RYUI? 2007-08-01 (水) 19:19:25

3次元グラフで座標値に文字列ラベルをつける方法

niko? (2007-07-27 (金) 15:41:34)

scatterplot3dを使って3次元グラフを作成しています.
グラフ中に表示される座標値に文字列ラベルをつける方法を探しています.
2次元グラフの場合は

plot(ai.i$points, type="n")
text(ai.i$points, labels=rownames(ai.i$points))

としていたのですが,3次元グラフ(scatterplot3d)での方法を見つけられずにいます.
どなたかご教示いただけないでしょうか.

  • scatterplot3d が返す,xyz.convert function which converts coordinates from 3D (x, y, z) to 2D-projection (x, y) of scatterplot3d. Useful to plot objects into existing plot. にヒントがあるのでは? -- 2007-07-27 (金) 16:08:42
    d <- structure(c(-4L, 1L, -8L, 4L, -3L, 10L, 15L, 11L, -21L, -5L, 
                     -3L, -11L, -9L, 11L, 4L, 0L, 6L, 14L, -18L, 11L,
                     0L, -8L, -9L, 12L, -5L, -3L, 16L, 5L, -18L, 10L),
                     .Dim = c(10L, 3L))
    obj <- scatterplot3d(d)
    text(obj$xyz.convert(d), labels=LETTERS[1:10], pos=4)
    scatterplot3d-example.png
  • ありがとございます!xyz.convert() を使うんですね.希望通りのグラフが作成できました.ご教示ありがとうございます. -- niko? 2007-07-30 (月) 16:52:30

barplotの凡例の順序について

K? (2007-07-24 (火) 16:15:10)

以下のようなスクリプトで、下図のヒストグラムを作りました。

  x <-    c( 4, 6, 8,10,12,14,16,18,20,22,24,26,28,30,32)    #直径階
remain<-  c( 8, 3, 7, 5, 4, 3, 5, 3, 4, 0, 2, 2, 3, 1, 1)      #スギ間伐前本数
fell  <-  c( 5, 2, 5, 3, 3, 2, 4, 2, 3, 0, 1, 1, 2, 1, 1)      #スギ間伐後本数
dead  <-  c( 3, 1, 1, 1, 2, 3, 1, 1, 1, 0, 1, 1, 2, 1, 0)      #スギ(Cryptomeria japonica D.Don)枯損木
DBH_dist <- function(x, remain, fell, dead, title, xlabel, ylabel){
  par(ps=20)
  par(mar=c(5,5,3,3))
  dummy <- x * 0
  sum   <-  matrix(rbind(remain, fell, dummy),nrow=3, length(x))
  yrange <- c(-max(dead)-10, max(remain + fell) + 10)
  par(new=F)
  barplot(-dead, names.arg=x,  ylim=yrange, col = "gray", axis.lty = 1, xlab=xlabel, ylab=ylabel, main=title)
  par(new=T)
  barplot(sum, ylim = yrange,  col = c("black","white","gray"),
           axes=F, ann=F, legend=c( "残存木","伐採木", "枯損木" ))
}
title = "テスト用のダミーデータ"
xlabel="DBH"; ylabel="直径階別本数分布(本/plot)"
DBH_dist(x, remain, fell, dead, title, xlabel, ylabel)


Kuni_070724.png

できれば凡例を上から順に「伐採木」「残存木」「枯損木」のように配置したいのですが、
legendのヘルプを見てもよく分かりません。どのようにすればよいのでしょうか。よろしくお願いします。

  • どこに問題点が存在するのか分からないのですが,2つめのbarplotを以下のように2つに分ければいいということではないのでしょうか?
    barplot 関数の legend 引数は,そのグラフに書かれた以外の凡例項目を描くには適していないと思います。 -- 2007-07-24 (火) 16:32:45
      barplot(sum, ylim = yrange, col = c("black","white"), axes=F, ann=F)
      legend("topright", c("伐採木", "残存木", "枯損木" ), fill=c("white", "black", "gray"))
  • 早速のコメント、ありがとうございます。barplotのlegend引数とlegend関数の区別ができていませんでした。勉強になりました。 -- K? 2007-07-24 (火) 17:16:03

モード/最頻値

とくめい? (2007-07-23 (月) 17:43:33)

Rはsummary関数でmedianとmeanを出してくれますが、modeは出してくれません。
modeを計算する関数はあるのでしょうか?

  • google などで検索しましたか? http://cran.r-project.org/doc/contrib/Lemon-kickstart/kr_dstat.html の他にも色々出てきますが,探すより,自分でプログラムする方が速いかも。計算対象変数が連続変数なのか離散変数なのかとかの問題もあるし。 -- 2007-07-23 (月) 17:59:37
  • 離散値データなら例えば次のようにする。 -- 2007-07-23 (月) 22:02:06
    > x                  # 数値データの場合
      [1]  1  1  2  5  6  3  5 10 10 10  1  2  3  3  5  5  7  8  9 10  0  1  3  2  7
     [26]  8  6  6  9 10  1  3  2  7  5  4  7  5  8 10  4  1  3  1  6  5  8  7  9 10
     [51]  1  4  3  5  4  5  9  6  8 10  1  1  2  0  4  5  5  8  8 10  1  3  4  5  8
     [76]  7  8  8 10 10  1  1  1  5  5  7  6  9 10 10  1  3  3  5  5  7  6  8  7 10
    > (y <- table(x))           # 度数分布表
    x
     0  1  2  3  4  5  6  7  8  9 10 
     2 15  5 10  6 16  7  9 11  5 14 
    > names(y)
     [1] "0"  "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10"
    > names(y)[y == max(y)] # モード
    [1] "5"
    > as.numeric(names(y)[y == max(y)]) # 数値が欲しければ
    [1] 5
    
    > x                     # 文字列データの場合
      [1] "x" "w" "w" "j" "c" "n" "x" "r" "r" "e" "n" "i" "v" "s" "c" "g" "h" "u"
     [19] "v" "v" "p" "q" "s" "h" "v" "o" "n" "w" "d" "o" "l" "k" "c" "s" "g" "e"
     [37] "j" "f" "l" "m" "o" "g" "n" "e" "x" "o" "d" "j" "s" "t" "p" "v" "s" "h"
     [55] "c" "v" "i" "q" "n" "x" "m" "r" "f" "r" "d" "h" "u" "k" "n" "m" "v" "x"
     [73] "o" "c" "j" "d" "j" "n" "p" "f" "m" "e" "g" "y" "c" "l" "y" "c" "k" "h"
     [91] "v" "g" "s" "n" "b" "s" "h" "x" "z" "y"
    > ( y <- table(x) )          # 度数分布表
    x
    b c d e f g h i j k l m n o p q r s t u v w x y z 
    1 7 4 4 3 5 6 2 5 3 3 4 8 5 3 2 4 7 1 2 8 3 6 3 1 
    > names(y)
     [1] "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q" "r" "s" "t"
    [20] "u" "v" "w" "x" "y" "z"
    > names(y)[y == max(y)]        # モード
    [1] "n" "v"
  • 定義できない(前述の定義では不定になる)場合があることをわすれないように。 -- 2007-07-23 (月) 23:01:33

Rmapのインストール

内田 直樹? (2007-07-20 (金) 17:25:44)

Rmap_1.1.0zipをインストールしましたがRのversion2.4.1では”有効なパッケージではありません”となってしまいます。ためしにversion1.9にインストールするとlibraryコマンドでロードできます。version2.4.1が日本語表示なのでこれを使いたいのですが、現在のversionにインストールする方法はあるでしょうか?

  • RTools を使ってソースよりビルドかな?maptools や sp パッケージなんかどうですか?谷村先生の解説文があるし。 -- 2007-07-20 (金) 17:37:31
  • ありがとうございました。谷村先生の論文を読んでみます。 -- 内田? 2007-07-21 (土) 12:35:33
  • 初級Q&A アーカイブ(6)「Rmapのインストールについて」「MacOS XでのRmapの使用」 -- 2007-07-21 (土) 18:13:10

barplotによるヒストグラムの作成方法について

K? (2007-07-19 (木) 17:07:04)

Rで以下のようなヒストグラムを作りたいのですが、思うようにいかずに困っています。

Kuni_070719.png

barplotを使ってやろうとしているのですが、helpをみても見当が付きません。
barの幅を0.4にして、間伐前と後のbarの位置をそれぞれ元の位置から
±0.25にしてやればできそうな気がするのですが、barplotのパラメータに
width = rep(0.4, length(x)) と書き加えても何も変わらないし、offsetは縦方向へのシフトだし...。
こんなイレギュラーなグラフは無理かなぁ、と思う一方、Rならきっとできるに違いないとは思うのですが。
どなたかお知恵を貸して下さい。よろしくお願いします。
以下が現状のスクリプトです。

# barplot03.R
#直径階別本数分布を、間伐前と間伐後を上向きに並べて、枯損木を下向きにして表示したいのだが...
x <-  c(4,6,8,10,12,14,16,18,20,22,24,26,28,30,32)         #直径階
N0 <- c(3,0,1, 0, 2, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0)    #スギ枯損木
L0 <- c(2,0,3, 2, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1)   #広葉樹枯損木
N1 <- c(8,0,7, 5, 4, 0, 0, 3, 0, 1, 0, 2, 0, 0, 1)         #スギ間伐前本数
L1 <- c(2,0,4, 2, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1)         #広葉樹間伐前本数
N2 <- c(6,0,5, 3, 3, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0)         #スギ間伐後本数
L2 <- c(1,0,3, 2, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1)         #広葉樹間伐後本数
sum0 <- N0+L0; sum1 <- N1+L1; sum2 <- N2+L2
before <- matrix(rbind(N1, L1),2, length(x))
after  <- matrix(rbind(N2, L2),2, length(x))
koson <- matrix(rbind(N0, L0),2, length(x))
par(ps=20)
par(mar=c(5,5,3,3))
xlabel="DBH"; ylabel="直径階別本数分布(本/plot)"
title = "テスト"
yrange = c(-max(sum0)-1, max(sum1)+1)
par(new=F)
barplot(before, width = rep(0.5, length(x)), space = NULL, names.arg=x, ylim = yrange,
       axis.lty=1, xlab="", ylab="", main="")
par(new=T)
barplot(-koson, names.arg=x,  ylim=yrange,
      axes = F,  xlab=xlabel, ylab=ylabel, main=title)
  • こんな感じで。 -- takahashi? 2007-07-19 (木) 17:53:03
    upper<-matrix(rbind(before,after),nrow=2)
    yrange<-c(-max(sum0)-1, max(sum1)+1)
    xlim<-c(0,7*ncol(upper)/2)
    barplot(upper, width = 2, space = c(1,0.5),  xlim=xlim ,ylim = yrange,       axis.lty=1, ann=FALSE)
    par(new=T)
    barplot(-koson, names.arg=x,  width=5, space=2/5, xlim=xlim,ylim=yrange, axes = F)
  • 期待していたグラフが描けました。どうもありがとうございます。 -- K? 2007-07-20 (金) 14:07:00

データフレームの列の(範囲)指定

PC初心者? (2007-07-19 (木) 04:00:45)

初心者です.(WindowsXP,R-2.5.1)
以下のサンプルデータ("sitsumon")があるとします.

q1 <- c("F","M","M")
q2_1 <- c(60,70,75)
q2_2 <- c(70,75,70)
q2_3 <- c(90,95,80)
q3 <- c(1,0,0)
q4 <- c(3,4,5)
sitsumon <- data.frame(q1=q1, q2_1=q2_1, q2_2=q2_2, q2_3=q2_3, q3=q3, q4=q4)

(*ここで,実際には,より多くのデータ数(行数)と質問の項目数(列数)があります)

実際に分析する際は,"q2_1"と"q2_2","q2_3","q3"のみを使用するので,扱いやすいように

UFO <- sitsumon[,2:5]

とすればよいのですが,実際には,列数が多い(100くらい)ので,"q2_1","q2_2","q2_3","q3"が何列目か非常にわかりにくいです.
また,抜き出す列数も実際には,もっと多いですので,

UFO <- sitsumon[ ,c("q2_1","q2_2","q2_3","q3")]

というのも避けたいです.

したがって,

UFO <- sitsumon[,"q2_1":"q3"]

のような感じで,2列目("q2_1")〜5列目("q3")を抜き出したいのですが,できませんでした.

そこで,皆様にご教授お願いしたいのが,「列名」からそれが何列目になるのかを数値で返す関数や,もしくは他のよい解決方法はございますでしょうか?
(できればOS(プラットフォーム)などに依存しないようなやり方がよいです)
愚直な質問ではございますが,ご教授くだされば幸いです.

  • help(subset) してみる(特にその example を見る)。 -- 2007-07-19 (木) 07:16:22
    > x <- data.frame(a=1:3, b=1:3, c=1:3)
    > colnames(x)
    [1] "a" "b" "c"
    > which(colnames(x)=="a")
    [1] 1
    > which(colnames(x)=="b")
    [1] 2
  • subset(sitsumon,se=q2_1:q3) -- takahashi? 2007-07-19 (木) 09:07:42
  • 普通はsubsetです。ラベル名の場所を知りたいならcharmatch("q4",names(sitsumon))も可 -- okinawa 2007-07-19 (木) 10:02:40
  • grep("q4", colnames(sitsumon))とか正規表現を使える関数はいろいろありますよ -- 2007-07-19 (木) 15:01:50
  • お返事が遅れてしまい,申し訳ございません.丁寧に回答していただきありがとうございました.皆様の回答により求めていた操作ができるようになりました. -- PC初心者? 2007-07-19 (木) 20:08:38

ベイズ−BUGSでのM−Hアルゴリズムの実行

azu? (2007-07-18 (水) 23:31:51)

ベイズでM-Hアルゴリズムを実行したい場合、
BUGSで動かすことは可能なのでしょうか?
BUGSの名前から考えるとできなさそうなのですが・・・・
M-Hアルゴリズムを動かすソフトウェアがあれば、
教えてください。
お願いします。

  • 何を質問しているのか、訳が分かりません。もっと、具体的にしゃべればいかが。それと、どうせ、Win環境なんでしょうが、それに関する情報も提供しておく方がいいでしょう。要するに、あなたの質問に答えてやっても良いかなと思わせるとか、将来自分もそのような問題に遭遇するかも知れないなとか、同じWin仲間なんだから救ってやろうかなとか、そういう「何か」がないと、回答がもらいにくいとは思いませんか?( WIn仲間なら、そんなことにはお構いなく、仲間意識で回答してやろうと思う人が大多数なんだろうか(^_^;)) -- 2007-07-18 (水) 23:38:49
  • 申し訳ないです。私の環境はWinですが、どんな環境をお使いでも、「M−Hアルゴリズムはこのソフトウェアでまわしてる」って方がいれば教えて頂きたいです。質問について説明が足りなかったので、詳細を書きます。ベイズ統計学を用いるためにMCMC法の実装を行っています。その中でBUGS(私の場合、windowsなので、R2WinBUGSなのですが)というパッケージソフトを利用しています。これでギブスサンプラーの実装は問題なく可能なのですが、M−Hアルゴリズムを用いる場合も同じBUGSで実装できるのか、もしくは別のソフトウェアが必要なのか、もしくはパッケージソフトに頼らず自分でプログラムを書くべきなのか、わからないままでいます。もし、MCMC法の実装を行われている方がいれば教えてください。お願いします。 -- azu? 2007-07-19 (木) 03:01:36
  • google でキーワード "Bayes BUGS M-H" で検索してみれば。 -- 2007-07-19 (木) 07:22:53
  • すいぶん遅いレスですが、Metropolis-Hastings法はたしかwinbugsでは使えないから自分で書いたと人づてに聞いたことがあります。BUGSはgibbsサンプラだけみたいです。英語でググレば出てくるかも? -- たけ? 2007-09-20 (木) 13:28:38

spgrass6のインストールが出来ません_ubuntu7.04

しまだ? (2007-07-18 (水) 20:40:33)

どなたか教えてください。
linux超初心者です。grassとRを使いたくて、ubuntu-7.04をセットアップしました。
ubuntu:7.04
grassのバージョン:6.2.1
Rのバージョン:2.4.1
R version 2.4.1 (2006-12-18)~ i486-pc-linux-gnu
ところが、install.packages("spgrass6")を実行すると、
下記のようなメッセージがでてspgrass6とrgdalをインストールすることができません。
spやmaptoolsは、r-base-devを入れたら、無事インストールできました。
どなたかアドバイスをいただけないでしょうか。

install.packages("spgrass6") 中で警告がありました:
   argument 'lib' is missing: using /usr/local/lib/R/site-library
--- Please select a CRAN mirror for use in this session ---
Loading Tcl/Tk interface ... done
* Installing *source* package 'spgrass6' ...
** R
** inst
** preparing package for lazy loading
Loading required package: sp
Loading required package: maptools
Loading required package: foreign
Error: package 'rgdal' required by 'spgrass6' could not be found
Execution halted
ERROR: lazy loading failed for package 'spgrass6'
** Removing '/usr/local/lib/R/site-library/spgrass6'
The downloaded packages are in
/tmp/RtmpQLX6fo/downloaded_packages
Warning message:
installation of package 'spgrass6' had non-zero exit status 
   in: install.packages ("spgrass6")
  • spgrass6 をインストールする前に rgdal をインストールせよ、といっているように読めますが。 -- 2007-07-18 (水) 21:21:50
  • ありがとうございます。早速やってみましたが、同様のメッセージがでました。いくつかネットもみて真似てみたのですが... -- しまだ? 2007-07-18 (水) 23:19:18
    install.packages("rgdal_0.5-13.tar.gz","/usr/lib/R/lib",repos=NULL)
    WARNING: invalid package 'rgdal_0.5-13.tar.gz'
    ERROR: no packages specified
    Warning message:
    installation of package 'rgdal_0.5-13.tar.gz' had non-zero exit status in:
      install.packages("rgdal_0.5-13.tar.gz", "/usr/lib/R/lib", repos = NULL)
  • rgdal パッケージは R から GDAL ライブラリにアクセスするためのツールのようです。GDAL, PROJ.4 ライブラリは既にインストール済みですか?その上で tar.gz ファイルをダウンロードして端末から sudo R CMD INSTALL rgdal_0.5-13.tar.gz したらうまくいきませんか。 -- 2007-07-19 (木) 10:36:48
    rgdal: Bindings for the Geospatial Data Abstraction Library
    
    Provides bindings to Frank Warmerdam's Geospatial Data Abstraction Library
    (GDAL) (>= 1.3.1) and access to projection/transformation operations from 
    the PROJ.4 library. The GDAL and PROJ.4 libraries are external to the package,
    and, when installing the package from source, must be correctly installed first. 
    Both GDAL raster and OGR vector map data can be imported into R, and GDAL raster
    data and OGR vector data exported. Use is made of classes defined in the sp package.
  • ありがとうございました.助言のおかげでバラバラに得ていた情報が,ようやく整理されました.きちんと解説を読むべきだったと反省しております.今は出先ですので,戻り次第動作確認したいと思います.本当にありがとうございました. -- しまだ? 2007-07-19 (木) 13:45:54
  • 確認しました。ありがとうございました。さきほど無事インストールできました。
    長年(本当に長かった)の懸案がようやく解決しました。
    Rjpwikiに投稿して本当によかったです。-- しまだ? 2007-07-19 (木) 22:08:17
    # apt-get install gdal-bin
    # apt-get install libgdal1-dev
    # apt-get install proj
    # (rgdal_0.5-13.tar.gzをダウンロード)
    # sudo R CMD INSTALL rgdal_0.5-13.tar.gz
    # ---Rを起動---
    # install.packages("spgrass6")
  • よかったですね。お役に立てて何よりです(自分で確かめもしない無責任アドバイスでしたが)。 -- 2007-07-20 (金) 00:00:37

boxplotの目盛りと外れ値

sh? (2007-07-15 (日) 11:10:40)

boxplot(c(0.1, 1:100, 1000), log="y") で、boxplotのy軸目盛りを「1e-01,1e+00,1e+01,1e+02,1e+03」ではなく「0.1,1,10,100,1000」と表示する方法を教えてください。また、○印(外れ値または異常値。ここでは1000)を決定するアルゴリズムを教えてください(標準偏差の何倍以上?)

  • 以下のようになさったらいかが。
    外れ値は、第1四分位数-1.5×四分範囲」より小さいデータと,「第3四分位数+1.5×四分範囲」より大きいデータ -- 2007-07-15 (日) 11:56:44
    boxplot(c(0.1, 1:100, 1000), log="y", axes=FALSE, frame.plot=TRUE)
    y <- 10^(-1:3) # 適切に設定(データから決定できるようにするとよい)
    axis(2, at=y, label=as.character(y))
  • こちらも検索不足。初心者コースの履歴にどちらも質問&回答されていますよ。質問をするより、検索した方が早い場合もある。 -- 2007-07-15 (日) 12:17:33
  • 軸の目盛りをどうするかということで、axis の使い方に関するものは何回も出ている。しかし、boxplot での axis については、他の場合と同じようにはいかない。確かに boxplot の説明を読めばよいのだが、他の関数と同じ感じで yaxt="n" とかやってもだめだし、初心者には難しいのかも知れない。 -- 2007-07-15 (日) 12:25:15
  • 解決しました。ありがとうございました。 -- sh? 2007-07-15 (日) 20:31:42
  • 単に、「解決しました」の報告だけではなく、色々付いたコメントに対するコメントもあればいいのになあと思いますが。 -- 2007-07-15 (日) 21:59:04
  • コメントです。『外れ値は、第1四分位数-1.5×四分範囲」より小さいデータと,「第3四分位数+1.5×四分範囲」より大きいデータ』という情報の出典(参考文献またはURL)を明記していただけると幸いに存じます。『# 適切に設定(データから決定できるようにするとよい)』の具体例を示していただけると幸いに存じます。関数plotやbarplotのx軸、y軸の目盛りの値を(関数histのcounts、midsのように)参照することは可能なのでしょうか。『初心者コースの履歴にどちらも質問&回答されていますよ。』の具体的な場所(タイトルまたはURL)または検索ワードを明記していただけると幸いに存じます。 -- sh? 2007-07-16 (月) 02:31:25
  • 何でも教えてもらえばいいというものではない.箱ヒゲ図の基本的な部分については統計の教科書に出ているし(中澤先生の「Rによる統計解析の基礎」には出ていたような?),自分の探したい情報の検索ワードはご自分が最もよくわかるはずでしょう.ここは自助努力が肝要,という場所です. --yuta? 2007-07-16 (月) 13:26:05
  • 追加の情報.青木繁伸先生のページが箱ヒゲ図について参考になります.私も詳しい経緯ははじめて知りました.勉強になった. -- yuta? 2007-07-16 (月) 13:35:55
  • 「縦軸に変数値をとって,第1四分位を下に,第3四分位を上にした箱を書き,中央値の位置にも線を引いて,さらに第1四分位と第3四分位の差(四分位範囲)を 1.5 倍した線分をヒゲとして第1四分位の下と第3四分位の上に伸ばし,ヒゲの先より外れた値を外れ値として○をプロットした図である」http://phi.ypu.jp/statlib/stat.pdf -- 2007-07-16 (月) 14:53:45
  • http://aoki2.si.gunma-u.ac.jp/lecture/Dosuu/box-and-whisker.html をどうぞ -- 2007-07-16 (月) 18:49:00
  • > barplotのx軸、y軸の目盛りの値を(関数histのcounts、midsのように)参照することは可能なのでしょうか
    online help の value の項を見れば分かるように(あるいは,簡単に str(boxplot(rnorm(100))) をしてみれば分かるように),boxplot が返すオブジェクトの中には,ご希望の情報はないですね。 -- 2007-07-17 (火) 17:53:23
  • > 『# 適切に設定(データから決定できるようにするとよい)』の具体例
    まず,何をもって適切とするかを決めないといけないがそれは,あなたが決めないとね。私が決めて良いとして,データがこんな範囲の値になるなら(つまり,データを常用対数で表したときの範囲が数個の整数を含む範囲であるような場合なら),以下のようなものも。 -- &new{2007-07-17 (火) 18:14:05};
    x <- c(0.1, 1:100, 1000)
    boxplot(x, log="y", axes=FALSE, frame.plot=TRUE)
    minx  <- floor(log10(min(x))) # 与えられた x だと,minx = -1
    maxx <- ceiling(log10(max(x))) # 与えられた x だと,maxx = 3
    y <- 10^(minx:maxx) # y は,c(0.1, 1, 10, 100, 1000)
    axis(2, at=y, label=as.character(y))
  • 更に一般的と言うことになると,ちょっと時間が必要だったなぁ -- 2007-07-17 (火) 18:59:52
    log.ax <- function(x) { # データの範囲を押さえる対数目盛り
    	minx  <- log10(min(x))
    	maxx <- log10(max(x))
    	y <- 10^(log10(1:9)+rep(floor(minx):ceiling(maxx), each=9))
    	l <- y >= 10^minx & y <= 10^maxx
    	minl <- min(which(l == TRUE))
    	if (minl > 1) l[minl-1] <- TRUE
    	return(y[l])
    }
    ############ 使用例
    x <- c(0.31, 2, 5, 40, 100, 1000, 5100)
    y <- log.ax(x)
    boxplot(x, log="y", axes=FALSE, frame.plot=TRUE)
    axis(2, at=y, label=as.character(y))
  • ああ・・・私自身も検索不足でしたね。「適切に設定」の具体例は私にも参考になりました。 -- yuta? 2007-07-17 (火) 19:59:26

バイナリデータの扱い方について

YURI? (2007-07-14 (土) 14:43:14)

初心者です。バイナリデータの扱い方が良く分りません。
バイナリファイルから読み込んだデータを4バイト単位にまとめるためにmatrixでncol=4を指定してみました。

> f <- file('/tmp/a.out', open='rb')
> a <- readBin(f, what='raw', n=40)
> a
[1] 7f 45 4c 46 01 01 01 00 00 00 00 00 00 00 00 00 02 00 03 00 01 00 00 00 90
[26] 82 04 08 34 00 00 00 50 1d 00 00 00 00 00 00
> matrix(a, ncol=4)
> dump("a", file="/tmp/a.R")
> system('cat /tmp/a.R')
"a" <-
c(7f, 45, 4c, 46, 01, 01, 01, 00, 00, 00, 00, 00, 00, 00, 00, 
00, 02, 00, 03, 00, 01, 00, 00, 00, 90, 82, 04, 08, 34, 00, 00, 
00, 50, 1d, 00, 00, 00, 00, 00, 00)
> source('/tmp/a.R')
Error in parse(file, n = -1, NULL, "?") : syntax error on line 2
> dput(a)
c(7f, 45, 4c, 46, 01, 01, 01, 00, 00, 00, 00, 00, 00, 00, 00, 
00, 02, 00, 03, 00, 01, 00, 00, 00, 90, 82, 04, 08, 34, 00, 00, 
00, 50, 1d, 00, 00, 00, 00, 00, 00)

matrixの結果は何も表示されません。what="raw"で読み込んだデータには matrixは使えない(1要素のサイズが分らないので?)ということでしょうか。
また、dumpしたものがsourceで読み込めないのも、こういうものなのでしょうか?
結局、"rawデータ"というものは、どのように使えるデータなのかを知りたいです。

  • すみません、整形に失敗しましたが、編集の仕方が分からず直せません。 -- YURI? 2007-07-14 (土) 14:57:47
  • どのように使いたいのかわからないが,readBin に対応するのは writeBin では?オンラインヘルプは参照しましたか?利用例も載っているようですよ。 -- 2007-07-14 (土) 17:06:51
  • もう少し勉強してみます。~私がmatrixで結果が表示されなかった原因が(多分)分かりました。sessionInfo()の結果を載せるべきでしたが、私のR環境はversion 2.1.0(と2.2.0)で、このバージョンでは動作しない模様です。試しに2.4.0をインストールしてみたら期待どおりの結果が表示されました。
    ありがとうございました。 -- YURI? 2007-07-14 (土) 21:36:09
  • 説明不足だったかもしれませんが、readBinを使ったのは、簡単にバイナリデータのサンプルを作りたかったためですので、writeBinはあまり関係ないと思います。質問の発端はmatrixが期待どおり動作しなかったことなのですが、それについては上記のとおり解決しました。~ -- 2007-07-14 (土) 22:10:01
  • 「バージョンが古かった」ってことですか。解決したなら良かったですね。バージョンは最新にしておくべきです。で、その結果としては、問題は全部解決したわけですね。これからは、全部自分でできるようになったわけですね。よかったよかった。 -- 2007-07-14 (土) 22:13:50
  • dump して source で読もうとしたわけですわな。確かに、普通のオブジェクトならそれでよかったんだけど、raw はそれではだめだったと。。。 -- 2007-07-14 (土) 22:19:41
  • すみません、まだ投稿法に慣れず、コメント入力の途中で投稿してしまいました(~の後にリターンを入力したら投稿されてしまった)。次の疑問は、dputが生成する外部表現と同じ入力(例えば、c(7f)など)がシンタクスエラーになってしまうことなのですが、もう少し自力で勉強してからにしたいと思います。長い道のりになりそうですが。。。 -- YURI? 2007-07-14 (土) 22:20:41

散布図を右90度回転

ssk? (2007-07-13 (金) 07:02:13)

plotで右90度回転させた散布図を書くためにはどうすればよろしいでしょうか?
boxplotやbarplotのようにオプションを指定して簡単にできないものでしょうか?

boxplot(1:10, horizontal=T)
barplot(1:10, horiz=T)
plot(1:10, horiz=T) # Warning messages:

環境:R version 2.5.0 (2007-04-23) i386-apple-darwin8.9.1
参考:「2変量散布図とヒストグラム図」「クラスター図を横向きで描画したい」

  • 形式的に回転するとy軸(元のx軸)の目盛が降順になってしまいますが?plot(x,y) と plot(y,x) では駄目ですか。最終的にepsファイル等の画像ファイルに落とすなら、画像ファイルを原稿に取り込む際に回転できます。 -- 2007-07-13 (金) 08:22:22
  • 「boxplotやbarplotのようにオプションを指定して簡単にできないものでしょうか?」って
    存在しないオプションを指定しても,エラーになるだけ
    どういう意図があって,どんな散布図が欲しいのか,上の人のコメントにあるが,plot(y, x) ではだめなのか,どこがだめなのか。
    plot(y, x) なんてことには考えがおよばなかったと言うことなのか。
    速やかなる返答を求む。 -- 2007-07-13 (金) 17:04:43
  • >速やかなる返答を求む。...たかだか8時間程度で,そんなことを言わなくても...別に,みんなが,ハードにこのページをチェックしているわけないですから。 -- 2007-07-14 (土) 00:04:50
  • 普段は一週間に一回見るか見ないかのような人でも、自分が何か質問したときには、回答がもらえるかどうか、それこそ1時間ごとにでも見るんじゃないの?そんなに急がない質問もあるとは思うけど。 -- 2007-07-14 (土) 09:33:02
  • こんな感じでしょうか?
    x <- 1:10
    y <- 20:11
    
    plot(x,y) #普通のplot()
    
    plot(y,x) #xとyが入れ替わっている->これだと、y軸の場所が違って嫌なのでしょうか?
    
    plot(y,x,xaxt="n", xlab="") #y軸の表示をOFFにして
    axis(side=3, at=y) #y軸を上に持ってきて
    mtext(text="y", outer=F, line=3) #ylabもmtextで追記する
    mainの引数もmtext()で自由になります -- 2007-07-14 (土) 11:21:04
  • at=y は,へんじゃない? -- 2007-07-14 (土) 12:11:14
  • > たかだか8時間程度で,そんなことを言わなくても
    さらに時間が経ちました。答えは必要じゃないということか、あるいは、既に回答は得られたということか。 -- 2007-07-14 (土) 21:28:13
  • アドバイスありがとうございます。返事が遅くなり申し訳ございません。離散データ(x, y)の両対数散布図、xの分布(対数目盛)図、yの分布(対数目盛)を横にした図、の3点を描きたいと考えています。右90度回転というのは私の勘違いでした。-- ssk? 2007-07-15 (日) 03:20:04
    def.par <;- par(no.readonly = TRUE) # save default, for resetting...
    x <- round(rlnorm(100)+1)
    y <- round(rlnorm(100)+1)
    xhist <- hist(x, breaks=seq(0.5,max(x)+0.5,1), plot=FALSE)
    yhist <- hist(y, breaks=seq(0.5,max(x)+0.5,1), plot=FALSE)
    nf <- layout(matrix(c(2,0,1,3),2,2,byrow=TRUE), c(3,1), c(1,3), TRUE)
    layout.show(nf)
    par(mar=c(3,3,1,1))
    plot(x, y, log="xy")
    par(mar=c(0,3,1,1))
    plot(xhist$mids, xhist$counts,type="h",log="x",xaxt="n")
    par(mar=c(3,0,1,1))
    plot(yhist$mids, yhist$counts,type="h",log="x",xaxt="n")#を横にしたい
    par(def.par)#- reset to default
  • 以下のようにすればいかが -- 2007-07-15 (日) 11:19:32
    # 最初の plot
    plot(x, y, log="xy", ylim=exp(range(log(y))))
    # 三番目の plot を以下の2行に
    plot(yhist$counts,yhist$mids, type="n",log="y",yaxt="n", ylim=exp(range(log(y))))
    junk <- sapply(1:length(yhist$counts), function(i) lines(c(0,yhist$counts[i]),
        rep(yhist$mids[i], 2)))
    plot999.png
  • だとしたら、sskさんの検索不足、このページの下に同じ「2変量散布図とヒストグラム」がありますね。ちゃんと調べてますか? -- 2007-07-15 (日) 12:13:46
  • 全ての可能性を検索するのも大変だし、下の方のは barplot を使っている訳で水平にするオプションを持つが、今回は plot なんだから。また、離散変数でしかも log を取るから、barplot は不適切なんでしょうね。結局今回の質問の問題点は、自分のやりたいことをちゃんと説明しなかったこと(追加説明するのにも時間が掛かったこと)でしょう。 -- 2007-07-15 (日) 12:20:42
  • 三番目のplotの凡例を描画するにはどうすればよろしいでしょうか? -- ssk? 2007-07-17 (火) 09:11:03
    legend(max(x), max(y), legend="(1)", box.lty=0) # 一番目のplotの凡例 - 表示される
    legend(max(xhist$mids), max(xhist$counts), legend="(2)", box.lty=0) # 二番目のplotの凡例 - 表示される
    legend(max(yhist$mids), max(yhist$counts), legend="(3)", box.lty=0) # 三番目のplotの凡例 - 表示されない
    legend(max(yhist$counts), max(yhist$mids), legend="(3)", box.lty=0) # 三番目のplotの凡例 - 表示されない
  • legend(locator(1), legend="(3)", box.lty=0) # 関数 locator() により対話的に凡例の位置を決定する -- 2007-07-17 (火) 12:36:50
  • 各 legend は,plot の直後に書くんですね。legend の第1引数と第2引数が凡例のどこの座標(右上,右下,左上,左下のどれ)を指しているかを確認すべきでしょうね。
    第一引数に "topright" などと書いて描いてみるのも一法かと。
    そのほかに指定できるのは,"bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right", "center"
    この例もまた,「ヘルプを良く読むべし」ですね。 -- 2007-07-17 (火) 15:54:54

メタファイルからビットマップファイルへの変更って?

チョメスケ? (2007-07-12 (木) 15:47:06)

最近,Rの勉強を始めたばかりの素人です。質問内容はごく簡単な内容だと思いますが,お分かりになる方がいらっしゃれば,よろしくご指導のほどお願いいたします。

質問内容:MCMCpackのMCMCregressを実行しました。独立変数が4つあるのですが,最初にグラフ化される(全2ページ中の1ページ目)のが切片を含め,x1,x2の3つの時系列プロットとposteriorが表示されます。ここで,次のページへ移行する前に,1ページ目の画像を保存で,どうしてもメタファイル形式しか選ぶことができません(2ページ目は保存形式が多数あるので,問題なありません)。WSwordに張り付けることを考えると,ビットマップ形式の方がストレスがなくて良いと考えているのですが,どのようにすれば1ページ目の画像をビットマップ形式での保存ができるのでしょうか?
いろいろと試行錯誤してみたのですが,できませんでした。もし,できないとしてもメタファイルで保存した画像をビットマップファイルで保存する方法を知っていらっしゃいましたら,ぜひご教唆いただけないでしょうか?
よろしくお願いいたします。

  • メニューから保存するのではなくて,直接ファイルへ書き出すようにプログラムすればよろしいのではないでしょうか。その方が手間が掛からないし便利だと思いますが。Winじゃできないのでしょうか(そんなことないですよね) -- 2007-07-12 (木) 16:16:08
  • ?bmpしてみればいいと思いますが、例えばこんな感じのを応用してはどうでしょう?
    x <- matrix(runif(10),nrow=5)
    png()
    image(x)
    dev.off()
    winなら普通にワーキングディレクトリに出ますよ。 あえて、メタファイルの画像ならばirfanviewやVIXなどフリーのソフトで見たり変換できるものがいくらでもあると思います。数枚ならパワーポイントとかOooとかでも行けますし。-- たけ? 2007-07-12 (木) 17:06:04
  • あなたが言っている「MCMCregressを実行しました」はMCMCpackのマニュアルに載っているMCMCregressのExampleのコードのことでしょうか?私の環境(WinXPsp2,R2.4.1)では問題なくメタファイルでもビットマップでもコピーできましたよ。あなたの実行環境が問題なのかもしれませんね。投稿するときには必ず実行環境と実行したコードを明示してください。 -- okinawa 2007-07-12 (木) 17:23:14
  • ご返答いただいた皆様,ありがとうございます。 -- チョメスケ? 2007-07-16 (月) 21:59:44
  • 実際に,実行したコードならびに実行環境をご説明させていただきます。まず最初にlibraryからMCMCpackを読み込ませます。その次に,Exampleのデータを自分のデータに編集しなおし,次のプログラムを実行いたしました(非常に見づらく申し訳ございません
    line<-list(df$y,df$x1,df$x2,df$x3,df$x4)
    posterior<-MCMCregress(y~x1+x2+x3+x4,data=df2,verbose=2000)
    plot(posterior)
    raftery.diag(posterior)
    summary(posterior)
    なお,Rは2.5です。 -- チョメスケ? 2007-07-16 (月) 22:14:02
  • 解決いたしました!実はとっても単純な話でございました。最初の図が表示されている状態でショートカットキーのCTRL+Cでビットマップファイルで保存ができました!(ちなみにCTRL+Wでメタファイル保存)本当に,いろいろとお世話になりました。これで,ストレスなく,wordに張り付けた状態でタイプできます。重ね重ね,いろいろな御助言ありがとうございました。 -- チョメスケ? 2007-07-16 (月) 22:25:33

barplotの横軸目盛り

K? (2007-07-12 (木) 06:37:14)

 barplotについての質問です。

x <-    c(4,6,8,10,12,14,16,18,20,22,24,26,28,30,32)
Grp1 <- c(8,0,7, 5, 4, 0, 0, 3, 0, 1, 0, 2, 0, 0, 1)
Grp2 <- c(2,0,4, 2, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1)
Grp3 <- c(3,0,2, 0, 2, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0)
y <- matrix(rbind(Grp1, Grp2, Grp3),3,length(x))
barplot(y, xlab="DBH")

として、下図を作図しました。このbarplotの横軸にベクトルxの値を目盛りとして付けたいのですが、どのようにしたらよいのでしょうか。よろしくお願いします。

Kuni_barplot.png
  • 回答は axis(1,at=1:length(x)*1.2-0.5, labels=x)
    ただし、1.2 と 0.5 と言う数値については、!is.matrix(height) || !beside のときの値である(詳しくは、barplot.default のソースを読むこと)-- 2007-07-12 (木) 07:35:15
    barplot2.png
  • ax<-barplot(y, xlab="DBH");axis(1,ax,x) -- takahashi? 2007-07-12 (木) 08:02:19
  • takahashi様の回答の方がシンプルに同じ結果が得られますね。助かりました。どうもありがとうございました。 -- K? 2007-07-12 (木) 08:37:18
  • barplotの引数namesを使って、一行で書くこともできます。> barplot(y,names=x) -- MMK? 2007-07-12 (木) 13:06:06
  • 原質問文を読むと確かに「横軸にベクトルxの値を目盛りとして付けたい」とあり,軸を描きたいとは書いていなかったなあ。ヘルプは良く読むべし。引数名はなるべく省略しない方がよいかも。 -- 2007-07-12 (木) 14:01:09
  • で,ヘルプを隅から隅まで読むと,barplot(y, names.arg=x, axis.lty=1) だけで軸も書けることがわかる。 -- 2007-07-12 (木) 14:10:20
  • 本当だ、helpを読むとnames.argやaxis.ltyのところにちゃんと書いてありますね。もっとヘルプを良く読むように努力します。軸も書きたかったので、axis.lty=1が必要ですね。ax<-barplot(y, xlab="DBH");axis(1,ax,x) よりも、barplot(y, names.arg=x, axis.lty=1)の方が(恐らく)引数名を省略していないので何をやっているのかがよく分かりました。ありがとうございました。 -- K? 2007-07-12 (木) 14:57:56

時系列データ処理の高速化の仕方について

? (2007-07-11 (水) 22:33:30)

アーカイブの(6)にも似たような内容がありまして恐縮ですが、どうしても自分の力で解決できないため、お力添えいただければと思います。

データ:2分間隔で測定した1年分のデータ、50箇所分。
目的 :週毎に平均した1日の変動の値を算出する。
(0:00の週の平均値、0:02の週の平均値…23:58の週の平均値の720のデータからなるセットを52週分算出する。)

#結果収納用データフレームを作成
y <- x[1,] #xは測定値が入っている262800行50列の行列

#52週間分平均データ計算
for (i in 1:52){
 for (j in 1:720){
  x <- rbind(y,
       apply(x[seq((i-1)*720*7+j, 720*7*i, by=720),], 2, mean))
 }
}

これだと、1回計算するのにかなりの時間(6〜7時間)がかかってしまうなので、なんとか高速化できないかと思っています。
行列を作ってみようともしたのですが、結局、平均するデータの行を指定する行列を繰り返し発生させるプログラムを書くことしか思いつきませんでした。

作業を効率的に行う方法をご教示いただけませんでしょうか。
よろしくお願い申し上げます。

  • 高速化が目的ならデータフレームやapply関数は使わず、ベクトル、行列、配列そして関連高速処理関数を使う、のが鉄則です。あなたの書いている「:00の週の平均値、0:02の週の平均値…23:58の週の平均値の720のデータからなるセットを52週分」という意味がよくわかりませんが。 -- 2007-07-11 (水) 23:11:46
  • 今あなたは、データを2次元配列で扱っているようだが、データは時刻×曜日×週×箇所の4次元データ。apply は何次元の配列についても計算できる。apply の第二引数は 1,2のようなスカラーばかりではない。時刻・曜日・箇所ごとに平均値を出すのだから、このように考えるのがぴったり。apply も、内部的には for を使うのと変わりないという話もあるが。
    テストデータとして、12ヶ月それぞれにおいて1週間のみ毎時0分の2カ所のデータを考える。つまり、24×7×12×2これを週の平均を取って24×12×2個の平均値を求めるには以下のようになる。 -- 2007-07-11 (水) 22:54:25
    x <- sample(100, 24*7*12*2, replace=TRUE)
    dim(x) <- c(24, 7, 12, 2)
    system.time(ans <- apply(x, c(1,3,4), mean))
    dim(ans)
    私のポンコツコンピュータでは以下のようになった
    > x <- sample(100, 24*7*12*2, replace=TRUE)
    > dim(x) <- c(24, 7, 12, 2)
    > system.time(ans <- apply(x, c(1,3,4), mean))
       ユーザ  システム      経過 
        0.057     0.002     0.071 
    > dim(ans)
    [1] 24 12  2
    あなたのデータ規模だと、どれくらいの時間が掛かるだろうか??
    単に比例計算してみて驚いたね。時間が掛からないのに驚いたのだ。
    これなら、実際に計算してみることができる。
    > x <- sample(100, 720*7*52*50, replace=TRUE)
    > dim(x) <- c(720,7,52,50)
    > system.time(ans <- apply(x, c(1,3,4), mean))
       ユーザ  システム      経過 
      318.294     6.466   367.592 
    > dim(ans)
    [1] 720  52  50
    5分くらいで終わった。たぶん、これであっていると思うんだが。
    検算してみよう。
    > ans[2,5,12] # 目的とする平均値の入っている要素を取り出し
    [1] 59.42857
    > mean(x[2,,5,12]) # 該当する7個のデータの平均値
    [1] 59.42857
    > ans[321,15,17] # 目的とする平均値の入っている要素を取り出し
    [1] 47.14286
    > mean(x[321,,15,17]) # 該当する7個のデータの平均値
    [1] 47.14286
    いいみたいだ。
  • ご教示、どうもありがとうございます。
    今まで配列というものを使ったことがなかったので、そのようなやり方があることを理解していませんでした。データを配列に組み替えて実行して、結果をご報告いたします。-- -- ? 2007-07-11 (水) 23:36:42
  • apply は使わざるを得ないと思うんだが。。。。 -- 2007-07-11 (水) 23:38:24
  • 第1の場所について、第1週すなわち1月1日0時0分~23時58分、2日0時0分~23時58分、…、7日について同じく、第2週について同じく、…、第52週まで、次いで第2の場所について、…、第50の場所についてという順でしょうから、単に dim 属性を変えるだけでよいと思いますよ。 -- 2007-07-11 (水) 23:44:38
    > a <- matrix(1:24, 12)
    > a
          [,1] [,2]
     [1,]    1   13
     [2,]    2   14
     [3,]    3   15
     [4,]    4   16
     [5,]    5   17
     [6,]    6   18
     [7,]    7   19
     [8,]    8   20
     [9,]    9   21
    [10,]   10   22
    [11,]   11   23
    [12,]   12   24
    > dim(a) <- c(3,4,2) # 3×4×2 配列にする
    > a
    , , 1
    
         [,1] [,2] [,3] [,4]
    [1,]    1    4    7   10
    [2,]    2    5    8   11
    [3,]    3    6    9   12
    
    , , 2
    
         [,1] [,2] [,3] [,4]
    [1,]   13   16   19   22
    [2,]   14   17   20   23
    [3,]   15   18   21   24
  • あなたのコードのまず大きな問題は rbind 関数の多用です。こうしたオブジェクトのサイズを後から変える操作は時間をくいます。さらに rbind 関数の結果はリストになってしまうので更にハンディを課すことになります。以下の例が参考になるでしょうか。まず結果を入れるのに必要なサイズの配列 res[52, 720, ] を用意しておいてそれに代入すべきでしょう。 -- 2007-07-11 (水) 23:53:35
    > x <- y <- 1:100
    > system.time(for (i in 1:9999) y <- rbind(y, x))
       ユーザ  システム      経過 
      106.122    38.282   146.565 
    > z <- matrix(0, 10000, 100)
    > system.time(for (i in 1:10000) z[i,] <- x)
       ユーザ  システム      経過 
        0.168     0.036     0.205 
    > str(z)
     num [1:10000, 1:100] 1 1 1 1 1 1 1 1 1 1 ...
    > str(y)
     int [1:10000, 1:100] 1 1 1 1 1 1 1 1 1 1 ...
     - attr(*, "dimnames")=List of 2
      ..$ : chr [1:10000] "y" "x" "x" "x" ...
      ..$ : NULL
  • 先ほど、ようやくもろもろ理解できました。 あまりに早くて(2分足らずで出来ました)涙が出そうです。どうもありがとうございました。配列恐るべしですね、勉強します。 -- ? 2007-07-12 (木) 03:16:44
  • 参考まで -- 2007-07-13 (金) 01:03:34
    > x <- sample(100, 720*7*52*50, replace=TRUE)
    > dim(x) <- c(720,7,52,50)
    > system.time(ans <- apply(x, c(1,3,4), mean))
       ユーザ  システム      経過 
      101.654     2.796   112.885 
    > ans2 <- array(0, c(720,52,50))
    > system.time(for(i in 1:720) for(j in 1:52) ans2[i,j,] <- colSums(x[i,,j,]) )
       ユーザ  システム      経過 
        3.636     0.624    23.429 
    > ans2 <- ans2/7  # 総和を平均に変換
    > identical(ans, ans2)
    [1] TRUE
    ついでに気づいた変なこと。
    > 720*52*7*50
    [1] 13104000
    > 262800*50
    [1] 13140000
    つまり一年は7の倍数(52週)ではないので、次元操作だけではまずい(なによりも、適正に配列化されているかどうか確認が面倒)。次のようなやりかたではどうでしょう。
    > x <- matrix(1:(262800*50), 262800, 50) # テスト行列
    > xx <- array(0, c(52,720,7,50))     # 作業用の配列の初期化
    ## 適切に配列化(これも少し時間がかかる)
    s <- seq(1, length=7, by=720)      # 作業用の長さ7の数列(毎回配列を生成しない)
    > for (i in 1:52) for (j in 1:720) xx[i,j,,] <- x[(i-1)*720*7+(j-1)+s,]
    > xxx <- array(0, c(52,720,50))  # 最終結果を入れる配列の初期化
    ## 高速な行和関数を使う
    > for (i in 1:52) for (j in 1:720) xxx[i,j,] <- colSums(xx[i,j,,])
    > xxx <- xxx/7     # 総和を平均に一括変換(一々割らない)して最終結果を得る
  • たぶん自動測定器のデータ等だろうから,週平均というのは必ずしも1月1日から始まらないし12月21日で終わらない(要するに52週分はあるだろうということ) -- 2007-07-20 (金) 11:15:11

データフレームからルックアップ

でみあん? (2007-07-11 (水) 10:01:43)

いつまで経っても,編集されないので,代理で編集しておく
以下のような二つのデータフレームがある場合、

df1<- data.frame(cbind(c("a","a","a","b","b"),c("abc","def","ghi","jkl","mno")))
df2<- data.frame(cbind(c("a","b"),c(3,1)))

df2の各レコードに対して、それぞれ一列目の値をキーに一定の基準(対応する最後のの行、など)に基づき、df1のレコードを選択してMergeする良い方法はありませんでしょうか?現在は、for ループを用いていますが実行にかかる時間を短縮する良い方法があればと悩んでいます。

最終的には、

  X1 X2  X3
1  a  3 ghi
2  b  1 mno

のようなデータフレームを作成したいのです。

  • もう少しわかりやすい(?)方法なら。いずれにしても簡潔に書くことはできても時間を短縮することはあまり期待出来ないでしょう。 -- 2007-07-11 (水) 13:48:35
    > df1 <- data.frame(x=c("a","a","a","b","b"), y=c("abc","def","ghi","jkl","mno"))
    > df2 <- data.frame(xx=c("a","b"), yy=c(3,1))
    > (temp <- df2$yy[df1$x])
    [1] 3 3 3 1 1
    > cbind(df1, data.frame(z=temp))
      x   y z
    1 a abc 3
    2 a def 3
    3 a ghi 3
    4 b jkl 1
    5 b mno 1
  • 「対応する最後の行」など という条件が付いていますねえ -- 2007-07-11 (水) 14:20:54

文字列を分解してベクトルに変換する方法

学生? (2007-07-10 (火) 18:01:03)

はじめまして。文字列を分割してベクトルに変換したいのですが、strsplit関数で

str <- "1.1 2.3 3.4"
str.v <- strsplit(  "1.1 2.3 3.4", " " )

などを自分では試してみたいのですが、うまくいきません。どのようにすればよいでしょうか?
よろしくお願いします。

  • どのような結果を期待していて,どのような結果しか得られないのかを書いてくださいな。上の str.v はちゃんとできているように思われますが?
    数値ベクトルにしたいというのでしょうか?だったら,str.v <- as.numeric(unlist(strsplit( "1.1 2.3 3.4", " " ))) とでも
    str.v <- scan(textConnection(str)) とでも,お好きに -- 2007-07-10 (火) 18:07:10
  • 属性を知ることから始めた方が近道と思いますよ -> R-Tips -- 2007-07-10 (火) 20:53:55
  • 解説ありがとうございます。やりたいことは、こうして作成したベクトルstr.vの数字の合計や平均値などを算出したいのです。ご指摘のとうり以下のようなプログラムを続けてかいています。
    str.v <-s.numeric(unlist(strsplit( "1.1 2.3 3.4", " " ))) 
    ave(str.v)
    [1] 2.266667 2.266667 2.266667
    という結果になります。これを
    [1] 2.266667
    のようになると思うのですが、どうして三つ平均値が表示されてしまうのでしょうか?
    • 学生? 2007-07-11 (水) 15:08:25
  • ?ave -- takahashi? 2007-07-11 (水) 15:32:58
  • すみません。aveじゃないですねmean(str.v)ですね。ちゃんと実行できました。 -- 学生? 2007-07-11 (水) 15:43:48

平方ユークリッド距離

KO? (2007-07-10 (火) 16:27:11)

非類似度として平方ユークリッド距離を使用したいのですが、helpで調べましたが、ユークリッド距離はありますが、平方ユークリッド距離については記載がありません。
計算可能でしょうか?

  • 平方ユークリッド距離は,ユークリッド距離を二乗しただけ(というか,平方根をとってユークリッド距離にする前の状態のものが平方ユークリッド距離) -- 2007-07-10 (火) 16:29:55
  • x <- matrix(rnorm(100), nrow=5); dist(x)^2; dist(x, method = "euclidean")^2 -- 2007-07-12 (木) 17:19:46
  • できました。ありがとうございます。 -- KO? 2007-07-13 (金) 16:34:36

アポストロフィーを含む文字列の読み込み

K? (2007-07-10 (火) 09:42:04)

read.table("ImportFile?", header=F, sep=";")のようにファイルを読み込む際、以下の例のようにデータの一部にアポストロフィーが入っていると、次にアポストロフィーが出てくる場所までがひとかたまりに読み込まれてしまいます。

2094;0;MACY'S;3;5;7
2095;0;MACY'S;3;4;8

read.tableのオプションなどでアポストロフィーを無視するような設定はできますか?

  • quote="" -- takahashi? 2007-07-10 (火) 09:47:34
  • ? read.table ってやらないんですか? -- 2007-07-10 (火) 10:30:09
  • すみません。help(read.table)やって読んだつもりだったんですが、すっかり見落としていました。おかげさまで見つけることができました。失礼しました。 -- K? 2007-07-10 (火) 10:32:30

Rcmd SHLIBでDLLファイルを作成する際にR2.5だとうまくいかない?

たけ? (2007-07-05 (木) 17:45:38)

仮に解決済みなので恐縮ですが、Rの他言語利用の、windowsの場合 にしたがって、dllファイルを作成してみたところ、Rの2.5.1もしくは2.5.0では以下のエラーになり・・・

D:?test?r>Rcmd SHLIB test.c
windres -I C:/PROGRA~1/R/R-25~1.1/include  -i test_res.rc -o test_res.o
c:?MinGW?bin?windres.exe: unknown format type `C:/PROGRA~1/R/R-25~1.1/include'
c:?MinGW?bin?windres.exe: supported formats: rc res coff
make: *** [test_res.o] Error 1

Rの2.3.1では以下のようにうまく行きました。

D:?test?r>Rcmd SHLIB test.c
latex: not found
gcc  -shared -s  -o test.dll test.def test.o -LC:/PROGRA~1/R/R-23~1.1/bin   -lR

これはバグでしょうか?私の環境はWinXPです。ちなみにtest.cとは上記サイトの例のmysum.c と同じです。ここで適切か分かりませんが、結構悩んだのでご報告まで。

  • 2.5.1 for Winですが、自家製Cプログラムは今まで通りコンパイルできてます。コメントを求める前に提示する情報が不足してますね。『仮に解決済み』って何ですか? -- 2007-07-05 (木) 20:30:04
  • Rの2.3.1 ならうまくいくので,バージョンダウンしたと言うことでは? -- 2007-07-05 (木) 21:44:08
  • はい、ダウンバージョンで解決済みという意味です。分かりにくい表現ですみません。あまり中身のことが分かっていないので必要な情報が分かっていないですが、他の方で問題がないようでしたら固有の問題かも知れません。失礼しました。 -- たけ? 2007-07-06 (金) 10:39:31
  • こちらの環境(2.5.1 for Win)ではmysum.cがコンパイルできました。RではなくMSYS側の問題なのでは? -- 2007-07-06 (金) 12:27:19

任意の相関係数を持つポアソン乱数列を作成したいのですが

taro? (2007-07-02 (月) 00:51:12)

与えられたlambdaを持つポアソン乱数列に対して、任意の相関係数を持つポアソン乱数列を生成したいと思っております。
こちらのサイトの「初級Q&A アーカイブ(1) - 層別の相関値」や、下記のサイトを拝見して、任意の相関係数を持つ乱数列を生成する関数 gendat2() と、ポアソン分布への適合度検定関数 poissondist() を知りました。

そこで、gendat2() 関数の最初のデータ作成部分を変更して、任意の相関係数を持つポアソン乱数列を生成するようにしてみて、その結果が果たしてポアソン分布に従っているのかを poissondist() 関数で確認しようとしましたが、うまくいきませんでした。

質問させていただきたいのは、
1. あるlambdaを持つポアソン乱数列に対して、任意の相関係数を持つポアソン乱数列の作成方法と、
2. poissondist() 関数の出力結果で、lambda が設定した値とかなりずれてしまっている原因について
の2点です。

下記に実行結果を記載させていただきましたが、gendat2() 関数を真似して実行した変換では、元の2つのポアソン乱数列を、異なるlambdaを持つ2つのポアソン乱数列へとまとめて変換してしまうこととなってしまいました。(確かに、変換後の2つのポアソン乱数列の相関係数は指定した通りなのですが、あるlambdaを持つポアソン乱数列に対して、指定した相関係数を持つポアソン乱数列を作成する、という本来の目的とは違ってしまいました)

また、下記以外にも、いろいろなlambdaを指定してrpois()関数で生成したポアソン乱数列を、poissondist()関数で検定してみたのですが、出力結果のlambdaが、設定したlambdaとかなり異なってしまっていて、原因が分からず困っております。

どうか、上記2点に関しまして、原因と対策をご教授いただけませんでしょうか。よろしくお願いします。
長文失礼しました。

nc <- 1000 # 乱数列の個数
r  <- 0.8  # 相関係数
z <- matrix(rpois(2*nc, lambda=1/10), ncol=2)
# 最初に生成したポアソン乱数列の確認
poissondist(hist(z[,1])$counts)
   $result
          n   lambda
   1000.000    0.419
   $table
       x   d            p            e
   c-0 0 899 6.577042e-01 6.577042e+02
   c-1 1   0 2.755781e-01 2.755781e+02
   c-2 2   0 5.773360e-02 5.773360e+01
   c-3 3   0 8.063460e-03 8.063460e+00
   c-4 4  98 8.446474e-04 8.446474e-01
   c-5 5   0 7.078145e-05 7.078145e-02
   c-6 6   0 4.942905e-06 4.942905e-03
   c-7 7   0 2.958682e-07 2.958682e-04
   c-8 8   0 1.549610e-08 1.549610e-05
   c-9 9   3 7.528501e-10 7.528501e-07
   $result2
    chi sq.     d.f.  P value
   2993.188    2.000    0.000
   Warning messages:
   1:  長いオブジェクトの長さが
            短いオブジェクトの長さの倍数になっていません in: x - e
   2:  長いオブジェクトの長さが
            短いオブジェクトの長さの倍数になっていません in: (x - e)^2/e
res <- eigen(r2 <- cor(z))
coeff <-  solve(r2) %*% (sqrt(matrix(res$values, 2, 2, byrow=T))*res$vectors)
z <- t((t(z)-apply(z, 2, mean))/sqrt(apply(z, 2, var)*(nc-1)/nc)) %*% coeff
x <- z %*% chol(matrix(c(1, r, r, 1), ncol=2))
# 変換して生成したポアソン乱数列の確認
poissondist(hist(x[,1])$counts)
   $result
          n   lambda
   1000.000    0.817
   $table
         x   d            p            e
   c-0   0 820 4.417549e-01 4.417549e+02
   c-1   1   0 3.609138e-01 3.609138e+02
   c-2   2   0 1.474333e-01 1.474333e+02
   c-3   3   0 4.015100e-02 4.015100e+01
   c-4   4 161 8.200841e-03 8.200841e+00
   c-5   5   0 1.340017e-03 1.340017e+00
   c-6   6   0 1.824657e-04 1.824657e-01
   c-7   7   0 2.129635e-05 2.129635e-02
   c-8   8   2 2.174890e-06 2.174890e-03
   c-9   9  16 1.974317e-07 1.974317e-04
   c-10 10   0 1.613017e-08 1.613017e-05
   c-11 11   0 1.198032e-09 1.198032e-06
   c-12 12   0 8.156599e-11 8.156599e-08
   c-13 13   1 5.442424e-12 5.442424e-09
   $result2
    chi sq.     d.f.  P value
   5244.693    4.000    0.000
   Warning messages:
   1:  長いオブジェクトの長さが
            短いオブジェクトの長さの倍数になっていません in: x - e
   2:  長いオブジェクトの長さが
            短いオブジェクトの長さの倍数になっていません in: (x - e)^2/e

下記の関数を参考にさせていただきました。

## http://cat.zero.ad.jp/~zak52549/R/haebara/chap04.txt
## 任意の相関係数を作成する関数を定義 gendat2(標本サイズ,相関係数)
gendat2 <- function(nc, r)
{
    # 仮のデータ行列を作る。この時点では変数間の相関は近似的に0。
    z <- matrix(rnorm(2*nc), ncol=2)
    # 主成分分析を行い,主成分得点を求める。この時点で変数間の相関は完全に0。
    res <- eigen(r2 <- cor(z))
    coeff <-  solve(r2) %*% (sqrt(matrix(res$values, 2, 2, byrow=T))*res$vectors)
    z <- t((t(z)-apply(z, 2, mean))/sqrt(apply(z, 2, var)*(nc-1)/nc)) %*% coeff
    # コレスキー分解の結果をもとに,指定された相関係数行列 を持つように主成分得点を変換。
    z %*% chol(matrix(c(1, r, r, 1), ncol=2))
}

## 青木先生のサイト http://aoki2.si.gunma-u.ac.jp/R/poissondist.html
## ポアソン分布への適合度の検定
poissondist <- function(d)                              # 度数分布ベクトル
{
        k <- length(d)                                  # 階級数
        n <- sum(d)                                     # データ数
        k1 <- k-1
        x <- 0:k1
        lambda <- sum(d*x)/sum(d)                       # 平均値(=λ)
        result <- c("n"=n, "lambda"=lambda)
        p <- exp(-lambda)*lambda^x/factorial(x)         # 確率
        p[k] <- 1-sum(p[-k])                            # 最後の確率は合計が 1 になるように
        e <- n*p                                        # 期待値
        table <- data.frame(x, d, p, e)                 # 結果をデータフレームにまとめる
        rownames(table) <- paste("c-", x, sep="")
        while (e[1] < 1) {                              # 1 未満のカテゴリーの併合
                d[2] <- d[2]+d[1]
                e[2] <- e[2]+e[1]
                x <- d[-1]
                e <- e[-1]
                k <- k-1
        }
        while (e[k] < 1) {                              # 1 未満のカテゴリーの併合
                d[k-1] <- d[k-1]+d[k]
                e[k-1] <- e[k-1]+e[k]
                x <- d[-k]
                e <- e[-k]
                k <- k-1
        }
        chisq <- sum((x-e)^2/e)                         # カイ二乗統計量
        k <- k-2                                        # 自由度
        p <- pchisq(chisq, k, lower.tail=FALSE)         # P 値
        result2 <- c("chi sq."=chisq, "d.f."=k, "P value"=p)
        return(list(result=result, table=table, result2=result2))
}
  • gendat2 関数の返す値はそもそも非負整数値になるのですか? 二変量ポアソン分布については、例えばこれを見る -- 2007-07-02 (月) 09:09:30
  • hist(rpois(n, lambda=bar))$count というのはまずいなあ。 -- 2007-07-02 (月) 10:32:56
  • hist() 関数で無理矢理整数化していたのですか。最初のコメントを注釈すれば、X0, X1, X2 をそれぞれ平均 c0, c1, c2 の互いに独立なポアソン乱数とすると、 Y1=X0+X1, Y2=X0+X2 はそれぞれ平均 c0+c1, c0+c2 のポアソン乱数になり、その共分散は c0 になります。したがって Y1, Y2 の相関係数は
    c0/sqrt((c0+c1)*(c0+c2))
    になります。c0 を適当に選べば、少なくとも (0, 1) 区間内の値を持つ相関係数は表現できます。負の相関はこの方法では、実現出来ません(そもそも可能かどうかも知りません)。 -- 2007-07-02 (月) 20:00:04
  • 二変量ポアソン分布の資料を教えていただいて、大変ありがとうございました。任意の相関係数を持つポアソン乱数列は、同時確率密度から計算できました。 -- taro? 2007-07-09 (月) 04:29:04

添付ファイル: fileKuni_070724.png 1522件 [詳細] filehistgram_k2.png 1801件 [詳細] filelegend_K1.png 534件 [詳細] fileplot999.png 1638件 [詳細] fileKuni_070719.png 1676件 [詳細] fileKuni_barplot.png 1504件 [詳細] filescatterplot3d-example.png 1642件 [詳細] fileggplot99.png 1627件 [詳細] filecif.gif 1907件 [詳細] fileiris78975839789.png 1741件 [詳細] filegroupedData2.jpg 1636件 [詳細] filebarplot2.png 1566件 [詳細] filenlme-plot.png 1588件 [詳細] fileggplot98.png 1582件 [詳細] fileconsolegamen.png 484件 [詳細] filegroupedData1.jpg 1630件 [詳細] filesp-grid-20080112.png 1961件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Google
WWW を検索 OKADAJP.ORG を検索
Last-modified: 2015-03-01 (日) 01:15:59 (1571d)