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

新規投稿はできません




凡例中に数式と数値変数を繰り返し書きたいのですが

TAG? (2007-06-27 (水) 21:50:40)

使用環境:R version 2.5.0 (2007-04-23), i386-apple-darwin8.9.1 です。

legendmiss.png

凡例として
「1 ミュー g/2 ミュー g/3 ミュー g」
のように数字と数式(ギリシャ文字)を混ぜて書きたいと思い

cols <- rainbow(3)
plot(sin, xlim = c(0, pi), ylim = c(-1, 1), ann=F, col=cols[1])
a <- NULL
b <- NULL
a <- expression(paste(1, mu, "g"))
b <- c(b, a)
a <- expression(paste(2, mu, "g"))
b <- c(b, a)
a <- expression(paste(3, mu , "g"))
b <- c(b, a)
legend(0.3, -0.2, b, col = cols, lwd=1)

上のようにコードを書いてみたのですが(プロットの部分はダミーです)
面倒なので凡例に関する部分を

a <- NULL
b <- NULL
for(i in 1:3)
{
	a <- expression(paste(i, mu, "g"))
	b <- c(b, a)
}
legend(0.3, -0.2, b, col = cols, lwd=1)

とすると、「iミューg/iミューg/iミューg」のように変数iが変数として認識されません。
paste関数の中なので当然なのかもしれませんが何か良い方法はないでしょうか?~

RjpWiki内のグラフィックス参考実例集:数式のプロットの例2に
「数式と数値変数の結合」とあるので参考にしてexpressionの代わりにsubstituteを使ってみたのですが

a <- NULL
b <- NULL
for(i in 1:3)
{
	a <- substitute(paste(x, mu, "g"), list(x=i))
	b <- c(b, a)
}
legend(0.3, -0.2, b, col = cols, lwd=1)

このようにすると「paste(1, mu, "g")/paste(2, mu, "g")/paste(3, mu, "g")」となり、
数値変数は反映されたのですがその他が駄目になるという結果になってしまいました。
解決策はあるのでしょうか?
どなたか対策をご存知でしたらご教授願えますでしょうか。

  • 今の場合なら,以下のようにしたらいかが? -- 2007-06-27 (水) 22:37:28
    cols <- rainbow(3)
    plot(sin, xlim = c(0, pi), ylim = c(-1, 1), ann=F, col=cols[1])
    b <- paste(1:3, "μg")
    legend(0.3, -0.2, b, col = cols, lwd=1)
  • なるほど!全角文字を使ってしまえば良いのですね。盲点でした。ありがとうございます。 -- TAG? 2007-06-28 (木) 08:55:34
  • legend(0.3, -0.2, sapply(1:3,function(i)(eval(substitute(expression(paste(mu,x,"g")),list(x=i))))), col = cols, lwd=1) #多分こう。 -- takahashi? 2007-06-28 (木) 09:06:08
  • legend(0.3, -0.2, parse(text=sapply(1:3,function(i)paste("paste(mu,",i,",g)",sep=""))), col = cols, lwd=1) #似たようなもん。 -- takahashi? 2007-06-28 (木) 09:14:31
  • さっそく試してみました。見事に上手くいきました!ありがとうございます。これならギリシャ文字に限らず分数や根号のようなものも混ぜられますね。 -- TAG? 2007-06-28 (木) 11:34:16
  • 一応付け加えておくと、こんなにややこしいのは、legendにvecto引数rとして渡さないといけないからで、そうでない場合(text, mtextなど)は、?plotmathに例があるとおり、もっと簡単です。 -- takahashi? 2007-06-28 (木) 13:03:38
  • 確かにtext, mtextだと何とかできたのですが、今回は自分の技量ではお手上げでした。再度お礼申し上げます。 -- TAG? 2007-06-28 (木) 17:14:47

Condor & R

HPC? (2007-06-27 (水) 00:31:04)

 Condor というPCクラスター上での Rの動作例ってありますか?

NAを含む行列を関数distで計算

su? (2007-06-25 (月) 13:34:43)

使用環境:R version 2.5.0 (2007-04-23) i386-apple-darwin8.9.1
NAを含む行列xyで行xと行yとの間の市街距離を計算しました。

x = c(1,1,NA)
y = c(1,2,NA)
xy = rbind(x,y)
dist(xy[,-3], method="manhattan") 	# 3列目を除いた距離は、1
dist(xy[,3], method="manhattan") 	# 3列目だけの距離は、NA
dist(xy, method="manhattan")		# 3列目を含めた距離は、1.5

dist(xy, method="manhattan")の値が、1でもなく、NAでもなく、1.5となるのは、何故でしょうか?

  • If some columns are excluded in calculating a Euclidean, Manhattan, Canberra or Minkowski distance, the sum is scaled up proportionally to the number of columns used. -- takahashi? 2007-06-25 (月) 14:38:33
  • 意訳すれば、NAを含む列は計算に使われないけど、その分距離の合計を上乗せするよ、ってことですね。 -- takahashi? 2007-06-25 (月) 14:47:44
  • ありがとうございました -- su? 2007-06-25 (月) 23:56:55

データフレームからデータが取り出せない

K? (2007-06-22 (金) 10:03:40)

度数分布表の値をもとに処理をしようとしているのですが、データフレームからデータが取り出せずに困っています。

> dbh_table <- table(data$DBH.Cls0)
> dbh_df <- data.frame(dbh_table)
> dbh_df
   DBH.Cls0 Freq
1        18    5
2        22   10
3        24   35
4        26   45
5        28   30
(中略)
20       58    5
21       60    5
22       72    5

としてデータフレームdbh_dfを作成したのですが、2列目のデータについては

> dbh_df[1,2]
[1]5

として欲しい"5"だけが得られるのですが、1列目のデータについては

> dbh_df[1,1]
[1] 18
22 Levels: 18 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 ... 72

となってしまい、必要な"18"の他に余分な1行が表示されてしまい、

x <- dbh_df[1,1]; y <- dbh_df[1,2]

としても、xに余分な1行が入ってしまい、処理できずに困っています。 どなたか、原因と対策をご教授いただけますでしょうか。 よろしくお願いします。

  • 処理が出来ないとは具体的にはどういうことでしょう? 1列目はfactor(因子)として扱われているので、処理によってはそのままのほうがよいし、数値に直す必要がある場合もある。as.numeric(x)とかas.integer(x)とかしてみると「余分な1行」は消えます。 -- takahashi? 2007-06-22 (金) 10:28:05
  • コメント、ありがとうございます。「処理ができない」と書いたのは、x <- dbh_df[1,1] としたときに、"18"だけを取り出すことができない、という意味です。また、as.numeric(dbh_df[1,1]) とすると 18 ではなく 1 が帰ってきます。どうやればdbh_df[1,1]を指定して"18"だけを得ることができるのでしょうか。(2行目以降についても同様。) -- K? 2007-06-22 (金) 11:09:27
  • 自己レスです。as.vector(x)とすることで、期待する値を得ることができました。どうもありがとうございました。 -- K? 2007-06-22 (金) 11:20:03
  • 18だけを取り出せているんですよ。ただし,18は数値ではなくてfactor ですよと。factor のリストは以下のようです,といっているだけです。もっとも,取り出されたのは factor なので,18という数値に変換することが必要なだけです。また,as.vector として取り出せるのは,character です。数値処理をしたければ,更に as.numeric する必要があるでしょう。
    data.frame を作った後に, dbh_df[,1] <- as.numeric(as.vector(dbh_df[,1])) としておくとよいでしょう。 -- 2007-06-22 (金) 11:48:09
  • なるほど、Rはfactorのリストを知らせてくれていたわけですね。as.numeric(as.vector())のご指摘、大変参考になりました。どうもありがとうございました。 -- K? 2007-06-22 (金) 13:41:48
  • ちなみに ?as.factor によれば、因子化された数値ベクトルを本来の表現に戻すお薦めの方法は as.numeric(levels(f))[f] です。 -- 2007-06-22 (金) 15:50:14

平行座標プロットの各ラインの選択

EDA? (2007-06-21 (木) 19:01:15)

 以前、マウスによるポイントの選択方法が掲載されていましたが、データマイニングで試用されている平行座標プロット(パラレルグラフ)をマウスで選択して、選択ラインの情報を見る方法はあるのでしょうか?

  • どのようなグラフで,どのように選択して,どのような情報を得たいのでしょうか?具体性がないので答えようがないかも。
    parcoord の中を見て,マウスでクリックされた点から一番近いデータ点を探し,そのy座標値から元の単位系へデータを逆算して,それと最も近いデータを持つケース(行)を選択して行を表示する(と同時にまあ,選択された線の色を変えてやる)。データ点でなくて線をクリックしても選択できるようにするには,クリックされた点の一番近くを通る直線を探し(クリックされた点の左右にある変数とそのy座標から直線の方程式をもとめ,一番近い線を探索する),あとは,前と同じように元のデータで一番近いものを探索して表示する。
    結構なプログラム量になりそう。 -- 2007-06-21 (木) 19:32:49

sub()で文字列を置き換えると,特定の文字で文字化けする

奥村泰之? (2007-06-21 (木) 13:38:26)

御世話になっております。 下記のコードに示したように,sub()で文字列を置き換えようとすると, 文字化けする現象が出て困っております。 対応策を御存じの方,お教え頂けると幸いです。

> sessionInfo()
R version 2.5.0 Patched (2007-04-27 r41346) 
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
> sub("AA", "山本", "AA氏")
[1] "山本氏"
> sub("AA", "五十嵐", "AA氏")
[1] "五女虫?x81"
  • sub("AA", "五十??嵐", "AA氏"); でいいんですが、内部エンコーディングって変えられなかったですかね? > 識者の方 -- takahashi? 2007-06-21 (木) 15:39:51
  • ただのバグです. CP932環境で, sub("A","?u8868","A") がうまく機能(2byte目が0x5cの文字)しません. 古い2.0.1日本語化までは動きます. UTF8の最適化の影響で結構前からCP932の扱いがうまくいかなくなりました. 私も直そうと思ってたり、直しかけのがあったりで手が回ってないのです.(^_^; バグレポート書くといいかも. -- なかま 2007-06-21 (木) 18:19:48
  • takahashi様,なかま様,ご回答を頂きありがとうございました。おかげさまで,古いバージョンのRを使うことで解決できました。 -- 奥村泰之? 2007-06-21 (木) 18:55:05
  • 先程パッチ付きレポートを送りました. -- なかま 2007-06-24 (日) 23:41:59
  • 2.5.1以降は大丈夫. -- なかま 2007-06-27 (水) 13:08:37
  • なかま様。R-develにレポートまでして頂き,誠に感謝しております。 -- 奥村泰之? 2007-06-27 (水) 14:52:16
  • version 2.5.1 で対応されていることを確認致しました。まことに,ありがとうございました。 -- 奥村泰之? 2007-06-29 (金) 15:54:46

散布図:同じような値を重ならることなく、横並びで表示したいのですが

ちょろ? (2007-06-20 (水) 14:29:14)

散布図を加工としております。
しかし、plot関数ではデータが一直線に並んでしまい、同じような値を重ねることなく、全ての個々の値の分布を直感的にとらえることは困難です。
一つ一つのドットが重ならないようにしたいのですが、どのように行えばよいか分からず困っております。
イメージ的には以下のURLの一番上にある画像の左上部のようなPlotをイメージしております。
http://www.ne.jp/asahi/room/kuro/ScatterMakerReadMe.html

ご教授いただければ幸いです。

  • library(gplots); ?space -- takahashi? 2007-06-20 (水) 14:37:56
  • 宣伝しておきます。http://aoki2.si.gunma-u.ac.jp/R/dot_plot.html 等はいかがでしょ。 -- 青木繁伸 2007-06-20 (水) 16:27:23
  • takahashi様、青木先生、アドバイスありがとうございます。アドバイスいただいた、gplotsパッケージのspace関数、ないしは青木先生の公開関数を利用し、チャレンジしてみます。本当にありがとうございますm(_ _)m -- ちょろ? 2007-06-20 (水) 16:49:36
  • 似たようなものとしてはstripchart(Y~X,vert=T,method="stack",pch=1,offset=1/2)という手もあります。 -- 中澤? 2007-06-21 (木) 15:43:46

3次元のデータをX-Yと色分けでplotする方法

ぼう? (2007-06-19 (火) 17:59:04)

お世話になります。XYZの3次元のデータがあり、Zを色分けしてx-y2次元のプロットに表示したい場合に以下のようにifを使って、一個づつ指定して離散的に色分けして表示することができますが、この部分を連続的にZの値にしたがったグラデーションで指定する方法をご存知の方がおりましたらお教え頂ければ幸いです。

mat <- matrix(rnorm(30),nrow=10)
plot(mat[,1],mat[,1],col=
 ifelse(mat[,3]>0.6,2,
   ifelse(mat[,3]>0.3,3,
     ifelse(mat[,3]>0,4,
       ifelse(mat[,3]>-0.3,5,
         ifelse(mat[,3]>-0.6,6,
         7
         )
       )
     )
   )
  )
)
  • データ数=色数でよいのなら -- g? 2007-06-19 (火) 18:46:11
    mat <- matrix(rnorm(30),nrow=10)
    plot(mat[order(mat[,3]),][,1],mat[order(mat[,3]),][,2],col=rainbow(10))
  • ありがとうございます。orderでzの順番にするとは妙技ですね。ですが、できれば データのレンジ=rainbowの色 のように値でグラデーションにできるのが理想です。 -- ぼう? 2007-06-19 (火) 19:17:27
  • 値を階級に分けて,その階級に色を適用するのでしょうから,データ数=色数ではちょっと困るんでしょうね。ということで,以下のような案を。 -- 2007-06-19 (火) 19:22:22
    mat <- matrix(rnorm(30), nrow=10)
    m = 5 # 階級の個数
    # suf は データの値により 1〜m になる
    # 例ではいい加減に分類したけど findInterval の第二引数で適切に指定のこと
    suf <- findInterval(mat[,3],  seq(min(mat[,3]), max(mat[,3]), length.out=m))
    plot(mat[,1], mat[,2], col=rainbow(m)[suf]) # rainbow でなくても良いから
  • じゃあ,無駄に大げさにして -- g? 2007-06-19 (火) 20:19:22
    plot(mat[order(mat[,3]),][,1],mat[order(mat[,3]),][,2],
    col=rep(c("red","orange","yellow","green","blue","violet"),
    hist(mat[,3],breaks=c(min(mat[,3]),-.6,-.3,0,.3,.6,max(mat[,3])),plot=F)
    $counts))
  • お二方ありがとうございます。分野によっては頻繁に使いそうなプロットですが意外とシンプルには行かないものなのですね。上の例はfindIntervalなんてソートに使えそうな関数がbaseにあったのは初めて知りました。gさんのはhistのcounts値を使ってcolに入れ込む大技ですね。大量だと処理時間がかかりそうですが、当面これをいじって解決しそうです。ありがとうございます。 -- ぼう? 2007-06-20 (水) 07:05:01

plotでx軸とy軸を任意の座標で交わらせることはできますか?

くに? (2007-06-19 (火) 10:24:13)

plotでx軸、y軸を(0,0)などの任意の座標で交わらせることはできますか?

kuni.png

のように。よろしくお願いします。

  • axis 関数のヘルプを読むべし -- 2007-06-19 (火) 10:58:24
    > x <- runif(20, min=40, max=100)
    > y <- runif(20, min=10, max=40)
    > plot(x, y, axes=FALSE, xlim=c(0,100), ylim=c(0, 80))
    > axis(1, pos=0)
    > axis(2, pos=0)
  • ほんとだ、helpにちゃんと書いてありますね。質問する前に、これからはもっとちゃんとヘルプを読むようにします(反省)。にもかかわらず、サンプルスクリプト付きでのコメント、ありがとうございました。 -- くに? 2007-06-19 (火) 11:20:04

等高線グラフィックのスケールについて

V9? (2007-06-13 (水) 14:31:14)

たまたま等高線の話題なので、続けさせていただきます。
contourで描いた等高線に、二次元の点をプロットしたく思っています。
単なるcontourではスケールをぴったり合わせられるのですが、filled.contourでは、凡例表示分だけ微妙にずれてしまうようです。
スケール(軸のメモリ)を合わせるのに、何かうまい方法はありますでしょうか。

  • 大変すみませんでした。具体的には、グラフィックス参考実例集filled.contour関数の例3のx軸、y軸のスケールが異なるので、それぞれ個別に設定したく思っておりますが、可能なのでしょうか? -- V9? 2007-06-13 (水) 16:00:59
  • 説明不足で大変申し訳ございませんでした。実は、例3の中で1点ではなく、20個くらいのデータをプロットしたく思っておりますが、エラーが出てしまう状況なのです・・以下にエラーceiling(length.out) : 数学関数に非数値引数が渡されました -- 2007-06-13 (水) 16:25:35
  • 大変すみません。削除の方法を教えて頂けないでしょうか。 -- V9? 2007-06-13 (水) 17:02:44
  • 実は、以下のようなものを書いて実行したら、エラーが出てしまいました。 -- V9? 2007-06-13 (水) 17:05:56
  • 上の方の「編集」をクリックして,編集画面の下の方に「テキスト整形のルールを表示する」というリンクがあるでしょう。それを良く読んで!
    ところで,何回も聞くけど ? filled.contour は読んだの? -- 2007-06-13 (水) 17:07:17
  • 以下のようなものを書いて実行したら、エラーが出てしまいました。 引数の長さが0とは、どういう意味なのでしょうか?
    以下にエラーif (axes) { : 引数の長さが0です
    
    filled.contour(
    	x,y,z,
    	levels=c(0,10,30,50,70,90,100),
    	xlim=c(-1,1),ylim=c(-100,100),key.axes=axis(4,seq(0,100,by=10)),
    	plot.title=title(main="main-title",xlab="x",ylab="y"),
    	key.title=title("ex"),
    	axes={axis(1);axis(2);points(exdata,pch=16,col="black")},
    	)
  • exdataは、以下のような型式です。 
         x  y
     1        -0.2 20
     2        0.8  40
     3        0.1  30
    ・・・
  • 描画できました。どうもありがとうございました。大変お手数かけまして申し訳ありませんでした。 -- V9? 2007-06-13 (水) 17:57:00
  • ちなみに、プロットした点は、値によって色を変えることは可能なのでしょうか。その場合、どのようなコマンドになるのでしょうか。 -- V9? 2007-06-13 (水) 17:58:08
  • ちなみに、プロットした点は、値によって色を変えることは可能なのでしょうか。その場合、どのようなコマンドになるのでしょうか。 -- V9? 2007-06-13 (水) 18:02:33
  • plot.axis の { } の中を,目的に従って書けば良いだけ。 -- 2007-06-13 (水) 18:03:58
  • ありがとうございました。点の色分けはできたのですが、判例がグラフ中に書き込めないようです。filled.contourの場合は、グラフ中に書き込むことはできない設定なのでしょうか。 -- V9? 2007-06-13 (水) 19:10:45
  • 判例の件、以下のようなものでうまく表示されませんでした。ちなみにRのバージョンは2.2.0でした。どうぞ宜しくお願い申し上げます。 -- V9? 2007-06-14 (木) 10:20:37
    filled.contour(
    	x,y,z,
    	levels=c(0,10,30,50,70,90,100),
    	xlim=c(-1,1),ylim=c(-100,100),key.axes=axis(4,seq(0,100,by=10)),
    	plot.axes={axis(1); axis(2); points(exdata, pch=16, col=c("red", 
                      "deeppink", "black", "deepskyblue", "blue"))},
           legend(0, 0, paste("value", c("0-20%", "20-40%",
                      "40-60%", "60-80%", "80-100%")), bg="gray",
                      col=c("blue", "deepskyblue", "black", "deeppink",
                      "red"), pch=16),
                      plot.title=title(main="main-title",xlab="x",ylab="y"),
           key.title=title("ex")
           )
  • そこで、plot.axexの中で修文したところ、うまく表示できました。お手数かけてご迷惑となり、大変申し訳ありませんでした。 -- V9? 2007-06-14 (木) 10:28:19

等高線用の外積計算

(2007-06-13 (水) 13:43:54)

等高線を描くために外積を計算しようとしていますが、lengthの長さがうまく設定できないためか、以下のようなエラーメッセージが出てしまいます。。
何かうまい設定法がございましたら、ご教示頂けると幸いです。

> x<-seq(-4,4,length=20)
> y<-seq(-10000,4000,length=20)
> z<-outer(x,y,function(xi,yj) 100/(1+exp(-(0.279-1.0344*x-0.0001*y)))
+ )
以下にエラーouter(x, y, function(xi, yj) 100/(1 + exp(-(0.279 - 1.0344 *  : 
        dim<- : dims [product 400] は object [20] の長さに整合しません
  • 基本的な間違いで、恥ずかしい限りです、、申し訳ありません。 -- 2007-06-13 (水) 14:22:07

無題

? (2007-06-10 (日) 20:00:15)

sourceで時系列オブジェクトを読み込もうとしたら、下記のエラーがでました。

以下にエラーattributes(.Data) <- c(attributes(.Data), attrib) : 
       不正な時系列パラメータが指定されました

このエラーはどのような意味なのでしょうか。
読み込もうとしているデータは以前に読み込みを成功した事があるデータを多少中身を変更しただけのものです。

  • 無題とはいやはや・・・ 読み込もうとしているデータも,変更の内容もわからないのにコメントのしようがない. 使用環境,問題の再現可能なデータセット,ロードしているライブラリなどの情報が ないと,誰も答えようがないと思いますが. -- 2007-06-10 (日) 20:52:44
  • 理由?不正な時系列パラメータを指定したからです!
    あ さん とは。如何に匿名とはいえ,もう少し考えたら? -- 2007-06-10 (日) 20:59:08

ヒストグラムのy軸目盛(スケール?)の変更は可能でしょうか。

くに? (2007-06-07 (木) 12:14:30)

 ヒストグラムのy軸目盛(スケール?)の変更は可能でしょうか。
 面積が0.3ヘクタールの調査区で樹木の胸高直径を測定しました。このデータから、直径階毎の本数分布をhist関数で作図したいと思います。ただし、そのとき縦軸の目盛を試験区(0.3ha)あたりの本数ではなく、1ヘクタールあたりの本数にしたいのです。
 そのようなことは可能でしょうか。よろしくお願いします。

  • 元のデータを 1/0.3 倍して、そのヒストグラムを描けばいいんじゃないっすか?それがいやなら,axis 関数を使って,自分で軸を描く。 -- 2007-06-07 (木) 12:41:40
  • 早速のコメントありがとうございます。元のデータをxとすると、barplotなら、barplot(table(x)/0.3)で作図できますが、hist関数ではhist(x/0.3)ではだめですよね。どうやったら「元のデータを1/0.3倍」できるのでしょうか。よろしくお願いします。 -- くに? 2007-06-07 (木) 13:25:28
  • 余計な御世話かもしれませんが、0.3ha毎の観測値が100個あれば、対象地域は30haですが、これを1haあたりの観測値100個とすれば、総計100ha観測したかのように見えなくもありませんね。それで構わなければ良いのですが。 -- 2007-06-07 (木) 15:31:23
  • ahaha,上のコメントで目が覚めた。私の提案していたのは,横軸の目盛りの変更だった。もうしわけない。申し訳ないついでに,ちょんぼのコメントも削除させてください。換算のための第二軸を設定する。 -- 2007-06-07 (木) 15:47:50
     set.seed(77789)
     x <- rnorm(100, mean=20, sd=5)
     par(mar=c(4,4,1,4))
     hist(x)
     axis(4, at=0:4*10, labels=round(0:4*10*(1/0.3),1))
  • 色々試行錯誤の末、少し求めているものに近づいてきました。  2cm刻みの直径階別本数を先に求めておいて、そのベクトルをtable関数に放り込んでからplot関数に渡したところ、求めるものに近いグラフが得られました(下図の下側)。
    upload.png
     でも、これではもう一息と言うところ。table関数で得られた表をデータフレームにしてから処理することで、なんとか達成できそうです。(でも、もっとスマートな方法はないものかしら。)アドバイスなどいただけましたら幸いです。よろしくお願いします。
    >area=0.3
    >layout(matrix(1:2),2)
    >DBHclass <- round((data$x+0.05)/2,0)*2 
    >#15cm以上17cm未満は18cmという、2cm刻みの直径階にするための処理
    >brks <- seq(14,74,by=2)
    >hist(DBHclass, breaks=brks,xlab="胸高直径階(cm)", ylab="立木本数(本/plot)")
    >
    >DBH_table <- table(DBHclass)
    >plot(DBH_table/area,type="h",xlab="胸高直径階(cm)",  ylab="立木本数(本/ha)")
    >
    >dbh<-data.frame(DBH_table/area)
    >dbh
      DBHclass Freq
    1        18    5
    2        22   10
    3        24   35
    4        26   45
    5        28   30
    6        30   55
    7        32   45
    8        34   60
    9        36   55
    10       38   90
    11       40   75
    12       42   65
    13       44   85
    14       46   30
    15       48   10
    16       50   10
    17       52   10
    18       54   15
    19       56    5
    20       58    5
    21       60    5
    22       72    5
    このデータフレームから、forループか何かで c(18,18,18,18,18,20,....,60,60,72)というようなベクトルをつくって、それをhist関数に渡す。力業?もっとエレガントな方法があるのでしょうか? -- くに? 2007-06-07 (木) 15:53:37
  • x<-rnorm(1000);hist(x,axes=F);axis(2,0:4*50,0:4*500) -- takahashi? 2007-06-07 (木) 16:27:16
  • x軸も一旦退避しておいてaxis()で書き込む -- 2007-06-08 (金) 08:59:12
     set.seed(77789)
     x <- rnorm(100, mean=20, sd=5)
     temp <- hist(x, breaks=seq(5,40,2), axes=FALSE, xlab="胸高直径階 (cm)", ylab="立木本数 (本/ha)", main="")
     axis(2, 0:5*3, round(0:5*3*10/3,0))
     axis(1, temp$breaks)
    hist.nag.png
  • parを活用すれば, -- なかま 2007-06-08 (金) 10:38:43
    set.seed(77789)
    hist(x,yaxt="n")
    axis(2,seq(par()$yaxp[1],par()$yaxp[2],5))
    固定値は限りなく少なくなる筈.
  • なるほど、
    axis(2, 0:5*3, round(0:5*3*10/3,0))
    としたところ、欲しいグラフが得られました。皆さん、ありがとうございました。 -- くに? 2007-06-11 (月) 15:57:12

和や部分リストを得る方法は

kd? (2007-06-06 (水) 10:16:48)

リストを与えて,和や部分リストを得たいのですが方法はあるでしょうか? はじめに下記のようにしました.

> dwt<-dwt(x,n.levels=8, boundary="reflection")
> str(dwt@W)
List of 8
$ W1: num [1:300, 1] -0.188  0.166 -0.150  0.170  0.115 ...
$ W2: num [1:300, 1]  0.6416  0.5408  0.0322 -0.3626 -0.2186 ...
  :
  :

次に,この dwt@W に対して和や部分リストを得たい思っています.例えば,

dwt@W[1:4]

に対して

dwt@W[[1]]+dwt@W[[2]]+dwt@W[[3]]+dwt@W[[4]]

を得たり,

dwt@W[1:4]

に対して

list(dwt@W[[1]][1:150],dwt@W[[2]][1:150],dwt@W[[3]][1:150],dwt@W[[4]][1:150]) 

を得たいと思っています.

  • dwt ってのは,どのパッケージに入っているのか?x, n にはどんなデータがどのように入っているのか?
    そんなこと教えなくても答えられる人だけ答えてくれというのなら,それで結構だが。
    上の方に長々と書いてある注意事項は読んでないんだろうね。「あなたの試みたことを回答者が再現できるような、データ、コード、エラー出力(そして、R のバージョン、使用 OS、使用非標準パッケージ名)を添えることをお勧めします」とあるよね。 -- 2007-06-06 (水) 10:32:46
  • 一個目はx<-list(1:3,4:6,9:12,13:15);do.call("+",x[1:2])、二個目はlapply(x[1:2],function(d)d[1:2])を応用する。 -- takahashi? 2007-06-06 (水) 10:34:06
  • ただし一個目は"+"が二項関数なので3項以上あると問題だったりする。多項版"+"ってないのかな? -- takahashi? 2007-06-06 (水) 10:43:00
  • そもそも,
    dwt@W[[1]]+dwt@W[[2]]+dwt@W[[3]]+dwt@W[[4]]
    ってやって正しい答えが出るならそれでいいじゃんか -- 2007-06-06 (水) 11:01:17
  • ご教示いただきありがとうございます.forループで回す関数で処理することにします.function(a) {b <- 0;for(i in c(1:length(a))) {b<-b+ai?};b} -- kd? 2007-06-08 (金) 11:38:17

指定データを除いたプロット

kd? (2007-06-03 (日) 15:05:17)

下記のように,指定データを除いてプロットすることは可能でしょうか?

最初に下記を行いました.

> newx<-modwt(x,n.level=8,boundary="reflection")
> str(newx2)
Formal class 'modwt' [package "wavelets"] with 11 slots
  ..@ W         :List of 8
  .. ..$ W1: num [1:3061, 1] -1.82e-06  3.58e-05 -1.12e-04  9.52e-05  4.66e-05 ...
     :
     :
  ..@ V         :List of 8
  .. ..$ V1: num [1:3061, 1] -2.14e-06  3.91e-05  1.48e-04  1.21e-04  7.28e-06 ...
     :
     :
> plot(newx)

このplot(newx) を実行すると,x,newx@W[1],newx@W[2],...,newx@W[8],newx@V[8] をプロットしてくれました.
次にnewx@V[8] だけを除いたプロットを行いたいので,下記の ???? に当たるような操作(?)を行うのかなと思っているのですが可能でしょうか?

> newx1<-????(newx)
> plot(newx1)

またnewx@W[8] だけをプロットすることも可能でしょうか?

> newx2<-????(newx)
> plot(newx2)
  • まず str(newx2) は str(newx) の、x,newx@W[1],... は newx@W[1],... の間違いですよね? newx はS4クラスオブジェクトらしく、W,V はそのスロットのようで、それぞれ8個のリストを含んでいるようです。クラス 'modwt' に対する plot 関数が何をプロットしてくれるのか(専用のメソッド関数があり特殊なプロットをしてくれる可能性があるので)がわからないので、正確な返答はできかねます。もし、単に複数のベクトルに対する散布図行列を書くだけでよいのなら、必要なベクトルだけを plot(newx@W[1], newx@W[2],...,newx@W[7],new@V[8]) もしくは plot(newx@W[8]) のように並べれば良いのではないですか。 -- 2007-06-03 (日) 17:16:33
  • どうもありがとうございます.計算結果 y1, y2, ... を単にプロットしたいだけなので,下記のようにしました.( plot(modwt(x,n.level=8,boundary="reflection")) もこれと似た感じのプロットを行ってくれます.) -- kd? 2007-06-05 (火) 15:44:25
  • par(mar=c(1, 1, 1, 1)) -- 2007-06-05 (火) 15:49:42

locfit パッケージのロードに失敗

kd? (2007-06-02 (土) 21:23:34)

locfit を v. 2.5 for MacOSX 上でロードしようとすると,下記のようなエラーが出て,ロードできないようです.

要求されたパッケージ akima をロード中です
エラー: パッケージ 'akima' をロードできませんでした
追加情報: Warning message:
'akima' という名前のパッケージはありません in: 
 library(pkg, character.only = TRUE, logical = TRUE, lib.loc = lib.loc) 

locfit は下記からダウンロードしたものです.
ftp://ftp.u-aizu.ac.jp/pub/lang/R/CRAN/bin/macosx/universal/contrib/2.5/locfit_1.5-3.tgz
ちなみにwavelets等の他のパッケージは正常にロードできて動いています.
ftp://ftp.u-aizu.ac.jp/pub/lang/R/CRAN/bin/macosx/universal/contrib/2.5/wavelets_0.2-2.tgz

  • 見たところ locfit パッケージは akima パッケージに依存するようですが、あらかじめインストールしてありますか。 -- 2007-06-02 (土) 21:33:46
  • この下の方にある「Macでパッケージvcdがロードできません 」を参考にされたし -- 2007-06-02 (土) 22:36:14
  • どうもありがとうございます.「パッケージの依存関係」なのですね.ただし当方では,パッケージインストーラでは,「一覧を取得」等を行うとエラーになってしまいますので,他の方法を調べる必要があるようです.(proxy設定済みで, curl -O ftp://ftp.u-aizu.ac.jp/pub/lang/R/CRAN/... 等も動いているので,本当はエラーの理由は理解していないのですが.) -- kd? 2007-06-03 (日) 14:58:43
  • ずぼらしないで、まずは akima  をインストールしたらどうですか。まだ他にも必要なら、また文句をいわれるでしょうが。 -- 2007-06-03 (日) 21:15:40
  • 重ね重ねご教示をいただき本当にありがとうございます.akimaをインストールでき,locfitがうまくうごきました. -- kd? 2007-06-04 (月) 01:52:58

カレントディレクトリのファイル名一覧を取得

mori? (2007-06-01 (金) 10:41:10)

カレントディレクトリに存在するファイル名の一覧をベクトル形式で取得したいのですがそのような関数はないでしょうか?たとえば、カレントディレクトリに file1 file2 file3 があったら c("file1", "file2", "file3") を返すような関数です。

  • dir <- dir(path=getwd(), pattern="file") のこと? -- 2007-06-01 (金) 10:52:37
  • ありがとうございます。dir()でできたんですね・・・ -- mori? 2007-06-01 (金) 11:08:43
  • list.files() http://takenaka-akio.cool.ne.jp/doc/r_auto/chapter_11.html -- 2007-06-04 (月) 14:23:58
  • dir と list.files は同じものです。 -- 2007-06-04 (月) 15:22:31
    > identical(dir, list.files)
    [1] TRUE

TIFFファイルをreadLinesで読み込みたい

akira? (2007-05-31 (木) 20:22:00)

測定器が出力するTiffファイルの1行目に測定条件が記入されています。エディター(xyzzy)では確認できるのですが、Rで読み込むとエラーになって読み込めません。

> x <- readLines("test.tiff", n=1)
Warning message:
 'test.tif' に関する readLines で不完全な最終行が見つかりました

画像ファイルをテキストで読み込む方法を教えてください。

OS:WinXPSP2
R version 2.5.0 Patched (2007-05-04 r41429) i386-pc-mingw32 です。

  • そうですか。MacOS X 10.4 では読めました(何の追加情報も提供できませんが) -- 2007-05-31 (木) 21:54:59
    こんなファイルを
    foo [9] > hexdump -c yellow.tiff
    0000000   M   M  ?0   *  ?0  ?0 002 220 301  ?0 301  ?0 301  ?0 356  ?0
    0000010 030 001   9 023  ?0   ]   ] 037  ?0 253   z   )  ?0 346   z   )
    読んでみました。
    R での結果は以下のようになりました。
    > x<-readLines("yellow.tiff", n=1)
    > x
    [1] "MM"
    これがあなたの望む結果なんでしょうか。あなたが読んだ tiff ファイルがどうなのかにもよるのかなぁ。
    > sessionInfo()
    R version 2.5.0 Patched (2007-05-18 r41630) 
    powerpc-apple-darwin8.9.1
  • exifならそんなパッケージがあったような気がしますし, バイナリーなら?fileとか?readBinあたりかな. 無理にテキストで読んでも, 画像のTIFFならあれげなデータになりますよ. -- なかま 2007-05-31 (木) 22:37:46
  • ありがとうございます。LINUXでもうまくいきませんでした。
    ubuntu Linux
    R version 2.4.1 (2006-12-18) i486-pc-linux-gnu 本当なら、実際に使用している画像をお見せしたいのですが、サイズが大きいので止めておきます。変わりにjpegでご説明します。
    デジカメ画像のjpegをEmacsでひらくと1行目に"Canon DIGITAL IXUS 950 IS"や"2007:03:20 13:08:11"といった文字列が見えます。この情報だけを取り出したいと思っています
    rgdalやrtiffのパッケージを用いるとtiffの数値情報をとれますが、この文字列がとれません(jpegにはrimageのread.jpegを使いましたがダメでした。)ただ、エラーが違うので参考にならないでしょうか。
    > x <- readLines(dir()[1])
    Warning message:
    '810IS_l0.jpg' に関する readLines で不完全な最終行が見つかりました
    > x
        [1]エラー:マルチバイト文字列が不正です
    なお、測定器のTiffファイルにも同じように測定条件の文字列が記入されています。
    なかま様、これはバイナリなのでしょうか?コーディングの問題なのでしょうか? -- akira? 2007-06-01 (金) 08:46:37
  • んと, JPEGやTIFFにはEXIF(ぐぐれば沢山あるでしょう)があります. RにはRexifというのもあります. これで読むか, バイナリのフォーマットデータをバイナリで扱って自分で読むかしないとダメでしょう. これはバイナリの問題でもコーディングの問題でも無く、「規格」の問題です. もちろん規格の中がUTF-16で書いてあったりISO-8859-1だったり色々はしますが. (昔,libjpegがこれでこけた覚えがある...) -- なかま 2007-06-01 (金) 09:30:49
  • なかまさま、ありがとうございます。バッチ処理をしたいので、Rexifを使ってみたいと思います。また、readBinは1行目の始めの文字しかとってこないので、引数をかえて試しています。
    > x <- readBin("test.tiff", what="character", endian="swap")
    > x
    [1] "II*"
    で、scan()を試してみたら、"エラー:マルチバイト文字列が不正です"はでていますが、最後まで読み込めたみたいです。
    > x <- scan("test.tiff", what="character")
    Read 26373 items
    > x
        [1] "II*"   
        [2] ""   
        [3]エラー:マルチバイト文字列が不正です
    これならいけそうです。
    これって、「ファイルTips大全の特定の文字コードのファイルを read」のこと? -- akira? 2007-06-01 (金) 10:43:07
  • readBin("test.tiff", what="", n=1) でどうなりますか?わたしんところの Win でやってみたところ,読めたようですが。 -- 2007-06-01 (金) 11:07:27
  • ありがとうございます。今はOSをubuntuで計算させているので、終わり次第、winでrebootして確認します。ちなみにubuntuではダメでした。 -- akira? 2007-06-01 (金) 11:26:36
  • バイナリをバイナリで読むには, -- なかま 2007-06-01 (金) 11:54:59
    con<-file("binary.file",open="rb")
    raw<-readBin(con,"raw",8)
    raw # Hexで出る
    # もし, rawの中がEUCで自分の環境がUTF-8なら"
    iconv(rawToChar(raw),"EUC-JP","utf-8")
    # めでたし,めでたし
    close(con)
    などとしないとダメでしょう.
  • なかまさまありがとうございます。教えていただいた通り実行してみたのですが、ダメでした(ubuntuで実行しています)。
    test.tiffの1行目の"II*"から目的とする文字列の間にある制御文字が後半を読めなくしているようです。jpegの場合も同じ状況のようです。 -- akira? 2007-06-04 (月) 06:45:47
  • IIはTiffのフォーマットがLittleEndian?と申しておりますから, そのまま読んで変換されているのですね. 制御部分を読み飛ばしさえすれば本来よいのですが... ubuntuならexifパッケージはありませんか? 無ければdebianからもってくればよいでしょう. 使いかたは man なり info なりでお調べいただき Rからは pipe で使えば良いでしょう. raw[4:10]とか切り出して使っていただけるものかと思ったのですが, 言葉足らずでお手数をおかけします. # 宇宙の言葉は話さないでと, 昔から言われてますが, 年々酷くなります. -- なかま 2007-06-04 (月) 11:22:24
  • ありがとうございます。生まれがwinなものでfileとかpipeとは縁遠い世界で生きております。すこし時間がかかってしまいますが、進展しましたらご連絡いたします。 -- akira? 2007-06-04 (月) 12:20:10
  • libtiffライブラリtiffinfoが望む結果を出すことがわかりました。
    libtiffはrtiffパッケージが使用しているのですが、tiffinfoは実装されていません。systemで呼び出すことにします。 -- akira? 2007-06-07 (木) 20:06:48
  • とりあえず、fileとscanだけでもほぼ望む結果が得られることが分かりました。
    con <- file(description=i, open="rb", blocking = TRUE)
    x <- scan(con, what="logical", n=20, sep=NULL, comment.char="")
    どの辺が「ほぼ」かというと、^@のコードの直後の文字が取得できていません。とりあえず満足です。 -- akira? 2007-06-19 (火) 15:54:23

オンボードRAIDは使ってますか?

通りすがりの佐々木? (2007-05-31 (木) 16:22:29)

R用PC(Athlon64,3200,2G,WXPSP2)をMBオンボードのRAIDを使っていたのですが、システム領域(RAIDなし)とデータ領域を独立したHDD(RAID1)に分けてしました。先月のはじめ、起動中のリソース使用率を見てみたら、いままで100%くらいだったのが70%になっていました。そろそろ再インストールかなと思っていたら、システムのHDDが突然死されました(RAIDなし!!)。BIOSのSMARTでディスクを交換せよとの警告が表示されてしまいました。データはデータ用のHDDに保存し、スクリプトはシステムにいれていたのですが、大部分のデータは幸いにも生きていましたが、1年分のスクリプトが消失しました。けっこうキツイ・・・。やはり、CPUを100%24時間*7日なんてことをずっとやっていたからでしょうか、アクセスが多くなるとHDDの熱がすごいんですね。CPUばかり気を使っていたので、HDDは対策が疎かになっていたのでした。反省。RAID2系統+ファン増設を行ったのはいうまでもありません。

  • 昔家のパソコンで、シミュレーションを夏休み中ぶっつづけでやったことがあります。冷房もないので家庭用扇風機を回しつづけました(笑)。ところで、これ質問ですか? -- 2007-05-31 (木) 16:50:58
  • みなさんのRAID環境をしりたかったりする。Rを使うためには、そんなに気合が必要なのだろうか、いや必要なのだ・・・。 -- 佐々木? 2007-05-31 (木) 16:54:14
  • 私は,RAID なんていうのは,今の今まで知らなかった。私は R を使う資格がないのだろうか。。。と悩んだりする。 -- 2007-05-31 (木) 22:02:40
  • RAIDなんてものを信用してはいけない。外付けHDに、こまめなbackupこそ命。 -- okinawa 2007-05-31 (木) 22:26:36
  • 性能的な物だと, 廉価版RAIDはOS+CPUに処理を丸無げするのでお薦めしません. 障害予防的な面ではRAIDより扇風機の方が効果が高いといわれています. ミラーリング(RAID1)とバックアップは等価では無いので, 別途行うのが良いです. 廉価なオンボードRAIDならソフトウェアRAIDでも大差はありません.(故障時の復帰時に操作ミスやバグでいや〜んと言う事もあります) 最近は少し高い「ニアラインストレージ」と言うディスクもあるので, 費用効果(いずれにせよバックアップは別に必要)を考えればオンボードRAIDの効果などはたかが知れています.(2000年から2005年頃のディスクは良く壊れたような気がします) -- なかま 2007-06-01 (金) 11:10:46
  • やはりオンボードRAIDは負荷がかかるのですね。実はRAIDにしてからRAIDアレーにアクセスするとCPUが40%も使われてしまうのです。IDEでは変化にきずきませんでしたので、ひょっとしてCPUに負荷をかけているのではないかと思っていたのですが、疑いは確信にちかづいてきました。扇風機もいいのですが、研究室内に扇風機を数台おいて使用すると竜巻が生じて・・・ついでに置き場所が・・・。ニアラインストレージ、初めて知りました。NASのようなものなのですね、頻繁にアクセスするものデータはカレントに置いて、終了したものは外部に移動するという手もありますね。シリコンディスクにファイルを置いて、処理が終わったらまとめて外部書き出し、次のデータ処理に移動という手もありですね。RAIDをやめて、外部にファイルサーバーを置いてNASのように使用し、シリコンディスクとのやりとりを置けば、一度にアクセスが集中することを少し軽減できそうですね。 -- 佐々木? 2007-06-01 (金) 12:36:57
  • その構成ならeSATAが使えればいいですね。NASだとEther自体がネックになると思う。(2次バックアップにNASを使うのはいいと思う。) -- okinawa 2007-06-01 (金) 16:11:18
  • GentooのHPにSOFTraidについての翻訳がありました。Linux RaidというHPの訳でしたが、これによれば、近年のSouth bridgeのRAIDはBIOSによるソフトRAIDとなっています。なるほど、安くて良い物などないのですね。 -- 佐々木? 2007-06-03 (日) 20:40:40
  • 実は、先日のクラッシュでHUBもお亡くなりになったので、ギガスイッチにしてみました。外部が遅いのですが、ハプ内では結構いいです。700Mが5分で移動できます。でもSATAの3G3にはかなわないですね、eSATAが使えるMBも検討する価値がありそうですね。計算中は反応がほぼなくなるので、あやしい臭いに注意しつつやってみます。 -- 佐々木? 2007-06-03 (日) 20:45:48
  • GigaEther?の欠点はjumbo packetにNIC&HUBが対応していないと転送スピードが上がらないところです。eSATA&JumboPacketReadyNIC&HUBだと言うことなし。しかし、クラッシュでHUBも死んだとは。落雷?漏電?どちらにしても電源系統が危険そうですね。 -- okinawa 2007-06-04 (月) 09:01:45

2変数関数の等高線プロットについて

GG? (2007-05-31 (木) 11:55:06)

2変数(x、y)関数の値(z)を求めることなく、
式の形から等高線を直接描ける方法はあるのでしょうか?
例えば、z=1/(1+exp(-(10x-3y)))のようなグラフを直接3次元で描画し、
適当に等高線で色分けした後に、
さらに重ねがきで幾つか点もプロットしたいと思っています。

  • 例示のようなきれいな関数なら,z=constant として,y=f(x, constant) を導けば描けないことはないでしょう。おやりになってみますか? -- 2007-05-31 (木) 12:38:38
  • 実際には例示よりももう少し複雑な式を考えておりますが、y=・・・を導かずに描くことは難しいでしょうか。 -- GG? 2007-05-31 (木) 13:24:44
  • y=もなくて,z も計算しなくて,コンピュータの中にこびとさんでもいれば,寝ている内に描いてくれるかも知れませんね(^_^) -- 2007-05-31 (木) 14:14:29

日本語入力してもグラフが表示されない

のり? (2007-05-27 (日) 20:39:34)

Rインストール後、コマンド(R console)で日本語入力しても全くグラフが表示されません。
(例)

x=0:50
y=dbinom(x,50,0.25)
plot(x,y,type=’h’,xlab=’x’,ylab=’y’,main=’二項分布’)
+

(注、type=hはヒストグラム表示、xlab、ylabはx、y軸の名前、mainはグラフのタイトルを表す。日本語を入力しない他のコマンドでは、グラフが表示されます。)
(使用環境 windowsXP)
この解決方法をご存じの方がおりましたら、ご教示いただけないでしょうか?
(尚、この質問は最初、初級者コースが閉鎖されていると認識し、中級者コースに投稿しました。ご指摘いただきましたokinawa様にはお手数おかけしましたことをお詫びいたします)

  • 張り付けられた文字と「日本語を入力しない他のコマンドでは」と言う点を加味して, おそらく入力後になにかメッセージは出てませんかね?. たぶんエラーだと言ってるんだと思いますが. 'と’の違いは見えませんか?. 上のとおりとすれば貴方が入力したのはRhightSingleQuotationMark?(?u2019)でApostorophe(?u0027)ではないんでないかと. -- なかま 2007-05-27 (日) 21:28:47
  • type='h' の h までもが,全角じゃ〜〜ないですか。なんで??半角にすれば,うごくでしょ。以下と,よく,見比べれば?まだ強情に = を使っているし。少しは,アドバイスを参考に再考してみればよいと思うのだけど。(初心者だからといって,なんでも許容されるものではない) -- 2007-05-27 (日) 21:37:04
    x=0:50
    y=dbinom(x,50,0.25)
    plot(x,y,type='h',xlab='x',ylab='y',main='二項分布')
    baz.png
  • 有り難うございました。日本語に注意を取られすぎ、いろいろやっておりました。)の部分が全角であったのに気づきませんでした。 -- のり? 2007-05-27 (日) 22:49:33
  • 本当の正解はこれなのでは? -- okinawa 2007-05-28 (月) 00:19:36
    x <- c(0:50)
    y <- dbinom(x,50,0.25)
    plot(x,y,type='h',xlab='x',ylab='y',main='二項分布')
  • あと、sessionInfo( )とコンソールから実行してみてくださいね。あなたのRの環境情報が出てきますから。実行しているOSによって回答が若干違ってくることもあるので・・・。 -- okinawa 2007-05-28 (月) 01:04:52
  • c(0:50) は冗長でしょう。0:50 でよい。= を使って悪いわけではない。慣例で,<- の方が優勢なだけでは?そう言う観点からは,' よりは " の法が優勢かも。 -- 2007-05-28 (月) 08:20:49
  • c(0:50)は冗長だとは思いますが、R言語の作法に沿ってコードを書いたつもりでした。(間違ってます?)あと、=は「Rの基礎とプログラミング技法」においても利用は勧められないと書いてあります。個人的には、初級者には冗長性や同義性の議論はかえって混乱を招くと思っていますので、難しいところです。 -- okinawa 2007-05-28 (月) 09:01:04
    x <- 0:50
    y <- dbinom(x,50,0.25)
    plot(x,y,type="h",xlab="x",ylab="y",main="二項分布")
  • 前にも投稿記事がありましたが、付値演算子として = を使うのは、やめたほうが賢明です。理由は foo(x <- 1:10) などというしばしばコードの簡略化に便利な構文が使えないこと。これは関数内では = が関数引数指定と紛らわしいからでしょうね。 -- 2007-05-28 (月) 09:31:41
    > order(x <- 1:10)
     [1]  1  2  3  4  5  6  7  8  9 10
    > x
     [1]  1  2  3  4  5  6  7  8  9 10
    > order(X = 1:10)
     [1]  1  2  3  4  5  6  7  8  9 10
    > X
     エラー: オブジェクト "X" は存在しません
  • order(x <- 1:10) がコードの簡略化に役立つとは思いませんが,議論はやめておきます(^_^)
    ? c で,This is a generic function which combines its arguments. とあり,更に c(..., recursive=FALSE) の Arguments は ... objects to be concatenated で,arguments と objects と複数になっています。一つでも悪くはない c(3) みたいに。 -- 2007-05-28 (月) 09:59:06
    > identical(c(1:10), 1:10)
    [1] TRUE

Macでパッケージvcdがロードできません

ピーター? (2007-05-25 (金) 18:46:30)

Mac OS X 10.4.9, R 2.5.0 GUI 1.19(4308)を使っています
「パッケージとデータ」のCRANからvcdをインストールし、「Rパッケージマネージャー」でvcdの状態のところをクリックしてもロードできません。コンソールには

要求されたパッケージ colorspace をロード中です
エラー: パッケージ 'colorspace' をロードできませんでした
追加情報: Warning message:
'colorspace' という名前のパッケージはありません 
    in: library(pkg, character.only = TRUE, logical = TRUE,  lib.loc = lib.loc) 

とあります。また、library(vcd)でも同様でした。
一体なぜののでしょうか。教えてください。お願いします。

  • colorspaceというパッケージをインストールする必要があると思うのですが、入っていますか? -- 2007-05-25 (金) 19:16:54
  • パッケージインストーラでインストールするときに,いつも,「パッケージの依存関係を解決する」をクリックしておいてから「選択をインストール」するようにしておけば間違いないですね。面倒だから,いっそのこと,全部のパッケージをインストールしておけばいちいち引っかかることもないので,心が平安に保てるでしょう。- 2007-05-25 (金) 19:41:49
  • なるほど!colorspaceをインストールしたら、vcdもロードできました。ありがとうございます。~これからそうします。ありがとうございます。 -- ピーター? 2007-05-25 (金) 23:36:04

三次元プロットについて

dmc? (2007-05-23 (水) 13:42:22)

scattaerplot3d関数で三次元プロットした点に、persp関数で二次元曲面を描きいれようと思っているのですが、軸やグラフの角度がうまく合いません。
点と曲面をうまくひとつのグラフ内に描く方法がございましたら、教えて頂けないでしょうか。

  • perspだけでやっちゃう。 -- takahashi? 2007-05-23 (水) 15:10:43
    x <- seq(-10, 10, length= 30)
    y <- x
    f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
    z <- outer(x, y, f)
    z[is.na(z)] <- 1
    pmat<- persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
                 ltheta = 120, shade = 0.75, ticktype = "detailed",
                 xlab = "X", ylab = "Y", zlab = "Sinc( r )")
    
    p<-list(x=runif(100,-10,10),
            y=runif(100,-10,10),
            z=runif(100,0,5))
    points(trans3d(p$x,p$y,p$z, pmat), col = 2, pch =16)
  • ご回答頂き、どうもありがとうございました。無事に描画できました。お礼申し上げます。 -- dmc? 2007-05-24 (木) 11:10:49
  • 追伸です。曲面を塗らず、透過させるにはどのようなコマンドがありますか?曲面の下にも点があるので・・ -- dmc? 2007-05-24 (木) 11:20:18
  • perspの引数colを0とかNULLとか"transparent"とかってことかな。 -- takahashi? 2007-05-24 (木) 11:48:00
  • 引数を0やtransparentにしてみたら、うまく描画できました。どうもありがとうございました。 -- dmc? 2007-05-24 (木) 15:14:31

データがマイナスとなる個数

port? (2007-05-23 (水) 12:02:52)

Rの超初心者です。
10000個くらいのデータがあって,そのデータの中でマイナスとなるデータ個数を知りたいのですが,どのうようにすればよいでしょうか?
色々調べたのですが,よくわかりませんでした。
どのたかご教授ください。

  • dat<-rnorm(10^4); sum(dat<0) -- takahashi? 2007-05-23 (水) 12:40:52
  • どうもありがとうございます。助かりました。 -- port? 2007-05-23 (水) 13:09:27
  • dat<-rnorm(10^4); length(which(dat<0))、でしょう。 -- とも? 2007-07-19 (木) 12:55:27

Excelからのコピペが上手くいきません

みかん? (2007-05-14 (月) 04:58:23)

「Tips紹介」→「R にファイルを読み込む tips 集(暫定版)」→
「Excel からコピーしたセル範囲を R に読み込む (何でも掲示版の青木繁伸さんの記事を転載)」
に記載されております方法が上手くいきません。
記載事項の通り、

> from.excel <- function(nc)
+ {
+        matrix(scan(""), byrow=TRUE, nc=nc)
+ }
> from.excel(3)

としたのち、 エクセルで

1 F 4
2 F 3
3 M 2

という3×3のセル範囲をコピーし、貼り付けると、

> from.excel(3)
1: 1F4
1: 2M3
以下にエラーscan("") : scan 関数は '$s' を期待したのに、得られたのは '$s' でした
> 3F2
エラー:"3F2" に構文エラーがありました
> 

というメッセージが出てきます。
一方同じExcelのコピー範囲を一度テキストエディタに貼り付け、"?t"を" "に置換してからコピーし、

> x <- read.table(stdin())

の後に貼り付けると、上手くいきます。
しかしこの方法でも、Excelでコピーしたそのままを貼り付けると、

> x <- read.table(stdin())
0: 1F4
1: 2M3
2: 3F2
3: 
> 

となってしまいます。
どうもExcelでコピーした列セル間の区切り子が、R上で認識されていない様子です。
どうしてこうなってしまうのか、ご教授をお願いいたします。
当方WindowsXP上で、ver.2.2.1を使用しております。Excelのバージョンは2000です。
よろしくお願いいたします

  • エクセルで範囲選択してコピー後、read.table("clipboard")じゃだめですか? -- takahashi? 2007-05-14 (月) 08:53:45
  • Windowsね。x をプリントしてみたらどうなっているの?mac だと以下のようにちゃんと行くが。 -- 2007-05-14 (月) 11:38:44
    > x <- read.table(stdin(),header=FALSE)
    0: 1	F	4
    1: 2	F	3
    2: 3	M	2
    3: 
    > x
      V1 V2 V3
    1  1  F  4
    2  2  F  3
    3  3  M  2
  • 青木先生の方法ではmatrix関数を使用しているため
    すべての値がすべて数値あるいは文字列になっていないとうまくいかないようです。

    そこで直接の回答ではありませんが、僕自身は
    from.excel2 <- function(h = 0)
    {
      dat <- (read.delim("clipboard", header=h))
    }
    のような関数を別に定義してExcelのデータを読み込んでいます。
    タイトル行がある場合は header=T とします。
    この場合にはデータフレームとして読み込まれるのでうまくいきました。

    また、僕自身はstdin()を知らないです。-- あら? 2007-05-14 (月) 08:22:32
  • dat に代入するのではなく(dat って,関数外のオブジェクトですよね),そのあとで return(dat) とした方がよいですよ。return(read.delim("clipboard", header=h)) でもよいけど。また,h=0 ではなく,h=FALSE とする方がよい。 -- 2007-05-14 (月) 17:26:38
  • 教えるつもりが逆に教えていただく型になってしまいました。ありがとうございました。 -- あら? 2007-05-14 (月) 18:18:20

直交実験を計画してくれる関数

TH? (2007-05-13 (日) 09:21:44)

直交実験を計画してくれる関数はどこかのRパッケージに含まれていないでしょうか?"Orthogonal Array Design"を中心にいろいろ検索をかけても見つかりません.S-PLUSには直交実験を計画してくれるとても便利な 関数oa.design() があります.もちろんoa.design() と打てば,その定義が画面にでてきます.その中に出てくる関数も含めて(途中で追えなくなるので,追えるものだけ)まとめて,Rの画面で流すとエラーがでてしまい,移植に失敗します.中身を少し読んでも私の実力ではとうてい書き換えられません.直交配列による実験は我が国では,大変多用されているので,どこかにないだろうかとお聞きする次第です.

  • CRAN にあるパッケージ crossdes あたりは使えませんか。 -- 2007-05-13 (日) 12:29:53
  • Rseek -- 2007-05-13 (日) 12:29:54
  • ありがとうございます.でもcrossdesはラテン方格と釣り合い型不完備ブロックのデザインでした.Rseek を見させていただきましたがどうやらないみたいですね. -- TH? 2007-05-13 (日) 13:16:02
  • oa.design とそれに関する関数を用意して,その段階ではどのようなエラーメッセージが出たのですか?こちとら,Sが使えませんので。 -- 2007-05-13 (日) 21:59:40
  • 移植しようと昨日からいろいろやっているうちに,全体の移植は無理だけど自分の研究に必要な機能の部分だけ,エラーを見ながら,スクリプトを読んで一部移植することに成功しました.おさわがせしました.有難うございます. -- TH? 2007-05-14 (月) 16:28:43

確率微分方程式

SDE? (2007-05-10 (木) 19:02:49)

確率微分方程式で検索してもヒットしなかったのですが、
解を求めるなど、確率微分方程式を取り扱ったパッケージは無いのでしょうか?
Rの守備範囲だと思うのですが…

  • どうやって検索したのか??site: r-project.org "Stochastic Differential Equations" でググっただけで,sde Simulation and Inference for Stochastic Differential Equations とか出てくるんだけど(もちろん,中身は見てないので適切なのかどうか知らない) -- 2007-05-10 (木) 19:32:09
  • ありがとうございます。検索が雑でした。CRAN.packagesで一覧を取って来ると出てこないんですよね…。 -- SDE? 2007-05-11 (金) 11:28:52

Fizz-Buzz問題

X Jr.? (2007-05-09 (水) 22:20:18)

どのカテゴリに書こうか迷ったのですが、とりあえず初級Q&Aにします。
どうしてプログラマに・・・プログラムが書けないのか?」という翻訳記事の中の「Fizz-Buzz問題」が一部のネット上で話題になっています(例えばはてなブックマーク)。その内容は、「1から100までの数をプリントするプログラムを書け。ただし3の倍数のときは数の代わりに「Fizz」と、5の倍数のときは「Buzz」とプリントし、3と5両方の倍数の場合には「FizzBuzz?」とプリントすること。」というものです。Rならどのようなプログラムができますか? お遊び感覚でお願いします。

  • [どうしてプログラマに・・・プログラムが書けないのか?]それは,そいつはプログラマじゃないからだ。 -- 2007-05-09 (水) 22:43:42
    > x <- 1:100
    > print(ifelse(x%%15==0,"FizzBuzz", ifelse(x%%3==0,"Fizz", ifelse(x%%5==0, "Buzz", x))))
      [1] "1"        "2"        "Fizz"     "4" 
      [5] "Buzz"     "Fizz"     "7"        "8"    
      [9] "Fizz"     "Buzz"     "11"       "Fizz" 
     [13] "13"       "14"       "FizzBuzz" "16"   
     [17] "17"       "Fizz"     "19"       "Buzz"
          :
     [93] "Fizz"     "94"       "Buzz"     "Fizz" 
     [97] "97"       "98"       "Fizz"     "Buzz"
  • やり方は幾つもあるけど,どれが一番「プログラマが書いたプログラム」というのだろうか? -- 2007-05-09 (水) 22:51:03
    x <- 1:100
    x[1:33*3] <- "Fizz"
    x[1:20*5] <- "Buzz"
    x[1:6*15] <- "FizzBuzz"
    print(x)
  • せめて,これを拡張して,1000までの素数をエラトステネスの篩のやり方で求めよとかにしてほしいな。 -- 2007-05-09 (水) 22:56:04
  • cat の方がいいかな。それと ==0 の代わりに ! を使う -- 2007-05-09 (水) 23:03:04
    > x<-1:100;cat(ifelse(!x%%15,"FizzBuzz",ifelse(!x%%3,"Fizz",ifelse(x%%5,x,"Buzz"))))
    1 2 Fizz 4 Buzz Fizz 7 8 Fizz Buzz 11 Fizz 13 14 FizzBuzz 16 17 
    Fizz 19 Buzz Fizz 22 23 Fizz Buzz 26 Fizz 28 29 FizzBuzz 31 32 
    Fizz 34 Buzz Fizz 37 38 Fizz Buzz 41 Fizz 43 44 FizzBuzz 46 47 
    Fizz 49 Buzz Fizz 52 53 Fizz Buzz 56 Fizz 58 59 FizzBuzz 61 62 
    Fizz 64 Buzz Fizz 67 68 Fizz Buzz 71 Fizz 73 74 FizzBuzz 76 77 
    Fizz 79 Buzz Fizz 82 83 Fizz Buzz 86 Fizz 88 89 FizzBuzz 91 92 
    Fizz 94 Buzz Fizz 97 98 Fizz Buzz
    2文字だけ短いもの(!?)
    x<-1:100;cat(ifelse(x%%15,ifelse(x%%3,ifelse(x%%5,x,"Buzz"),"Fizz"),"FizzBuzz"))
  • データを考える言語の立場から, 条件分岐はどうでしょうね. お題は決まってますが, いじくり回すと言うスタイルで. -- なかま 2007-05-09 (水) 23:24:59
    M<-expression(1:100)
    x<-eval(M)
    x[!eval(M)%%3]<-"Fizz"
    x[!eval(M)%%5]<-"Buzz"
    x[!eval(M)%%3&!eval(M)%%5]<-"FizzBuzz"
  • Rでのベクトルと行列の扱いなどの特徴を活かして(^_^;) -- 2007-05-09 (水) 23:54:15
    x<-1:105
    dim(x)<-c(3,105/3)
    x[3,]<-"Fizz"
    dim(x)<-c(5,105/5)
    x[5,]<-"Buzz"
    dim(x)<-c(15,105/15)
    x[15,]<-"FizzBuzz"
    cat(x[1:100])
    これも,ほんの少しだけ短くなる
    a<-function(i,s){dim(x)<<-c(i,105/i);x[i,]<<-s}
    x<-1:105
    a(3,"Fizz")
    a(5,"Buzz")
    a(15,"FizzBuzz")
    cat(x[1:100])
  • いろいろと面白い例をありがとうございました。Rはループを使わなくてもよいので結構短いものもできますね。なかまさんの例は私にはとても思いつきそうにありません。「ほかにもこんなふうにできるぞ」という方がおりましたら是非。 -- X Jr.? 2007-05-10 (木) 22:14:47
  • 人に聞く前に,あなたのプログラム例をまな板に載せてから聞いた方がよかったかもね(^_^;) -- 2007-05-10 (木) 22:51:54
  • 初心者コーナーとあれば一番Rらしくない(笑)のもあげておくべきでしょう。-- 2007-05-10 (木) 23:04:52
    > for (i in 1:100) {
    + if(i%%15==0) cat("FizzBuz,")
    + else if (i%%3==0) cat("Fizz,")
    + else if(i%%5==0) cat("Buz,")
    + else cat(i,",")}
  • そうですね。順序が逆になりましたが、初心者が書くとこんなふうになるものが、上の皆さんのようにできますという例です。 、、、と書いていたらどなたかが似たものを載せてくれました。すみません(^^; -- X Jr.? 2007-05-10 (木) 23:10:49
  • いや,実際,挙げられた一番単純なものが一番確実に動く,一番時間を掛けずに(それこそ2分で)書けるプログラムだと思いますよ。私が採用担当者なら,このプログラムを書いた人を採用するかも知れませんね。(こったプログラムとかは,最初のアルゴリズムを用いて,文字数を最小にするには,とか,工夫しているだけ。場合によっては,処理系の持っている特殊性を活用するために,処理系を選択しているだけでしょう。ある言語は,この問題を解決するための関数を提供しているかも知れませんよね。というか,Rだって,この問題を解く関数パッケージがあれば,FIzzBuzz?(3,5,100) ですむかもしれませんものね(^_^;) -- 2007-05-10 (木) 23:31:37
  • 私が採用係なら直前のコードを書いた人物(=> 実は私本人)は採用しませんね。15,3,5 の順で条件判断すべきではなく、3,5,15の順で(つまり条件が満たされやすい順で)判断すべきです。参考までに3,5,15を順に変えたコードの千回の実行時間(コンソール出力は抑制)を書いておきます。一番時間が短いのが3,5,15の順で条件判断したものです。-- 2007-05-11 (金) 00:09:35
    [1] 2.744 0.000 2.757 0.000 0.000
    [1] 2.737 0.000 2.748 0.000 0.000
    [1] 2.712 0.000 2.723 0.000 0.000
    [1] 2.676 0.004 2.698 0.000 0.000
    [1] 2.628 0.000 2.633 0.000 0.000
    [1] 2.632 0.000 2.649 0.000 0.000
  • 私が採用担当なら, プログラマーになりたいと言う人は, 思い留めさせます... -- なかま 2007-05-11 (金) 09:48:07
  • > 15,3,5 の順で条件判断すべきではなく、3,5,15の順で
    注意しないと,15 の倍数を3の倍数だと判断して(間違いじゃないが),Fizz と表示してしまう(これは間違い)ことにならないように注意が必要。 -- 2007-05-11 (金) 15:24:36
  • 実は,一番起こりやすいのは 3,5いずれの倍数でもないということ。しかし,これを真っ先に判定する人はいないと思われる。
    3,5,15の順で正しく判定するためには,かえって手間が掛かる。実際,所要時間は無視できないことになる。
    ちょっと計算時間が短くなっても,バグ持ちプログラムでは意味がない。 -- 2007-05-11 (金) 16:14:00
    > a <- function(){
    + 	for (i in 1:1e7) {
    + 	if(i%%15==0) 15
    + 	else if (i%%3==0) 3
    + 	else if (i%%5==0) 5
    + 	else 0 }
    + }
    > system.time(a())
       ユーザ  システム      経過 
       34.522     0.116    34.522 
    
    > b <- function() {
    + 	for (i in 1:1e7) {
    + 	if (i%%3==0 && i%%5!=0) 3
    + 	else if(i%%3!=0 && i%%5==0) 5
    + 	else if(i%%15==0) 15
    + 	else 0 }
    + }
    > system.time(b())
       ユーザ  システム      経過 
       46.421     0.133    46.393
    
    > bb <- function() { # こうすれば,無駄な比較をしないで,かつ「正しい」
    + 	for (i in 1:1e7) {
    + 	if(i%%3==0) {
    + 		if (i%%5==0) 15 else 3
    + 	}
    + 	else {
    + 		if (i%%5==0) 5 else 0}
    + 	}
    + }
    > system.time(bb())
       ユーザ  システム      経過 
       29.060     0.138    29.072
  • one-liner -- 2007-05-12 (土) 19:23:10
    replace(replace(replace(1:100,1:33*3,"Fizz"),1:20*5,"Buzz"),1:6*15,"FizzBuzz")

.RData の怪

昔の名前で出ています? (2007-05-09 (水) 16:13:25)

R 2.5.0 をインストールして、使おうとしている途中です。~
ディレクトリの設定をして、q() で「作業スペースを保存」しても、実際には何も保存されません。~
それでは、と、「作業スペースの保存」をして、できた .RData をダブルクリックしたら、R 2.4.0 が立ち上がります。なぜでしょう。

> sessionInfo()
R version 2.5.0 (2007-04-23) 
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"  "methods"  
[7] "base" 
> # この後、ワーキングディレクトリを設定    
> getwd() # 設定の確認
[1] "C:/Documents and Settings/Owner/デスクトップ/r"
> # 作業スペースの保存
> save.image("C:??Documents and Settings??Owner??デスクトップ??r??.RData")
> #ここで、いったん終了
> q()

デスクトップ??r??.RData をダブルクリックすると、なんと、R 2.4.0 が立ち上がる!

> sessionInfo()
R version 2.4.0 (2006-10-03) 
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] "methods"   "stats"     "graphics"  "grDevices" "utils"     "datasets" 
[7] "base"
  • んと. explorerを起動して「ツール:フォルダオプション」の「ファイルの種別」の中からRdataを探して目的のRgui.exeに関連付けして下さい. インストール時に関連付けをしなかっただけでしょう. -- なかま 2007-05-09 (水) 17:05:05
  • お返事ありがとうございました。インストール時の「Rを.RDataに関連づけ」はちゃんとしてあった(というか、変えなかった)のですけど。再度インストールしてもう一度やっても同じでしたが、あちこちいじくって(どこをどういじったかわからなくなるほど)いたら、できるようになりました。しかし、.RDataをクリックしてR2.5.0が立ち上がるようにはなったのですが、起動後、ディレクトリを変更して(作業スペースを保存するを直接選ばず)、q()で終了するときに作業スペースを保存するにしても、ディレクトリの変更は保存されていません(もちろんそれ以外のいろいろな変更も)。 -- 2007-05-09 (水) 17:37:24
  • q()で終了したとき、イメージはそのときのワーキングディレクトリに保存されます。次回起動時は、デフォルトのワーキングディレクトリで起動します。そもそもディレクトリの変更は保存されません。 -- takahashi? 2007-05-09 (水) 21:52:54
  • q()で終了したとき イメージは ワーキングディレクトリの .RData に保存されるんですね。その .RData には,ディレクトリの変更は保存されていないのでしょうか?明示的に「作業スペースの保存」をして,q() では,保存をしないときは,.RData には,ワーキングディレクトリは保存されるのでしょうかされないのでしょうか。前の場合と違った結果になるんでしょうか。q() で作業スペースを保存する方を選択して,その後,そのワーキングディレクトリにある,.RData をクリックしても作業ディレクトリの変更が保存されていないのです。(わかりにくいですか)
    単に作業スペースの変更の保存だけでなく,関数の定義なども保存されないのです。
    私の Windows が,おかしいのかな。 -- 2007-05-10 (木) 00:07:47
  • .RDataにはワーキングディレクトリは保存されません。.RDataがあるディレクトリがワーキングディレクトリになります(.RDataをクリックしてRを起動した場合)。たとえば、保存した.RDataをどっか別のディレクトリにコピーしてクリックしてRを起動したら、ワーキングディレクトリはコピー先のディレクトリになるわけです。関数の定義などはもちろんその中に保存されています。あなたのwindowsはおかしくないと思いますよ。 -- takahashi? 2007-05-10 (木) 00:33:57
  • コメントありがとうございました。以下のようにまとめて宜しいでしょうか。
    (1) 現在の作業ディレクトリに .RData がない場合には,q() で「作業スペースを保存」するを選んでも,新たに .RData が作られることはない。
    (2) .RData がないときには,明示的に「作業スペースの保存」を選ぶことにより,新たに .RData が作成され,同時に関数定義なども保存される。
    (3) 作業ディレクトリは .RData には保存されない。
    (4) .RData をダブルクリックして R を立ち上げると,立ち上げに使った .RData のあるディレクトリが作業ディレクトリになる。
    (5) 既に .RData がある場合には,q() で「作業スペースを保存」を選ぶと,それまでの関数定義などは,作業ディレクトリにある .RData に保存される。 -- 2007-05-10 (木) 16:21:28
  • (1)は違います。(2)は(1)が違うので、別に明示的に作らなくてもいいです。(3)、(4)はそうです。(5)は既に.RDataがなくても新たに作られて保存されます。ひとつの例外は、作業ディレクトリに.RDataがなく、かつ保存するべきオブジェクトが何もないときで、この時は新しく.RDataは作られません。 -- takahashi? 2007-05-10 (木) 19:24:23
  • ああ,「一つの例外」に陥っていたのですね。 -- 2007-05-10 (木) 19:34:17

hist で「'breaks' の数が無効です」とは?

norion? (2007-05-06 (日) 21:50:20)

hist(5,x)というカキコミをして、(xにはあらかじめ数値列がつくってある)

'breaks' の数が無効です というエラーメッセージがでます。

5のところを500とかに変えても同じです。データの数は全部で2000足らずです。
よろしくお願いします。

  • ?hist してみましょう。正しい引数を正しい順で与えていますか? -- 2007-05-06 (日) 22:00:13
  • ?histは見ましたが、breaksというものの意味は(英語が得意でない事もあり)良くわかりませんでした。しかしデータは1つ目に入れなくてはいけないことがわかったので hist(x) と入れてみたところ 'x' は数値でなければなりません というエラーメッセージが出ます。作ってあるxの中のデータは数値のみなんですが・・・。何がいけないのか良くわかりません。
    • norion? 2007-05-06 (日) 22:17:27
  • hist(5,x) というようにしても,仰るようなエラーにはなりませんね。x にはどのようなデータが入っているのか。といっても,たとえば2000個もあるxを列挙されても困りますが。エラーを再現できる最小限のデータと共に質問なさるようにお勧めします。 -- 2007-05-06 (日) 22:24:35
  • 1行目に文字があったことに気づかなかっただけで、これを取り除いたら出来ました。ありがとうございました。 -- norion? 2007-05-06 (日) 22:24:41
  • breaks とはデータ範囲をどのように分割(簡単にいえば棒とその幅を幾つにするか)するかを指定する引数です。 -- 2007-05-07 (月) 09:12:26

複数ウインドウを開いて、描き出すウインドウを選択したい

ななし? (2007-05-05 (土) 13:35:36)

x11() #ウインドウ一つめ
plot(1:10)
x11() #ウインドウ二つめ
plot(2:11)

このとき一つ目のウインドウにグラフを追記したいのですが可能でしょうか?
ひとつのウインドウを分割する場合、split.screenで分割して、screenで描き出す場所を指定するというのがありますがそれのウインドウ版のようなのを探しています。
御存じの方教えてください。
R version 2.4.0
os linux-gnu vine/linux4.0

  • こんな感じの回答がお望みでしょうか。 -- aa? 2007-05-05 (土) 16:15:35
    x11() #ウインドウ一つめ
    plot(1:10, main="1つめ")
    x11() #ウインドウ二つめ
    plot(2:11, main="2つめ")
    dev.set(dev.prev()) ##1つめに移動
    lines(1:10) #1つめに追記
    dev.set(dev.next()) ##2つめに移動
    lines(2:11) #2つめに追記
    graphics.off()
  • まさしくそれです。ありがとうございました。 -- ななし? 2007-05-05 (土) 17:33:17

garchFit で「引数 "formula" がありません」というエラー

たく? (2007-05-01 (火) 12:46:50)

gaarchFitで解析したいものがあり、試しにUKgasでやってみたらこういうエラーが出ました・・・。
どうすれば使えるようになるのでしょうか。 バージョンは2.4.1、OSはXPです。

> data(UKgas)
> library(fSeries)
> UKg.d<-diff(UKgas)
> UKg.m <-garchFit(formula.mean =~arma(2,0),formula.var=~garch(1, 1),series = UKg.d)
以下にエラー.modelSeries(fake = FALSE, lhs = TRUE) : 
       引数 "formula" がありませんし、省略時既定値もありません
  • その使い方はどこから探し出してきたのか?UKg.m <-garchFit( formula = ~arma(2,0)+aparch(1,1),data = UKg.d) とかじゃないの? -- 2007-05-01 (火) 15:14:36
  • http://www1.doshisha.ac.jp/~mjin/R/0606_35.pdf このサイトを参考にしました。ちなみにおっしゃるようなやり方でやっても同じエラーが発生しました。-- たく? 2007-05-01 (火) 15:42:00
  • どこかのpdfは見るけど,オンラインヘルプは見ていないのですか?
    「おっしゃるようなやり方でやっても同じエラー」おかしいなぁ。formulaを与えたんだから,formula がないなんて,言われるはずないんだけど。Rのバージョンを上げてみたらいかが?当方では,以下のように途中エラーは出るが一応結果は出ているみたいだけど(初期値かなんかが必要なんでしょうかね) -- 2007-05-01 (火) 15:50:25
    > UKg.d<-diff(UKgas)
    > UKg.m <-garchFit( formula = ~arma(2,0)+aparch(1,1),data = UKg.d) 
    [1] "arma(2, 0)"   "aparch(1, 1)"
    
    Series Initialization:
     ARMA model:                arma
     Formula mean:              ~ arma(2, 0)
     GARCH model:               aparch
     Formula var:               ~ aparch(1, 1)
     ARMA Order:                2 0
     Max ARMA Order:            2
     途中省略
    Hessian Matrix:
                    mu         ar1         ar2        omega
    mu       0.4038139  -14.661207   -4.759985  0.005485700
    ar1    -14.6612069 3721.675424  502.178082 -0.430102036
    ar2     -4.7599849  502.178082 4486.891725 -1.201257840
     途中省略
    gamma1  -9.5579802   57.61092815   -4.6013946   -8.7055680
    beta1   85.2086937   -4.60139464  240.7460409 -101.1148148
    delta  -22.7605912   -8.70556795 -101.1148148   58.0273266
    
    --- END OF TRACE ---
    
    Warning message:
     計算結果が NaN になりました in: sqrt(diag(fit$cvar)) 
    > summary(UKg.m)
    
    Title:
     GARCH Modelling 
    
    Call:
     garchFit(formula = ~arma(2, 0) + aparch(1, 1), data = UKg.d) 
     途中省略
    Information Criterion Statistics:
          AIC       BIC       SIC      HQIC 
    -9.893399 -9.693561 -9.903576 -9.812388 
    
    Description:
     Tue May  1 15:54:32 2007 by user:
  • 上のをコピペしてみたらできました・・・。さっきはどこかタイプミスしてたのかもしれません。お手数おかけして申し訳ないです。ありがとうございました。 -- たく? 2007-05-01 (火) 16:19:05

INDSCAL(個人差MDS,3元データのMDS)の関数

TH? (2007-04-29 (日) 08:13:46)

私が探した範囲ではRでは,多次元尺度法を行う関数として,以下の関数が用意されています.

 計量的多次元尺度法のための関数
      パッケージ'stats'の関数 cmdscale
 非計量多次元尺度法のための関数
      パッケージ'MASS'の関数 isoMSD ,関数 sammon 
      パッケージ'vegan'の関数 metaMSD 

しかし個人差を扱う,3元データ用のMDSモデルのINDSCALを実行する関数がみあたりません.MDSのなかでは有名な分析方法なので,私の探し方が悪いようにも思うのですが,どなたかINDSCALを実行する関数をご存知の方がいらっしぃませんか.

  • SensoMineR に indscal : Construct the Indscal model for Napping data というのがあるようですが,あなたの求めているものかどうか? -- 2007-04-29 (日) 08:30:56
  • 私の知識不足で,Napping data typeというのがいまひとつ分かりませんが -- TH? 2007-04-29 (日) 11:49:15
  • 上切れました(失礼)私の知識不足で,Napping data typeというのがいまひとつ分かりませんが参照に以下の基本文献がありましたので,ずばりこれです. Carroll, J.D. & J.J. Chang (1970). Analysis of individual differences in multidimensional scaling via an N-way generalization of "Eckart-Young" decomposition. _Psychometrika_, 35, 283-319. すばやいレスに心より感謝申し上げます.help.search("indscal")を試してみなかった自分を反省! -- TH? 2007-04-29 (日) 11:51:25
  • 投稿は編集できるんですけどねえ。下のR/qtl も,自分でちゃんとやれば良かったんだけど -- 2007-04-30 (月) 00:13:36
  • Rと多次元尺度法 -- 2007-04-30 (月) 04:42:52
  • 便乗ですみません。MDSのALSCALって、Rでサポートされているのでしょうか? -- Others? 2007-05-07 (月) 15:22:55
  • こういう時は RsiteSearch? でキーワード検索します。結果は 検索結果 で、どうもないようですね。既存の関数のオプションで代用できる可能性はありますが。 -- 2007-05-07 (月) 15:36:57
  • ALSCAL MDS R でググって,出てくるものを調べれば,http://www.cuddyvalley.org/psychoR/ にあるのではないですか。http://www.cuddyvalley.org/psychoR/alscal/alscal.R その解説は http://www.cuddyvalley.org/psychoR/alscal/paper/alscal.pdf -- 2007-05-07 (月) 16:04:55

R/qtl が load できません

saoki? (2007-04-27 (金) 23:29:18)

R/qtlがloadできません。初めてRを使う初心者です。どなたかご教授いただけましたら幸いです。
Mac OSX バージョン10.4.7のintel MacにR-2.50をインストールしqtl_1.05-2.tarでパッケージqtlを出してlibrary中に入れ、library(qtl)とtypeしたところ以下のように出てきます。

> library(qtl)
Error in dyn.load(x, as.logical(local), as.logical(now)) : 
共有ライブラリ 
'/Library/Frameworks/R.framework/Versions/2.5/Resources/library
/qtl/libs/i386/qtl.so' を読み込めません
  dlopen(/Library/Frameworks/R.framework/Versions/2.5/Resources/library/qtl/libs/ i386/qtl.so, 6): 
Library not loaded: /Library/Frameworks/R.framework/Versions/2.4/ Resources/lib/libRlapack.dylib
  Referenced from: /Library/Frameworks/R.framework/Versions/2.5/Resources/library/qtl/libs/i386/qtl.so
 Reason: image not found
 以下にエラーlibrary(qtl) :  .First.lib は 'qtl' に対して失敗しました

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

  • intel Mac が今ここにないので確かめられないのですが,iBook G4などで(その他のマシンでも同じ方法でやっていますが)いつもやっている方法ですが,「パッケージとデータ」-->「パッケージインストーラ」経由でやっても同じですか? -- 2007-04-27 (金) 23:45:57
  • ありがととうございます。iBook G4があったので(OSXバージョン10.4.6)おっしゃる通りやってみたところ同じコメントが出て参りました。 -- saoki? 2007-04-27 (金) 23:52:44
  • 以下で解決しました。現在のR/qtl (qtl_1.05-2.tar)はR-2.4.0で初めて動く(R-2.4.1や2.5.0では動かない)ことを確かめました。皆様ありがとうございます。 -- saoki? 2007-04-28 (土) 11:22:14
  • おかしいですね。少なくとも iBookG4 の R2.5.0 で qtl_1.05-2 は動きますけど?以下参照 -- 2007-04-28 (土) 14:11:36
    > sessionInfo()
    R version 2.5.0 Patched (2007-04-26 r41343) 
    powerpc-apple-darwin8.9.1 
    
    locale:
    ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8
    
    attached base packages:
    [1] "stats"     "graphics"  "grDevices" "utils"     "datasets" 
    [6] "methods"   "base"     
    
    other attached packages:
         qtl 
    "1.05-2" 
    
    > library(qtl)
    > example(bayesint)
    
    baysnt> data(hyper)
    
    baysnt> ## Don't show: 
    baysnt> hyper <- subset(hyper, chr=c(1,4))
    
    baysnt> ## End Don't show
    baysnt> hyper <- calc.genoprob(hyper, step=0.5)
    
    baysnt> out <- scanone(hyper, method="hk")
    
    baysnt> bayesint(out, chr=1)
               chr  pos      lod
    c1.loc33     1 36.3 2.454921
    c1.loc44.5   1 47.8 3.562391
    c1.loc80.5   1 83.8 2.511073
    以下略
  • 再度行ったところ前述のエラーがやはり出て参りました。sessionInfo()では教えていただいた参照と違いother attached packages: qtl "1.0.5-2"の所が出て参りませんでした。libraryにqtlを入れただけなのが問題なのでしょうか。なおR-2.4.0ではlibrary(qtl)の後example(bayesint)が動くことを確かめました。 -- saoki? 2007-04-28 (土) 21:54:19
  • おかしいですね. R2.6.0でも試しましたが,動きましたよ。「パッケージとデータ」-->「パッケージインストーラ」経由でやったのですか?場合によっては何らかの状況下で開発環境がダメージを受けていてちゃんとインストールできていないのかも知れませんが。まあ,ちゃんと動かない方が異常であることは確かです。R2.4.0でしか動かないというのは,いつの時点の,どの文書によるものでしょうか。 -- 2007-04-28 (土) 22:24:20
    R version 2.6.0 Under development (unstable) (2007-04-26 r41343)
    Copyright (C) 2007 The R Foundation for Statistical Computing
    ISBN 3-900051-07-0
       中略
    'q()'と入力すればRを終了します。
    
    [Workspace restored from /Users/***/***/***/.RData]
    
    > library(qtl)
    > library(help=qtl)
    > 要するに,ここまでで何のエラーもないし,この後もちゃんと動く
  • intel Mac でも,ok でしたよ。 -- 2007-04-29 (日) 07:55:39
  • 手動で入れた、[/Library/Frameworks/R.framework/Versions/2.5/Resources/library/qtl]を手動で消す。[ユーザー名/ライブラリ/R/library/]もしくは[ユーザー名/ライブラリ/R/2.5/library/]、にqtlがあれば手動で消す。再度「パッケージとデータ」-->「パッケージインストーラ」経由で(パッケージの依存関係を解決)をチェックしてインストール。 -- okinawa 2007-04-29 (日) 08:18:18

サンプルについて

goto? (2007-04-25 (水) 17:57:25)

どなたか教えてください。
以下のサンプルでは要素(球体)やノードを大きくしたり加工しています。
どのようにしたらこのようになるのでしょうか。
gplot関数のパラメータを変更するのでしょうか。

http://erzuli.ss.uci.edu/R.stuff/sna/screenshots/sna_sample3.png
sna_sample3.png
  • タイトルが不適切。引用で画像のURLを書いたら表示されるとはいうものの,大きすぎだし,どこから引用したのか(出典)が表面には現れないし。不適切なんだと思いますが。それに,説明もまるっきりないし。回答を望んでいないとしか思えない。 -- 2007-04-25 (水) 22:43:41
  • 説明不足の件は大変失礼致しました。snaパッケージを利用したものです。サンプルは ttp://erzuli.ss.uci.edu/R.stuff/ ←にありました。availableというところにあります。他はどのような情報が必要でしょうか。3Dでランダムの要素とノードを矢印で描くことはできますが、要素やノードの大きさを変えたり、もしくは画像全体の縮尺を変更したりする方法がわかりません。Rを初めて3日目なのでRに関しての知識は薄いです。上記の画像を見てこのような関数を利用しているのではと仮説でもいただける方よろしくお願い致します。 -- goto? 2007-04-26 (木) 10:14:14
  • ?gplotとして、ヘルプファイルを熟読する。vertex.cex: expansion factor for vertices; may be given as a vector, if vertices are to be of different sizes. あたりの引数を弄るとか。 -- takahashi? 2007-04-26 (木) 11:29:42
  • 引用される図があまりにも大きすぎて不都合なので,ダウンロード後,縮小して,ここへ再アップ。 -- 2007-04-26 (木) 14:40:02
  • 図の縮小ありがとうございます。gplotのパラメータを弄ってみたのですが、3Dの映像にできませんね。他の関数を利用しているのでしょうか。。 -- goto? 2007-04-26 (木) 16:30:55
  • 自己レスです。gplot3という関数が用意されていました。解決です。 -- goto? 2007-04-26 (木) 17:13:02

Rでお絵かきツール

ロトリング? (2007-04-24 (火) 21:24:31)

 Rで出力したグラフに、ポリゴンやテキストをマウスで描画できるお絵かきツールのようなものはないでしょうか?

  • > ポリゴンやテキストをマウスで描画
    その程度のことなら,テキストを指定して,描画位置をマウスクリックで指定してとか,マウスクリックした点を順次線で結ぶとかで実現可能でしょう。どの程度の性能をお望みなのかな? -- 2007-04-24 (火) 21:34:10
  • ご返答ありがとうございます。返事が遅れて申し訳ございません。利用例としては、散布図を描いて、マウスで対話的にポリゴンを描いた後、ポリゴン内のポイントを列挙できないものかと考えております。 -- ロトリング? 2007-04-26 (木) 14:30:15
  • ポリゴンにもよるけど,ある点がポリゴンの内か外かを判定するのはなかなか難しかったような(凸ポリゴンの場合は比較的簡単)。 -- 2007-04-26 (木) 15:18:03
  • splancsパッケージにinpip()を使えば楽だと思います。マルチポリゴンには使えないのでGIS的使用は難しいですが、この場合には十分な機能と思います。 -- 2007-04-26 (木) 20:52:05
  • みなさん、情報ありがとうございます。ところで、最初の質問のどのような関数を使ってマウスをクリックしながらポリゴンを描けばいいのでしょうか? -- ロトリング? 2007-04-27 (金) 16:19:41
  • 単にこういうことを知りたかったとか? -- 2007-04-27 (金) 22:38:31
    poly.xy <- locator(type = "l")

95%信頼楕円

suz? (2007-04-21 (土) 11:23:24)

95%信頼区間楕円を描く方法として以下の操作は正しいでしょうか?

x = rnorm(100, sd=1)
y = rnorm(100, sd=3)
plot(x, y)
library(ellipse)
r = cor(x, y)
stdev = c(sd(x), sd(y))
centre = c(mean(x), mean(y))
polygon(ellipse(r, scale=stdev, centre=centre, level=0.95))
  • (元記事質問者と二番目のコメント投稿者は同じだと仮定して)次の例を見れば想像がつくのでは。 -- 2007-04-22 (日) 00:28:18
    > V <- rbind(c(var(x), var(x,y)),c(var(x,y),var(y)))
    > a <- ellipse(V, centre=centre, level=0.95)
    > b <- ellipse(r, scale=stdev, centre=centre, level=0.95)
    > identical(a,b)
    [1] TRUE
  • 注意:コメントがすでにあった場合、元記事を書き直すのは違反です。コメントとして再質問してください。 初級者用Q&Aということで注釈すれば、2変量 x,y に対し cor(x,y), var(x,y) は相関・共分散(スカラー値)を計算しますが、相関行列・共分散行列を計算するわけではありません。共分散行列(上の例の V)からは cor(x,y), sd(x), sd(y) は計算できますから、scale 引数は不要です。しかし相関係数だけでは sd(x),sd(y) の値は決められないので scale 引数が必要になるわけです。-- 2007-04-22 (日) 07:10:32
  • 再質問します。-- 2007-04-22 (日) 12:47:38 95%信頼区間楕円を描く方法として以下の操作は正しいでしょうか?
    set.seed(101)
    x = rnorm(100, sd=1)
    y = rnorm(100, sd=3)
    xy = cbind(x,y)
    library(ellipse)
    stdev = c(sd(x), sd(y))
    centre = c(mean(x), mean(y))
    plot(x,y)
    polygon(ellipse(var(xy),  scale=stdev,  centre=centre, level=0.95)) #1
    polygon(ellipse(var(xy),  scale=c(1,1), centre=centre, level=0.95)) #2
    polygon(ellipse(var(xy),                centre=centre, level=0.95)) #3
    polygon(ellipse(cor(x,y), scale=stdev,  centre=centre, level=0.95)) #4
    polygon(ellipse(cor(x,y), scale=c(1,1), centre=centre, level=0.95)) #5
    polygon(ellipse(cor(x,y),               centre=centre, level=0.95)) #6
    95%信頼区間楕円を描く方法として#1から#6のどれが正しいのでしょうか?

    #3と#4は同じ楕円を描きます。 scaleのデフォルトはc(1,1)ですが、#2と#3は同じ楕円を描きません。(#5と#6は同じ楕円を描きます)

  • まあ,基本的には,その関数のオンラインヘルプをちゃんと読む&exampleを実行してみる。さらには,答えの分かっている例題を解いてみる。さらには,まあ,多数決とはいわないまでも,同じことをやっている色々なサイトで,実際のテストデータを実行してみて,どれが同じ結果を得るのかを見てみるのも判断の基準となるでしょう。それらを総合して,どのような点が疑問であるかを明示して質問すべきでしょう。質問のVersion 2でも,何が知りたいのか・問題と思っているのかはまるっきり明らかではありませんね。質問するにおいても,作法というものがあることをよく理解されることを期待したいと思います。 -- 2007-04-23 (月) 22:04:06
  • > ご指摘ありがとうございます。元記事です
    厳しく言えば,Wiki においては編集は無限の自由度があるのだから,不要な修正をしなかった場合にどのようなスレッドができあがったかを再現することだってできるわけで,その再現工事は,原質問者の義務でもあると思うんですが,そこまで望むのは酷ですか?
    というか,直せば直すほど修正箇所が増えて,どうしようもなくなることは,めにみえているわけですが。だからこそ,へんな修正はしないほうがよいわけで。。。 -- 2007-04-23 (月) 22:32:28
  • ellipseのhelpを読み、exampleを実行してみました。答えの分かっている例題を解いてみたり、実際のテストデータを実行してみたいのですが、同じことをやっているサイトを発見できませんでした。そこで、こちらのサイトに質問させていただいた次第です。#3と#4が正しい95%信頼区間楕円を描いているのかどうか判断することができませんでした。scaleのデフォルトはc(1,1)なのに、#2と#3が同じ楕円を描かないのが理解できませんでした。 -- 2007-04-24 (火) 02:04:14
  • 確率楕円で検索してみたのでしょうか?http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc027/01846.html などもありますね -- 2007-04-24 (火) 10:08:34
  • 「統計学関連なんでもあり」の過去ログ // 教えて下さいDensity Ellipse http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc008/083.html // 等確率楕円の長軸 http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc015/290.html // これは,確率楕円ですか? http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc027/01846.html // 二次元正規分布の確率楕円 http://aoki2.si.gunma-u.ac.jp/lecture/stats-by-excel/vba/html/kakuritudaen.html // 散布図(各種の描画機能付き)散布図を描き,棄却楕円(確率楕円),回帰直線,回帰直線の信頼限界帯,MA,RMA による回帰直線を描く。 http://aoki2.si.gunma-u.ac.jp/R/scatter.html -- 2007-04-24 (火) 10:32:10

mfrowの値とplotmeansの動作

suzu? (2007-04-20 (金) 16:32:43)

グラフィックパラメータpar(mfrow=c(行数, 列数))の行数によってplotmeansの動作が異なります。

library(vegan)
data(BCI)
sp.n = specnumber(BCI)
group = gl(5, 10)
library(gplots)
par(mfrow=c(7,1)) # 行数によってplotmeansの動作が異なる
plotmeans(sp.n ~ group)

において

par(mfrow=c(6,1))

の場合、正常に動作しますが、

par(mfrow=c(7,1))

の場合、group 4のエラーバーが表示されず、以下のWarning messagesが出力されます。

Warning messages:
1: 長さゼロの arrow は角度が不定ですので,スキップされました 
2: 長さゼロの arrow は角度が不定ですので,スキップされました 

エラーバーを正常に表示させるためには、どうすればよろしいでしょうか?

使用環境 R version 2.4.0 (2006-10-03) i386-apple-darwin8.8.1 参考文献 http://cse.naro.affrc.go.jp/takezawa/r-tips/r/68.html

  • 描画領域を大きくする。 -- takahashi? 2007-04-20 (金) 17:22:49
  • アドバイスありがとうございます。「描画領域を大きくする」というのは「par(mfrow=c(行数, 列数))の行数を小さくする」という意味でしょうか?「描画領域」とは、「プロット領域(plot region)」または「作図領域(figure region)」のことでしょうか? http://cse.naro.affrc.go.jp/takezawa/r-tips/r/52.html -- suzu? 2007-04-24 (火) 16:20:57

表形式でないファイルの読み込み

mori? (2007-04-19 (木) 15:03:13)

ファイルからデータを読み込む方法に関して質問があります。

[test.txt]
milk bread
bread butter
beer
milk bread butter
bread butter

上のようなtest.txtからデータを読み込んで下のように作ったリストと同じリストを返す関数がほしいのですがどのようにファイルを読み込めばよいでしょうか。read.table()は列数が同じものしか読み込めないらしくうまく読み込めません。アドバイスお願いします。

> x <- list(c("milk","bread"),
c("bread", "butter"),
c("beer"),
c("milk", "bread", "butter"),
c("bread", "butter"))
  • read.table の fill パラメータを調べましょう。
    スマートではないが,以下のようなものも -- 2007-04-19 (木) 15:55:11
    x <- as.matrix(read.table("test.txt", header=FALSE, fill=TRUE))
    dimnames(x) <- NULL
    apply(x, 1, function(a) a[a != ""])
    こんな別解も
    a <- sapply(readLines("test.txt"), function(s) unlist(strsplit(s, " ")))
    names(a) <- NULL
  • もしくは類似解には -- 2007-04-19 (木) 16:28:59
    > sapply(as.list(readLines("test.txt")), function(i) strsplit(i, split=" ")) 
    [[1]]
    [1] "milk"  "bread"
    
    [[2]]
    [1] "bread"  "butter"
    
    [[3]]
    [1] "beer"
    
    [[4]]
    [1] "milk"   "bread"  "butter"
    
    [[5]]
    [1] "bread"  "butter"
    
    Warning message:
    'test.txt' に関する readLines で不完全な最終行が見つかりました
    注:上の warning は,ファイルの最後に EOF がないために出る。
      ちゃんとしたファイルなら warning は出ない。
  • 類似解ですが,文字数が少ない(^_^;) -- 2007-04-19 (木) 16:56:40
    mapply(strsplit, as.list(readLines("test.txt")), split=" ")
  • readLines()というのを知らなかったので大変勉強になりました。ありがとうございます。 -- mori? 2007-04-20 (金) 09:42:10

garchFit(パッケージfSeries)のエラーについて

ケンゴ? (2007-04-14 (土) 16:05:43)

こんにちは、garchFit(パッケージfSeries)関数を使用し時系列データの計算を
おこないたいと思い以下のような処理を致しました。
結果「エラー:!missing(data) is not TRUE」と表示されうまく計算できません。
どこか間違っているのでしょうか?教えて下さい。
使っているのはR2.4.1 OSはXPです。

> library(fSeries)
> UKgas
       Qtr1   Qtr2   Qtr3   Qtr4
1960  160.1  129.7   84.8  120.1
1961  160.1  124.9   84.8  116.9
1962  169.7  140.9   89.7  123.3
1963  187.3  144.1   92.9  120.1
1964  176.1  147.3   89.7  123.3
1965  185.7  155.3   99.3  131.3
1966  200.1  161.7  102.5  136.1
1967  204.9  176.1  112.1  140.9
1968  227.3  195.3  115.3  142.5
1969  244.9  214.5  118.5  153.7
1970  244.9  216.1  188.9  142.5
1971  301.0  196.9  136.1  267.3
1972  317.0  230.5  152.1  336.2
1973  371.4  240.1  158.5  355.4
1974  449.9  286.6  179.3  403.4
1975  491.5  321.8  177.7  409.8
1976  593.9  329.8  176.1  483.5
1977  584.3  395.4  187.3  485.1
1978  669.2  421.0  216.1  509.1
1979  827.7  467.5  209.7  542.7
1980  840.5  414.6  217.7  670.8
1981  848.5  437.0  209.7  701.2
1982  925.3  443.4  214.5  683.6
1983  917.3  515.5  224.1  694.8
1984  989.4  477.1  233.7  730.0
1985 1087.0  534.7  281.8  787.6
1986 1163.9  613.1  347.4  782.8
> UKg.d<-diff(UKgas)
> UKg.m<-garchFit(formula=~arma(2,0)+garch(1,1),series=UKg.d)
エラー:!missing(data) is not TRUE

※条件付き平均(mean)=arma(2,0)
※条件付き分散(variance)=garch(1,1)

  • 使ったことがないから適当なコメントですが、引数 formula の与え方は適切ですか。 -- 2007-04-14 (土) 17:36:24
    ?garchFit によれば:
        garchFit(formula.mean = ~arma(0, 0), formula.var = ~garch(1, 1), 
            series = x, init.rec = c("mci", "uev"), delta = 2, skew = 1, shape = 4,
            cond.dist = c("dnorm", "dsnorm", "dged", "dsged", "dstd", "dsstd"), 
            include.mean = TRUE, include.delta = NULL, include.skew = NULL,
            include.shape = NULL, leverage = NULL, trace = TRUE,  
            algorithm = c("sqp", "nlminb", "lbfgsb", "nlminb+nm", "lbfgsb+nm"), 
            control = list(), title = NULL, description = NULL, ...)
  • アドバイスありがとうございます。 -- ケンゴ? 2007-04-14 (土) 19:37:31
  • で、アドバイスは役に立ったんですか、的はずれだったんですか。 -- 2007-04-14 (土) 22:57:05
  • どうも,質問の出所というか参考にしているところは同じらしい。正しくは,
    UKg.m<-garchFit(formula=~arma(2,0)+garch(1,1),data=UKg.d)
    のようですよ(?) -- 2007-05-01 (火) 18:10:52

Rにインストールは必要ない?!

豊田秀樹? (2007-04-14 (土) 07:03:17)

信じられないようなことですが,ちゃんとインストールした
Rのフォルダを単純にUSBフラッシュメモリにコピーして
おくと,そのUSBフラッシュメモリを持ち歩くだけで,Rをインストール
していない計算機でRを使えます.(そのフォルダは,全ての
パッケージがインストールされており1ギガくらいの容量です)
更にそのUSBフラッシュメモリのRののフォルダを単純に
別の計算機にコピーするとその計算機でもRが使えるようになります.
ここで「コピーで済むなら何のためのインストールだったん
だろう?」という疑問が沸いてきました.きっとこの方法で
コピーを繰り返すと「何か困ったことがそのうち起きる」あるいは
「何か使えない高度な機能がある」としか思えません.
文科系の学部なのでインストールはハードルが高い学生もいますが,
コピーだけならだれでもでき,各学生に任せられます.
とても便利なので,つい学生にも薦めて見たくなっていますが,
あとで変なことが起きると困るので,どなたか「インストール
しないと,こういう点で後で困るよ」というアドバイス
がありましたら教えてください.それともインストールなんて
必要ないのでしょうか?そんなことないですよね.

  • すみません環境・バージョンを書き忘れました.OSは Windows XP と 2000 で確認しましたバージョンはR2.2.0 と最新のR2.4.1です. -- 豊田秀樹? 2007-04-14 (土) 07:17:09
  • 面白いですね。いわれてみればなるほどな、ですね。最近もVistaで管理者権限を持たない人がパッケージを使うために、パッケージをUSBにインストールしているという話を聞きました。調べてみたら、既に Windows FAQ に記事がありました。それによるとCD ROMでも良い(当然DVD ROM でもよいのでしょうね)そうです。いくつか注意すべき点はあるようですが、基本的にはよさそうです。Vista では?だれか検証おねがいします。関連する r-help 記事。考えられる欠点は、アクセス速度が遅くなりがち、他のRが利用するプログラム、ライブラリ等が使えなくなる可能性がある、等でしょうか。 (ついでですが、Windows R FAQ に日本語を使う方法という項目があるのを発見しました。)-- 間瀬茂 2007-04-14 (土) 09:19:28
    2.6 RをCDやUSBから起動できますか?
    
    注意すればできます。基本的なRインストールは再配置可能(どこに置いても起動可能)
    ですから、インストールしたRのイメージをハードディスクにコピーしたり、直接
    フラッシュメモリUSBドライブ等のリムーバブル記憶装置にインストールすることが
    できます。(もしパッケージを私的なライブラリ(フォルダー)にインストールして
    いれば、それらの絶対パスがHTMLパッケージリストに記録されます。)
    
    Rの実行には書き込み可能な一時ディレクトリとホームディレクトリへのアクセスが必要
    で、もしそうしたものがなければ現在のディレクトリが使われます。これは、適正に
    編成されたNTベースのWIndowsのバージョンでは何の問題も起こさないはずですが、
    もし問題があれば書き込み可能なフォルダーへのショートカット無しにはRを実行する
    ことはできないかも知れません。
    
  • R は, 全てのOSで再配置可能です. 自分のトップ(R_HOME)があればそれに従い実行可能です. Windowsの場合はDCOMの為にレジストリにもR_HOMEを書き込みするだけです. 10Gぐらいあれば(リムーバブルHDしか無いか?),CRANとBioConductorのパッケージを全部入れられますね. -- なかま 2007-04-14 (土) 11:18:26
  • 早速のアドバイスを,どうも有難うございました.了解いたしました.早速1つの実験をしてみました.もともと正式にハードディスクにインストールしたRの入っている計算機(WindowsXP,R2.4.1)に,Rのフォルダを有するUSBフラッシュメモリを更に挿し,USBフラッシュメモリからRを起動して,筑波からパッケージをインストールしてみたところ,ハードディスクではなくUSBフラッシュメモリのRにインストールされました. -- 豊田秀樹? 2007-04-14 (土) 19:21:03
  • Rguiの引数に, R_LIBS=c:/temp 等を加えれば, パッケージのインストール先を更に制御出来ます. *nix系だと, `R_LIBS=/tmp R' 等としますが, Windows版だと引数を環境変数の変わりに代用するようになっているので Rgui なら, `Rgui.exe R_LIBS=c:/temp' 等(ショートカットのプロパティとかでも設定出来ます)とする事が出来ます. 環境変数R_DEFAULT_PACKAGES(R_LIBS同様の方法)なんかでRcmdrなど加えてしまえば, ポータブルメニュー付きR解析キットにも出来るかと.(USB1個で出張解析...嬉しいのか悲しいのか私はプログラム屋なので存じませんが...) -- なかま 2007-04-15 (日) 22:39:57
  • 貴重な情報をありがとうございます.勉強になりました.計算機の苦手な文科系の学生にも手間をかけずに持たせることができることがわかりました. -- 豊田秀樹? 2007-04-17 (火) 09:18:26
  • なかまさんが言っているような「R解析キット」が分野別にあるとuserの裾野が広がっていいかもしれませんね。インストールやパッケージの追加が結構面倒だったりするので。(最近のパッケージの増加についていけてない・・・) -- okinawa 2007-04-17 (火) 11:50:08
  • うーん, "「R解析キット」が分野別=>パッケージの追加が面倒" と言う考えはたぶん間違っています. 裾野を広げる意味であたらしいバイナリ配布形態を考えるなら, それらのソースバージョンを保持しないといけないでしょう. Rは「R(統計)を学んだり改変したりする自由を担保している」フリーソフトウェアに過ぎません. 作るとすれば, 「R解析キット作成キット」で自分のコピーを依存関係を調べて自由に作れるものなら問題ないでしょう(それらは単なる出力だし, 作成はユーザの自由になります). こっそり(たとえば身内だけ)そういった物を配るのは自由です. これが自由を担保する仕組の上に成り立っている事を忘れてはなりません. -- なかま 2007-04-18 (水) 11:46:05
  • そうそう、私が言いたかったのは「R解析キット分野別作成キット」ですね。(言葉の内容的に誤解させてたらご免なさい) -- okinawa 2007-04-19 (木) 09:39:39
  • よいことを聞いたと、うれしくなり、とりあえずデスクトップの任意のディレクトリにインストールしてみました。しかし、ライブラリが見つけられないということです。そして、その状態ではほとんど何もできません。sessionInfo() することもできません。R2.5.0 を、Windows XP にインストールしているところです。でてくるエラーメッセージは以下の通りです。
    Warning messages:
    1:  'lib.loc' 中に如何なるライブラリー木も見つかりません 
      in: library(package, lib.loc = lib.loc, character.only = TRUE, logical = TRUE,  
    2: package "methods" in options("defaultPackages") was not found 
     起動準備中です ー 12 件の警告がありました (警告を見るには warnings() を使って下さい) 
    何かが足りないか(レジストリとかをいじらないといけないようですが、やりかたがわかりません)、設定を間違えているかでしょうけど、それがわかりません。 -- 2007-05-09 (水) 14:04:21
  • .Rdataをダブルクリックしたまたは, コピーが不完全じゃないのかな? Rgui.exeをダブルクリックしても同じ? getwd();.libPaths() あたりはどうですか? -- なかま 2007-05-09 (水) 17:02:11
  • お返事ありがとうございます。binディレクトリのRgui.exeをダブルクリックしています。コピーが不完全とは思われません。なぜならば、R-2.5.0ディレクトリを含むRをもとのProgram Files ディレクトリに戻すと、ちゃんと動きますから。 -- 2007-05-09 (水) 17:41:19
  • 私の思い付く限りは, せいぜいみょーにフォルダ名の長いところに移したとか, profileで凄いことをしない限り, ちゃんと動くと思うんですが.? -- なかま 2007-05-09 (水) 23:46:03
  • デスクトップの r というディレクトリに移しただけで,profile では何もしていません(途中に日本語を含むパス名が不当だとか?) -- 2007-05-09 (水) 23:58:26
  • 日本語ディレクトリはなるべくなら避けてほしいなぁ(希望). -- なかま 2007-05-11 (金) 10:56:50
  • 避けて欲しいというか,それが原因のようですね。Cドライブのトップにディレクトリを作って,その名前に日本語を使ったりやめたりの実験で,日本語が使われるとだめになることが明らかになりました。一番わかりやすい場所がデスクトップということもあるし,日本語のディレクトリ名を許してくれるようになると良いなあ(希望)。 -- 2007-05-17 (木) 11:23:21

Windows版コンソールで使えるフォント

mori? (2007-04-12 (木) 14:32:35)

Windows版のR2.4.1なのですが、GUIプリファレンスで指定できる日本語フォントがCourierとMSゴシックとMS明朝しかありません。他の日本語フォントを指定することはできないのでしょうか?

  • たとえば,Terminal とか FixedFont? とか,これで,日本語が表示できるのかと思われるようなものを選んで反映をクリックすると,ちゃんと字体が変わるものが何種類かあるようですが?あなたがインストールしているフォント(たとえば草書体とかヒラギノ明朝など)が他にあって,それを使うのにはどうしたらよいかという質問なんでしょうか? -- 2007-04-12 (木) 16:01:11
  • あれはリストボックスじゃなくって, コンボボックスだったような記憶があるんですが? -- なかま 2007-04-12 (木) 18:35:43
  • コンボボックスに直接指定すればダイジョブ。これは「HGS行書体」。なんだか趣きのあるRですね。 -- takahashi? 2007-04-13 (金) 08:08:05
    ss2.png
  • アドバイスありがとうございます。リストにないフォントでも入力したら使えるようになりました。コンボボックスだとは気がつきませんでした・・・ -- mori? 2007-04-16 (月) 08:11:28

2変量散布図とヒストグラム図

teru? (2007-04-11 (水) 19:56:49)

xとyの2変量散布図と,xおよびyのヒストグラムを書きたいのですが,
Rではどのようにすればいいのでしょうか?

  **|      oo
  ***|     ooo
  ****|   oooo
  ***|   oo
   **|  oo
   *| o
    -------------
     *********
      *******
      ****
       **
  • 出来合いのものはないと思われるので,自分で関数を書くしかないのでは?自分で書くとすれば,おのおのの流儀があるでしょう。
    pairs 関数の example を(実行して)見れば,望むものにちょっと近いものがあるのに気づくでしょう。それを直せばなんとかなるかなぁと。 -- 2007-04-11 (水) 21:21:05
    pairs.png
  • ヒストグラムではないですが、chplotパッケージがレイアウト的には近いですかねぇ。 -- 2007-04-11 (水) 22:15:07
  • グラフィックス参考実例集:箱型図の「箱型図を軸の装飾に使う」という記事あたりが参考になるかと.他にも役に立つスクリプトはいろいろ探せると思うのでひととおり見回ってみると良いでしょう -- 2007-04-11 (水) 22:24:21
  • ヒストグラムと密度の推定の「二変量データのヒストグラム:marginal」はご覧になりましたでしょうか?似ていると思います。 -- 2007-04-12 (木) 09:47:47
  • par()でmfrowでやるより,layout(matrix(c(1,1,4,2,2,3,2,2,3),3,3))とでもして,#3の描画をやめると,バランスがいいんじゃないかと思います。 -- 中澤? 2007-04-12 (木) 10:33:18
  • 大変助かりました.みなさん,どうもありがとうございます -- teru? 2007-04-12 (木) 17:39:13

ハッシュ

JR? (2007-04-11 (水) 14:24:09)

Rにはperlなどのスクリプト言語でいうハッシュに相当するデータ型はないのでしょうか。例えば、

a => 1
b => 2
…
z => 26

というハッシュと

v=c("a","b",…"z")

というベクトルがあったとき、vの各要素をハッシュのキーとして、値に置換するような処理をしたいです。データフレームとsubset関数を使えば強引にできるのですが、遅いです。よい方法はありますでしょうか?

  • こういうことでしょうか -- 2007-04-11 (水) 14:48:33
    > a <- 1:26
    > names(a) <- letters
    > a
     a  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  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 
    > a['b']
    b 
    2
  • a[names(a) == "b"] でないと最初に見付かったものしか選ばれません。ちなみにデータフレーム経由で取り出すやり方と比べるとそんなには早くないですね。 -- 2007-04-11 (水) 17:54:31
    > x <- sample(1:100, 1e6, rep=TRUE)
    > y <- sample(letters, 1e6, rep=TRUE)
    > xx <- x
    > names(xx) <- y
    > xxx <- data.frame(x=x, y=y)
    > xx["a"]
     a 
    81  
    > length(xx[names(xx) == 'a'])
    [1] 38482
    > system.time(xx[names(xx)=="a"])
    [1] 0.164 0.000 0.166 0.000 0.000
    > system.time(xxx$x[xxx$y =="a"])
    [1] 0.248 0.000 0.254 0.000 0.000
  • ありがとうございます。ドンピシャなものは無いけど、ベクトルで何とか頑張れば似たようなことはできそうですね。速度や重複の問題を解決するには本物のハッシュ実装するしかなさそうですね…。 -- JR? 2007-04-12 (木) 13:23:28

区切り文字が1文字以上の空白で、#を含むNAに当たる文字が要素に含まれるテーブルの読み込み

TT? (2007-04-07 (土) 17:30:52)

こんにちは、以下のような、テキストファイルに入ったテーブルを読み込もうとしています。区切り文字は1文字以上の空白で、-1.#INDは不定値なのでNAと読ませたいと考えています。

1  2 3 4 5
11 12 13 -1.#IND 15

これを、

tt<-read.table("hoge.txt")

で読み込もうとすると

'2' 行目には,5 個の要素がありません

というエラーが出ます。どうも2行目の4個目の要素の-1.#INDの#をコメントの始まりと解釈してそのあとを読まないようです。そこで、今度は、

tt<-read.table("hoge.txt",sep=" ")

とすると、今度は、1行目の1個目と2個目の要素の間の2個の空白が、2個区切り文字がありその間に、空の要素があると判断されて、1行目は6個の要素があると判断され、

'2' 行目には,6 個の要素がありません

とエラーが出ます。何とか、1文字以上の連続した空白を1個の区切り文字と判断させて、-1.#INDをNAと判断して読み込む方法はありませんでしょうか?

  • hoge というのはどうも。以下のようにしたらいかが?さらなる極意は read.csv のオンラインヘルプを読みましょう。 -- 2007-04-07 (土) 20:11:36
    > x <- read.csv("foo.dat", header=FALSE, sep=" ", na.strings="-1.#IND")
    > x
      V1 V2 V3 V4 V5
    1  1  2  3  4  5
    2 11 12 13 NA 15
  • read.csv は,read.table にシュガー・コートを着せたものですから,read.table でも当然同じことができるんですよ。read.table にあるたくさんのオプションの設定のデフォルトを特定のものにしたものが read.csv, read.csv2, read.delim, read.delim2 なんですから。read.table で,先の read.csv と同じことをするためには,以下のように指定するのです。 -- 2007-04-07 (土) 20:57:56
    > x <- read.table("foo.dat", sep=" ", na.string="-1.#IND", comment.char="")
  • 回答ありがとうございます。read.tableのオンラインマニュアルは読んでいたのですが、read.csvとread.tableとで、sep以外にも既定のオプションが違うのですね。ただ今の場合は一行目の1個目と2個目の要素の間に2個半角の空白が入っているので
    > x <- read.csv("foo.dat", header=FALSE, sep=" ", na.strings="-1.#IND")
    > x
      V1 V2 V3 V4 V5 V6
    1  1 NA  2  3  4  5
    2 11 12 13 NA 15 NA
    となってしまい、意図した結果になりません。もう少しご教授をお願いします。 -- tt? 2007-04-07 (土) 20:47:13
  • 見れども見えず。read help. -- 2007-04-07 (土) 21:15:01
    sep	 the field separator character. Values on each line of the file are 
             separated by this character. If sep = "" (the default for read.table) 
             the separator is “white space”, that is one or more spaces, tabs, 
             newlines or carriage returns.
    ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
    > read.csv("foo.dat", header=FALSE, sep="", na.strings="-1.#IND")
      V1 V2 V3 V4 V5
    1  1  2  3  4  5
    2 11 12 13 NA 15
  • もしNA値に置き換えるべき項目いつも"-1.#IND" ではなく、「単に # を含むもの」ということならば(多少面倒ですが)次のように段階的にやるしかないでしょう。ただしこれは列数が5ということがわかっていることが前提です。 -- 2007-04-07 (土) 21:27:31
    > x0 <- scan("temp.txt", what=character(1)) Read 10 items
    > x0
     [1] "1"       "2"       "3"       "4"       "5"       "11"      "12"     
     [8] "13"      "-1.#IND" "15"     
    > x1 <- gsub(".+#.+", "NA", x0)
    > x1
     [1] "1"  "2"  "3"  "4"  "5"  "11" "12" "13" "NA" "15"
    > x2 <- as.numeric(x1) # 直後のコメントとからもわかるようにここは単に as.numeric(x0) で可
    > x2
     [1]  1  2  3  4  5 11 12 13 NA 15
    > x3 <- matrix(x2, ncol=5, by=TRUE)
    > x3
         [,1] [,2] [,3] [,4] [,5]
    [1,]    1    2    3    4    5
    [2,]   11   12   13   NA   15
    > x4 <- as.data.frame(x3)
    > x4
      V1 V2 V3 V4 V5
    1  1  2  3  4  5
    2 11 12 13 NA 15
  • そこまでやる必要はないのでは?
    要するに,数値を表す文字列以外を NA にするなら,以下のようにすればよい(適用範囲は,より広い)。この場合の warning は無視すればよい(望むべき処理をしてくれたことを報告してくれているに過ぎないのだから)。何行・何列のデータであってもかまわない。 -- 2007-04-07 (土) 21:39:04
    > x <- read.csv("foo.dat", header=FALSE, sep="", as.is=TRUE)
    > y <- as.data.frame(matrix(as.numeric(unlist(x)), nrow(x)))
    Warning message:
    強制変換により NA が生成されました 
    > y
      V1 V2 V3 V4 V5
    1  1  2  3  4  5
    2 11 12 13 NA 15
    無意味な空行はインデントを狂わすのでご注意
  • ありがとうございました。オンラインマニュアル、まさしく見れども見えずで、read.tableの場合、タブも1文字以上の空白文字も、1個の区切り文字と解釈しているのには気づいていたのですが、そこから先で躓きました。今度からもう少し注意します。さらに数値を表す文字列以外をNAとするやり方も大変勉強になりました。ありがとうございました。 -- tt? 2007-04-08 (日) 17:58:33

パッケージのインストール

tyamada? (2007-04-06 (金) 18:13:53)

教えてください。Rのパッケージ(たとえば、vegan)をインストールしたいと思っています。
使っているのはR2.4.1 OSはビスタです。

chooseCRANmirror()
と打って、適当なミラーを選択し(日本の3つとも試しました)
install.packages(c("vegan"))
と打ちますと、はじめは調子よくいっているように見えるのですが、少しすると

URL 'ftp://ftp.u-aizu.ac.jp/pub/lang/R/CRAN/bin/windows/contrib/2.4/vegan_1.8-5.zip' を試しています
ftp data connection made, file length 1176434 bytes
開かれた URL
downloaded 1148Kb

以下にエラーzip.unpack(pkg, tmpDir) : ファイル 'C:/Program Files/R/R-2.4.1/library/file7bb9902/vegan/chtml/vegan.chm' を開くことができません

などというメッセージが返されて、止まってしまいます。

また、ローカルにあるZIPファイルからのインストールでも、同じエラーメッセージが返されてしまいます。
どこが間違っているのでしょうか?教えてください。

  • Vistaではユーザの権限管理がunix系並みに厳しくなっているので、その関係のような感じがします(追試はまだしていませんが)。Rを起動するときに右クリックをして「管理者として実行」すると改善するのでは。 -- 2007-04-06 (金) 22:43:13
  • ありがとうございます。ご教示していただいた方法で、解決いたしました。 -- tyamada? 2007-04-09 (月) 09:10:10

ファイルタイプの識別

cam? (2007-04-02 (月) 11:33:16)

Rではsave関数でXDR形式のバイナリデータを作ることができますが、データをこのXDRファイルで保存しておく場合と、csvのようなテキストで保存しておく場合と、色々あると思います。
データファイルを読み込む時に、まずはファイル名からこのファイルがXDRなのか、テキストなのかを判別して処理を切り替えたいのですが、ファイルタイプを識別できるような方法はありますでしょうか?

  • ?data しかしこれはファイル拡張子でファイルタイプを識別しているようですから、保存する際拡張子を標準的なものにしておくこころがけが肝要でしょう。 -- 2007-04-02 (月) 13:27:04
  • system 関数を使えばよろしいでしょう -- 2007-04-02 (月) 13:52:36
    > system("file test.XDR")
    test.XDR: gzip compressed data, from Unix
    > system("file test.data")
    test.data: UTF-8 Unicode text
  • ありがとうございました。参考になります。 -- cam? 2007-04-02 (月) 15:16:32

作業スペースの日本語が文字化けしてしまいます

akiba? (2007-03-31 (土) 21:40:55)

http://cran.md.tsukuba.ac.jp/から、R-2.4.1-win32.exeをインストールしたのですが、作業スペース上の日本語が文字化けしてしまっています。
文字化け部分をテキストに貼り付けると日本語としてちゃんと表示されます。
どうすれば、文字化けを直すことができるでしょうか。
ご教示をお願いします。

  • まず 日本語化掲示板 の過去記事を読みましょう。その上でまだ解決できなければ、ただ「できません」ではなく、もう少し情報をあたえないと他人にはわかりません。 -- 2007-04-01 (日) 00:07:12
  • とりあえず直すだけなら「編集」メニューの「GUIプリファレンス」ででる設定画面で Font のところを MS Gothicとかにすれば読めるようになります。FAQページが必要か? -- 2007-04-02 (月) 16:12:01
  • GUIプリファレンスのFontを変えたら、直りました。ありがとうございます。 -- akiba? 2007-04-09 (月) 13:05:39

パッケージのインストールができません

奥井正俊? (2007-03-29 (木) 13:02:38)

Rをインストールしたのちに、パッケージのインストールを試みましたが、つながりません。「以下にエラーopen.connection(file,"r"):コネクションを開くことができません 追加情報:Warning message:'cran.r-project.org'をポート80でコネクトできません」と出ます。どうすればよいのかご教示をよろしくお願いします。

  • 現在、Aizu&Tukubaにうまくつながっていないようです。CERNミラーサイトの選択で、Tokyoを選んでみてください。 -- okinawa 2007-03-29 (木) 13:17:44
  • たぶん質問された方の環境だとプロキシ設定をしないとインターネットにつながらないとかではないでしょうか。「ヘルプ」の「Windows版RのFAQ」のページにある2.19番の回答通り, Rの起動時にダブルクリックするショートカットのプロパティを開いて、「リンク先(T)」の....Rgui.exe" を ...RGui.exe" --internet2 と追記するといいかもしれません。 -- 2007-04-02 (月) 16:17:54

Rの2軸プロット

nick? (2007-03-28 (水) 19:52:19)

Rのプロットで,y軸の主軸,2軸の設定は可能でしょうか?
色々調べたのですが,わかりません。
どなたかご存知の方お願いします。

  • 例えば、plot(0,yaxt="n");axis(2,-2:2*0.1);axis(4,-4:4*0.2) -- takahashi? 2007-03-28 (水) 20:04:07
  • お礼が遅くなりました。どうも有難うございます。 -- nick? 2007-03-31 (土) 14:27:01

プロットしたい行を読みとばしたい

内藤? (2007-03-28 (水) 15:41:16)

Rの初心者です。
プロットしたいデータが10000行あって,行を100行おきとかに読み飛ばしたいのですが,何かコマンドみたいなものはあるのでしょうか?
gnuplotでいうとevery 100みたいな。
以上,どうぞよろしくお願いします。

  • 例えば、 dat<-rnorm(10005); plot(dat[1:100*100]) -- takahashi? 2007-03-28 (水) 17:43:59
  • どうもありがとうございます。助かりました。 -- 内藤? 2007-03-28 (水) 19:25:10

行列の1列目の中からある項目だけを取り出す

norion? (2007-03-27 (火) 23:56:46)

他で作った表をRに読み込んで作業しようとしています。表全体の読み込みはできました。
行列の1列目には色々な名前が書いてあるんですが、この中からひとつの名前を指定し
指定した名前に1列目が該当している行の全ての列を取り出したいのです。
どのような関数を書けばそれが可能でしょうか。

不適切な質問だったらすみません。

  • なんとなく行列ではなくて、データフレームではないかという気がしますが。いずれにしても列ラベル(一行目ではなく)が付いているようですから次のようにすればよいのでは。 -- 2007-03-28 (水) 05:41:43
    > M     # 行、列ラベルつき行列
       C1 C2 C3
    R1  1  4  7
    R2  2  5  8
    R3  3  6  9
    > M[,"C2"]  # 列ラベルで第二列を取り出す
    R1 R2 R3 
     4  5  6 
    > M["R2",]  # 行ラベルで第二行を取り出す
    C1 C2 C3 
     2  5  8 
    > str(M)     # 行列なら次のように表示されるはず
     int [1:3, 1:3] 1 2 3 4 5 6 7 8 9
     - attr(*, "dimnames")=List of 2
      ..$ : chr [1:3] "R1" "R2" "R3"
      ..$ : chr [1:3] "C1" "C2" "C3"
    > str(M)     # データフレームなら次のように表示されるはず
    `data.frame':   3 obs. of  3 variables:
     $ C1: int  1 2 3
     $ C2: int  4 5 6
     $ C3: int  7 8 9
  • 例えば,
    > A
         [,1] [,2] [,3]
    [1,] "a"  "1"  "7" 
    [2,] "b"  "2"  "8" 
    [3,] "c"  "3"  "9" 
    [4,] "b"  "4"  "10"
    [5,] "c"  "5"  "11"
    [6,] "a"  "6"  "12"
    > A[A[,1]=="b",]
         [,1] [,2] [,3]
    [1,] "b"  "2"  "8" 
    [2,] "b"  "4"  "10"
    こんな感じでいいのでは? -- g? 2007-03-28 (水) 10:46:12
  • ありがとうございました。 -- norion? 2007-03-28 (水) 22:58:11

関数が効かない

sunny? (2007-03-26 (月) 23:01:23)

初めて書き込みします。
Rを使い始めてわずかの、全くの初心者です。
barplotと入れても、表示をしてくれなくなってしまいました。
どうやら、何かの弾みにbarplot関数が壊れてしまったようです。
壊れた関数を再設定するにはどうしたらよいのでしょう。

  • ls() としてみてください。もし barplot と表示されたら、関数をうっかり再定義してしまったわけです。もしそうなら rm(barplot) とすれば本来の定義に戻ります。 -- 2007-03-26 (月) 23:06:33
  • ls()をやってみましたが、その中にはbarplotはありませんでした。すると、何か他の原因で効かなくなっているということでしょうか。 -- sunny? 2007-03-26 (月) 23:36:48
  • > 表示をしてくれなくなってしまいました
    とありますが,エラーメッセージも何も表示されないのですか?
    引数なしで,barplot とだけ入力したら何か表示されますか? -- 2007-03-27 (火) 06:43:47
  • Rをもう一度インストールしてみたらいかがですか。 -- 2007-03-27 (火) 07:44:21
  • エラーメッセージも何も表示されなかったんですが、インストールしなおしたら直りました。ありがとうございました。 -- sunny? 2007-03-27 (火) 23:10:26

グラフのX,Yの刻み指定と目盛線の引き方

のっち? (2007-03-23 (金) 00:45:59)

目的としては、CPUの負荷傾向を日々のCPU負荷状態(%)およびJOB稼働係数からグラフ化をしたいと思っております。
Rを入れて色々と調べているのですが、10, 20,30... といったY軸の刻み指定と任意の刻みでの目盛線の引き方が分かりません。
教えて君で申し訳ありませんが、ヒントなどいただけませんでしょうか。

  • 目盛線については、plot(...,tck=1)で引くことが出来ました。ただし、実線になるためどう変更するのか色々試していますがなかなか見つかりません。。 -- のっち? 2007-03-23 (金) 01:53:03
  • abline 関数,および線種ならばその引数として渡せる lty を調べてみましょう。lty については,par 関数を調べます。 -- 2007-03-23 (金) 07:21:17
  • 任意の刻みはaxisで指定しましょう。 -- 2007-03-23 (金) 09:21:55
  • 例として、
    plot(1:10, runif(10)*100, ylim = c(0, 100));abline(h=c(10,20,30,40,50,60,70,80,90), lty="dotted")
    で意図する刻みでの線が引けました。ありがとうございます。ただ、axisでX軸の刻みは指定出来るようですが、Y軸にも使えるんでしょうか? -- のっち? 2007-03-26 (月) 00:51:39
  • axis のヘルプ見ました?第一引数は何だったか分かりませんでしたか?
    c(10,20,30,40,50,60,70,80,90)は1:9*10でよいね。 -- 2007-03-26 (月) 08:29:25
  • 線を引くのに悩みましたが、axis( 2, 1:9*10, lty="dotted", tck=1 )で刻み指定と線が引けました。ありがとうございました。 -- のっち? 2007-03-26 (月) 23:09:12

小数以下が0を明示したいとき

X Jr.? (2007-03-21 (水) 21:22:39)

例えば、

y <- 1:10/2; plot(1:10, y); text(1:10, y, y, pos = 1)

とすると、小数以下が0のものは整数部分のみが表示されますが、".0" も表示させたいときはどのようにすれば良いのでしょうか? よろしくお願いいたします。

  • ?format として, nsmall をみませう. -- なかま 2007-03-21 (水) 22:02:17
  • y <- 1:10/2; plot(1:10, y); text(1:10, y, format(y, nsmall = 1), pos = 1) 
    でできました。ありがとうございました。-- X Jr.? 2007-03-21 (水) 22:13:05
  • y <- 1:10/2; plot(1:10, y); text(1:10, y, sprintf("%.1f", y), pos = 1) も,考えてみましょう。なお,呈示例では,文字が欠けてしまう場合もあるので,par(xpd=TRUE) も調べてみましょう。 -- 2007-03-21 (水) 22:28:33
  • なるほど。上の例で y <- 1:10/3 とした場合を試してみました。 ありがとうございました。 -- X Jr.? 2007-03-21 (水) 22:38:03
  • ?formatとして, nsmallを見ると, digitsというのも(略 -- なかま 2007-03-21 (水) 23:37:01

nlsで出力される回帰係数の読み込みについて

myuhe? (2007-03-14 (水) 18:55:45)

 はじめまして。Rのド素人です。次のことについて、ご教示いただければ幸いです。

 nls()を実行すると、例えば次のような結果が得られます

Nonlinear regression model
 model:  inverseHeight ~ 1/(A * DBH^B) + 1/C
  data:  testtrees
        A          B          C
0.2689986  3.1630292 24.6279807
residual sum-of-squares:  0.0003014474

 この中の回帰係数(上の例では、0.2689986 3.1630292 24.6279807)を利用して、続けて処理を行いたいと考えています。
 続けて行いたい処理は、出力された回帰係数で任意の独立変数の値の時の従属変数の値を連続して取得していくというものです。独立変数は、あくまで任意です。

 最初は、この結果をデータフレーム等に変換して、そこから回帰係数を取得しようと考えていたのですが、出力結果は独自のオブジェクトのようで変換できませんでした。
 
 predict()を使った方法も検討しましたが、任意の独立変数の値をうまく指定することができませんでした。

 初歩的なものかもしれませんが、よろしくお願いします。

  • 例えば、nls()が出力したnlsオブジェクトがhogeだとすると、-- 2007-03-14 (水) 19:11:42
    as.vector(hoge$m$getAllPars())
  • coef(hoge) -- takahashi? 2007-03-14 (水) 19:25:41
  • 何事によらず R オブジェクトの構造を見たければ str() 関数を使うのが定石です.nls() 関数の結果 xxx の構造を見たければ str(xxx) としてください.係数だけ取り出したければ総称的関数 coef() を使って coef(xxx) とする。予測をしたければ predict(xxx, yyy) とする。ここで yyy は当てはめに使ったデータと同じ構造のデータフレームです. -- 2007-03-14 (水) 19:30:25
  • coef()を使えば良かったのですね。str()もこれまでデータの要約程度の認識しかありませんでしたが、本来はデータ構造を見るものなのですね。  -- myuhe? 2007-03-14 (水) 19:53:47
  • coef()を使った方法でできました。多くのアプローチまで教えていただき、ありがとうございました。 -- myuhe? 2007-03-14 (水) 19:57:08

作成したモデル(予測式)の保存・呼び出し方法

ADMET modeler? (2007-03-09 (金) 15:37:11)

こんにちは。初めて投稿します。超初心者で非常に基本的な質問で恐縮なのですが、教えてもらえるとうれしいです。

仕事の関係で判別予測モデル(手法としてはBayesian, Random Forest, Support Vector Machineなどです)を作成して利用しようと考えているのですが、一旦作成したモデルを保存する方法がわかりません。

たとえば、Random Forestの場合、

hlm <- read.csv('C:/work/data.csv',head=TRUE)
trainset <- hlm[1:1952,]
testset <- hlm[1953:2439,]
model <-  randomForest(Class ~ ., data = trainset)

というふうにモデルを作成して、
そのままpredictすると別のデータを予測できますよね。

pred <- predict(model, testset)

しかし、できればここで作成したモデルを保存しておいて、後日呼び出して再び使いたいのです。もちろん、その際にもう一度モデルを作成しなおしてもいいのですが、データが多いのでモデル作成に時間がかかってしまいます。

この例では、Random Forestですが、一般的な回帰モデルなども含めて一度作成した「モデル」を保存・読み込みをするにはどのように行うか教えてください。

  • Rのオブジェクトをそのままディスク上に保存しておくには、?saveと?loadを読んでみる。 -- takahashi? 2007-03-09 (金) 15:52:31
  • あるいは、終了時にそのまま「ワークスペースを保存」しておけば次回起動したときにも残っているはず -- 2007-03-09 (金) 18:32:50
  • ファイルを読み込む tips 集(暫定版) で save.image() 関数の使いかたを見てください. -- 2007-03-13 (火) 17:27:44

連続する整数をまとめる方法

きゃらっと? (2007-03-08 (木) 08:31:34)

初めて投稿させていただきます.初心者です.
(1:3)が(1,2,3)やとなる変換の逆を行わせる方法を教えてください.

(1,2,...100, 250,251,...500, 10001, 10002,... 10005)
変換後
(1:100, 250:500,10001:10005)
よろしくお願いします.

  • 関数を書くしかないのではないでしょうか?その関数は,以下の通り(と書いておくと誰かが付け加えてくれるだろう) -- 2007-03-08 (木) 10:30:41
    vecpack<-function(x) {
            S<-x[-1]
            E<-x[-length(x)]
            TF<-E-S!=-1
            s<-c(x[1],S[TF])
            e<-c(E[TF],x[length(x)])
            myfun<-function(x) {
                    if(s[x]==e[x])
                            as.character(s[x])
                    else
                            paste(as.character(s[x]),as.character(e[x]),sep=":")
            }
            vlist<-lapply(1:length(s),myfun)
            paste("c(",paste(vlist,collapse=","),")")
    }
    早速付け加えてくれたようですね
    以下のように考えたのですが,ちょっと美しくないけど,せっかく書いたので。
    > set.seed(836478)
    > ( x <- unique(sort(sample(100, 90, replace=TRUE))) ) # テストデータ
     [1]  2  3  5  6  7  8 12 14 15 16 18 19 20 24 28 29 30 31 33 35 36 37 38 39
    [25] 40 42 45 49 50 52 53 54 55 59 61 62 63 65 66 68 69 70 71 72 73 74 77 80
    [49] 82 83 84 86 88 90 91 92 93 97 98 99
    > n <- length(x) # ベクトルの長さ
    > s <- c(TRUE, x[-n]+1 != x[-1]) # 新しいシークエンスが始まると TRUE を持つベクトル
    > res <- "" # 結果を得る文字列
    > for (i in 1:n) { # いやらしいループで
    + 	if (s[i]) { # 新しいシークエンスが始まると文字列に「, 整数値」を加える
    + 		res <- paste(res, x[i], sep=", ")
    + 		last <- x[i] # 最後に加えた整数値を覚えておく
    + 	}
    + 	if (i < n) { # シークエンスが終わるときは文字列に「:整数値」を加える
    + 		if (s[i+1] && x[i] != last) { # ただし,k:k のようなのは避ける
    + 			res <- paste(res, x[i], sep=":")
    + 		}
    + 	}
    + 	else if (x[i] != last) { # 最後の要素の処理
    + 		res <- paste(res, x[i], sep=":")
    + 	}
    + }
    + res <- paste(res, ")") # 最後に「)」を付けて
    > ( res <- sub(",", "c(", res) ) # 先頭が「,」で始まるので,「c(」で置き換える
    [1] "c( 2:3, 5:8, 12, 14:16, 18:20, 24, 28:31, 33, 35:40, 42, 45, 49:50, 52:55,
    59, 61:63, 65:66, 68:74, 77, 80, 82:84, 86, 88, 90:93, 97:99)"
    vecpack の vlist lapply を使わなくてもできるので,より分かりやすくなる
    vecpack2 <- function(x) {
    	S <- x[-1]
    	E <- x[-length(x)]
    	TF <- E-S != -1
    	s <- c(x[1], S[TF])
    	e <- c(E[TF], x[length(x)])
    	vlist <- ifelse(s==e, s, paste(s, e, sep=":")) # これだけで O.K.
    	paste("c(",paste(vlist,collapse=","),")")
    }
    巧い!
    vecpack3 <- function(x)
    {
     s <- x[c(TRUE, x[-n]+1 != x[-1])]
     e <- x[c(x[-n]+1 != x[-1], TRUE)]
     return(paste("c(",paste(ifelse(s == e, s, paste(s, e, sep=":")),
                             collapse=","),")"))
    }
    マイナーな修正ですが,未定義の n とスペースを除いて
    vecpack4 <- function(x)
    {
     s <- x[c(TRUE, x[-length(x)]+1 != x[-1])]
     e <- x[c(x[-length(x)]+1 != x[-1], TRUE)]
     return(paste("c(",paste(ifelse(s == e, s, paste(s, e, sep=":")),
                             collapse=","),")", sep=""))
    }
  • 早速のご返答ありがとうございました.早く皆さんのようにスマートな関数がかけるようになりたいです. -- きゃらっと? 2007-03-08 (木) 17:28:12

実験計画法

ちゃい? (2007-03-05 (月) 21:56:05)

Rで実験計画法ができるのでしょうか

  • 「実験計画法」という言葉で,あなたは何を意味していますか?実験計画法から得られたデータの分析なら,Rであろうと,Excelであろうと,自作プログラムであろうと(ちゃんと作られていれば)分析できないわけはありませんね。
    実験計画を立ててくれというなら,それはたぶん無理なんじゃないでしょうか?(できるのかも知れないけど,知らない) -- 2007-03-05 (月) 22:09:25
  • CRAN の貢献パッケージリストを最初から読む、または RSiteSearch? でキーワード(例えば "design of experiments")検索. -- 2007-03-05 (月) 23:25:32
  • 有難うございました -- ちゃい? 2007-03-06 (火) 21:48:47

integrate()で一様関数を積分できますか?

藤巻十三? (2007-02-28 (水) 17:52:22)

こんにちは。
表題の一様関数の積分は例えで、具体的にやりたいことは他にあるのですが、例として簡単なので取り上げさせて頂きたくお願い申し上げます。

f<-function(x) 1
integrate(f=f,lower=0,upper=1)

とすると答えとして「1」を期待したいところでしたが、

以下にエラーintegrate(f = f, lower = 0, upper = 1) : 
       evaluation of function gave a result of wrong length

となります。この原因と回避する方法を教えて頂けませんか?
尚、「f<-function(x) ceiling(x)」なら答えは1が返ってきますが、
「f<-function(x) max(x,1)」では駄目です。

  • すみません。わかりました。 -- 藤巻十三? 2007-02-28 (水) 17:57:42
  • なぜエラーになるかはおわかり(ヘルプ文章に書いてあるように定数関数は駄目)のようですが、真の理由は恐らく積分値を離散近似 S_n から計算する際の収束判定(相対誤差 |S_{n+1}-S_n|/|S_n-S_{n-1}| < epsilon)で分母がゼロになるためと思われます。ceiling(x) では端点 x=0 での値がたまたま 0 という 1 以外の値をとるため相対誤差の計算ができたからだろうと思われます。 -- 2007-02-28 (水) 21:02:08
  • そうじゃないです。fが返すベクトルの長さの問題です。f<-function(x)1+0*xとかしたら吉。helpではVectorizeを使えと書いてあります。 -- takahashi? 2007-03-01 (木) 17:30:42
  • 失礼しました.なるほど、離散点での関数値を同時にサンプリングしているようですね.すると max(x, 1) の例もベクトル化最大値 pmax(x, 1) を使えば良いわけですね. -- 2007-03-01 (木) 19:47:39

リスト内の各ベクトルから特定の要素のみを抽出する方法

生物系大学院生? (2007-02-26 (月) 06:11:57)

こんにちは、初めて質問させていただきます。
下記のようにして作成したブートストラップデータのリスト(z)内の各ベクトルから、特定の要素のみを抽出するいい方法をご存知ないでしょうか?

> x<-c(1:100)
> y<-function(){sample(x,100,replace=TRUE)}
> z<-lapply(1:100,function(i)try(y()))

塩基配列データに基づくパラメータの推定値の分散をブートストラップ法で求めたいと考えています(式は100塩基座の100回のブートストラップ試行の場合です)。推定するパラメータの性質上、情報をもつ特定の塩基座(informative site;例えば1〜10)がそれぞれ何回各ブートストラップの試行で利用されたかが重要になってくるため、それだけを抽出したいと考えています。

各ベクトルについては、

> z[[1]][z[[1]]<11]

として11未満の数字のみを抜き出す方法や、

> sort(z[[1]])

として先頭にinformative siteを集める方法にはたどり着いたのですが、これをリスト全体に適用する方法が分かりませんでした。
どなたかいい方法をご存知ないでしょうか?

  • もう一歩。初心者ということならまずは for ループを使う。その上で、すでに lapply 関数の存在を知っているのだから応用を考える。 -- 2007-02-26 (月) 08:24:55
    > z <- lapply(1:100, function(i) sample(1:100, 100, rep=TRUE))
    > A <- list(rep(NULL,100))  # 空リスト
    > for (i in 1:100) A[[i]] <- z[[i]][z[[i]] < 11]
    > B <- lapply(1:100, function(i) z[[i]][z[[i]] < 11])
    > identical(A, B)
    [1] TRUE
  • 早速解法を教えていただきどうもありがとうございました。 新しいリストを構築するということまで考えがおよびませんでした。 教えていただいた方法で思い描いていた結果が得られました。 どうもありがとうございました。 -- 生物系大学院生? 2007-02-26 (月) 10:06:17

陰関数の表示について

もも? (2007-02-24 (土) 17:27:16)

単発の質問で申し訳ないのですが
陰関数を表示する関数、あるいはパッケージをご存知の方教えてくださいませ。

  • 「陰関数を表示する」とはどういう意味ですか。 -- 2007-02-24 (土) 20:34:52
  • すみません。自己解決しました。等高線=0に気づいていませんでした。 -- もも? 2007-02-25 (日) 10:43:29
  • そのノウハウを知っているといつか役にたつ(?)関数達(2)とかにまとめられると良いと思います。 -- 2007-02-25 (日) 12:30:28

他言語で作成したプログラムのコンパイル方法について。

Y.Takenaka? (2007-02-18 (日) 16:52:55)

mt19937の乱数の改良版が発表されていましたので、ラッパーを作った。そこでわからない事にぶち当たりました。altivecやsse2といったものにも対応しているんですが、対応させる方法がわかりません(コンパイルオプションを付ける方法です。R cmd SHLIB ファイル オプション?)もしよろしければ教えていただけませんか?情報が載ってる先がわかるだけでも助かります。

altivecやsse2と行ったものに対応していないものでは使えるようになったのです。公開してます。
http://sun.s167.xrea.com/Archives/HomePage.html (r_sfmt19937.cです。)

  • もう少し状況が分かるように説明した方がよろしいのでは? -- 2007-02-18 (日) 17:57:12
  • コンパイルオプションを付けるなら,
    R CMD COMPILE CFLAGS=-msse2 a.c b.c c.c
    R CMD SHLIB -o abc.so a.o b.o c.o
    こんな感じになります.
    R  CMD COMPILE -h
    R  CMD SHLIB -h
    を見てください.-- なかま 2007-02-18 (日) 19:05:28
  • 情報ありがとうございます。2つめの方法でうまくいったです。(なかまさんのかかれてるR CMD SHLIB -o .... のもの)あっ、それと説明ですが、mt19937ってのはこのwikiでも名前がでてた事がある乱数発生方法で高速かつ周期が長いというアルゴリズムですね。それが改良されて高速化されたということが公開先(C言語で公開されている。)に書かれていたのでそれをRで使えるようにしようと思ったんです。そこでC言語でラッパー(Rで使えるようにするプログラム)を作ってみたのです。 -- Y.Takenaka? 2007-02-18 (日) 20:33:37
  • Rで使えるようにと言う意味ならば, help(Random.user)するといいですよ. runifは元より*norm,sampleなどもその乱数を使ってくれます. ただ乱数をベクタライズしている訳では無いので, 関数呼び出しのオーバヘッドが大きく, 高速化のメリットはあまり無いかも. -- なかま 2007-02-18 (日) 21:11:11
  • なるほど、おっしゃるように関数のオーバヘッドは大きい事は影響が強くでそうですね。これは少し残念ですけど。help(Random.user)みました。改良を試みてみます。なかまさん!適切なアドバイスありがとうございました。 -- Y.Takenaka? 2007-02-18 (日) 22:23:26

自己回帰係数の算出方法

千成? (2007-02-16 (金) 13:17:45)

ARモデルの自己回帰係数の算出方法として、
Yule-Walker法、最小二乗法、最尤法、Burg法などがありますが、何故デフォルトではYule-Walker法が選択されているのでしょうか?

Yule-Walker法が他の方法に比べてどの点で秀でているのでしょうか?
助言をお願い致します。

  • ここで尋ねるような話では無いですね。時系列の本を(何冊か)あたってみてください。ただ、色々あるということは、どれも一長一短ということに他ならないのでしょう。 -- 2007-02-16 (金) 13:59:15
  • やはりそうですか…。ありがとうございます。 -- 千成? 2007-02-17 (土) 11:10:41

整数の割り算について

廣瀬_敏之? (2007-02-15 (木) 02:16:56)
次のような結果をどう考えればいいのでしょう? 98-BASICの整数の割り算の結果とあまりに違うのでとまどっています。 どちらが正しいか、という問題でない事はわかっていますし、整数部をとったり、切捨てなら専用の関数があることも知っていますが、それでも初心者にはとまどう『出来事』です。 各々のround(x,0)をとってから割り算をしていると考えても少しおかしいし・・・? この計算の定義までマニュアルに書いてあったでしょうか? いわゆる『R言語のマニュアル』と、過去の『初心者向けQ&A』に対する検索で『整数』『割り算』『除算』では、目的とする内容がヒットしないようなので、せめてキーワードだけでも教えてください。多分、私の探し方が足りないと思うので・・・(礼)。 使用環境はDOS/V系のパソコンで、『XP』と『Vista』で同じでした。それと

R version 2.2.0, 2005-10-06, i386-pc-mingw32 
attached base packages:
[1] "methods"   "stats"     "graphics"  "grDevices" "utils"     "datasets" 
[7] "base"     

です。本当は複数のバージョンでチェックすべきなのでしょうが、『R』の使用環境が変わっては困るので、そこまではご勘弁を。
なお『R』自体のソースコードをどうぞ、と言われましても、『C』は読めません。

> -2.3%/%1
[1] -3
> test<-c(-5.99999,-5.00001,5.00001,5.99999)
> test%/%1
[1] -6 -6  5  5
> test%/%2
[1] -3 -3  2  2
> test%/%2.5
[1] -3 -3  2  2
> test%/%2.2
[1] -3 -3  2  2
> test%/%2.7
[1] -3 -2  1  2
> test%/%-1
[1]  5  5 -6 -6
> test%/%-2
[1]  2  2 -3 -3
> test%/%-2.5
[1]  2  2 -3 -3
> test%/%-2.2
[1]  2  2 -3 -3
> test%/%-2.7
[1]  2  1 -2 -2
  • もうしわけありませんでした。helpで分かるそうです。 -- 廣瀬_敏之? 2007-02-15 (木) 05:42:32
  • 演算記号に対してまでhelpが使えるとは思いませんでした。パソコンについては初期世代なので、適切なhelpの使い方がいまだにわからない場合があります。ご迷惑をおかけしました。 -- 廣瀬_敏之? 2007-02-15 (木) 05:47:46
  • %/% はいわゆる整数商ですが、98-basic (懐かしい響き!)の対応する結果とはどんなものなのでしょうか。逆に気になりました。 -- 2007-02-15 (木) 07:15:18

plot()のy軸ラベルと目盛りの表示について

kokko? (2007-02-14 (水) 18:22:11)

x <- 1:5
y <- c(100,1000,10000,100000,1000000)
plot(x,y, log="y", xlab="X_num", ylab="Y_num", las=1)

こんにちは。このようにy軸を対数目盛にしたグラフを作製したいのですが、y軸のラベル "Y_num"を時計回りに45°回転させてx軸と平行にし、さらに目盛りの数値と重ならないようにしたいです。
また目盛りの値が "1 e+05" のように1とeの間にスペースが入ってしまい、このスペースを消したいのですがどのようにすればよいでしょうか?
どなたかよろしくお願いします。

  • > y軸のラベル "Y_num"を時計回りに45°回転させてx軸と平行
    時計回りに 90 度回転ではないですか?それはともかく,text を使って描けば,どこへでも描けると思いますがどうでしょう?
    > "1 e+05" のように1とeの間にスペースが
    お使いの R のバージョンと,OS の種類およびバージョンは?
    私は R 2.5.0 で Mac ですが,スペースが入っているようには見えないんですけど。
    もし,どうしても気になるということなら,axis 関数で,文字列として自分で描いてやればよろしいかと思います。
    図や表の体裁も,自分ならこのように描く(書く)のにどうしてできるようになっていないんだろうと思うかもしれませんが,そのようになっているのはその世界では慣用的(常識的)なことで,あえて自分流を主張するとかえって異様に思われたりすることもあるかもしれません。-- 2007-02-14 (水) 18:30:22
  • ご回答ありがとうございます。余白に書き込むためにはmtext()という関数を使えばよかったのですね。windowsでR2.3.1を使っていましたが、2.4.1にしたら1とeの間のスペースがなくなりました。 -- kokko? 2007-02-14 (水) 22:27:04

Windows版 help() の表示方法の変更

X Jr.? (2007-02-03 (土) 21:31:51)

R-2.4.0以降のWindows版で、インストール時にヘルプの表示形式(テキスト, chm, HTML)を選択できますが、これをインストール後に変更する方法を教えて頂けますでしょうか。

  • Rprofile.siteと言うファイルがR_HOMEのetc以下にあると思います. 漢はテキストなのでそれらをコメントにします. -- なかま 2007-02-04 (日) 14:12:35
  • 遅くなりましたが、ありがとうございました。私はヘルプをHTML表示にしました。Firefox + greasemonkey で、英単語を選択すると自動的に和訳が出てくるようにして使ってます。ちょっと便利ですので、ご参考までに。 -- X Jr.? 2007-02-18 (日) 23:12:30

wilcox.exact で,P 値が 1 より大

ウィリアム (2007-02-02 (金) 07:57:57)

Wilcoxon signed rank test を行おうと思いましたが、同じ値がありましたので Exact Wilcoxon signed rank test を行いました。 小さなデータセットで試してみたところ P 値が 1 より大きくなりました。

> wilcox.exact(c(8,8,-5,-10))

       Exact Wilcoxon signed rank test

data:  c(8, 8, -5, -10) 
V = 5, p-value = 1.125
alternative hypothesis: true mu is not equal to 0

P 値が 1 より大きくなってよいのですか?

  • 貢献パッケージ中の関数(ですよね)に言及する時はパッケージ名もそえること。それと help(wilcox.exact) は読みましたね。 -- 2007-02-02 (金) 08:16:27
  • すみません。書くのを忘れました。 >library(exactRankTests?) でした。 -- ウィリアム 2007-02-02 (金) 08:28:24
  • 質問の意図が、バグ報告なのか、どうしても計算したいのに困った、なのか判然としませんが、help(wilcox.exact) を読めば、タイがある場合は「正確」ではなく正規近似を使うと書いてありますから、こうした結果が出てもおかしくない(近似式の精度の問題)し、バグでもないと思います。基本パッケージのそれを例えば c(8, 8+1e-10,-5,-10) , オプション exact=TRUE でやれば済むことと思いますが。 -- 2007-02-02 (金) 12:22:52

指定されたベクトルのサイズが長すぎます

szk? (2007-02-01 (木) 09:03:07)

以下のエラーを解決する方法をお尋ね申し上げます。

> dist(rnorm(560327))
以下にエラーvector("double", length) : 指定されたベクトルのサイズが長すぎます

以下は正常に計算されます。

> dist(rnorm(560327)[1:10])
  • dist(1:10) が何を返すかは理解できていますか?
    dist((1:10)[1:3]) とはどこが違うか分かりますか?
    あなたが要求した操作には(実数を格納するには)560327×560327×8バイトくらい必要だということはご存じですか?
    そして,あなたのコンピュータはそれだけのメモリを使えますか? -- 2007-02-01 (木) 09:14:30
  • 2.3テラバイト! 東工大の誇るグリッドスパコン TSUBAME の総合メモリ容量は21.4テラバイトだそうだから、実行できるかも。 -- 2007-02-01 (木) 09:25:01
  • 1.2Tぐらいかと. TSUBAMEは共有型クラスタではなかったと思うので,無理だったかと. たとえ64bit Rでも, dist(rnorm(2^16))付近が限界(作れそうな範囲と言う意味)で, ちゃんと使える(いろいろ加工と言う意味)レベルだと, dist(rnorm(2^14))ぐらいだと思います(たとえ64bitでも). -- なかま 2007-02-01 (木) 13:19:17

添付ファイル: filesna_sample3.png 1806件 [詳細] filess2.png 1698件 [詳細] filebaz.png 1683件 [詳細] fileupload.png 2158件 [詳細] fileK_barplot.png 894件 [詳細] filehist.nag.png 1843件 [詳細] filesna_sample4.png 969件 [詳細] filepairs.png 1689件 [詳細] filekuni.png 1914件 [詳細] filelegendmiss.png 2005件 [詳細] filefjok.PNG 1009件 [詳細] filehistgram.png 961件 [詳細] file3dplot.png 876件 [詳細]

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