R および RjpWiki に関する質問コーナー

注意:新規記事用の入力欄は以下の目次の直後にあります。各記事への追加コメントの入力欄は、各記事の後にあります。質問が長くなる場合は、まず出だしだけ投稿し、次に編集ボタンでそれを修正、追加するのが便利でしょう。

もし新規に質問をされる場合は、効率的にコメントを得るためにも、善意のコメント者の時間を無駄にしないためにも、以下に特にご注意下さい。R のメイリングリストへの投稿記事のガイド も参考になるでしょう。

  • 広い意味で R に関係する話題に限定!
  • 学校の宿題は自分で考えるべきです
  • 具体的な状況がわかる様に背景説明をけちらない
  • CRAN の検索エンジンマニュアル を調べた上での質問かどうかを書き添えれば、空振りコメントが減るでしょう。また必須ではありませんが r-help の投稿ガイダンス R のメイリングリストへの投稿記事のガイド も参考にして下さい。
  • 自分には馴染みでも、他の人には馴染みの無い手法・概念もあります。日本語は得てしてわかりにくい(英語名をそえると案外わかるものです)ものです。
  • もし問題が解決した(結局解決しない場合を含め)場合はその旨報告していただくとコメント者や他の人に参考になります。また後で関連情報をこのWikiに要約記事としてご投稿いただくと今後の参考になります。

    注意:このコーナーは重量オーバーで表示に時間がかかり過ぎるようになりました。今後は過去記事へのコメントを除き、新設の Q&A コーナーを利用下さい。



群散布図(Group Scatter)を作成するパッケージはありますか?

(2004-05-13 (木) 22:15:48)

ここにあるような群散布図(Group Scatter)を作成するパッケージはありませんか。

  • パッケージ ade4 にそれらしき関数が多数用意されていますね.
> library(ade4)
> data(lascaux) 
> scatter(dudi.acm(lascaux$ornem, sca = FALSE), csub = 3)
plot-00.png

例えばこんな関数がありました.

> x <- c(0.5,0.2,-0.5,-0.2) ; y <- c(0.2,0.5,-0.2,-0.5)
> eti <- c("toto", "kjbk", "gdgiglgl", "sdfg")
> plot(x, y, xlim = c(-1,1), ylim = c(-1,1))
> scatterutil.eti.circ(x, y, eti, 2.5)
> abline(0, 1, lty = 2) ; abline(0, -1, lty = 2) 

> x <- c(0.5,0.2,-0.5,-0.2) ; y <- c(0.2,0.5,-0.2,-0.5)
> eti <- c("toto", "kjbk", "gdgiglgl", "sdfg")
> plot(x, y, xlim = c(-1,1), ylim = c(-1,1))
> scatterutil.eti(x, y, eti, 1.5)

> plot(runif(10,-3,5), runif(10,-1,1), asp = 1)
> scatterutil.grid(2)
> abline(h = 0, v = 0, lwd = 3) 

> x <- runif(10,0,1) ; y <- rnorm(10) ; z <- rep(1,10)
> plot(x,y) ; scatterutil.star(x, y, z, 0.5)
> plot(x,y) ; scatterutil.star(x, y, z, 1) 

> x <- c(runif(10,0,0.5), runif(10,0.5,1))
> y <- runif(20)
> plot(x, y, asp = 1) # asp=1 is essential to have perpendicular axes
> scatterutil.ellipse(x, y, rep(c(1,0), c(10,10)), cell = 1.5, ax = TRUE)
> scatterutil.ellipse(x, y, rep(c(0,1), c(10,10)), cell = 1.5, ax = TRUE)

> x <- c(runif(100,0,0.75), runif(100,0.25,1))
> y <- c(runif(100,0,0.75), runif(100,0.25,1))
> z <- factor(rep(c(1,2), c(100,100)))
> plot(x, y, pch = rep(c(1,20), c(100,100)))
> scatterutil.chull(x, y, z, opt = c(0.25,0.50,0.75,1))
> par(mfrow = c(1,1))
plot-01.png

関数が多数あるので全部は調べ切れませんでした.すみません. -- 舟尾 2004-05-13 (木) 22:52:40

  • いろいろ調べるのが面倒なので,作る方が楽。 -- 青木繁伸 2004-05-14 (金) 00:06:03
    # テストデータの作成
    mean <- c(3,5,6,8,7)
    sd <-c(0.3, 0.4, 0.5, 0.2, 0.4)
    
    x <- rep(1:5, each=20) # 5群(20ケースずつ)のデータ,平均値・標準偏差は上述の通り
    y <- NULL
    for (i in 1:5) y <- as.integer(100*c(y, rnorm(20, mean=mean[i], sd=sd[i])))/100
    
    #accu 測定値(縦軸)の精度(階級幅)
    #stp  横軸方向のずらし量
    # 以上二つは,試行錯誤で決定(^_^;)
    
    graph <- function(x, y, accu, stp)
    {
      y <- round(y/accu)*accu
      x1 <- unique(x)
      for (i in 1:length(x1)) {
        freq <- table(y[x==x1[i]])
        for (j in 1:length(freq)) {
          if (freq[j] >= 2) {
            for (k in 1:length(y)) {
              if (abs(y[k]-as.numeric(names(freq)[j])) < 1e-10 && abs(x[k]-x1[i]) < 1e-10) {
                x[k] <- x[k]+(freq[j] <- freq[j]-1)*stp
              }
            }
          }
        }
      }
      plot(x, y)
    }
    
    graph(x, y, 0.2, 0.05)
graph.png
  • 質問した者です。ありがとうございました。青木先生の関数を使わせていただきます。 -- 2004-05-14 (金) 11:31:01
  • 凄いっ!いつもながら鮮やかなソースコードですね.簡潔さと分かり易さを兼ね備えたコードはそう簡単に書けるものではないですよ!勉強させてもらえました. -- 舟尾? 2004-05-15 (土) 18:07:32

巡回セールスマン問題

(2004-05-06 (木) 17:15:50)

巡回セールスマン問題用の関数がRのあるパッケージにあったようですが失念してしまいました。どなたか、ご存知ありませんか?

  • ?optim です。なんでも掲示板で以前教えていただきました。 -- なかま 2004-05-06 (木) 17:34:15
  • example(optim) を実行すると、組み込みデータセット eurodist を使った解析例が見られます。SANN (simulated annealing) 法を使った例です。 -- 2004-05-06 (木) 18:24:16

matplotの軸

さち? (2004-05-03 (月) 00:00:58)

matplotで軸の正の方向を左向きにしたいと思い、

matplot(c(30,-80),c(50,-50))

のようにしたのですが、グラフにしてみると軸は正の方向が右側のままでした。
matplotで軸の方向を変える方法を知っている方がいたらアドバイスをお願いします。

  • x 軸の正の方向を左向きにする場合は matplot() に引数 xlim を与えて, xlim で調節してください.この方法は他の高水準作図関数 (例:plot()) でも使えます.
sines <- outer(1:20, 1:4, function(x, y) sin(x / 20 * pi * y))
matplot(sines, xlim=c(20,1), type="o", pch=1:4, col = rainbow(ncol(sines)))
matplot-00.png

y 軸側も同様にして 引数 ylim で調節してください.-- 舟尾? 2004-05-03 (月) 14:35:36

sines <- outer(1:20, 1:4, function(x, y) sin(x / 20 * pi * y))
matplot(sines, pch = 1:4, xlim=c(20,1), ylim=c(1,-1), 
        type = "o", col = rainbow(ncol(sines)))
  • xlimを求めたい場合
    matplot(c(30,-80),c(50,-50))
    matplot(c(30,-80),c(50,-50),xlim=par("xaxp")[2:1]
    とか、2度描画させるのが一番楽なのでしょうか? -- 2004-05-04 (火) 05:05:43
  • matplot() で描いた直後は座標系の値がグラフィックスパラメータ "xaxp" に保存されていますので,matplot() を実行した後に par("xaxp") を実行して下さい.
par(xaxp=c(0,1,5))
par("xaxp")        # xlim は [0 1]

[1] 0 1 5

plot(1:10)         # xlim ( c(xaxp[1],xaxp[2]) ) がこの時点で変更される
par("xaxp")        # xlim は [2 10] に

[1]  2 10  4

par(xaxp=c(0,1,5)) # 自分好みの座標系にする場合はxaxpを設定した後,xlimで指定する
plot(1:10, xlim=c(xaxp[1],xaxp[2]))
par("xaxp")

[1] 0 1 5

余談ですが,par() で設定して (永続的に) 自分好みの座標を設定することは出来る のでしょうか?"usr" や "xaxp" は値がコロコロ変わってしまうので, 仕方なく xlim , ylim で一時的に変更しているのですが.どうも私の使い方が悪い せいか,"usr" や "xaxp" の適当な使い方が思いつきません.-- 舟尾? 2004-05-04 (火) 11:14:34

# 原点を通る座標軸を描く
plot(rnorm(30), rnorm(30), axes = F, xlab = "", ylab = "")
axis(1, pos = 0, at = pretty(par("usr")[1:2]), adj = 0)
axis(2, pos = 0, at = pretty(par("usr")[3:4]), las = 2, adj = 1)
box()

[質問]ベクトル

さち? (2004-04-30 (金) 05:19:49)

ある変数に値を代入していき、その変数に値をベクトルのように
記憶させていくプログラムを作りたいです。例えば、
x<-1
x<-3
x<-9
x
1 3 9
のようにしたいです。このとき、
x[3]
のように長さを与えれればよいのですが、
いま代入する変数の数がいくつになるのかわからない状態です。
このような時はどうしたらよいのでしょうか?
初歩的な質問ですが、アドバイスなどお願いします。

  • R のベクトルは任意に長さを延長できますから気にせずに代入していけばOKです(このへんが慣れれば C なぞの低級言語を使えなくなる理由の一つです)。ただし効率性の観点(巨大なベクトルを扱う際実行時間が問題になる)からは、可能な最大サイズを見積もり x <- numeric(maxsize) としておく方がベターでしょう。-- 2004-04-30 (金) 08:14:46
> x <- numeric(0)  # x <- NULL でもOK (この時点では x は空の変数)
> x
numeric(0)         # この時点では x は空の変数
> x[1] <- 1
> x[2] <- 2
> x[3] <- 3
> x
[1] 1 2 3 
> x[5] <- 5 # 途中飛ばしてもOK
> x
[1]  1  2  3 NA  5  # 未定義の x[4] には自動的に NA 値が代入される
  • append という関数を使うと吉-- 青木繁伸 2004-04-30 (金) 10:33:09
    > x <- NULL
    > x <- append(x, 1)
    > x <- append(x, 3)
    > x
    [1] 1 3
    > x <- append(x, 2, after=1)
    > x
    [1] 1 2 3
  • ありがとうございます。 -- さち? 2004-04-30 (金) 16:15:16
  • ありがとうございました。早速ためしてみます! -- さち? 2004-04-30 (金) 16:16:24

princompのloadingsは固有ベクトルそのもの?

石原茂和? (2004-04-27 (火) 07:42:09)

(すみません,同内容を旧Rwikiの掲示板にまちがえてかいてしまいました)

こんにちは.いま気がついたのですが,princompのloadingsは
固有ベクトルそのものですね.
helpでも

loadings: the matrix of variable loadings (i.e., a matrix whose columns
   contain the eigenvectors). 


と,そのように書いてあります.
しかし,主成分負荷量は,加工していないもとの変数での値と,主成分スコアの相関係数
のはずですが,
確かに固有ベクトルとたいがいの場合,同じような傾向になりますが
固有ベクトルそのものを負荷量というのは誤解のもとでは?

  • example(Harman74.cor) を実行してみてください。おっしゃる意味の負荷量と、固有ベクトルの両方が表示されませんか。 -- 2004-04-29 (木) 07:34:04
  • 出ますけど,原質問者の質問の答えになっているのかな?
    それと,必ず回転された答えが出るんでしょうか。何回転になっているのかな?(教えてちゃんモード)
    二変数の分析してみると,教科書と違った解になるようですが。
    	> x
    	            [,1]       [,2]
    	[1,] -0.24451455  0.5090637
    	[2,]  0.04563039 -0.5012150
    	[3,]  0.97775313 -1.1832550
    	[4,]  1.26989942  1.4814664
    	[5,] -0.14870401 -1.3809846
    	> loadings(princomp(x))
    	
    	Loadings:
    	     Comp.1 Comp.2
    	[1,]  0.273  0.962
    	[2,]  0.962 -0.273
    	
    	               Comp.1 Comp.2
    	SS loadings       1.0    1.0
    	Proportion Var    0.5    0.5
    	Cumulative Var    0.5    1.0
    • 2004-04-29 (木) 19:22:43
  • 失礼、Harman74.cor の example は princomp ではなく、factanal を使った例でした。たまたま loadings という単語がでていたので早とちりしました。 -- 2004-04-29 (木) 21:42:41

数量化1類で消費されるメモリ量について。

koichi? (2004-04-19 (月) 23:48:32)

はじめまして。Rの初心者です。
Rで数量化1類を行おうとしているのですが、メモリの問題で手を焼いています。
具体的には、
1列目(V1)を『0から1までの8桁の乱数』、2列目(V2)を『1から1000までの整数のどれか』としたレコード50000件(テキストファイルで820KB)をd0に読み込み、V2を因子として水準1000の数量化1類を行います。

d0$V2 <- factor(d0$V2)
result <- summary(lm(V1 ~ V2,data = d0))


すると

Error:cannot allocate vector of size 390625 Kb
Reached total allocation of 1022Mb


と出て、タスクマネージャーを見るとRは「応答なし」で止まっています。

このデータ量で390625 Kbもの行列を生成してしまうのでしょうか。それにしても2行目を見ると1022Mbにはまだ余裕があるように思えます…。

メモリ関連の操作や、計算法の工夫を通してなんとか計算させる方法をご存知の方がおられましたら、是非ご教授お願いしたいと思います。

当方の環境はWindowsXP,Pentium4 2.8GHz,1GBRAMです。
よろしくお願いいたします。

  • 変数の長さが n ならば、回帰分析で登場するハット行列は nxn ですから、そのメモリー量は n^2*8 です。n=50000 なら 2.33GB ですから失敗するのは当然では。5万件のデータをそのまま解析する必要が本当にあるのでしょうか。 -- 2004-04-20 (火) 00:31:09
  • 一つのやりかたは、必要な情報だけを公式通りに自前で計算することでしょうか。R の lm() 関数は、二次的解析・出力用に解析結果に対する詳しい附属情報を隠れて保持しています。当然、大きなデータに対しては、サイズが莫大になりがちです。そんなに遠くない(?)将来、メモリーたっぷりの64 bit マシンに R が本格的に対応すれば、こうした問題も解消するのでしょうが。 -- 2004-04-20 (火) 07:34:23
    # 参考 10055007 は文字列、数値の項目数ですから、
    # メモリー量はざっとその 8 倍程度になると考えるべきでしょう
    > n=100; x=rnorm(n); y=factor(floor(1000*runif(n))); z=lm(x~y);  length(unlist(z))
    [1] 10902
    > n=1000; x=rnorm(n); y=factor(floor(1000*runif(n))); z=lm(x~y);  length(unlist(z))
    [1] 661272
    > n=10000; x=rnorm(n); y=factor(floor(1000*runif(n))); z=lm(x~y); length(unlist(z))
    [1] 10055007 # z に含まれる項目総数

  • 提示された条件に従ってデータを生成して数量化I類を(定義により)やってみましたが,なぜか非常に時間がかかります。スワップが頻発しているのでしょうか。
    そもそも,このような設定の数量化I類が必要なのか適切なのか疑問です。 -- 2004-04-20 (火) 17:18:05
  • 解析対象のデータがわかりませんから、野次馬コメントですが、一般に5万ものデータになると、はたして皆同質とみなして良いかどうか心配になりますね。一度ご検討ください。スピードが遅いというのも詳細がわかりませんから無責任なコメントになりますが、ひょっとして愚直な for loop なぞ使っていませんか、R はベクトル演算に最適化されていますから、コードを今一度点検されることをお勧めします。最も、早いコードを工夫するよりも愚直に計算する方が結果的には早いということもありそうですが。 -- 2004-04-20 (火) 21:36:25
  • みなさん、早速のご回答ありがとうございます。また、つたない編集によってご迷惑をおかけしています(^^;
    参考のために提示された式、とてもためになりました。Rに習熟するとそんなスマートに示せるんですね。n^2*8のメモリが必要になるという件ですが、水準数を100にすると50000件のデータでも問題なく計算されます。消費されるメモリ量には件数だけでなく、水準数も(むしろ件数より大きく)関与しているようなのですが、その関与の仕方は当方にはわかりません。length(unlist(z))を色々なnに対して試してみることによって(結果が出るかどうかによって)、メモリ的に扱えるデータ量が大体わかってきました。length(unlist(z))によってわかるのは、結果を保持するためのメモリですが、計算途中で消費されるメモリ量も大体これに比例すると考えればよいでしょうか。
    提示された条件に従ってデータを生成して数量化I類を(定義により)やってみたとのことですが、後学の参考にコードを示していただけませんか?Rプログラミングの超初心者で、とてもためになると思いますので。 -- koichi? 2004-04-20 (火) 21:58:23
  • 数量化I類のプログラムは,ここ においてあります。
    for を使わないとかえって複雑怪奇になるところもあるので(と思っているので)あんなプログラムです。ダミー変数を使ったデータ行列に展開する所なんかも,for の方がわかりやすい(^_^)。
    実行したのは Mac G5 です。
    必要メモリー数ですが,原データが50000×2×8バイト,原データをカテゴリーデータに展開したものが50000×1000×8バイトでしょう(最初のコメントの方は変数の個数を n と書いていながら,その後でn×nの計算のとき50000^2としてしまっています),分散共分散行列(相関係数行列)は1000×1000×8バイト
    cov 関数で分散共分散行列を計算するようにしたので,展開したダミーデータ行列を残しておく必要がありましたので(あとで,予測値を計算するときにも必要になるし。。。)
    プログラミング言語で書くときには,こんなにメモリを贅沢に使うコードは,もったいないお化けが出るので書けません(^_^) -- 青木繁伸 2004-04-21 (水) 00:34:19

.Call/.Externalで使用するDLLの作成方法について

どうむ (2004-04-19 (月) 15:20:39)

現在R(Windows版)から呼び出せるCで書いたDLLを作成中です。環境は
■R1.8.1
■コンパイラー:「Borland C++ Compiler 5.5」
■OS:WindowsXP
です。いままでRやSPlusでCで書いたルーティンを作ってきたのですが、今回はCからRの関数を呼び出すことに初めてチャレンジせねばならず、この点に関してアドバイスをいただけると幸いです。「Writing R Extensions」の第4章には一応目を通しましたが。。。
まず、簡単なCルーティンですが

extern "C" __declspec(dllexport) void R_TEST001(int *x, int *ret);

__declspec(dllexport) void R_TEST001(int x[], int *ret)
{
  *ret = x[0]+x[1]+x[2]+x[3];
}

整数型のベクトルを引数にとってその合計を「ret」で返している簡単な関数です。
DLLを作成するのにMakefileを作っておりコンパイルとリンクは以下のとおりです。

bcc32 -u- -I"***" -O2 /c ***.cpp
ilink32 /w /m /Tpd /L"***?Lib"
    "C:?borland?bcc55?Lib?C0D32.OBJ" ***.obj,
    ***.dll, 
    CW32.LIB IMPORT32.LIB
(mapファイルの指定が欠如していますが気になさらぬよう。。。)

「***.cpp」が該当のCルーティンです。
これを呼び出すRのプログラムは

> dyn.load("***.dll")
> is.loaded("R_TEST001") #確認
[1] TRUE

と、なります。結果は

> ret <- .C("R_TEST001",as.integer(c(1,2,3,4)),ret=integer(1))
> ret$ret
[1] 10

です。
ここまでは、良いのですがRの「.Call」関数を使用できるようにまずヘッダーファイルを加えてみます。

#include <R.h>
#include <Rinternals.h>

extern "C" __declspec(dllexport) void R_TEST001(int *x, int *ret);

__declspec(dllexport) void R_TEST001(int x[], int *ret)
{
  *ret = x[0]+x[1]+x[2]+x[3];
}

たった、これだけで上のコンパイルが通らなくなりました。
が以下のように「-A」オプション(ANSI 準拠の予約語のみを使用する)を加えることでなんとかコンパイルは通るようになりました。

bcc32 -u- -I"***" -O2 -A /c ***.cpp
ilink32 /w /m /Tpd /L"***?Lib" 
    "C:?borland?bcc55?Lib?C0D32.OBJ" ***.obj, ・・・

今度は、「Writing R Extensions」の第4章にあるように関数を書いてみました。

#include <R.h>
#include <Rinternals.h>

extern "C" __declspec(dllexport) void R_TEST001(int *x, int *ret);
extern "C" __declspec(dllexport) SEXP R_TEST003(SEXP x, SEXP ret);

__declspec(dllexport) void R_TEST001(int x[], int *ret)
{
  *ret = x[0]+x[1]+x[2]+x[3];
}
__declspec(dllexport) SEXP R_TEST003(SEXP x, SEXP ret)
{
  int i, j, nx, ny;
  SEXP ans;
  
  nx = length(x);
  PROTECT(ans = allocMatrix(REALSXP, nx, 1));
  UNPROTECT(1);
  return(ans);
}

このようにすると、コンパイルは通るのですがリンクができません。

ilink32 /w /m /Tpd /L"***?Lib" "C:?borland?bcc55?Lib?C0D32.OBJ" ***.obj,
   ***.dll, CW32.LIB IMPORT32.LIB
Turbo Incremental Link 5.00 Copyright (c) 1997, 2000 Borland
Error: 外部シンボル 'Rf_length(SEXPREC *)' が未解決(***.OBJ が参照)
Error: 外部シンボル 'Rf_allocMatrix(unsigned int, int, int)' が未解決(***.OBJ が参照)
Error: 外部シンボル 'Rf_protect(SEXPREC *)' が未解決(***.OBJ が参照)
Error: 外部シンボル 'Rf_unprotect(int)' が未解決(***.OBJ が参照)

** error 2 ** deleting ***.dll

R.DLLから「implib.exe」でR.libを作成してそれを指定しても結果は同じでした。
このようなエラーが出るということは、致命的な間違いを、犯しているのでしょう。が、マニュアルを読んでも良くわかりません。
解決につながるような情報がございましたら、よろしくお願いいたします。

  • 関数が「Rprintf("Hello!?n");」のみだったら通りますね。 -- どうむ 2004-04-20 (火) 09:53:23
    > ret <- .Call("R_TEST003",as.double(c(1,2,3,4)),ret=double(1))
    Hello!
    > ret
    [1] 1 2 3 4
    ということは、
    length
    allocMatrix
    PROTECT
    UNPROTECT
    が、いけないことになりますね。適切なヘッダーを指定していないのでしょうか?
    なにか、大きな過ちを犯しているような。。。
    何を見落としているのでしょう???
  • 自己レスばかりですがなんとなくわかりました。 -- どうむ 2004-04-21 (水) 09:28:23

    Rprintfが定義してある「Print.h」には

    #ifdef  __cplusplus
    extern "C" {
    #endif
    とあって、C++でコンパイルされることを想定しているようですが「length」等が定義してある「Rinternals.h」にはそれらがありません。そのために、C++形式でエクスポートされていたようです。
    ですので、C++ではなくCだと思わせてコンパイルをかけることで解決しました。
    「Rinternals.h」に「extern "C"」を書けばC++でもコンパイルできるとは思うのですが、どこに影響が出てくるかわからないのでそのままにしておきます。。。
    同一ファイルでC++のクラスが使えないのは不便かなぁ。
  • こうした使い方ができることは噂には聞いていましたが、実際試みられた人は日本では初めてではないでしょうか。後進の参考になると思います。 -- 間瀬茂 2004-04-21 (水) 11:01:59
  • これにあたってはRから他言語利用を、参考にさせていただいておりました。今回のポイントは、拡張子が「cpp」だったのでC++でコンパイルされた、ということなのでしょう。CとC++をそんなに厳密に考慮していなかったのが原因です。拡張子が「cpp」でもCだと思わせてコンパイルするか、そもそも拡張子を「.c」にしてしまえばよかったようです。その場合もちろん「extern "C"」は必要ありませんが。それにしても、提供されているヘッダーファイルが一方はC++対応でもう一方は考慮されていないというのは仕様なのでしょうか?時間のあるときに調べてみます。このままではC++でクラスを用いたコードはRで使用できますが(Writing R Extensions,Chapter4,「4.6Interfacing C++ code」,P35)、そのクラスの中でRのAPIは使えないということになりますね。あれ、別のファイルに書けばいいのかな? -- どうむ 2004-04-21 (水) 14:35:00

Mac OS X の作図デバイス

舟尾? (2004-04-18 (日) 15:49:17)

R を起動しますと,UNIX では x11 が,Windows では windows が,Mac OS 8/9 では macintosh がそれぞれ自動的に起動します (この情報も古い?) が,Mac OS X でも macintosh が起動するのでしょうか?あと,Linux 系の R でも x11 が起動するのでしょうか?興味本位の質問ですが,お答え頂けたら幸いです.

  • Mac OS X では、x11 が起動することも、Aqua (Mac OS X のグラフィックデバイス)が起動することもできます。R-1.9.0/src/modules/aqua/AquaReadMe?.txt に記載されているように、configure のオプションでどちらを使うか、両方使うかを選択できます。 -- なかお? 2004-04-18 (日) 21:53:05
  • 回答ありがとうございます.なるほど,Aqua を起動させることも出来るのですねぇ. -- 舟尾? 2004-04-19 (月) 00:35:09
  • Linux系の場合、Rの起動コマンドオプションに --gui=X11(デフォルト)か--gui=gnomeまたは--gui=noneを指定します。 -- 谷村 2004-04-20 (火) 08:26:55
  • 回答ありがとうございます.勉強になりました.基本はどのOSでも x11 がデフォルトなのですね. -- 舟尾? 2004-04-20 (火) 20:25:23
  • バイナリのレベルでは、RAqua.dmg (1.9.0ではR.dmg)からインストールしたバージョンではquartzデバイスがデフォルト、finkからインストールすればX11がデフォルトです -- 2004-04-30 (金) 13:22:07

GLMM()でできるオブジェクトの参照方法

林 (2004-04-14 (水) 13:38:20)

lme4パッケージのGLMM()でできるオブジェクトは
is.list()で調べるとリストなのですが、names()では名前はなく、[[1]]とリストの添え字で参照しようとしても参照できません。str()の出だしは

list()

となります。通常のリストでは

List of 5 

などといくつ要素があるとか返ってくるので、GLMM()でできるオブジェクトはなにか変わっているのでしょうか?

1.9.0がmakeできない

谷村 (2004-04-13 (火) 00:42:16)

VineSeed?でR-1.9.0がmakeできません。

making dataentry.d from dataentry.c
making devX11.d from devX11.c
making rotated.d from rotated.c
making rbitmap.d from rbitmap.c
make[4]: ディレクトリ `/home/hoge/tmp/R-1.9.0/src/modules/X11' から  出ます
make[4]: ディレクトリ `/home/hoge/tmp/R-1.9.0/src/modules/X11' に入ります
gcc -I. -I../../../src/include -I../../../src/include  -I/usr/X11R6/include
  -I/usr/local/include -DHAVE_CONFIG_H  -D__NO_MATH_INLINES -mieee-fp -fPIC  
  -g -O2 -c dataentry.c -o  dataentry.lo
In file included from dataentry.c:31:
/usr/X11R6/include/X11/Xlib.h:1400: error: 文法エラー before "_Xconst"
/usr/X11R6/include/X11/Xlib.h:1488: error: 文法エラー before "char"
/usr/X11R6/include/X11/Xlib.h:1516: error: 文法エラー before "_Xconst"
/usr/X11R6/include/X11/Xlib.h:1520: error: 文法エラー before "char"
/usr/X11R6/include/X11/Xlib.h:1542: error: 文法エラー before "_Xconst"
[snip]

誰かVineSeed?でmakeに成功している方はいらっしゃいますか?
XはXOrg-6.7.0-0vl2です。

  • コンパイルが通りました。dataentry.cの#include <X11/X.h> の上方に#include <stdio.h>を書いて、#define NeedFunctionPrototypes? 0 をコメントアウトします。 -- 谷村 2004-04-21 (水) 00:10:46

RjpWiki 頁のアーカイブのうまい作り方?

間瀬茂 (2004-04-04 (日) 23:06:46)

例えばこのQ&Aコーナーがかなり長くなってきたので、アーカイブを作り、古いものを移したいのですが、添付画像ファイルが多いとどうしたら一緒に引越ししたら良いのか?いっそ「Q&A(2)」という名前で新しい頁を立てれば簡単ですが、それもかっこ悪い。うまい方法をご存知の方はお教え下さい。そもそも既存の頁の名前をリンク関係を保存したまま変えることはできないのですよね?

SJava

高階知巳? (2004-03-31 (水) 10:56:07)

S/RとJavaのインタフェース SJava を、1年半ほど前に試したときは、動作する環境が非常に限定されている感じで、ちゃんと動作しなかったのですが、最近、どうなっているんだろうと思って omegahat のウェブサイトを見に行こうとしたところ、昨日も今日も繋がらないですね... どなたか、SJavaをうまく使えたという方はいらっしゃいますか?

  • SJavaうまくいきました。omegahatは今は大丈夫なようで、しかも、SJavaが0.67-3にアップデートされています。次のようなJavaプログラムでRをエンジンとして呼び出せます。便利でしょ? -- 高階知巳? 2004-04-05 (月) 17:05:20
import org.omegahat.R.Java.*;

public class REvalSample {
   public static void main(String [] args) {
       String [] rargs = {"--slave", "--vanilla"};

       System.out.println("Sample program to call R engine from Java");
       ROmegahatInterpreter interp
           = new ROmegahatInterpreter(ROmegahatInterpreter.fixArgs(rargs),
                                      false);

       REvaluator e = new REvaluator();
       
       Object val = e.eval("x <- sin(seq(0, 2*pi, length=30))");
       val = e.eval("x * 2.0");
            
       if (val != null) {
           double[] objects = (double[])val;
           for (int i=0; i<objects.length; i++) {
               System.err.println("("+i+") " + objects[i]);
           }
       }
   }
}
  • R-1.9.0+SJava-0.67.3でコンパイルは成功しました.が,実行時に,libRSNativeJava?.soのロードでエラーをはきます.>undefined symbol:STRING_ELT っておこられます.原因がわかる人はいらしゃらないですか. -- riki? 2004-04-28 (水) 16:55:22

pairs(data.frm, upper.panel = points, lower.panel = 残差プロット)

みやむら? (2004-03-26 (金) 13:42:27)

pairs() の描画のようにして,上三角部分に通常の散布図を,下三角部分に重回帰による残差の散布図を描きたいのですが,適当な方法はあるでしょうか?

  • pairs用 回帰直線つきpanelを参考に対応する lower.panel 用の関数を定義して、指定すれば良いのでは。 -- 2004-03-26 (金) 14:52:04
  • 単回帰の残差であればよいのですが,この場合重回帰ですので,そんなに単純ではないと思います.lower.panle に渡される引数についての情報があればよいのですが -- 2004-03-27 (土) 14:29:51

help.startができなくなりました。

(2004-03-26 (金) 13:15:11)

help.startでブラウザが起動しなくなりました。どこをチェックしたら良いのか
ご存知の方は教えてもらえないでしょうか?
こんなコメントが出ます。Making links in per-session dir ...If /usr/bin/open is already running, it is *not* restarted, and you must switch to its window. Otherwise, be patient ...
使用ソフトRAqua ver1.8です。

  • Macintosh HD/private/tmp/Rtmp988/.R/doc/html/index.html を開くアプリケーションがブラウザ以外のものに設定されていませんか。ただし,このファイルはGUI ではパスをたどって開けないので,やっかいです。本来ブラウザが開くファイルの「情報を見る」で,「このアプリケーションで開く」がそのブラウザになっていることを確認してください。もし違うものになっていたら,ブラウザが開くように設定して,「全てを変更」をクリックしてみてください。うまくいくといいですが。 -- 2004-03-26 (金) 19:36:46
  • GUI で見えないのでどうでもいいのですが,Rtmp988 みたいなディレクトリ名で,数値の部分は固定ではないみたいです。 -- 2004-03-26 (金) 19:42:35
  • 私の場合だと,その index.html ファイルを jedit で開くように設定を変えているので,help.start() すると,Jedit で開かれてしまいます。Jedit は,command キーを押しながらタイトルバーをクリックすると,そのファイルを含むフォルダを選択して開くことができるので,フォルダ中のファイルを選択して「情報を見る」で,開くアプリケーションを変更できます。ほかの場合はどうなっているのか分かりません。 -- 2004-03-26 (金) 22:15:02
  • コメントどうもです。アドバイスどおりやってみましたが、、うまくいきません。open /private/../index.html -- しん? 2004-03-26 (金) 22:46:56
  • open index.htmlだとちゃんとひらくのですが、、この前OpenOffice?をいれたのがまずかったのかも、、 -- しん? 2004-03-26 (金) 22:48:54

2変量のヒストグラム

sak? (2004-03-24 (水) 23:53:51)

2変量のヒストグラムを描く事は可能でしょうか?

  • 何種類かパッケージがありますが,その一例はここで紹介されています. -- 舟尾? 2004-03-25 (木) 00:09:17
  • うーん、なんでもできそうですね。すごい。これからRの世界にはまっていく事にいたします。ありがとうございました。 -- sak? 2004-03-25 (木) 09:49:27

[質問]DOSの |more や /p にあたるコマンドは?

inagaki00? (2004-03-23 (火) 22:27:03)

 長ーいデータを見るとき、DOSでは

type hogehoge.txt | more とか

type hogehoge.txt /p のように打ちます(UNIXだとcatですね)
これと同じスイッチがRにはありますでしょうか。

ありそうで見つからないので困っています。
ご存知の方、どうぞよろしくお願い致します。

  • 普段は x[1:100] のように一部だけ表示しますのでおっしゃるような機能の必要をあまり感じませんでしたが、たしかにあると便利ですね。あると便利な機能が「探せば必ずある」のが R のいいところです。早速定石の help.search("pager") で調べると案の定あった! page(x) を試してみてください。おかげでまた一つ賢くなりました。ちなみに page() 関数で起動されるページャ(help 関数で使われている)は options()$pager で決められており、既定値は "/usr/lib/R/bin/pager" (Debian GNU/Linux の場合)です。これをお好みのページャに変えることも可能です。-- 2004-03-23 (火) 22:40:08
page(base)                      Invoke a Pager on an R object
ObjecttaskCallbackManager(base) Create an R-level task callback
managergrid.newpage(grid)       Move to a New Page on a Grid 
Devicetkpager(tcltk)            Page file using Tk text widget
  • なるほど!Pagerというんでしたか。 早速のご回答有り難うございました。helpはこう使うのですね、早速いろいろ試してみます。 -- inagaki00? 2004-03-23 (火) 23:47:41

ページのヘッダやフッタ描画、またはページの4隅のユーザ座標を得るには?

Я・・R? (2004-03-23 (火) 18:38:15)

レイアウト、グラフ、余白設定などに関係なく、ページの特定位置にヘッダやフッタを描画したいのですが、mtext(), title() などでは描画内容や余白設定によって位置が動いてしまいます。
ヘッダやフッタを text(ユーザ座標) で描こうと思ったのですが、4隅の座標を得る方法がわからず、少なくとも par() の usr, plt, fig の値から逆算してもレイアウトによって定まらないことがわかりました。
ところが

box(which="outer", lwd=4, lty=2, col="blue")

とやると必ず正しく4隅を囲むので、どうやって4隅の座標を得るのかを見ようと

?box

としましたがソースコードが見れませんでした。
ヘッダやフッタ描画、またはページの4隅のユーザ座標を得る方法をご存知でしたらよろしくお願いします。

  • identify、locator 関数でインタラクティブに指定位置の座標値を得たり、オブジェクトをおく場所を決めるというのではまずいですか? -- 2004-03-23 (火) 18:59:09
  • はい、ページ数が非常に多いので、マウス入力はちょっと無理です。逆にいうと、ページ数が多いので、ヘッダやフッタ・ページ数などの印字が必要になってくるともいえます。この方法がない場合ですが、1ページ分の作図後または作図前にレイアウトをいったん標準状態に初期化して見えない図を1つ描いてから標準の特定位置にヘッダやフッタを描くようにすれば、つねに同じ位置に描くことができそうです。でも、もっとちゃんとした方法があるような気がするので・・。あるいは、box関数のソースが見れればと思うのですが・・ -- Я・・R? 2004-03-24 (水) 08:48:35
  • フリーソフトですから、ソースは(その気にさえなれば)見るのは簡単です。とはいえ、私もそこまでやるのはひるみますが。たしかにグラフィックス画面を完全に自前でコントロールするやりかたを私もいつか知りたいものと思っていますが。 -- 2004-03-24 (水) 09:29:47
  • 一応解決しました。1ページ分を作図した後に par(new=T, cex=1, omd=c(0.05,0,0.05,0)) などとして上下にすき間を再設定し、title(sub="これはフッタです", cex.sub =.6, line=0, adj=1, xpd=NA, outer=T) などとやるととりあえずはフッタが描けました。ありがとうごさいました。 -- Я・・R? 2004-03-24 (水) 11:59:39

ESSの中の日本語

(2004-03-23 (火) 11:33:40)

gnome-terminal上では

> a <- "日本語"
> a
[1] "日本語"

と日本語が使えますが、ESS上では

ess-20040323.png

という感じになります。gnome-ternimal上では問題がないのでESSの設定で
何とかなればと思っています。ご助言などがございましたらよろしくお願いします。
[環境]
VineLinux2.6
ess-emacs-5.2.0beta1-1vl1
R-1.8.1-4vl1(本家版)

ファイル書き出し時の文字コードの指定方法

matsu? (2004-03-18 (木) 19:06:00)

はじめまして、matsuと申します。

早速ですが、write()やwrite.table()でオブジェクトをファイル
に書き出す場合、ファイルの文字コードを制御することはできるの
でしょうか。

私の環境( R ver1.6.1 Win2k )では sjis(改行コードCRLF)で
書き出されるのですが R側の設定でeuc(改行コードLF)に変更
はできるものなのでしょうか。

ポインタ等でもお教えいただければ幸いです。
よろしくお願い致します。

  • オリジナルの R は当然日本語コードなど知るはずもないからこれは無理な注文では。中間バージョンでも無理?SJISファイルを Unix で nkf -e foosjis.txt > fooeuc.txt とすれば済むことですね。 -- 2004-03-18 (木) 21:36:46
  • 中からperl(WindowsならActivePerl?)を呼び出すことにしてjcode.plとかJcode.pmとかEncode.pmとかと組み合わせた変換関数を書けば,原理的にはできると思います。この情報がヒントになるかもしれません。でも,逆にいえば,オリジナルのRは日本語について何も操作しないので,プログラムやデータ中の文字列をもともとEUCで作っておけば,出力もEUCのままになるはずです。改行コード制御をするにはwrite()やwrite.table()ではなくて,save()を使うしかないと思います。または出力そのものをperlでさせてしまうとか。 -- 中澤? 2004-03-18 (木) 23:30:03
  • 中澤コメントに補足。例えば write 関数には file="|cmd" の形のオプションがありますね。ここで cmd は任意の(R が理解できる)システムコンマンド名です。これをファイルのコード変換と、適当なファイルへの出力を組み合わせたスクリプト名にすればもしかしてうまくいくかも知れません。 Unix ならともかく、MSW ではこうしたやりかたがうまくいくのか私は知りません。 -- 2004-03-19 (金) 07:11:24
  • みなさまありがとうございます。ソースをEUCで記述したり、writeの書き出しをpipeでつなぐ方法は知りませんでした。もともとwin版Rで吐き出したhtmlをlinuxにアップする行程でのことでしたのでおとなしくlinux上でnkfで変換することに致しました。本当にありがとうございました。 -- matsu? 2004-03-22 (月) 10:37:12
  • 以前、save(file="hoge.txt",ascii=T,list=hoge)ですと、ascii以外はオクタルで出力されてたので日本語版ではそのまま出してみたんですが、実はバイト長もあったのでnkfで変換しても使えない事がわかり、オクタル出力に戻しました。(^^;英数字だけの流通ならうまくいくと思います。 -- なかま 2004-03-23 (火) 09:24:32

グラフィックウインドウの名前

(2004-03-18 (木) 11:20:43)

windows()などで表示されたウインドウ
に自動的にふられる名称
「R Graphics: Device x」
を変更することは可能でしょうか?

  • Rソースコードのsrc/modules/X11/devX11.c の中で、XライブラリのXChangeProperty?関数でハードコーディングされているので、ソースコードから修正しないと無理のようですね。逆にソースコードを修正すればどのようにも変更できそうです。 -- 岡田 2004-03-18 (木) 15:48:50
  • そうですか。作成グラフ数が多いときにできれば便利です。上位バージョンに期待します。 -- 2004-03-18 (木) 16:23:06

非線形回帰の初期値

ぶんしょ? (2004-03-17 (水) 17:40:55)

非線形回帰の初期値がすぐ入らないのですが、いい方法があれば教えてください。 今やっているのが

tpot<-nls(HAV~(a*SI+b)*(1-exp(k*AGE))^m,start=list(a=1,b=-10,k=-0.005,m=1.5),data=plot)

という式です。 見当をつけて初期値を入れれば大体はすぐ当てはまるのですが、かたよりが大きい場合、

Error in numericDeriv(form[[3]], names(ind), env) : 
Missing value or an Infinity produced when evaluating the model

Error in nls(HAV ~ (a * SI + b) * (1 - exp(k * AGE))^m, start = list(a = 0.5,  : 
singular gradient

Error in nls(HAV ~ (a * SI + b) * (1 - exp(k * AGE))^m, start = list(a = -50,  : 
step factor 0.000488281 reduced below `minFactor' of 0.000976563

といったエラー文が出ます。これらの意味もできれば教えてください。

  • 最適化の初期値の選び方についてはチップス 関数の最大・最小化 を見てください。全然参考にならないと不満かも知れませんが。試行錯誤(そして幸運)以外に王道はありません。エラーメッセージは簡単にいえば計算途中に「とんでもない値が現れてギブアップ」ということでしょうか。これらの意味を詮索するよりは、より良い初期値を試行錯誤する方が生産的かと思います。 -- 2004-03-17 (水) 18:15:49
  • 四つの変数を適当な範囲と、刻目で変えながら、非線形誤差自乗(二乗,2乗)和の値を計算して、一番小さな値を与える変数値を初期値にして実行してみるというのが、一番素直でしょうか。関数 optim で利用可能なさまざまなアルゴリズムを順に試してみるのもありかも知れません。初期値の選択にある程度頑健な Nelder-Mead 法や SA 法は、困難な非線形最適化問題の最初のステップにはお勧めです(これで見つかった値を初期値にして再試行してみる) -- 2004-03-17 (水) 18:21:01
  • ありがとうございます。やはり神頼み的な思考錯誤しかないようですね。optim 関数もなかなか当てはまりません。根気良くがんばってみます。 -- ぶんしょ? 2004-03-18 (木) 13:06:41
  • 誤差の入り方を無視し log(HAV) ~ log(a*SI+b)+m*log(1-exp(k*AGE)) で一次的当てはめというやりかたもあるかも知れません。 -- 2004-03-18 (木) 17:05:05
  • すみません勉強不足なんですが、これは対数で回帰させるのですか? -- ぶんしょ? 2004-03-19 (金) 02:53:06
  • そうです。単なる思いつきです。目的変数を log(HAV) にとります。(誤差の入り方を無視すれば)式そのものは同値ですが、数値的振舞はしばしば異なる(良くなるとは限りませんが)ことがあります。もし運良くそれらしい値が求まれば、それを初期値に本来のモデル式で再挑戦したら...ということですが、「藁をも掴みたい」時ならやってみる価値はあるかもしれません。こうした元の式を変形した式で試すというのも試行錯誤の定石の一つです。 -- 2004-03-19 (金) 07:03:36
  • 何度もすみません。対数で試したら
    Error in numericDeriv(form[[3]], names(ind), env) : 
           Missing value or an Infinity produced when evaluating the model
    In addition: Warning messages: 
    1: NaNs produced in: log(x) 
    2: NaNs produced in: log(x)
    といったエラー文が出るのですが何がいけないのでしょうか? -- ぶんしょ? 2004-03-19 (金) 08:19:17
  • これは当然 1-exp(k*AGE)) が試行中に非正値になったということでしょう。期待させたら申し訳ありませんでいたが、やはりうまくいかないということでしょうか。k の値を負値に制約する必要がありますが、思いつきですからそこまでやる価値があるかどうか。当然、この場合も適当な初期値の選択は不可欠ですが。元の式でも k は負の値という自然制約があると思いますが、制約付きの最適化は試されましたか。-- 2004-03-19 (金) 08:33:40
  • もし差し支えなければ、そしてデータがあまり大量でなければ、一例を紹介されるともう少し突っ込んだコメントが得られるかも知れません。ただし、もともと簡単なケースとは思われませんから、期待されても困りますが。 -- 2004-03-19 (金) 08:46:00
    AGE	HAV	     SI
    18 	197.98	 21.48
    20 	266.02	 23.78
    24 	695.26	 39.88
    26	429.21	 25.14
    27 	312.04	 22.79
    33 	484.55	 24.15
    35 	708.11	 30.79
    25 	743.29	 35.23
    28 	714.40	 32.50
    29	377.93	 19.36
    29 	684.48  29.50
    42 	557.06	 20.37
    47	483.70	 16.62
    50 	725.47	 20.88
    24 	400.58	 33.11
    25 	458.79	 32.69
    25 	296.57	 24.12
    27 	382.04	 27.50
    32 	343.38	 21.45
    33 	450.72	 26.33
    33 	375.25	 22.31
    34 	420.39	 21.82
    36 	384.44	 20.51
    が当てはまらない一例なんですが、初期値a,m=正、b,k=負の値を取ることが多いです。どんな感じでしょう? 
    • ぶんしょ? 2004-03-19 (金) 09:06:58
  • 少し解析をしてみましたが、やはり最初に注釈があったように、バラツキ(特に HAV の値)が大き過ぎるというのが諸悪の根源のような気がします。試しに HAV 値が600以下のデータだけで当てはめてみると、決して良いとはいえませんが SANN 法で収束はしますね。またその時単に HAV ~ AGE + SI で線形回帰してみると修正済み R 自乗値は0.7835となりました(おそらくモデル式は何らかの意味を持つのでしょうから蛇足でしょうが)。 待てよ、k*AGE の値が小さければ (1-exp(k*AGE))^m は (-k*AGE)^m で近似されるから全く阿呆な考えでもないかも知れません。実際元データでこの近似を行なうと Nelder-Mead 法では収束しました。ただしその値が役立つとは思えません。もう少しましなデータの場合、初期値の検討をつけるのに使えるかも知れません。-- 2004-03-19 (金) 12:23:15
    > optim(c(0.1, 0.1, -0.1, 0.1), SF)
    $par
    [1] 10.907251 -1.975936  -0.057869  1.013191
    $value
    [1] 403297
    $counts
    function gradient
         313       NA
    $convergence
    [1] 0
    $message
    NULL
  • わざわざありがとうございます。確かに粗悪なデータです。自分ももっと勉強したいと思います。 -- ぶんしょ? 2004-03-19 (金) 15:14:58

lm()関数の不思議

杉本? (2004-03-17 (水) 16:31:34)

x<-c(8,3,4,2,1,4,8,4,9)
y<-c(4,1,5,1,1,2,5,2,3)

として x=a*y+bの形で線形回帰したいときは

lm(x~y)

ですが x=a*y^2+bの形で線形回帰するときはどうすればよいのでしょうか?

lm(x~y^2)

とすると上と同じ結果になってしまいます。もちろん x=a*y+b*y^2+cの形で線形回帰するときも

lm(x~y+y^2)
Call:
lm(formula = x ~ y + y^2)
Coefficients:
(Intercept) y 
1.586       1.197  

となり、x=a*y+bの形で線形回帰してしまっています。新たにデータとして作り直せばよいのですがlmでなんとかなる方法をご存知でしたら教えていただけないでしょうか。

  • Rの統計解析関数Tips を見て下さい。R のモデル公式では "*", "+", "^" は特殊な意味を持つので、本来の意味で使いたい時は例えば I(y^2), I(x*y) 等と表します。I(.) は as.is 関数と呼ばれます。?I で確認してください。 例えば x ~ y + I(y^2) としてみたらどうなりますか。-- 2004-03-17 (水) 16:35:20
  • ありがとうございます。うまくいきました。 -- 杉本? 2004-04-08 (木) 15:01:00

無題

大熊 大郎? (2004-03-12 (金) 23:58:36)

関数 width.SJ はなんですか?おしえてください。

  • 悪い質問の典型ですね。標準パッケージ以外のパッケージ中の関数は、該当パッケージを引用しないと?です。help(width.SJ) でオンラインヘルプは御覧になりましたか。MASS パッケージは Venable & Ripley の本のコンパニオンパッケージですから、原著に当たられるのが一番確実で間違いがないです。日本語訳も出てます。リンク集を参照してください。 -- 2004-03-13 (土) 00:07:20

ヒストグラムを描くhist()

Rin? (2004-03-04 (木) 11:57:28)

hist()を使いヒストグラムを描くとき、matrixデータを与えると列データに対し、色分けして描いてくれますが色分けの中身(データのタイトル)はどのようにしたら表示できるのでしょうか?どなたか知っておられる方がおりましたら、教えてください。
よろしくおねがいします。

  • 意味がいまひとつ理解できません。具体例をあげていただけませんか。hist の引数は数値ベクトルで、行列を与えるとベクトルに強制変換されると思うのですが。もしかして barplot の間違い?-- 2004-03-04 (木) 15:50:00
  • 上の方が仰るとおり,hist()では多次元データはベクトルに強制変換されますね.barplot() ならば引数 col に色の種類を指定することが出来ます. -- 舟尾? 2004-03-04 (木) 21:47:41
    > x <- matrix(1:24,6,4)
    > barplot(x)
    > barplot(x,col=rainbow(24))
    > barplot(x,col=rainbow(12))
    > barplot(x,col=rainbow(6))
    > barplot(x,col=rainbow(3))
  • 理解しにくい質問の仕方をして申し訳ありません。凡例はどのようにしたら表示できるのでしょうか。色の種類の指定の仕方、参考になりました。ありがとうございます。 -- Rin? 2004-03-05 (金) 09:57:23
  • ここのような凡例でしょうか? -- 舟尾? 2004-03-05 (金) 22:48:14
  • はい。そうです。早速、legend()を使い試してみましたが、やはり凡例のつけ方がわかりませんでした。hist()とlegend()を使い、どのようなスクリプトを書けば凡例がつけられるのでしょうか。 -- Rin? 2004-03-08 (月) 17:30:26
  • histogram につける気の利いた凡例が思いつかなかったのが申し訳ないのですが,とりあえずこんな感じでしょうか?-- 舟尾? 2004-03-08 (月) 18:30:26
> data(islands)
> hist(islands)
> legend( 5000,y=30,paste("data",c(1:3)), col = c(1:3), lty=1, merge = TRUE, bg='gray90')
> legend(10000,y=30,paste("data",c(1:3)), col = c(1:3), lty=1, merge = TRUE, bg='gray90')
> legend(10000,y=40,paste("data",c(1:3)), col = c(1:3), lty=1, merge = TRUE, bg='gray90')
hist.png
  • 凡例をつけることができました!やっと解決できて嬉しいです。有難うございました。 -- Rin? 2004-03-09 (火) 17:27:34
  • いえいえ.?legend や help(legend) でさらなる情報が得られますので,良かったらご覧下さい. -- 舟尾? 2004-03-11 (木) 03:01:26
  • 初めまして。matlab初心者の者です。棒グラフで値を表示したいと思っています。 -- kuro? 2006-10-27 (金) 18:09:01
  • 続き)matlabでは棒1本1本の色をかえることができずその方法があれば教えてください!!お手数をおかけして申し訳ありませんが、どうぞよろしくお願いします。 -- kuro? 2006-10-27 (金) 18:11:37

帯グラフ

(2004-03-02 (火) 12:57:28)

グラフィックス参考実例集には、なぜか帯グラフがありません。帯グラフを描くときは皆さんどうしていらっしゃいますか?barplotをいじらないとだめですか?

  • リンク集、グラフィックス参考実例集にもリンクが張られている Statistique avec R に barplot 関数を使った例があります。コツはデータを列次元=1の行列(縦ベクトル)とすることのようです。barplot 関数はもともと行列データを与えられると、各列毎に帯グラフからなる棒を書くように設計されています。ただし、参考コードのままでは正方形になりますから、作図領域をあらかじめ長方形にしておかないと帯グラフらしくなりません。 -- 2004-03-02 (火) 13:57:34
    > data(HairEyeColor)
    > x <- apply(HairEyeColor, 2, sum)
    > par(fin=c(6,2))  # グラフィックス画面を 6x4 インチに指定
    > barplot(as.matrix(x), horiz=TRUE, col=rainbow(length(x)), legend.text=TRUE)
    obigraph.jpg
  • ありがとうございます。私のイメージとはちょっと異なるので、 Statistique avec R のリンクを拝見させていただきました。どうやらt()を使うとよいようでした。解決の糸口を与えてくださり、感謝します。 -- 2004-03-02 (火) 17:07:07
    > data(HairEyeColor)
    > a <- as.table( apply(HairEyeColor, c(1,2), sum) )
    > b <- a / apply(a, 1, sum)
    > barplot(t(b))
    obi.png
    > d.seg.xs <- rep((0.2 + 1) * 1:(ncol(t(b)) - 1), each=nrow(t(b)))
    > d.seg.xe <- rep((0.2 + 1) * 2:ncol(t(b)) - 1, each=nrow(t(b)))
    > d.seg.ys <- apply(t(b)[,1:(ncol(t(b)) - 1)], 2, cumsum)
    > d.seg.ye <- apply(t(b)[,2:ncol(t(b))], 2, cumsum)
    > segments(d.seg.xs, d.seg.ys, d.seg.xe, d.seg.ye)
    obi2.png
  • おぉ、horiz=TRUEなんてあるのですか。早速、使ってみました。 -- 2004-03-02 (火) 17:45:57
    > barplot(t(b),horiz=TRUE)
    > d.seg.ys <- rep((0.2 + 1) * 1:(ncol(t(b)) - 1), each=nrow(t(b)))
    > d.seg.ye <- rep((0.2 + 1) * 2:ncol(t(b)) - 1, each=nrow(t(b)))
    > d.seg.xs <- apply(t(b)[,1:(ncol(t(b)) - 1)], 2, cumsum)
    > d.seg.xe <- apply(t(b)[,2:ncol(t(b))], 2, cumsum)
    > segments(d.seg.xs, d.seg.ys, d.seg.xe, d.seg.ye)
    obi3.png
  • あとは、0.0-1.0ではなくてパーセント表示にしたいのですが、100をかけてtext()で(%)を追加するしかないでしょうか。 -- 2004-03-02 (火) 17:52:06
  • % 表示はやはり手作業技でしょうね。ところで元記事質問者の方の希望は、帯グラフを並べたグラフだったのでしょうか。 -- 2004-03-02 (火) 18:21:35
  • 名乗らなくてすみません。元記事質問者です。上の3つの帯グラフは、回答にあったリンクを参考に私自身が私の希望通りに作成しました。%の件は了解です。 -- 元記事質問者? 2004-03-02 (火) 19:06:37
  • barplot()の中でaxes=Fとしておいて,barplot()後にaxis(1,seq(0,1,by=0.2),seq(0,100,by=20))とすれば100をかけた値表示になります -- 2004-03-02 (火) 20:42:42
  • ありがとうございます!なるほど。%が欲しかったので、axis(1,seq(0,1,by=0.2),paste(seq(0,100,by=20),"%"))とすると望み通りになりました。多謝! -- 元記事質問者? 2004-03-03 (水) 21:57:17
  • 質問者の分際でなんなのですが、しかも自力ではないのですが、自作グラフィックスへこの帯グラフを書き込んでおきます。 -- 元記事質問者? 2004-03-03 (水) 22:01:44

Rが動作するOS

(2004-03-01 (月) 18:58:11)

Rの紹介記事にRが動作するOS一覧を載せたいのですがどこかにそういうものはありますか?
現在の把握状況は下記の通りです。

  • MacOS X、Windows、Linux、FreeBSD、NetBSD、OpenBSD、Solaris、AIXでは完全に動作する。
  • 古いMacOSは古いRなら動く?
  • Linux版のザウルスは、制限がありながらも動作する。

    上記以外にRが動作したOSがありますか?
    LindowsやBeOS上でRは動作しますか?超漢字はどうですか?
  • 2002年2月の r-help 記事に WindowsCE ではビルドできないという話がありました。動かない断定はしていませんが。 -- 2004-03-01 (月) 19:54:58
  • Mac OS 9 以下では,R 1.7.1 までしか動きません。 -- 2004-03-01 (月) 23:41:42

ESSでのコメント行の色

さとう? (2004-02-28 (土) 09:46:19)

Emacs+ESS+Rを利用してます。emacsでRファイルを編集するとき、コメントマーク「#」の数により、その行の色を変更できると便利ですよね。
とりあえずの対処療法として、

M-x hi-lock-mode
M-x highlight-regexp [RET] ### [RET] hi-red-b
M-x highlight-regexp [RET] ## [RET] hi-green-b

などとしてます。(順序が大事)
拡張子がR(またはr)のファイルを開いたときには、いつもこのように色づけするスマートな方法はないものでしょうか?

Rで並列処理計算(グリッド・コンピューティング?)

(2004-02-23 (月) 15:18:03)

 Rで利用できる並列処理のチュートリアルのようなものはあるのでしょうか?

  • 例えば R news 2002/2 号に "Rmpi: Parallel Statistical Computing in R" という記事があります。他にも CRAN の検索機能で grid, parallel 等のキーワードで調べればいくつかの記事が。 -- 2004-02-23 (月) 17:16:24

64ビット環境でのRの動作

(2004-02-23 (月) 15:16:40)

 最近、64ビットのPCが盛んに話題になっていますが、Rは32ビットよりも、64ビット上の方がパフォーマンスが段違いのでしょうか?

  • C言語のdoubleは現在64bitであることが多いですから、doubleの演算に64bit CPUの命令を効果的に使えるコンパイラと数値計算ライブラリの組み合わせのもとでコンパイルしたRはそれなりに早くなるはずじゃないかと思います。PowerMac? G5なら手元にあるけど、gccは64bit対応済んでたかな... -- 2004-02-23 (月) 16:46:18


階級間隔が異なるヒストグラム

さかもと? (2004-02-20 (金) 11:30:56)

Rにて階級区間が異なるヒストグラムを作成したいのですが,どうしてもうまくいきません・・・.間隔dに対してある階級からは間隔fにしたいのですが,
(度数)*d(基準とする標準間隔とでも言いますか・・・)/f
を間隔fの高さにすればいいのですが,どうすればRでそれが可能なのでしょうか?
ちなみに,初心者です・・・.

  • 引数 breaks (vector型)に、binの区切りとなる値を与えれば良いのではないでしょうか? -- テツヲ? 2004-02-21 (土) 09:20:47
  • ありがとうございます。さっそく試してみます。小生が試したところ、引数widthに間隔値を入れたのですが・・・(初心者)。 -- さかもと?? 2004-02-21 (土) 14:56:28
  • 確か引数 width は S の頃に(関数density()などで)使われていたものだったと思われます.あと,度数は「区切り幅を決めた後」に決まる値ですから,さかもとさんの様な処理を行うのであれば「区切り幅を決める」->「度数を求める&関数などを使って得る」->「自分の要望に合う様に区切り幅なり高さなりを決めて図を描く」という処理になるかと思われます.区切り幅が決まった後の「度数を求める&関数などを使って得る」処理は,データを1つずつチェックするか,次のような plot=F と $counts などを使うことで実現できます.
> x <- rnorm(50) # データ
> h <- 0.8
> bins <- seq(min(x)-0.1, max(x)+0.1+h, by=h)
> pointlist <- hist(x, breaks=bins, prob=T, plot=F)$counts
> hist(x, breaks=bins) # 普通にヒストグラムを描く場合(区切り幅はbinに従う)

リストpointlistには,長さ h ごとに区切られた場合の左から数えた度数が入っています.-- 舟尾? 2004-02-21 (土) 20:32:00

免疫アルゴリズム

(2004-02-17 (火) 18:01:23)

R や S-plus で免疫アルゴリズム(IA)を扱った例はあるのでしょうか?

ライブラリのバージョン表示

(2004-02-17 (火) 13:16:00)

ライブラリパッケージが最新バージョンかどうかを確認するため、インストール
しているライブラリのバージョンを表示させたいのですが、方法が分かりません。自分なりに検索してみましたが、方法を見つけることができませんでした。

  • update.packages() で解決しました。 -- 2004-02-17 (火) 14:02:31
  • read.dcf(file=system.file("DESCRIPTION",package="tree")) とかでパッケージの情報の取得が可能です。dcfはdebian control format の事ですが、Winでも使えます。 -- なかま@金沢市民? 2004-02-22 (日) 12:55:23

Rmapやmaptoolsでの縮尺バーや方位記号の表示

(2004-02-17 (火) 00:17:51)

Rmapやmaptoolsで表示した地図には縮尺バーや方位記号を入れられないのでしょうか。

  • arrows などを使って自作するしかないのでは -- 2004-02-18 (水) 01:05:40
  • 凡例はデモにハワイの火山表示のコードが参考になるかも -- 2004-02-18 (水) 14:58:18

plotで縦軸と横軸の単位幅を同じ長さにする方法

(2004-02-16 (月) 19:10:03)

FAQでしたら、削除しください。
x軸とy軸がどちらも単位がm(メートル)であるplotで、x軸の1mとy軸の1mを同じ
長さにしたいのですが、勝手に調整されてそろいません。x軸とy軸の単位長を
同じにする方法はありませんか。

  • xlim と ylim で座標軸範囲を同じにすれば良い?-- 2004-02-16 (月) 19:27:22
plot(runif(10), 2*runif(10),xlim=c(0,2),ylim=c(0,2))
  • ありがとうございました。確かに、座標範囲を指定すると軸とy軸の単位長さが等しくなりました。しかし、無用な空白が図の大部分を占め、何とも美しくありません。余分な空白を作らずに単位長を等しくする方法はありませんか。 -- 2004-02-16 (月) 21:09:23
  • サンプルコードを紹介して頂かないと、コメントのしようがありません。 -- 2004-02-16 (月) 21:41:20
  • 例えば、下のような場合では、上に大きな余白ができます。 -- 2004-02-16 (月) 22:43:37
plot(2*runif(10), runif(10,min=0,max=.1),xlim=c(0,2),ylim=c(0,2))
  • x,y 座標を同じスケールにする以上これは避け難い。グラフ領域の形状そのものを正方形以外にしたいのでしょうか? -- 2004-02-16 (月) 22:49:54
  • 申し訳ありません。最初にいうべきでしたが、x,y 座標が同じスケールの「長方形のグラフ」を描きたいのですが、無理なのでしょうか。 -- 2004-02-17 (火) 00:00:54
  • par(mai=c(1,2,6,2))ってしてみるとか,postscriptデバイスとかpngデバイスならwidthやheightを指定できるので長方形のグラフは書けます -- 2004-02-17 (火) 01:50:49
  • dev.copy(x11,whdth=値,height=値)で書き直させる? -- 2004-02-17 (火) 08:37:07
  • widthやheightを指定する方法は、作図領域の指定ではなく、ディバイス領域の指定なので、xyのスケールを正しく合わせることができませんでした。 -- 2004-02-17 (火) 10:58:26
  • library(grid)を使えばできそうな気がします -- 2004-02-17 (火) 18:32:58
  • 正確な指定の仕方はわかりませんが、次のようにするとそれらしい図が描けます。まず par(fin=c(x,y)) で作図領域のサイズをインチ単位で与えます。次に xlim, ylim で両軸の座標範囲を指定しますが、その際(この例では) ylim に比例定数 3/6 をかけておきます。-- 2004-02-17 (火) 18:53:23
    > par(fin=c(6,3))
    >  plot(6*runif(100),1.5*runif(100),xlim=c(0,6),ylim=c(0,3)/2)

比例定数 3/6 というのは本質的ではなく、実際の画面を見て適当な比率を試行錯誤するというのが現在思いつく解決策です。なお fin パラメータの既定値は fin=c(6.990803, 6.997416) でこれより大きなサイズは指定でき無いようです。par(fin=c(6,3)) を実行する際は

oldpar <- par(no.readonly = TRUE) # 現在のグラフィックスパラメータ退避
... 作業 ...
par(oldpar) # パラメータ復帰

を実行しないと後で困るかも知れません。

#ref(): The style ref(filename,pagename) is ambiguous and become obsolete. Please try ref(pagename/filename)

  • 追加。これでも実際は余計な余白がたっぷり付いている可能性が高いです。最終的にはポストスクリプトファイルに落し、余計な空白を後でクリップする必要がありそうです。 -- 2004-02-17 (火) 19:08:55
  • 大体様子がわかったような気がします。グラフィックス画面の大きさは par()$fin パラメータで決められます(インチ単位)。一方グラフィックス画面中の作図領域の大きさは par()$plt パラメータ(グラフィックス画面への割合で指定)で決められます。従って作図画面の実際の大きさはインチ単位で -- 2004-02-18 (水) 21:09:28
    par()$fin[1]*(par()$plt[2]-par()$plt[1]) # x 方向のサイズ(インチ)
    par()$fin[2]*(par()$plt[4]-par()$plt[3]) # y 方向のサイズ(インチ)
    と決まります。従って x,y 座標のスケールを同じにするには xlim と ylim の幅をこのサイズの等倍にすれば良いようです。参考コードを紹介しておきます。Xlim と Ylim の倍数 10 は任意ですが同じでないとスケールが異なってしまいます。つまり Xlim と Ylim の幅の比率が実際の比率に等しい必要があります。ただし、これでも正確にいうと座標隅に入る遊び分を考慮していませんから、微妙に異なっているはずです。
 > par(fin=c(6,4))  # グラフィックス画面を 6x4 インチに指定
> par(plt=c(0.1,0.9,0.1,0.9)) # 作図領域をグラフィックス画面への割合で指定
> Xlim <- 10*par()$fin[1]*c(0, par()$plt[2]-par()$plt[1])
> Ylim <- 10*par()$fin[2]*c(0, par()$plt[4]-par()$plt[3])
> plot(1:20,1:20, xlim=Xlim, ylim=Ylim)

#ref(): The style ref(filename,pagename) is ambiguous and become obsolete. Please try ref(pagename/filename)

参考までにメモすれば plt, fin の既定値は次のようになっています。

 > par()$plt
[1] 0.1172537 0.9399432 0.1457143 0.8828571
> par()$fin
[1] 6.990803 6.997416

従って既定値での作図画面のサイズはインチ単位でそれぞれ 5.751261, 5.158095 です。

  • お陰様でうまくいきました。ありがとうございました。 -- 元質問者? 2004-02-18 (水) 21:50:43
  • その後、ただ単にasp=1というオプションを付ければよいことが判明しました。
    例)
    > x <- 1:10
    > y <- x * 10
    > plot(x,y,asp=1)
    結果)

    #ref(): The style ref(filename,pagename) is ambiguous and become obsolete. Please try ref(pagename/filename)

正準相関分析での第 i 正準変量の求め方

(2004-02-14 (土) 20:55:31)

library(mva)の関数cancor()を用いれば正準相関係数行列が得られますが、
第 i 正準変量を求める場合はどうすればよいのでしょう?例えば このような分析 を R で行いたいのです.
(exampleでは x %*% xcoef と計算していますが何か違う気がして・・・)


latex形式のヘルプファイルの見方

(2004-02-06 (金) 23:25:05)

libraryごとにhtmlやhelpなどヘルプ文書がありますが、latex配下にあるTeX文書は、documentclassもなくいきなり?Headerで始まっていますが、これはどうすればコンパイルできるのでしょう。helpやhtmlでは数式に難があるので、できればxdviなどでヘルプを読んでみたいです。

  • ソースから R システムを構築するとき make install-dvi とすれば dvi ファイルが一挙にできるようです。また、やはりソースファイル中にあるヘルプファイルの原型の xxx.Rd ファイルを個別に R CMD Rd2dvi xxx.Rd と処理すればやはり dvi ファイルが得られます。しかし、たしかにご質問のように xxx.tex から dvi ファイルを作る方法は? -- 2004-02-07 (土) 07:41:46
  • もしかすると 森さんのWeb「R の日本語文章 (pdf 版)」 にある R-exts.jp.pdf の2章にある記事がお役に立てるかもしれません.-- 舟尾? 2004-02-07 (土) 18:45:19
  • 全く無知で申し訳ないのですが、R CMD Rd2dvi パッケージ名.Rd だと個々のコマンドのhelpが含まれていないようです。latexディレクトリ配下に個々のコマンドのTeXファイルがあるのにもどかしいです。 -- 2004-02-09 (月) 10:24:57
  • 例. library/nlme/latex/nlme.tex だけを make する手順 (R CMD Rd2dvi とまったく同じものができるわけではありませんが...): (1) 作業用のディレクトリー (たとえば, /tmp/r-dir) を作る; (2) $R_HOME/share/texmf/*.sty を /tmp/r-dir 下にコピーする; (3) /tmp/r-dir 下に空の tex ファイル (たとえば, tmp.tex) を作る; (4) tmp.tex の中身は: ?documentclass{book} ?usepackage{Rd} ?begin{document} ?chapter*{} [ここに, library/nlme/latex/nlme.tex の中身を挿入] ?end{document} ; (5) このディレクトリー上で (p)latex tmp を実行する -- yy? 2004-02-09 (月) 10:47:38
  • なお, 上記の手順は, bin/Rd2dvi の内部で実行されていることをマニュアルで再現したものです (bin/Rd2dvi では実際には, makeindex 等, もっといろいろなことをしていますが). -- yy? 2004-02-09 (月) 10:50:54
  • わざわざ latex 形式のヘルプファイルを一揃い用意しているからには、on the fly で必要なヘルプの dvi 化(さらにはおそらく xdvi での表示)が出来るメカニズムが用意されていると考えるのが、順当だと思うのですが。調べてみてもそれらしい記述は見つからないですね。 -- 2004-02-09 (月) 11:20:53
  • 謎が解けました。やはり help 関数には tex 形式のヘルプファイルを要求に応じてその場で処理し、ps ファイルを生成するオプションがありました。?help でオプション offline を見て下さい。必要に応じ環境変数を設定する必要があるようです。なお、ヘルプ関数は途中で出来た dvi ファイルは抹消する仕様になっているようです。これが欲しければ help 関数を適当に改変する必要がありそうです。本来プリントアウトを考えた機能のようです。-- 2004-02-09 (月) 14:00:04
    > options(latexcmd="/usr/bin/latex")
    > options(dvipscmd="/usr/bin/dvips") 
    > help(optim, offline=TRUE)
    Saving help page to 'optim.ps'  # optim.tex から optim.ps が作られた
    思いついて次のようにしたら dvi ファイルのオンライン表示が出来ました。
    > options(dvipscmd="/usr/bin/xdvi")
  • おぉ、すごいです。無事にxdviでコマンドのhelpが読めました。ありがとうございました。素晴らしい。 -- 元質問者? 2004-02-09 (月) 22:16:37
  • さらにdvips を dvipdfm に変えてみました。すると無事にpdfが作成されました。カレントではなく/tmp配下ですが。 -- 元質問者? 2004-02-09 (月) 22:21:25
  • パッケージのhelpと各々のコマンドのhelpをxdviで表示する方法は分かったのですが、それらをすべて1つまとめてdvi化またはpdf化する方法は用意されていないのでしょうか。パッケージに含まれるコマンド全てのhelpを一度に閲覧(印字)できると嬉しいのですが。。 -- 2004-02-09 (月) 22:32:54
  • 課題としては面白いですが、まず気になるのは、例えば base パッケージのオブジェクトすべてのヘルプを印刷するとおそらく数千ページになりそうだということですね。さらに頻繁にバージョンアップされますから悲鳴をあげること確実になりそうです。ですからそうした機能がもともと用意されているとは思えません。小さなパッケージなら問題はないかもしれませんが。それなら個別に印刷してもそれほど手間はとらない。関連して私もまだ把握していませんが、vignette という比較的新しいメカニズムが用意されていて、パッケージの簡略紹介が作成できるようです。 -- 2004-02-10 (火) 07:49:24
  • help() に offline という引数があったとは.勉強させてもらいました! -- 舟尾? 2004-02-10 (火) 21:53:02
  • dvipscmd をだましてやると、dvi ファイルも保存できます。ただし中間ファイルで名前が変。例えば -- 2004-02-11 (水) 18:45:32
    cp ./*.dvi /tmp 
    というスクリプトファイル(名前を仮に makedvi )に実行属性をつけた上で
    > options(dvipscmd="./makedvi")
    としてやる。もう少しスクリプトを工夫する余地?

文字列のオブジェクトへの変換方法 (なんでも掲示版より移動)

sakamoto? (2004-01-30 (金) 19:35:15)

まずは次のプログラムを見てください。

money<-1:10
>a<-readline()
money
>money[a>5]
 [1]  1  2  3  4  5  6  7  8  9 10
> a<-as.name(a)
> a
money
> money[a>5]
 [1]  1  2  3  4  5  6  7  8  9 10

ここでしたい事は、以下の事である。

>money[a>5]
[1]  6  7  8  9 10

つまり、文字列として取り込んだものをオブジェクトとして扱いたい。いろいろ探してみたのですが、見つかりませんでした。どなたか教えてもらえないでしょうか?

  • ご質問の趣旨はこういうことですね。すでに a という名前の長さ10の数値ベクトルが存在し、それをコンマンドラインから指定したい。 get 関数はいかがでしょう。-- 2004-01-30 (金) 21:32:35
> a <- runif(10)     # a という名前の数値ベクトルを作る 
> a           
 [1] 0.6313190 0.4538556 0.3934363 0.4643964 0.3375608 0.3453617 0.3785413
> aa <- readline()   # その名前をコンマンドラインから入力
a                    # 名前は a
> aa
[1] "a"              # aa は文字列 "a"
> get(aa)            # aa という文字列を名前に持つオブジェクトの値を得る
 [1] 0.6313190 0.4538556 0.3934363 0.4643964 0.3375608 0.3453617 0.3785413

まとめてみると

> a <- runif(10)
> b <- runif(10)
> x <- 1:10
> x[get(readline())<0.5]
a                          # コンマンドラインから a を入力    
[1]  3  4  5  6 10         # x[a<0.5] が得られた
> x[get(readline())<0.5]  
b                          # コマンドラインから b を入力 
[1] 1 3 4 5 6 7 8          # x[b<0.5] が得られた
  • money[eval(parse(text=a))>5] としたいと言ってるのかな? -- なかま@金沢市民? 2004-01-30 (金) 23:37:06
  • 中間さんの解を注釈すると -- 2004-01-31 (土) 00:05:55
> a 
[1] 0.557755694 0.756125389 0.062838319 0.116221907 0.287878570 0.490844365 
> y <- readline()
a
> y
[1] "a"
> eval(parse(text=y)) # 文字列 y をオブジェクト名 a に変換し評価した 
[1] 0.557755694 0.756125389 0.062838319 0.116221907 0.287878570 0.490844365
  • ご回答ありがとうございます。おかげさまで問題が解決しました。 -- sakamoto? 2004-01-31 (土) 15:20:57

hist()の引数breaksの決め方

(2004-01-29 (木) 17:57:01)

ヒストグラムを描くhist()の引数のbreaks=seq(min,max,by=bin)の minとmaxはどういった基準で決めればよいのでしょう?

  • データ数が大きい場合(50以上)はあまり気にしなくても良いはずですが,データ数が少ない場合はかなりのバラつきが出てしまいますね。 -- 2004-01-29 (木) 19:31:35
  • 代替案的なものですが,ASH という方法があります.ただ,これを使うならカーネル密度推定を行う方が良いでしょう.詳しくは「ヒストグラムと密度の推定」をご覧下さい. -- 舟尾 2004-02-10 (火) 21:54:48

ベクトルと行列

Soichi? (2004-01-29 (木) 05:39:56)

Matlabから乗り換えようかとRを試し始めました。
ベクトルと行列のどちらも引数として受け取ることができ、
行列を受け取ったときは各列に処理を、
ベクトルを受け取ったときは(n,1)行列に変換した上で、その後の処理を行いたいのですが、
Rでそのようなことは可能でしょうか。

  • as.matrix()をかませておけばいいのでは? -- 2004-01-29 (木) 07:04:03
  • 次のようにすれば良い。問題によりますがいちいち行列に変換しなくてもすむばあいもある? -- 2004-01-29 (木) 08:20:46
    if(is.vector(x)) x <- as.matrix(x)  #もしくは
    if(!is.matrix(x)) x <- as.matrix(x) 
    ちなみに縦ベクトルにしたければ
    if(is.vector(x)) x <- t(as.matrix(x))
    判断する必要がないのなら単に
    x <- as.matrix(x)   # x がすでに行列なら何もしない

R の sleep() 機能

φ(oo-)o? (2004-01-28 (水) 09:11:25)

このサイトの「R プログラミング Tips 大全」のページで

指定秒数プログラムをアイドリングする関数
は標準関数の Sys.sleep(sec) ではどうでしょうか。ヘルプには

This function may not be implemented on all systems.
と書かれていますが、Windows2000では動くようです。
(R 1.8.1を使用しています)

  • 本当ですね。知りませんでした。「R プログラミング Tips 大全」に追加しました。なお、勝手に追加・修正歓迎です。 -- 2004-01-28 (水) 10:17:39

作ったページの削除は?

mnyu? (2004-01-23 (金) 01:08:50)

作ったページの削除は どうしますか

  • 編集モードで内容をすべて削除し、更新すると一人でになくなります。 -- 2004-01-23 (金) 07:24:04

非線形には対応していますか

tia? (2004-01-22 (木) 17:22:56)

Rに関してまったくの初心者です.
Rで非線形問題を解くことは可能でしょうか.
ご教授いただければ幸いです.

  • 問題によりますが、汎用最適化関数 optim や nlm を使えばたいていのことはできます。tips 集を調べて下さい。CRAN パッケージリストを一覧すれば、個別の非線形問題用のパッケージも色々ありますので調べて下さい。 -- 2004-01-22 (木) 17:30:08
  • 早速のお返事ありがとうございます.いろいろと試してみます. -- tia? 2004-01-22 (木) 17:41:04
  • 非線形最小自乗法用関数 nls については Rの統計解析関数Tips 参照。 線形不等式制約付きの最適化関数 constrOptim も参照。-- 2004-02-10 (火) 23:30:59

CCDFについて。

sato? (2004-01-22 (木) 14:24:52)

Rを使って、1‐累積密度分布であるCCDFを描きたいのですが、パッケージ等では用意されていないのでしょうか?
またもしない場合は、パッケージstepfunの中のecdfを使って描く必要があるのでしょうか?

  • stepfun パッケージの ecdf では困る、もしくは避けたい理由があるのでしょうか。 -- 2004-01-22 (木) 16:21:54
  • ecdfで困るということではないのですが、それ用に用意されているのかな、と思いまして。。 -- sato? 2004-01-22 (木) 16:57:11
  • 無いようですね。累積和関数 cumsum を使えば自作は十分可能でしょうが、stepfun の ecdf の利用が簡便では。 -- 2004-01-22 (木) 17:27:48
  • わかりました。ありがとうございました。 -- sato? 2004-01-22 (木) 17:48:02

Rでビット演算

(2004-01-19 (月) 16:29:43)

Rでビット演算をするには、どうすれば、良いのでしょう

  • R にはそのものずばりのパッケージはないようです。具体的になにをなさりたいのでしょうか。バイナリーデータを扱いたいのか、それとも整数を二進数と考えてビット毎の演算をしたいのでしょうか。 -- 2004-01-19 (月) 17:50:40
  • Colorオブジェクト同士のビット演算をしたいのですが。 -- 2004-01-19 (月) 19:03:52
    > heat.colors(9)
    [1] "#FF0000" "#FF2A00" "#FF5500" "#FF8000" "#FFAA00" "#FFD500" "#FFFF00"
    [8] "#FFFF40" "#FFFFBF" 
    のようなカラーを与える16進数(を表す文字列)間の and, or, xor, negate 等の演算という意味でしょうか。固定長で大量に処理する必要がない限り、腕力で関数を書くことは難しくはないような気がします。まず二進ベクトルに変換し、桁毎に演算することになるでしょうか。ついでにいえば、同じ長さの二進数ベクトルに変換すれば、論理値ベクトルに対する組み込み演算 &,|,!,xor が使えます。-- 2004-01-19 (月) 21:11:25

直線回帰の信頼区間について

ぶんしょ? (2004-01-18 (日) 20:38:47)

直線回帰の95%信頼区間をRで計算したいのですがどうしたらよいでしょうか?

  • predict()関数を使うのが便利です。直線回帰でいいなら,example(predict.lm)とすると,使い方がわかると思います。 -- 中澤? 2004-01-18 (日) 23:45:15
  • ありがとうございます。さっそく試してみます。 -- ぶんしょ? 2004-01-19 (月) 12:58:34
  • プロットまでやる例をhttp://phi.ypu.jp/swtips/sumfunc.html#linear.regに作ってみました。ご参考まで。 -- 中澤? 2004-01-23 (金) 15:40:22
  • 久しぶりに見ました。丁寧にありがとうございます -- ぶんしょ? 2004-01-27 (火) 14:23:46

スレッド表示

(2004-01-13 (火) 21:19:06)

Q & Aではできないのでしょうか?

[質問] lm() 関数の返り値のF統計量のp-value <= なんでも掲示板から移動

舟尾? (2004-01-10 (土) 21:12:33)

lm()関数で重回帰分析を行った際,

result <- summary(lm(y~x))

として result[[k]] や result$names (k:整数,names:名前)
をすることで推定値やt値などが取り出せますが,F-statistics
のp-valueだけがどうやっても取り出せません.names()で調べても
F-statisticsのp-valueにあたるものが見当たらないのですが,
この件に関して情報をお持ちの方,レス頂ければ幸いです.

間瀬? 2004-01-10 (土) 22:52:45
たしかに summary.lm 関数は F 統計量の p 値自体は計算しないようですね。次のようにすれば計算可能ですが、summary.lm が表示する p 値はいつの時点で計算されるのか?

> x=1:10 
> y=c(rnorm(5),rnorm(5,mean=1)) 
> z=summary(lm(y~x)) 
> w=z$fstatistic 
> pf(w[1], w[2], w[3], lower.tail=FALSE)     
      value 
0.1643454  

中澤? 2004-1-12(月)11:42:00
print.summary.lmの中のようです。プロンプトでprint.summary.lmと打てば,

if (!is.null(x$fstatistic)) {
       cat("Multiple R-Squared:", formatC(x$r.squared, digits = digits))
       cat(",  Adjusted R-squared:", formatC(x$adj.r.squared, 
           digits = digits), "?nF-statistic:", formatC(x$fstatistic[1], 
           digits = digits), "on", x$fstatistic[2], "and", x$fstatistic[3], 
           "DF,  p-value:", format.pval(pf(x$fstatistic[1], 
               x$fstatistic[2], x$fstatistic[3], lower.tail = FALSE), 
               digits = digits), "?n")
   }

となっているのがわかります。

  • なるほど。気づきませんでした。R は部品化されているのでこの辺が気になるとわかりにくいですね。おさらいすると、lm 関数が返すオブジェクトには属性 "lm" がつく。これに summary 関数を適用すると暗黙のうちに summary.lm 関数が呼び出される。さらにそれを表示しようとすると print.summary.lm 関数が呼び出されるという具合でしょうか。> lmオブジェクト による内容表示の際も暗黙のうちに適当な print method が選ばれ、それに応じて整形されて出力されるわけですね。 -- 間瀬? 2004-01-12 (月) 12:09:54
  • 間瀬さん、中澤さん、非常に丁寧な回答をありがとうございます。私も関数pf()で計算しなければいけないのかと思っておりましたが、printする際に print.summary.lm 関数が呼び出されて出力が整形されていたのですね。以後、他の関数を用いた際に同様の症状が出た場合はprint.●●...を疑ってみます。 -- 舟尾? 2004-01-13 (火) 18:08:24

ファイルが読み込めません

haru? (2004-01-09 (金) 14:52:09)

初心者なのですが、R1.8.1をダウンロードして早速Excelでデータを入力し、customernumberという名前でテキスト形式で保存。Rで読み込もうと思って

customernumber <- read.table("costomernumber.txt")
としましたが
Error in file(file, "r") : unable to open connection
In addition: Warning message:
cannot open file `costomertable.txt'
というエラーメッセージがでてしまい読み込むことができません…。
read.delimなども使ってみましたが同じ結果でした。
http://phi.ypu.jp/swtips/Data2R.htmlに書いてある説明通りにやったはずなのですが…
Rに関する質問なのかもよくわかりませんが、よろしくお願いします。

  • 原因は,ファイル名の指定が正しくないことだと思われます。read.tableやread.delimに与えるファイル名は,カレントディレクトリにあるのでなければ,フルパスで指定しないとだめです。例えば,フロッピーにあるなら,read.delim("a:/customernumber.txt")としなくてはだめですし,My Documentsにあるなら,Windows2000でユーザusernameでログインしている場合,read.delim("C:/Documents and Settings/username/My Documents/customernumber.txt")としなくては読めません(こうするのは面倒なので,例えばCドライブにworkという名前のディレクトリ(フォルダ)を作っておいて,RのGUIを起動するアイコンのプロパティのショートカットタブの作業フォルダをC:?workと設定しておき,Excelで名前をつけて保存するときにC:?workに保存するようにすれば,いちいち振るパスで指定する必要がなくなります)。また,Excelで名前を付けて保存するとき,拡張子をつけなかったなら(通常は勝手につくはずですが,つかない場合もあります),読み込むときも拡張子なしです。 -- 中澤? 2004-01-09 (金) 15:14:20
  • できました!ご説明本当にありがとうございました! -- haru? 2004-01-11 (日) 12:48:26

R初心者

文字の読み込み? (2003-12-16 (火) 01:33:52)

文字列(例えば3030A6-0010など)の入ったファイルから各々の文字列を読み込みたいのですが、scan()関数を使ってcdata <- scan("worldstnid.dat")??
stnid <- cdataとしても、Error in scan("worldstnid.dat") : "scan" expected a real, got "3030A6-0010"と、怒られてしまいます。どうもscan()関数は数値データにしか対応しいないように思えて仕方ないのですがそうなのでしょうか?また、このような文字列を読み込む関数がほかにあるのであれば教えていただきたいです。よろしくお願いします。

  • scan("worldstnid.dat", what=character())とwhat=を引数で指定します。 -- かくれR狂? 2003-12-16 (火) 08:22:56
  • scan にデータがどういうタイプか教える必要があります。既定では実数値とみなします。?scan でヘルプをじっくり調べて下さい。

macでRのproxy設定

mac? (2003-12-14 (日) 05:55:22)

OSXのパンサーを使っています。
Rでproxyを利用したいのですが、どうすればいいのでしょうか。
教えてください。

  • すでに試されたかもしれませんが、googleでproxy site:r-project.orgとすると、http://www.r-project.org/nocvs/mail/r-help/2001/4560.htmlというのがでました。help(download.file)にproxyk -- detour? 2003-12-14 (日) 17:48:45
  • 切れてしました。続き。proxyの設定方法が書いてありました。ただ環境変数というのかいまひとつわかりません。Unix shellの環境変数のことのようです。私はMacではなく、proxyをまだ必要としていないので、これ以上はわかりません。何か参考になれば。 -- detour? 2003-12-14 (日) 17:50:55
  • ありがとうございました。出来ました。他にも困っている人がいるかもしれないので、やり方を書いておきます。 -- mac? 2003-12-14 (日) 19:16:10
  • まず、インストールした直後に、Sys.putenv("http_proxy"="http://xxxx.com:xxxx)と入れます, -- mac? 2003-12-14 (日) 19:17:21
  • ftpも同様に入れてください。Rをインストール直後にこれを行いませんと、中の環境変数が再設定されません。これでOKです(私は再インストールしました、しかしパッケージなどはそのまま残っているのでほぼ問題なかったです。) -- mac? 2003-12-14 (日) 19:19:06
  • 環境変数はプログラムの起動時ごとに設定されるので、インストール直後ではなく、毎回Rを起動後にSys.putenv()をする必要があると思います。起動直後に行う必要があるのはRの仕様で、インターネットからのダウンロードをするinstall.packages()等を起動後初めて使う際にしか環境変数を読まないようです。 -- 岡田 2003-12-14 (日) 19:59:23
  • また、毎回起動時に環境変数を設定する必要をなくすには、/Applications/StartR.app/RAqua.app/Contents/R というunix shellで書かれたスクリプトを編集し、環境変数を適切に設定するように変更する必要があります(現時点では)。本来はStartR.appの機能に含めるべきでしょうね。 -- 岡田 2003-12-14 (日) 20:02:06
  • あと、インストールしたパッケージとかは ~/Libaray/RAqua/ にあるので、ここを消さない限り残ります。一方、StartR.appの設定などは~/Libaray/Preferences/RAqua.plist です("~"はホームフォルダのことです)。このへんTipsページにまとめておいたほうがいいかもしれませんね。 -- 岡田 2003-12-14 (日) 20:07:00

同時実行

(2003-12-12 (金) 16:58:27)

Windows用のRを2つ立ち上げて、
それぞれを操作した場合に、
それぞれの実行結果は完全に独立だと思ってよいのでしょうか?

  • 独立とは統計的な独立(乱数を使用した場合に限り意味のある質問ですが)という意味ですか、それとも通常の計算結果の独立(例えば双方で同じ名前の変数を使用した場合それらが交じりあわない)という意味ですか。おそらくどちらの意味でも独立だと思います。ただし終了時問い合わせにともにyesと答えると困るでしょうね。MSW では事情がわかりませんが、Linux では違ったディレクトリで立ち上げれば、結果の保存もそれぞれのディレクトリの .Rdata .Rhistory に保存されますので全く問題ありません。実行速度は当然落ちるでしょうが。 -- 2003-12-12 (金) 22:16:32

グラフがかけない X11 is not available

TAKU? (2003-12-12 (金) 11:33:19)

超初心者で、はじめてLinuxをインストールし、Rをemacs上で立ち上げる
とこまでできたのですが、postscript以外のグラフィックコマンドが
使えないのです。png("hoge")やx11()どを入力すると
X11 is not available
とエラーメッセージがでてしまいます。Rの質問ではないのかもしれませんが、
助け舟をよろしくお願いします。
ディストリビューションはdebianです。

  • X Window Systemは起動しているんですね?それならば X11(display=":0") としてみてください。 -- 2003-12-12 (金) 14:07:16
  • Error in x11(display = ";0") : X11 is not available というメッセージが出てしまいました。gnuplotのほうでは書けるのですが・・・。kterm上でのRでも試しましたが、同じメッセージでした。お返事ありがとうございます。 -- TAKU? 2003-12-12 (金) 14:18:58
  • 私もdebianです。同じことがおこったことがあります。もしかして、 Rはsourceから../configure;makeしましたか? 私の場合はdebianのX11の開発環境がインストールされていなかったからでした。 -- detour? 2003-12-12 (金) 20:03:51
  • で、debianのX11クライアント開発環境のインストールですが、今、www.jp.debian.orgからpackageの一覧や説明が載っているホストがdownしていてどのpackageだったかはっきり確認できません。が、たぶん、xlibs-devだと思います。dpkg -s xlibs-devでインストールされているかどうか調べてみて下さい。 -- detour? 2003-12-12 (金) 20:33:51
  • もしインストールされていなければ、apt-get -s install xlibs-devで試してみて、本当にインストールするならapt-get install xlib-devでインストールできるはずです。 -- detour? 2003-12-12 (金) 20:35:51
  • 開発環境がインストールできたら、Rのsourceに戻って、./configureし直します。config.logをみて、checking for Xや#define HAVE_X11 1を確認します。X11が認識されていることを確認して、makeとmake installします。 -- detour? 2003-12-12 (金) 20:40:45
  • 私は試していませんが、make が面倒なら、www.r-project.orgのCRANのdownloadにはdebian用のパッケージが用意されているようです。日本のミラーサイトftp://ftp.u-aizu.ac.jp/pub/lang/R/CRAN/bin/linux/debian/ReadMeなどをよみ、設定を行えばapt-getでRが直接インストールできるでしょう。1.8.0ですが。たぶん、debianの必要なライブラリなどは字動的にインストールされるはずです。 -- detour? 2003-12-12 (金) 20:46:12
  • detourさんほんとうにいろいろありがとうございます!(;_;)。さっそくやってみます。 -- TAKU? 2003-12-12 (金) 21:00:56
  • ありがとうございました!解決しました。このことでずっと迷っていたので、ほんとうに助かりました。debianはパッケージが難しいですね。Rを使う目的でlinuxにしたので、使えてうれしいです。このこともメモしておきました。ありがとうございました! -- TAKU? 2003-12-12 (金) 21:12:03
  • いえいえ、Knoppix-jp CDROM から HD へインストールすれば簡単に Debian マシンに変身ですよ。 -- 2003-12-12 (金) 22:39:08
  • テスト中CRANミラーにもDebian用バイナリがあるのでご利用ください;-) 実はHDD容量のほとんどをLinux用バイナリがしめていることがミラーしてみてわかりました(^^;; -- 岡田 2003-12-14 (日) 20:15:46
  • ほんとにいろいろありがとうございました。まだまだ初心者でわからないことだらけですが、このページを見ながら研究を進めていきます。 -- TAKU? 2003-12-15 (月) 11:19:12

Rプログラミングで得られた結果を別のファイルへ出力するには?

ぽいと? (2003-12-11 (木) 18:21:16)

現在、R言語でプログラムを書いたのはいいのですが、その結果(様々な数値)を別のファイルへ格納したいのですが、その方法が分かりません。扱っているデータが膨大な上、outputも多いのでいちいちコピー&ペーストと言うわけにもいきません。どなたか知っておられる方がおりましたら、教えてください。よろしくおねがいします。

  • プログラムの冒頭でsink("c:/temp.txt")のようにすれば、temp.txtに出力されます。解除はsink()。-- 2003-12-11 (木) 19:50:02
  • Tips 集の R出力の記録 が参考になります。 -- 2003-12-11 (木) 21:48:09
  • 早速の回答をありがとうございます。投稿後、色々調べてみたところ、write()関数でも出力できることを知りました。sink()関数でも試してみようと思います。有難うございました。 -- ぽいと? 2003-12-12 (金) 00:52:39

Rパッケージのインストール&ロード方法

R初心者? (2003-12-11 (木) 16:03:04)

現在Linux上でパッケージ「evd」をinstall.package("evd")で確かにインストールしたのですが、search()でインストールされているパッケージ情報を参照しても、インストールされているとの表示がでません。またWindows版のRではLoad Packagesで簡単にロードできたのですが、Linux上ではそのロード方法が分かりません。どなたかご存知の方がおられましたら、ご教授願います。

  • 該当パッケージをダウンロードした上、ルートになって R CMD INSTALL .....tgz というのが私がいつも使っている方法ですが、install.packages("evd") でも無事にインストール出来ましたが? ルートで実行していないとか? ちなみにライブラリのロード方法は libarary(evd), ライブラリ中のオブジェクト一覧は library(help=evd) です。-- 2003-12-11 (木) 21:08:28

大規模データからランダムにデータをサンプリングする

msasalo? (2003-12-10 (水) 21:48:07)

はじめまして。最近Rを使いはじめるようになった者です。
200万件38Mbのデータを処理しようと思っていたのですが、読み込みできず、あきらめました。
そこで、データのサブセットを作って36万件6Mbのデータを作ったのですが、読み込みはできたものの、Strip spot anovaができません・・・・。1way-anovaが何とかできる程度です。
マシンのスペックは、P4-1.6GHz,512Mで、ページングメモリーは4G確保しました。実行すると15分ほどで、CPU使用率は100%、ページングファイルは0.8Gとなり、12時間経過すると、1.3G程度となり・・・何も変わりません・・・1週間待つべきでしょうか・・・。
データファイルを元にランダムにとる方法がRにあるとありがたいのですが・・・。どなたか、お手すきの方がいらっしゃいましたらアドバイスお願いします。

  • sample 関数が探しているものでは?-- 2003-12-10 (水) 23:18:23
    > x=rnorm(2000000) # 200万個の正規乱数発生 
    > length(x)
     [1] 2000000 
    > 2000000*8/1024/1024 
    [1] 15.25879   # 15 Mbyte
    > xx = sample(x, 10000)  # その中から一万個を無作為抽出
    > length(xx)
     [1] 10000  
  • Rに読み込んであるデータならばsample()という関数が使えます。Rに読み込んでないデータなら、データベースサーバを立てて、そちらの手を借りるのでしょうか。 -- 悠遊? 2003-12-13 (土) 01:22:34
  • そうなのです、以前はPostgreSQLを立ち上げてデータをODBCとかで読んでいたのですが、DHCPになってからは、アドレスをいちいち変更するのが大変で、ついつい、Accessを使うようになり・・・。悪いのではないのですが、どうもAccessと性格が合わず、仕事がはかどりません・・・。 -- msasaki? 2003-12-16 (火) 14:01:01
  • ありがとうございました!! なんとか解決しました!! Rで216万の数から21.6万の数を発生させ、それをテキストに落とし、Accessでよみこんで元のデータにリレーションをかけました。あとは該当レコードを抽出し、テキスト経由でRに読み込みました。本当にありがとうございました。助かりました!!! -- msasaki? 2003-12-16 (火) 13:57:07

[質問]2変量t乱数の作り方

(2003-12-10 (水) 16:59:27)

2変量乱数の発生でt乱数の作り方を教えて下さい。

  • CRAN にあるアドオンパッケージ mvtnorm は多変量の正規分布とT分布用のパッケージのようです。-- 2003-12-10 (水) 20:33:13

アルゴリズムを書く代わりに命題みたいなものを示します。

  • R:自由度mのカイ二乗分布に従う確率変数
  • Z:p変量正規分布N(0,Ip)に従う確率ベクトル(独立な標準正規乱数がズラッと縦に並んでいる)
  • V:正定値対称行列(ちらばりを表す行列で分散共分散行列とは少し異なるもの)
  • V=C%*%t(C):コレスキー分解
    • X = sqrt(m/R)*Z はp変量楕円t分布met(p,m,0,Ip)に従います
    • Y=μ+C%*%X はp変量楕円t分布met(p,m,μ,V)に従います

拡張性があるように2変量ではなく3変量楕円t乱数に従う乱数を生成する関数を以下に定義します(2変量の場合は met3 中のfor文の中身をいじって下さい)。

met3 <- function(m, mu, V, n) {
  # m  : 自由度
  # mu : 平均ベクトル
  # V  : 散らばり行列
  # n  : 乱数の個数
  U  <- svd(V)$u
  V1 <- svd(V)$v
  D  <- diag(sqrt(svd(V)$d))
  B  <- U %*% D %% t(V1)
  
  w <- c()
  for (i in 1:n) {
    R <- 0
    for (j in 1:m) {
      R <- R + rnorm(1)^2
    }
    w <- append(w, list(mu + B %*% (cbind(rnorm(3))*sqrt(m/R)))) 
  }
  return(w)
}
# 実行例
mu <- cbind(c(1,1,1))
V  <- array(c(2,1,1,1,2,1,1,1,2), dim=c(3,3))
n  <- 1000
w  <- met3(5, mu, V, n)
sm <- 0
for (i in 1:n) { sm <- sm + w[[i]] }
sm <- sm/n
sv <- 0
for (i in 1:n) { sv <- sv + (w[[i]]-sm) %*% t(w[[i]]-sm) }
sv <- sv/n

この分布の平均ベクトルはμ、分散共分散行列はm/(m-2)*V(m>2)なので、 多変量楕円t乱数がうまく出来ているかどうかはこれで確かめることが出来ます。-- 舟尾? 2004-01-16 (金) 23:08:24

[質問]hist()関数とnclass()関数の関係

(2003-12-06 (土) 11:09:38)

どちらの関数でもSturgesの方法によって階級が求められるようですが、nclass.Sturges()は、hist()による出力結果ではどの部分になるのでしょうか?

  • ご質問の意味が今一不明確ですが、?hist に書いてあるように、hist 関数のクラス分割指定オプションには4種類あり、それぞれ
  • (1) 単に hist(x) は hist.default(x, breaks = "Sturges") で nclass.Sturges を使う。他にも breaks = "Scott", breaks="FD" がある
  • (2) 明示的に分割点ベクトルを与えても良い、
  • (3) 単に整数 n を一つを与えれば、区間が n 等分される、
  • (4) 分割点を計算するユーザー定義関数を与える

ということでしょうか。 -- 2003-12-06 (土) 17:00:57

  • 具体的には(1)に関してのことになりますが、以下ではsは常に同じ返り値であり、bではそうはなりませんでした。hist.defaultはSturgesの公式で階級数を求めているのでしょうか?
x<-rnorm(80,mean=0,sd=1)
s<-nclass.Sturges(x)
b<-length(hist.default(x,breaks = "Sturges")$breaks)
s;b
  • たしかに変ですね。関係する関数は内部関数なので詳細はソースコードを覗いてみないと? -- 2003-12-08 (月) 23:59:51
> test <- function(x)  {c(nclass.Sturges(x),nclass.scott(x),nclass.FD(x),
+                                   length(hist(x)$breaks),length(hist.default(x)$breaks))}
> test(rnorm(1000))
[1] 11 21 26 16 16
> test(runif(1000))
[1] 11 10 10 11 11
> test(rnorm(1000))
[1] 11 19 25 15 15
> test(rexp(1000))
[1] 11 20 32 15 15
> test(rcauchy(1000))
[1]   11  109 3057   13   13  # ずいぶんと違うもんですね!グラフでは br="scott" が一番まし
  • R の ML に Sturges 方式を hist の既定にするほうが良いかどうかというコメントがありました。Sturges 方式は単なる惰性で、使用は止めたほうが良さそうですね。 -- 2003-12-10 (水) 12:03:37
  • おなじことは他の方式でもおきているようです。一つのありそうな回答は、hist は nclass 関数の計算した分割数を参考にするが、そのまま使うわけではないということでしょうか。-- 2003-12-10 (水) 12:16:36
    > nclass.scott(x)
     [1] 9 
    > length(hist(x,breaks="scott")$breaks) 
     [1] 8 
     > nclass.FD(x)
     [1] 18 
    > length(hist(x,breaks="FD")$breaks)
     [1] 14 
    > nclass.Sturges(x) 
     [1] 9 
    > length(hist(x)$breaks)
     [1] 8 

R(D)COMの使い方

f? (2003-12-06 (土) 00:37:19)

R(D)COMを使って、VC++で作ったプログラムにRを組み込みたいのですが、上手く出来ません。パスの設定はしたのですが・・。
R(D)COMをお使いの方がおられましたら、使用法等、ご助言お願いいたします。

[質問]条件分岐

さかい? (2003-12-04 (木) 18:59:26)

i<500のときx=1、500<=i<1000のときy=1、
それ以外のときz=1とするプログラムを作成したかったので、
以下のようなプログラムを作成したのですが、
iが1000以上のときx=0,y=1,Z=0となり、うまく値を返してくれませんでした。
どのようにしたら解決できるのか分かる方がいたらアドバイスをお願いします。

以下のプログラムはRプログラミングTips大全集の条件実行 if , if elseを参考にして作成しました。

#プログラム
x <- 0
y <- 0
z <- 0
i <- 2000
if (i < 500) {x <- 1}  else if (500 <= i < 1000) {y <- 1}  else z <- 1
#結果
> x
[1] 0
> y
[1] 1
> z
[1] 0
  • disc <- function(i) { x <- 0; y <- 0; z <- 0; if (i<500) { x <- 1 } else if (i < 1000) { y <- 1 } else { z <- 1 }; list(x,y,z) }として,disc(2000)とかすれば? -- 中澤? 2003-12-04 (木) 19:24:07
  • 2番目のifは1番目のifに対するelseの後なので,500<=iの部分は要らないということです。eval(500 <= 2000 < 1000)はTRUEを返します。 -- 中澤? 2003-12-04 (木) 19:29:56
  • ありがとうございます。やってみます。 -- さかい? 2003-12-04 (木) 19:43:32
  • eval(500 <= 2000 < 1000)はなぜTRUEを返すのですか? -- さかい? 2003-12-04 (木) 20:03:19
  • if (i < 500) {x <- 1} else if (500 <= i && i < 1000) {y <- 1} else {z <-1} とすべきだったのです。もちろん中澤さんの書かれたように 500 <= i の判断は冗長で無用です。ついでながら直前のご質問は面白いポイントをついていますね。評価 500 <= 2000 < 1000 はおそらく (500 <= 2000) < 1000 という順序で行なわれ、従って最終評価は TRUE < 1000 で真と評価されたということなのではないでしょうか。ちなみに R では TRUE は必要に応じて整数 1 に強制変換されます。 ご質問の計算は次のようにしても可能です。
x <- ifelse(i < 500, 1, 0)
y <- ifelse(i < 1000, 1, 0)
z <- ifelse(i >= 1000, 1, 0)

ついでですが、今チェックをしていて気づきましたが R 1.8 からは暗黙のリスト返り値 return(x,y,z) が非難されるようになったのですね。便利だと思うのですが。

> test <- function(x,y) return(x,y)
> test(1,2)
$x
[1] 1
$y
[1] 2
Warning message:
multi-argument returns are deprecated in: return(x, y)
  • 2003-12-04 (木) 21:38:49
  • もっと簡単に(+FALSE は論理値を整数に強制変換するトリック。+0 でも良いが実数になってしまう。別に論理値のままで後々困ることはないが気持ちがわるいかもしれない)
    > x = (i < 500)+FALSE; y = (i < 1000)+FALSE;  z=(i >=1000)+FALSE 
    > x;y;z 
    [1] 0 
    [1] 0 
    [1] 1  
  • 2003-12-04 (木) 22:02:08

Rdbi RdbiPgSQL

なかま (2003-12-04 (木) 10:42:45)

http://www.bioconductor.org/ のが 1.0.2 で cran.r-project.org の contrib/Devel のは RPgSQL_1.0-0.tar.gz http://sourceforge.net/projects/rdbi/ のは 0.1.2 とすると、
bioconductor のが最新?

  • ところで、RdbiPgsql? はCRAN にはバイナリがありませんが、これはMinGWでビルドできるのでしょうか? -- 2003-12-04 (木) 15:25:38
  • PostgreSQL Native Win32だから、libpq作ってからMingwで(素直にonfigure出来るかは謎。できなくても大していじらなくてもいいとおもう)ビルドかな?。 -- なかま 2003-12-04 (木) 17:56:11
  • やってみると、面倒くさかったのでVine2.6(PostgresSQL-7.2.3)に合わせてWin用バイナリを置いておきます。 -- なかま 2003-12-04 (木) 22:50:54
  • Rdbiはbioconductorの方にバイナリがあります。 -- なかま 2003-12-04 (木) 22:52:50
  • やっぱりbioconductorのが最新のようですね(他にもあったりして)。 -- なかま 2003-12-05 (金) 00:24:39
  • なかまさん、RdbiPgSQLのバイナリありがとうございます。わたしも、postgresql7.4で libpq のビルドを試みているのですが、pq.dll はできても、libpq.dll はできないでいます。ところで、なかまさんのサイトにこのバイナリについて「サーバのバージョンが異なると動きません」とありますが、この「サーバ」とはPostgreSQLのバージョンを指しているのでしょうか? -- 2003-12-05 (金) 16:33:10
  • libpqとpostmasterのバージョンが異なると通信出来なかったと思います。そういう意味でのサーバです。(解かりにくいかな?)ですから、7.2.3用のみです。 -- なかま 2003-12-05 (金) 16:43:52

[質問]棒グラフで対数軸

(2003-11-29 (土) 17:20:27)

棒グラフを作図するのには barplot(1:10) のようにすれば良いとわかりました。
そこで、Y 軸を対数軸として作図しようと log="Y" を指定したところ以下のようなエラーとなりました。

barplot(1:10, log="y")
Error in plot.window(xlim, ylim, log = "", ...) :
formal argument "log" matched by multiple actual arguments
対数軸の棒グラフはどのように指示すれば描くことが出来るのでしょうか?
よろしくお願いします。

  • 残念ながら base パッケージの barplot 関数は対数軸をサポートしていないようです。おそらくあまり需要がないと考えているのでしょうか。また par(ylog=TRUE) の指定も無効にされるようです。アドオンパッケージ gregmisc 中には両対数軸オプションを持つ barplot 関数があるそうです。ついでながら、こうした疑問は CRAN の過去 ML 記事検索を使うと殆んど誰かが過去に質問をしていますので一度試してみて下さい。 -- 2003-11-29 (土) 18:22:07
  • ありがとうございます。アドオンパッケージを調べてみます。CRAN(英語)は、私にとって敷居が高くて.....(申し訳ありません) -- 2003-11-29 (土) 20:27:13

subsetで値を取り除く

detour? (2003-11-28 (金) 14:49:55)

subsetで名義変数の特定の値を除いたサブセットを作成したとき、例えば、
a<-data.frame(choice=c("ab","ab","ab","c","c","d","d","d"))
summary(a)
summary(subset(a,choice!="c"))
"c"というのが一応残っているようです。これはどうやったら消せるのでしょうか。初心者ですが、よろしくお願いいたします。

  • 次のようにすれば一応ラベル c は消えますが、おそらくこれはご希望の方法ではないでしょうね。 -- 2003-11-28 (金) 16:14:04
    > b=data.frame(choice=a[a!="c"]) 
    > summary(b)  
    choice  
    ab:3  
    d :3
  • 使用の目的にもよるのでしょうが、気にしないのもありかと。summary 関数を使わなければ度数ゼロの項目 c は表示されませんね。 -- 2003-11-28 (金) 19:19:45
  • 消せないせいでうまく行かなかったことがありましたが、勘違いかもしれません。あまり気にしないことにします。 -- detour? 2003-11-29 (土) 08:28:44
  • おまけ:subsetを見付ける前はbquoteを使ってみようと思ったりしました。[choice!="a"]を各カラムに挿入するのに。 -- detour? 2003-11-29 (土) 08:31:39
  • 気にしない事にしましたがやはりうまく行かない事もありました。上の例ではchoiceが説明変数に入っている場合などは頻度0のレベルがあることからエラーになってしまうこともあります。 -- detour? 2003-12-10 (水) 00:57:42
  • で、できました。b<-subset(a,choice!="c")だとすると、levels(b$choice)<-c("ab","ab","d")で一応、望むようになったように見えます。 -- detour? 2003-12-10 (水) 01:00:23
  • おめでとうございます。苦労して得られたノウハウは是非然るべき Tips にメモを。誰かが悩んだことは、人も悩むところです。 -- 2003-12-10 (水) 08:23:55

Rのインストール方法(UNIX)

R初心者? (2003-11-28 (金) 14:21:53)

ただいまRをUNIXにインストールしようとしています。R-1.8.1のソースをダウンロードしてきて解凍し、./configureは成功したのですが、makeができません。以下のような感じで怒られます。

make
Make: Cannot open /share/make/vars.mk. Stop.
vars.mkファイルは確かに/share/make/に存在しましたし、原因が分かりません。どなたか解決方法をご存知でしたらご教授下さい。お願いいたします。

  • Debian GNU/Linux 3.0ですが、問題なくmakeできました。Unixとは? permissionの問題は? -- detour? 2003-11-28 (金) 14:53:42
  • 使っているUNIXマシンはCompaq Tru64です。 -- R初心者? 2003-11-28 (金) 14:59:28

添付ファイル: fileasp.png 1473件 [詳細] fileobi.png 2613件 [詳細] filehist.png 2709件 [詳細] fileobi2.png 2622件 [詳細] fileplot-00.png 2798件 [詳細] fileobigraph.jpg 2814件 [詳細] filegraph.png 2683件 [詳細] fileess-20040323.png 591件 [詳細] filepar.fin.jpg 1500件 [詳細] fileobi3.png 2694件 [詳細] fileplot-01.png 2620件 [詳細] filematplot-00.png 2534件 [詳細] filefinplt.jpg 1576件 [詳細]

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