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

新規投稿はできません




起動中の vignette 呼び出し

青木繁伸 (2005-11-07 (月) 12:40:49)

Macintosh OS 10.4.3 で,R 2.2.0 を使っていますが,R を起動すると vignette という関数(<environment: namespace:utils>)を呼び出して大変時間を食います(今までのバージョンではこのようなことはなかった。少なくともこんなに時間がかかる処理はなかった)。どこかの設定でこれをやらないようにできるのでしょうか。~
現在は,ダミーの vignette 関数を定義してバイパスするようにしております。

  • 私はPowerBook? G4(1GHz)、Mac OS 10.4.3でR 2.2.0を使っています。確かに、起動時にウィンドウ最下部のステータス表示を見る限りではvignette()$resultsを実行していますが、5秒程度で起動でき問題ありません。 -- ara? 2005-11-07 (月) 17:19:29
  • パッケージはどれくらいインストールなさってますか。
    私は全部インストールしているせいか(^_^;),dual 2.5GHz PowerPC G5 で vignette の実行が開始されて終わるまで27秒ほどもかかります(^_^;) -- 青木繁伸 2005-11-07 (月) 19:56:19
  • お、多いですね・・ 私の場合、追加インストールしたパッケージは5つです。初心者なので詳しいことはわからないのですが、パッケージの内容(構成)にも関係あるみたいですね。 -- ara? 2005-11-07 (月) 20:54:25
  • 起動中の vignette にどのような意味があるかはよく分からないのですが,起動後に vignette やってみると,ま,あまり役に立たない情報を集めているようだし,結局それはインストールしたパッケージを全部走査するようなので,たくさんインストールしているとそれだけ時間が掛かっても文句は言えないということのようですね。
    ということで,取りあえずは vignette <- function() return(0) という埋め草を定義してやっています。 -- 青木繁伸 2005-11-07 (月) 21:30:04

計量経済分析向けの商用ソフト

どうむ? (2005-10-26 (水) 11:08:06)

株式市場の(主に国内)分析等でRをバリバリ使っていますが、「比較的安価で広く出回っている計量経済分析向けの商用ソフト」って、なになのか、気になります。S-PlusとかSASとかSPSSとかSHAZAMとかTSPとか使ったことがありますが、どれも先の要件にあてはまらないような気がします。思い当たるのはS-Plusぐらいでしょうか。なので、とても、気になります。(計量経済分析関連のページ、作っていただいて、意味もなく喜んでいたりします)

  • よく考えたら、TSP も Eviews も RATS もそんなに安くはないですね。学生版は格安ですが。まあ、MATLAB よりはだいぶ安いはずなので "比較的"。私も、統計もそれ以外の数値計算もほとんど R でこなしています。エコノメのページ、近々書き足します。 -- 蓮見? 2005-10-28 (金) 00:48:47

独立成分分析のチュートリアル

 ?(2005-10-26 (水) 10:54:37)

Rを使ったICA(独立成分分析)のチュートリアルってどこかにないでしょうか?

時系列データの解析について

山形? (2005-10-26 (水) 10:52:53)

時系列データの解析を行なうにあたり、Box-Jenkins ARIMAを使いたいと考えています。
Rでもこの方法は出来ますか?出来るとしたらどのパッケージをインストールしたら良いのでしょうか?
それとも、普通のARIMAで出来るのでしょうか?

  • google で,site:r-project.org Box-Jenkins ARIMA を検索すると,
    CRAN Task View: Computational Econometrics - [ このページを訳す BETA ]
    Time series modelling : Classical time series modelling tools are contained in
    the stats package and include arima() for ARIMA modelling and Box-Jenkins-type
    analysis. Furthermore stats provides StructTS() for fitting structural time ...
    cran.r-project.org/src/contrib/Views/Econometrics.html - 12k - キャッシュ - 関連ページ
    みたいのが出てきますが -- 2005-10-26 (水) 11:13:28

バスエラー - コアダンプに悩まされております(T_T)

ちょろべぇ? (2005-10-25 (火) 13:34:32)

現在はじめてUNIXでのR実行にチャレンジしています。
ただ、そこでPDF形式で画像を出力しようとしますと

”バスエラー - コアダンプしました。”

というエラーが出てしまい、異常終了いたします。
またUNIXとLINUXの双方でプログラムをじっこうしてみた所、
UNIXではコアダンプしてしまいますが、なぜかLINUXでは問題なく
動作します。
本件、どのように解決すべき問題なのか、ウェブ等色々調べましたが
まったく見当がつかない状態です。お気づきの点等ありましたら
アドバイスいただければと思います。

====== 以下 詳細情報
・Rのバージョン 2.0.1(UNIXもLINUXも)
・実行したOS ソラリス8(UNIX)、スージーLINUX
・実行したコマンド

     TestChars <- function(encoding="ISOLatin1", ...)
    {
        pdf(encoding=encoding, ...)
        par(pty="s")
        plot(c(-1,16), c(-1,16), type="n", xlab="", ylab="", xaxs="i", yaxs="i")
        title(paste("Centred chars in encoding", encoding))
        grid(17, 17, lty=1)
        for(i in c(32:255)) {
            x <- i
            y <- i
            points(x, y, pch=i)
        }
        dev.off()
    }
    ## there will be many warnings.
    TestChars("ISOLatin2")

#↑?pdfで出てくるサンプルスクリプト
・読み込んでいるライブラリー Rの標準もの

  • TestChars?()も同様に失敗するのでしょうか。環境変数LANGには何が入っていますか?pdfのhelpにあるencoding欄の記述が気になります。 -- 谷村 2005-10-25 (火) 19:49:32
  • 自分で make した物なんでしょうか.だったら,自分で gdb 等しないと原因は第三者にはわかりにくいのでは? 2005-10-25 (火) 19:57:43
  • 谷村様、TestChars?も同様に失敗します。なお失敗するUNIXの環境変数LANGにはja_JP.PCKがはいっています。またアドバイスいただきましたPDFのencodingオプション指定ですが、試してみましたがうまくいきません。非常に初心者じみた質問で申し訳ないのですが、UNIXでRを使うときには環境変数R_HOMEをRのlibファイルが置かれている場所として定義しておかないと、フォント等よみこまれないのでしょうか? -- ちょろべぇ? 2005-10-28 (金) 16:42:19

グラフの大きさ

初心者K? (2005-10-23 (日) 18:50:45)

かなり初歩的な質問で申し訳ありません.
毎回同じ大きさグラフをつくりたいのですがデバイスの大きさを指定できません.
どなたかご教示お願いします.

  • どのグラフィックデバイスをお使いですか?それと,OS名とそのバージョン。それを書かないとだあれも答えられませんね。
    私は pdf を使うのですが,pdf(ファイル名, width=インチ, height=インチ) でちゃんと大きさ指定できますが。-- 2005-10-23 (日) 19:27:30
  • 初級Q&A アーカイブ(1)の「出力画像のサイズを変更したい」で解決しませんか?ただ、何も指定しなければデフォルト固定で毎回同じサイズになると思うので、サイズが固定されなくなる何らかの作業をされていると思います。または、「同じ大きさのグラフ」の定義が異なるとか。 -- 谷村 2005-10-24 (月) 10:43:48
  • 「毎回同じ大きさにならないグラフ」のスクリプトを掲載した方が解決が早いと思いますよ。
    私はxinXPproでRは1.90ぐらいから現バージョンまで使っていますが、このような障害はでていません。 -- Akira? 2005-10-24 (月) 13:43:59
  • ご回答ありがとうございます.初級Q&Aアーカイブ(1)等を参照しましたがうまくいかなかったためすべて削除後,再インストールをしましたところ問題なくできました.皆様ご迷惑をおかけしました.また何かありましたらよろしくお願いします. -- K? 2005-10-24 (月) 21:31:47

連続尺度の変数を名義尺度にコーディングしたい。

bob3? (2005-10-23 (日) 01:46:36)

連続尺度の変数を名義尺度の変数にコーディングしたいと思っています。
たとえば、以下の例では AREA が1〜5の場合は「P」という区分に、6〜9の場合は「Q」という区分にコーディングしたいのですが、全て「Q」になってしまいます。

x <- data.frame(AREA=1:9)
AreaCode <- c('P','Q')
ifelse(x$AREA>5, x$A2<-AreaCode[1], x$A2<-AreaCode[2])
# ここでは期待した結果が出力されるのですが、
x
# データフレームをみると全て「Q」になっています。


上手なやり方があれば、教えてください。
なお、使っているのは 闇R 2.2.0 です。
よろしくお願いいたします。

  • 不等号が逆なのでは?ちなみに x$AREA < 6-- 2005-10-23 (日) 09:36:13
  • 闇R2.2.0以外のほかのバージョンではどうなるんですか。それと,問題解決のためにも,OS名とバージョンもちゃんと報告のこと
    Mac R2.1.1では以下の通りになりましたとさ
    > x <- data.frame(AREA=1:9)
    >  AreaCode <- c('P','Q')
    >  ifelse(x$AREA>5, x$A2<-AreaCode[1], x$A2<-AreaCode[2])
    [1] "Q" "Q" "Q" "Q" "Q" "P" "P" "P" "P"
    > x
      AREA A2
    1    1  Q
    2    2  Q
    3    3  Q
    4    4  Q
    5    5  Q
    6    6  Q
    7    7  Q
    8    8  Q
    9    9  Q
    ifelse の第2,第3オペランドで代入しているから変なことになるようですね。 不等号の件はご愛敬か
    ifelse(x$AREA>5, x$A2<-AreaCode[1], x$A2<-AreaCode[2])
    は冗長なので
    x$A2<-ifelse(x$AREA>5, AreaCode[1], AreaCode[2])
    これで,問題解決。
    または
    x$A2<-AreaCode[1+(x$AREA>5)]
    などとも -- 2005-10-23 (日) 09:36:34
  • ご回答ありがとうございます。不等号の件、失礼いたしました。 OSは WindowsXP HOME と Windows98SE で確認しました。
    出力結果を載せなかったのもまずかったです。
    ひとまず、解決いたしました。
    当初の方法だと、以下のようにデータフレームには全て「Q」が入ってしまいます。
    > x <- data.frame(AREA=1:9)
    > AreaCode <- c('P','Q')
    > ifelse(x$AREA<=5, x$A2<-AreaCode[1], x$A2<-AreaCode[2])
    [1] "P" "P" "P" "P" "P" "Q" "Q" "Q" "Q"
    > x
      AREA A2
    1    1  Q
    2    2  Q
    3    3  Q
    4    4  Q
    5    5  Q
    6    6  Q
    7    7  Q
    8    8  Q
    9    9  Q
    そして、ご回答いただいた方法にするとうまく行きました。
    > x <- data.frame(AREA=1:9)
    > AreaCode <- c('P','Q')
    > x$A2<-ifelse(x$AREA<=5, AreaCode[1], AreaCode[2])
    > x
      AREA A2
    1    1  P
    2    2  P
    3    3  P
    4    4  P
    5    5  P
    6    6  Q
    7    7  Q
    8    8  Q
    9    9  Q
    
    >  x <- data.frame(AREA=1:9)
    >  AreaCode <- c('P','Q')
    >  x$A2<-AreaCode[1+(x$AREA>5)]
    >  x
      AREA A2
    1    1  P
    2    2  P
    3    3  P
    4    4  P
    5    5  P
    6    6  Q
    7    7  Q
    8    8  Q
    9    9  Q
    ところで、名義尺度の水準が三つ以上の場合は、以下のように ifelse を入れ子にする方法のほかに上手いやり方(冗長でない方法)はありますでしょうか。
    > x <- data.frame(AREA=1:9)
    > AreaCode <- c('P','Q','R')
    > x$A2<-ifelse(x$AREA<=3, AreaCode[1], ifelse(x$AREA>3 & x$AREA<=6, AreaCode[2], AreaCode[3]))
    > x
      AREA A2
    1    1  P
    2    2  P
    3    3  P
    4    4  Q
    5    5  Q
    6    6  Q
    7    7  R
    8    8  R
    9    9  R
    よろしくお願いいたします。 -- bob3? 2005-10-23 (日) 13:22:47
  • x$AREA をどのような規則性により変換するかにもよりますね
    例では1〜9の整数値ですが,最初の質問では連続値とありましたので,もし1以上10未満の値を取るx$AREAを3区分(例題のように等間隔)するとすれば,
    x$A2 <- AreaCode[as.integer((x$AREA-0.00001)/3)+1]
    のようなのも一つの方法ですね(0.00001 というのは本当は良くない)
    1以上4未満,4以上7未満,7以上10未満というのなら
    x$A2 <- AreaCode[as.integer((x$AREA-1)/3)+1]
    データに応じた,臨機応変なプログラミングが必要になるでしょうが,for と if-elseif-else で地道にコーディングするのが一番一般的ですね。使い捨てプログラムなら,あれこれ考えて時間をつぶす方がもったいないですから。
    ちなみに
    x$A2<-ifelse(x$AREA<=3, AreaCode[1], ifelse(x$AREA>3 & x$AREA<=6, AreaCode[2], AreaCode[3]))
    は冗長で,
    x$A2<-ifelse(x$AREA<=3, AreaCode[1], ifelse(x$AREA<=6, AreaCode[2], AreaCode[3]))
    とするだけでよい -- 2005-10-23 (日) 17:02:02
  • なるほど。実際に分析したいデータは必ずしも等間隔ではないので、こつこつif-elseで書いた方が良いのかもしれませんね。ありがとうございました。 -- bob3? 2005-10-23 (日) 17:26:56
  • 質問文をよく読んでませんでした。ifelse の第2,第3パラメータに付値がありましたので,ベクトル演算の結果,一番最後の結果がリサイクルされて全て同じ値になるのですね。冗長性を排除したら,必然的にバグもなくなり万々歳ですね。 -- 2005-10-23 (日) 19:42:10
  • cut(x$AREA, breaks=c(-Inf,4,7,Inf), right=FALSE, labels=AreaCode?) ではダメですか? -- 2005-10-24 (月) 11:56:11
  • もちろん,これが一番ですね。なにかあったはずと思いつつ,思い出せませんでした。 -- 2005-10-24 (月) 22:32:14

dynlmというパッケージについて

山形? (2005-10-21 (金) 18:32:18)

Rを用いて時系列解析を行ないたいと考えています。

Rは初心者です。
先輩から dynlm というパッケージを用いてやると良いと聞いたのですが、
help(dynlm)を見たときに、Examplesに

## multiplicative SARIMA(1,0,0)(1,0,0)_12 model fitted
## to UK seatbelt data
uk <- log10(UKDriverDeaths)
dfm <- dynlm(uk ~ L(uk, 1) + L(uk, 12))

というプログラムがありました。
L(uk,1)やL(uk,12)の1や12は何を表しているのでしょうか?
これと似たようなプログラムを実行したいので、よろしくお願いします。、

  • An example would be d(y) ~ L(y, 2), where d(x, k) is diff(x, lag = k) and L(x, k) is lag(x, lag = -k), note the difference in sign. -- 2005-10-21 (金) 23:46:33
  • dynlm のヘルプに書いてありましたね。よく読んでなく、すみませんでした。こんな質問にコメントを書き込んでくれた方、ありがとうございました。 -- 山形? 2005-10-24 (月) 13:44:54

メモリの上限に関する質問です

花子? (2005-10-20 (木) 14:55:17)

大規模データを用いて数量化理論を実行したいと思っている
Rの超初心者です。
掲示板の記事によると
「Rは実メモリ又は1024Mという上限がある」ようですが、
これは明確な情報としてどこに明記されているのでしょうか?
(疑っているわけではありませんが、私の中でホントに「1024Mの上限があるのか?」で止まっています。)
実際に、試しに1Gから2Gにメモリを増やしてRを起動し、
memory.limit()をみてみましたが、
1Gの時と変わらないのをみると、やはり1024Mが上限なのかなぁという気はしております。

「どこに明記してあるか?」というあいまいな質問ですが、
人からの情報ではなく、なにか説明書のようなものに上限が書かれているのでしょうか?という質問です。

  • 恐らく MSW 版で R を使っているのでは。Windows 版 R の FAQ (RjpWiki に和訳あり) を見て下さい。R 自体にはメモリー使用上の制限は無いのですが、ハードとOS由来の制限があります。それと、何もオプション無しで起動した場合の制限(コンパイル時のオプション)は区別する必要がありますね。 -- 2005-10-20 (木) 15:24:22
  • 返答ありがとうございました。ちなみに私の使用しているWindows版のRです。説明不足ですみません。また自分なりに解決していきたいと思います。 -- 花子? 2005-10-21 (金) 15:37:03

data.frame の番号の書き換え

近藤? (2005-10-20 (木) 08:59:58)

例えば

data(sleep)
sleep[sleep$g==2,]

とすると

>sleep[sleep$g==2,]
  extra group
11   1.9     2
12   0.8     2
13   1.1     2
14   0.1     2
15  -0.1     2
16   4.4     2
17   5.5     2
18   1.6     2
19   4.6     2
20   3.4     2

となります。
この左端の番号 11 - 20 を 1 からふり直すにはどうしたら良いのでしょうか?

  • 行名なので,例えば data.frame(sleep[sleep$g==2,],row.names=1:10) みたいにします. -- なかま 2005-10-20 (木) 10:01:44
  • row.names を表示しないというのは可能ですか? -- 近藤? 2005-10-20 (木) 11:31:02
  • cat を使って表示するような汎用関数を書くと良いでしょう。 -- 2005-10-20 (木) 12:01:54
  • がんばってみます。 -- 2005-10-20 (木) 12:34:46

さっぱり訳が分かりません・・・・・

リサ? (2005-10-14 (金) 08:03:35)

intvlest<-function(x,n,sd,a){
h<-abs(qnorm((1-a)/2)*sd/sqrt(n))
return(c(x-h,x+h))
}
conf.intvl<-function(n,k,a=0.95,mu=0,sd=1){
Llim<-rep(0,k)
Ulim<-rep(0,k)
for(i in 1:k){
x<-rnorm(n,mu,sd)
lim<-intvlest(mean(x),n,sd,a)
Llim[i]<-lim[1]
Ulim[i]<-lim[2]
}
matplot(cbind(Llim,Ulim))
abline(h=mu)
ng<-0
for(i in 1:k){
if(Llim[i]>mu|Ulim[i]<mu)ng<-ng+1
lines(c(i,i),c(Llim[i],Ulim[i]))
}
return(ng/k)
}

初めまして。

以上のプログラム、学校の授業で作成したのですが、黒板に書かれたものをただ打ち込んだだけで、さっぱり意味が分かりません・・・・

intvlestっていうプログラムに何かを代入して、グラフを作成しているのは分かるんです。
これが結果的に何を計算しているかも分かるんです。
中心極限定理ですよね?
でも、それぞれの言語の意味が分かりません。
かなり初歩的なことだとは思うのですが、このプログラム1行1行が何の作業しているのか軽くコメントしてくださると嬉しいです。
よろしくお願いします。

  • ご質問は「関数が分からない」のか、「文法が分からないのか」のどちらですか?
    関数が分からないなら、R上で「?(関数名)」とすれば使い方が分かります。このHPの「単語検索」を使ってもかなりのことが分かりますよ。文法が分からないなら、R-tipsが良いと思います。
    しかしそれ以前に、学校の授業ならプログラムが黒板に書かれるまでにいろいろと説明があったのではないでしょうか?それに授業なら直接先生に聞くのは一番早いと思いますが。 -- 2005-10-14 (金) 08:47:47
  • 文法が分かりません。 -- リサ? 2005-10-14 (金) 09:00:59
  • help("Syntax") とか help("{") とか help("if") とか help("for") とか help("function") とか. でも,一番いいのは help("先生!") と思われる. -- 2005-10-14 (金) 09:26:04
  • 注意:このコーナーの上に「学校の宿題は自分で考えるべきです」と書いてありますね。 -- QDU? 2005-10-14 (金) 12:18:14
  • その通り! -- 先生? 2005-10-14 (金) 14:09:24
  • プログラムが,あまりエレガントでないので書き換え
    conf.intvl <- function(n, k, a = 0.95, mu = 0, sd = 1) {
    x <- apply(matrix(rnorm(n*k, mu, sd), n, k), 2, mean)
    h <- qnorm((1-a)/2)*sd/sqrt(n)
    lim <- cbind(x+h, x-h)
    matplot(lim)
    abline(h = mu)
    sapply(1:k, function(i) lines(c(i, i), c(lim[i,1], lim[i,2])))
    return(sum(lim[,1] > mu | lim[,2] < mu)/k)
    }
    こんな風にしてみましたが,sapply が,「いやん」な感じ -- 2005-10-14 (金) 15:19:57

maptoolsのインストール

たろう? (2005-10-13 (木) 17:32:01)

R 2.2.0で,maptoolsのzipファイルをR上でインストールしたのですが,
以下のように表示されてしまいます.改めてzipファイルをXP上でlibraryに
展開しても同じでした.library()と入力すると,mapttolsも表示されるのですが,library(maptools)と入力するとエラーが表示されます.
 対策をご教示くださいませ.

utils:::menuInstallLocal?()
package 'maptools' successfully unpacked and MD5 sums checked
updating HTML package descriptions

library(maptools)
要求されたパッケージ foreign をロード中です
エラー:'%s' が要求したパッケージ '%s' は見つけられませんでした

  • そういうときは、まずwebででもパッケージのマニュアルを読んでみましょうね。foreignパッケージは入れていますか(あと,sp パッケージも必要です)? -- 2005-10-13 (木) 19:20:51
  • まあ,「エラー:'%s' が要求したパッケージ '%s' は見つけられませんでした」はエラーメッセージにはなってませんね。 -- 2005-10-13 (木) 21:19:52
  • ありがとうございます。packageのmanualに書いてあるDependsというのは、これらのパッケージを入れておくことが条件と解釈すればよいのですね。 -- たろう? 2005-10-13 (木) 22:07:36
  • 対処方法のめどがたったなら,そのようにして頂いて,うまくいきましたという報告があれば,後進の参考にもなると思われます。よろしく,ご配慮の程お願い申し上げる次第でござリマする。 -- 2005-10-13 (木) 22:12:47
  • うまくいきましたぁ! 英語の勉強にもなりました。御礼申し上げます。 -- たろう? 2005-10-13 (木) 22:59:21
  • どうやったらうまくいったかを書くべきですね。 -- 2005-10-13 (木) 23:49:09
  • 失礼しました.Dependsに書いてあるforeign, spパッケージをインストールするとmaptoolsが使えるようになりました. -- たろう? 2005-10-14 (金) 12:59:49

確率分布同士の演算の信頼区間

takahashi? (2005-10-11 (火) 20:54:48)

Rとは直接関係ないのでここに書いてみます(説明はRでしてみます).

> x<-c(-3,-2,-1,1,2,3)
> y1<-c(0.1,0.2,0.3,0.7,0.8,0.9)
> y2<-c(0.2,0.3,0.4,0.6,0.7,0.8) #データは適当です

で,
y1~N1=pnrom(x,0,sd1),y2~N2=pnorm(x,0,sd2) #これはRのコードではありません.
と正規分布に従うことを仮定して,N1,N2に対するsd1,sd2の推定値と信頼区間は得ることが出来ます.ここで,N3~N1+N2としたとき,N3に対するsd3(推定値はsd1+sd2)の信頼区間を得るにはどのような方法があるのでしょうか?

ちなみにあんまり関係ないんですが,↑の過程でpsignifit(http://www.bootstrap-software.com/psignifit/)のR用ラッパー(historyにはto doになっていたものの開発停止でリリースされる見込み無し)を書いたんですが,c側のコードを実行中にたまに落ちます.メモリ使いすぎなのかな?

  • >を書いたんですが,c側のコードを実行中にたまに落ちます.メモリ使いすぎなのかな?
    プログラムのミスなんでは?
    コードも何も示されないでそんなこと言われても,反応は限られていますね。 -- 2005-10-11 (火) 21:59:41
  • 自然なパラメータは sd1^2+sd2^2 もしくはそのルートだと思うのですが、sd1+sd2 というパラメータ値が特に意味を持つ場合があるのでしょうか。 -- MKR? 2005-10-12 (水) 09:06:37
  • ああ,すみません.勘違いです.そうですね.pnormの引数は標準偏差ですからsd3=(sd1^2+sd2^2)^(1/2)です. -- takahashi? 2005-10-12 (水) 09:13:32
  • とすれば(二種類のデータの数が同じという前提で)適当にペア毎の和(実際はペアの取り方には依存しませんが)をとれば、N(0,sd1^2+sd2^2) に従うデータが得られたわけですから、後は標準的な手法で sd3=(sd1^2+sd2^2)^(1/2) の信頼区間が求められるのでは。 -- MKR? 2005-10-12 (水) 09:28:18
  • すみません,ペア毎の和をとるという操作が良くわからないのですが・・・N3に従うデータが取れれば確かに信頼区間は求められます. -- takahashi? 2005-10-12 (水) 09:50:20
  • >ペア毎の和をとるという操作が良くわからないのですが
    MKRさんが言っておられるように,適当に2つずつとって和を求めるってことです。 -- 2005-10-12 (水) 10:30:56
  • N1,N2からサンプルを生成して和をとりモンテカルロ法でN3の擬似サンプルを得るということですか? -- takahashi? 2005-10-12 (水) 10:50:55
  • そういうことでしょう。ただし,今回の場合は,2つの独立な正規分布に従う確率変数の和の分布は平均値がμ1+μ2,標準偏差はsqrt(σ1^2+σ2^2)であることが理論的に言えるわけですから,モンテカルロ法も,疑似サンプルも無関係ということでは? -- 2005-10-12 (水) 11:38:32
  • 度々すみません.N3のパラメータについて点推定であれば理論的に平均値がμ1+μ2,標準偏差はsqrt(σ1^2+σ2^2)でよいのですが,μ1,μ2,σ1,σ2がそもそも真の値ではなく推定値なので,N3のパラメータについて区間推定を行いたいのです.このとき,σ1,σ2の区間推定は観測データから可能なので,これを使ってN3のパラメータの区間推定をどうやったらいいか,という問題です. -- takahashi? 2005-10-12 (水) 11:48:34
  • 母数がわかっていることなんか滅多にないことなので,サンプル値を使えばよいのだと思います。sd1,sd2がσ1,σ2の推定値なのだから,sqrt(σ1^2+σ2^2)は sqrt(sd1^2+sd2^2)でよいでしょう。
    > y1 <- rnorm(200, mean=5, sd=3)
    > y2 <- rnorm(200, mean=8, sd=2)
    > y3 <- y1+y2
    > mean(y1)+mean(y2)
    [1] 12.82907
    > mean(y3)
    [1] 12.82907
    > sqrt(var(y1)+var(y2))
    [1] 3.726669
    > sd(y3)
    [1] 3.84084
    となりますが,試行回数を多くするなりしてシミュレーションしてみれば? -- 2005-10-12 (水) 12:04:08
  • 上の結果は勿論,試行回数を無限に取れば理論値に収束します.でもここで問題なのは,y1,y2を生成するときのsd1,sd2が推定値である,ということです(何度もすみません).なのでsd1,sd2は標準誤差を持っています(だから区間推定が可能).でも上のシミュレーションにはその情報は含まれていないと思います(sd1,sd2を推定するためのデータのあてはまりが良くても悪くても一緒).sd3の推定値自体はデータのあてはまりに関わらず同じになりますが,区間推定を行えば結果は異なると思うんですが...このsd3の区間推定のやり方を探索しています. -- takahashi? 2005-10-12 (水) 12:25:01
  • 上記試行を繰り返し,sd3 の経験分布を求めれば良いだけでは?
    > sim <- function(n, trial=1000)
    + {
    + 	result <- numeric(trial)
    + 	for (i in 1:trial) {
    + 		y1 <- rnorm(n, mean=5, sd=3)
    + 		y2 <- rnorm(n, mean=8, sd=2)
    + 		result[i] <- sd(y1+y2)
    + 	}
    + 	result <- list(sd3 = mean(result), se = sd(result), result=result)
    + 	class(result) <- "sim"
    +  	return(result)
    + }
    > print.sim <- function(result)
    + {
    + 	cat("sd3=", result$sd3, "?tse=", result$se)
    + }
    > sim(10)
    sd3= 3.477083 	se= 0.8488288
    > sim(200)
    sd3= 3.601848 	se= 0.1748326
    > sim(400)
    sd3= 3.602084 	se= 0.1328714
    とか。。se は sd/sqrt(2n) かなあ。ちゃんと数学で解かないといけないカモね
    以下で図を描くと,ほぼあっていそう。
    > result1 <- result2 <- numeric(100)
    > for (i in 1:100) {
    + a <- sim(i*10)
    + result1[i] <- a$sd3
    + result2[i] <- a$se}
    > plot(1:100, result2)
    > lines(1:100, result1/sqrt(1:100*20), col="red")
    そんないい加減でいいのか。。 -- 2005-10-12 (水) 12:33:13
  • まだ混乱があるようですね。y1~N1=pnrom(x,0,sd1) と書いたら、y1 は理論平均,"理論"標準偏差 sd1 の正規分布に従うデータという意味になりますが。sd1 が推定値という意味がよくわからない。 もしおっしゃる意味が yi~Ni=pnrom(x,0,σi)、 sdi はσi の(yi から求めた)推定値という意味なら、sqrt(sd1^2+sd2^2) の信頼区間というのは意味不明で、sqrt(σ1^2+σ2^2) の信頼区間という意味なのではないでしょうか。 -- MKR? 2005-10-12 (水) 14:55:59
  • すみません.全くそういうことです.>MKRさん.N3の理論標準偏差(σ1^2+σ2^2) の信頼区間をN1,N2の実測データから推定したい,ということです. -- takahashi? 2005-10-12 (水) 15:24:19
  • おさらいしておきます。x1,x2,...,xn ~ N(0,σ1^2), y1,y2,...,yn ~ N(0,σ2^2) とします(データ数は同じとします)。x データは互いに独立、y データも互いに独立、さらに x,y データどうしも互いに独立、とします。この仮定の下では、zi=xi+yi はN(0,σ1^2+σ2^2) に従う互いに独立なデータになりますから、標準的な方法で sqrt(σ1^2+σ2^2) の信頼区間が求められます。もちろん x データと y データをどうくみあわせるかは便宜的で、どうでもよいわけです。もしデータ数がことなるなら、勿体ないがどちらか余分な方を捨てるのでしょうか。捨てない方法は知りません。Q & A コーナー向きの話題でしたね。
    -- [[MKR]] &new{2005-10-12 (水) 18:46:59};
  • よくわかりました.今回最初から私の説明が悪く,皆様に余計な誤解を与えてしまったようです.具体的にはxデータ,yデータは確率変数ではなく,それぞれの水準における比率です.以下に具体例を示します. -- takahashi? 2005-10-13 (木) 16:11:15
  • 比率が正規分布するんですか。 -- 2005-10-13 (木) 16:33:22
  • たとえ話とか,話を簡単にするとかは,論者がよっぽどの技量がないとうまく行かないものです。うまくいかないどころか,まるっきりあさっての方向へ突っ走ってしまい,あらぬ解が得られたと言うことに質問者も回答者も気づかないという,悲劇が往々にして訪れる。 -- 2005-10-13 (木) 21:43:36
  • Rottenbach -- Austria? 2005-10-22 (土) 03:03:24

Ruby-RMathlib

るびお? (2005-10-05 (水) 23:37:07)

ここのページに書いてあるRuby-RMathlibを使ってみたいのですが、
http://www.okada.jp.org/RWiki/index.php?cmd=read&page=Ruby

そのためには、libRmath.a standalone libraryをCRANからとってこなければならないようなのです。しかし、どうしても見つかりません。
どこでこのライブラリのwindows版もしくはcygwin版等を手に入れることができるのでしょうか?
また、Ruby-RMathlibを使うのに、R自体はインストールする必要はあるのでしょうか?

環境はWindows XP + cygwinです。

  • Rのソースをとってきて,Rを構築(Mingw)した後,src/nmath/standaloneに移って make -f Makefile.win(この時cygwinにすればいいのかも) で *たぶん* 出来ると思います.(旨く行かなきゃ多分ちょっとヘッダを弄る必要があるかも). あと, standalone libRmath の乱数生成は "Marsaglia-Multicarry" なので, あとでR(標準は"Mersenne-Twister")での結果と比べると頭をかかえるかもしれません.(このへんはあまり力は入ってないようで...) -- なかま 2005-10-06 (木) 10:27:16
  • VineLinuxだったらlibRmath-devel-2.1.1-0vl1.i386.rpmに/usr/lib/libRmath.aは入っている飲んですけどね:-p -- 2005-10-07 (金) 13:25:32
  • なかまさんのやり方でできました。ありがとうございました。しかし本物のRに比べるとできることはすごく限られますね。RRbというのができるともっといろんな機能が使えるようになるのかなあ。開発進んでないみたいだけど… -- るびお? 2005-10-08 (土) 12:58:48

パッケージの読み込みができない

hatena? (2005-10-03 (月) 16:56:52)

何故か、lme4のパッケージをインストールした後、読み込みができず、利用できません。

library(lme4)と入力すると、
エラー:'%s' が要求したパッケージ '%s' は見つけられませんでした
という表示が出ます。

何か解決法がございましたら、よろしくお願いいたします。

  • 使用している OS, R バージョン は最低述べるべきです。 インストールした、というのは本当ですか? エラーメッセージは(R から参照できるようには)インストールされていないといっているように見えますが。「インストールした」というのは具体的に何をしたのか説明して下さい。 -- QDU? 2005-10-03 (月) 22:38:23
  • lme4はパッケージMatrixとlatticeに依存しています。すでにインストール済でしょうか。 -- 2005-10-03 (月) 23:45:55
  • require(lme4) としてみるとか. -- なかま 2005-10-04 (火) 09:34:18
  • Matrixがインストールされておりませんでした。無事に、読み込みができました。ありがとうございました。初歩的なミスです。気を付けます。 -- hatena? 2005-10-04 (火) 15:50:06

The R Tipsの正誤

mori? (2005-09-28 (水) 16:54:19)

船尾さんのThe R Tipsの練習問題の解答(p.360)にt.test(data, mu=0)の結果として
p値が0.8417なので、母平均が0であるという帰無仮説は棄却される
と書いてあるのですが、p値>0.05の場合は帰無仮説は棄却されないと思うのですがこの記述は間違いでしょうか。アドバイスお願いします。

  • mori さんの解釈が正しいです -- 2005-09-28 (水) 17:37:48
  • すみません・・・。moriさんのおっしゃる通りです。拙著は多数の間違いがありますので、サポート頁 を参照いただけましたら幸いです。-- 舟尾? 2005-09-28 (水) 18:31:58
  • サポートページがあると思って探していたのですが見つけられませんでした。アドバイスありがとうございます。 -- mori? 2005-09-29 (木) 08:01:20
  • 一番最後のページに載ってましたね・・・ -- mori? 2005-09-29 (木) 08:03:00
  • ご参照いただきましてありがとうございます。 -- 舟尾? 2005-09-29 (木) 21:46:00

Rmapのインストールについて

ひろ?? (2005-09-28 (水) 14:17:41)

はじめまして。WINDOWSで、R2.1.1に、Rmapのインストールを試みています。Rmap_1.1.0.zipをダウンロードしてきて、ローカルからのインストールまではうまくいったみたいなのですが、library(Rmap)を実行すると、
「'Rmap' は有効なパッケージではありません。バージョン < 2.0.0 でインストールされたかも知れません」
というメッセージがでてきてしまいます。
ちなみにbinフォルダにprojとshapelibは入れています。
どのように対処すればよいかどなたかご教示いただけないでしょうか。
よろしくお願いします。

  • これはバイナリパッケージがR2.0以降には対応していないためで、Rmapのパッケージを再ビルドする必要があるようです。一番簡単なのは旧版の1.9x等を探して(The R BookのCD-ROMにある?).これをインストールして、Rmapをロードすることでしょう。 -- 2005-09-28 (水) 15:28:06
  • 1.9xでインストールして使うことができました。ありがとうございました。 -- ひろ? 2005-09-29 (木) 14:48:54
  • すでに解決済みとこのことですが、OGRがどうしても必要だと言う訳でもないなら、maptoolsはどうでしょう。(OGRサポートを除いて)Rmapと比べて見劣りがするパッケージではありませんし、もちろん最新のRで動作します。 -- 谷村 2005-09-29 (木) 18:43:41

Rのソースを自動インデントするエディターはありますでしょうか

ちょろちゃん? (2005-09-21 (水) 23:58:08)

R駆け出し中のちょろと申します。
この度はRのソースを自動インデントするエディターを
ご教授いただきたく投稿しました。

現在、勉強もかねてRのスクリプトを作成しております。
そのスクリプトですが、数百行を超えてきたあたりから、
for,ifなどのカッコの対応を追うのが非常に困難になってきました。
そこで、Javaで言うところのEclipseのソース→フォーマットのような
既存のソースをオートインデントするソフト、ないしはエディタが
何か無いかと思い、xyzzyのR-modeの設定を行いましたが、

URL http://plaza.umin.ac.jp/~takeshou/xyzzy/rmode.html

R-modeではカラーリングはばっちりされるのですが、いかんせん
インデンテーション機能が弱いように感じております。


何か、既存のRのソースを整形するよりよいエディターやソフトウェア
などありましたら、ご教授いただきたく考えております。
どうぞヨロシクお願いいたします。m(_ _)m

  • perl とか AWK とかで,自分好みのインデンてーションをするプログラムを書くとか。
    google で site:r-project.org indent などとすると,いくつかみつかるかも ess にも s のインデンてーションをする機能があるとかという情報が。。 -- 青木繁伸 2005-09-23 (金) 08:52:51
  • ESSを使うと、括弧の上にカーソルがあるとそれに対応する括弧をハイライトしてくれて、対応を追うのが困難になることはありません。そのため、後付のインデントを考えたことがないけど、ESSにもその機能があるかも。ちなみに、ESSのインデント機能に今のところ不満はありません。EclipseのRのプラグインは試してみられましたか? -- 谷村 2005-09-23 (金) 09:31:40
  • ESSでご要望のことができるようです。ファイル(バッファ)全体を選択してから、C-M-?でインデント整形されます。if文やfor文などブロックの中だけならC-M-qで、1行だけならTABです。自動インデント機能のおかげで後からインデント整形する必要がほとんどなく、いままで気が付きませんでしたが、インデント整形したい行のどこにカーソルがあってもTABキーで(TABがその場に挿入されるのではなく)ちゃんと整形できるのですね。これは楽だ。-- 谷村 2005-09-23 (金) 10:01:54
  • 先輩方いろいろとアドバイスありがとうございます。まずはアドバイスいただきましたESSを使ってみようと思います。本当にありがとうございます。大変感謝しています。m(_ _)m -- ちょろちゃん? 2005-09-24 (土) 14:08:12
  • R-mode、放置していてすみません。インデントの部分については、僕も弱いなぁと思っているのですが、今は手を出せない状況です。もし代替コードなどありましたら、いつでもご連絡ください。 -- 竹内? 2005-09-29 (木) 16:24:01
  • 昨日、ノムラさまにまたコードをいただきました。マージする時間がなかったので、アップロードさせていただいただけですが、よろしければどうぞ。 -- 竹内? 2005-09-30 (金) 20:02:20

multiv パッケージがない

lost? (2005-09-21 (水) 12:10:56)

 multiv パッケージがいつのまにかCRAN から消えています。いといろと検索しましたが行方がわかりません。

 どなたか、ご教示願います。

R Site Search によれば

multiv is ORPHANED.
If you want it you can go the orphaned subdirectory on CRAN 

メンテナーがいなくなったか、他に良いものがあるのか、どちらかでしょう。 -- MKR? 2005-09-21 (水) 12:18:32

  • 貴重な情報ありがとうございました。 -- lost? 2005-09-21 (水) 13:08:20

「Rで自己組織化マップを」で、library(som)が呼べない

K.Morgen? (2005-09-14 (水) 16:41:34)

「Rで自己組織化マップ」を参照し、前段のlibrary(class)から、掲載中のプロセスは、出力まで、確認されました。レーダチャートの数値に、若干の相違がありましたが。 次に、後半の例に従って、library(som)を実行しましたが、”som"のパッケージはない、とのエラーメッセージが出ます。利用可能なRパッケージリストを開いて見ますと、確かに、somの名称のパッケージは見当たりません。
現在、R2.1.1Patchedを使用中ですが、インストールに問題が有ったのでしょうか。

  • 同じバージョンを使用していますが、"som"は筑波のミラーからDL・インストールできましたが。 -- 2005-09-14 (水) 16:50:54
  • 早速のアドバイス恐縮です。私も、R2.1.1Patchedは、筑波のミラーから、DL。インストールしていますが、貴殿のご指摘の内容は、"SOM"についてのお話でしょうか。R2.1.1Patchedに"SOM"パッケージが常駐しているのではないのでしょうか。以上ご教示下さい。 -- k.Morgen? 2005-09-15 (木) 07:47:36
  • install.packages("som")は実行済みでしょうか? -- 2005-09-15 (木) 08:04:16
  • ご指摘感謝申し上げます。お手数をお掛けし申し訳有りません。今後ともよろしくお願い申し上げます。 -- k.Morgen? 2005-09-15 (木) 08:35:37

OS XのRcmdrのフォント

ショーゾー? (2005-09-12 (月) 10:09:31)

Mac OS X 10.4.2 + X11 1.1 + R 2.1.1 + Rcmdr 1.1.1の環境です。Rcmdrで日本語が出たのには驚きでしたが,表示されるフォントが汚いです。これをWin環境でのRcmdrのように奇麗なフォントで表示させるにはどうしたらいいでしょう?
Rcmdrの「ツール」->「オプション...」->「デフォルトのフォント」には「1.0」とありますが,ここの修正でどのように例えばOsakaフォントを指定したら良いのでしょうか?
生半可知識ですが,X11のfontsフォルダにKochi-substituteをコピーしてmkfontdirを実行したのですが,fonts.dirに反映されませんでした。よろしくご教示お願いします。

  • 一度Commanderを終了してから、options(Rcmdr=list(default.font="-misc-fixed-medium-r-normal-ja-18-*-100-100-c-*-*-*")) を実行して、Commander()で再起動してみると...少しはきれいでしょうか(^^; -- 2005-09-12 (月) 17:21:20
  • できました! しかし15インチディスプレイでははみ出てしまいます。で,18(多分,フォントサイズ)を13か14にしたら丁度いい具合になりました。R-intro(日本語版)のAppendix Bにしたがって,このoptions()をRを起動するディレクトリ内の.Rprofileに書いておきました。これで,やや奇麗なRcmdrが立ち上がりました。 -- ショーゾー? 2005-09-12 (月) 23:37:32

error message "x and y lengths differ"どうすれば?

yanagi? (2005-09-10 (土) 17:07:37)

g<-expression(2*x^2+3*x+5)
D(g,"x")
y<-D(g,"x")
plot(y,-10,10)
をしましだがerror message "x and y lengths differ"がだます
どうすればいいですか?

  • yanagiさんの実行したいことはyで定義した関数を-10から10の範囲でプロットしたいということでしょうか?
    curve(2*x^2+3*x+5, xlim=c(-10,10));
    ではだめですか? -- Akira? 2005-09-10 (土) 20:28:22
  • y(というかD()のreturn value)はfunctionではなくcallなので,plot(function(x){eval(y)},-10,10) かな -- takahashi? 2005-09-10 (土) 22:31:37

APLとRの関数の比較表

Fleischmann? (2005-09-08 (木) 13:11:16)

 昔のAPLのコードをRに変換したいのですが、APLとRの関数の比較表のようなもの、ないでしょうか?

大規模データからサンプリングせずにデータ抽出する

msasaki? (2005-09-08 (木) 11:23:54)

6000行×8列の時系列データから、ある列から100個のデータを抽出しようとしています。Samplingを使うとランダムになるため、うまく行きませんでした。関数を探していましたが無いようでして、何か良い方法はございませんでしょうか。皆様におかれましては、大変ご多忙中のこととは存じあげますが、どうぞご指導のほど、よろしくお願い申し上げます。

  • 質問の意味が良く分かりませんが。x[1:100, i] とすれば,i 列の 1〜100 行目の要素が取り出せるというのはご存じ?x[c(...), i] とすれば,c( ) の中に書いてある要素を取り出せますし。何がやりたいのかもう少しちゃんと書けばよろしいかと思います。 -- 青木繁伸 2005-09-08 (木) 11:32:30
  • もうしわけございませんでした。6000個のデータから等間隔で100個のデータを抽出したいのです。 -- msasaki? 2005-09-08 (木) 11:38:06
  • 上に書いたものの応用でしょう。x[seq(start, end, by=n), i] を試してみたらいかがでしょう。 -- 青木繁伸 2005-09-08 (木) 11:49:07
  • ありがとうございました!!おかげさまで1歩前進しました!!青木先生の御書きになられましたコメントはいつも拝読させていただいております。先生にコメントしていただいて感激です。本当にありがとうございました。 -- msasaki? 2005-09-08 (木) 11:56:02

RでKNN(K近傍法)を行うには?

nori? (2005-09-08 (木) 09:48:14)

はじめまして。1000*65のデータの欠損値をKNNで解析を行いたく、今回初めてRにふれている者です。パッケージEMVをインストール後の作業について参考になるページ等ございましたらよろしくお願いいたします。

  • 書くまいと思ったものの。。。どのようなことを聞きたいのかちゃんと買いたほうがよいでしょう。あなたの参考になるかならないかわからないページを薦めることもできないでしょう。どんなことでも参考になりますというわけではないでしょう? -- 青木繁伸 2005-09-08 (木) 11:35:08
  • 青木繁伸様、失礼いたしました。しかしながら、Rについて基本しか知らないためソースコードを見ながらどの部分を改変すれば良いか検討しようと思っていました。お気を悪くさせてしまいすみませんでした。 -- nori? 2005-09-08 (木) 13:02:48

c() の最大の要素数

初心者 ち? (2005-09-06 (火) 12:35:22)

c()の引数に16000個の引数を取ることは可能でしょうか?エラーが出ているのですが,要素数の最大値の調べ方が分かりません.google等でしらべたのですがわかりません(^ ^;)

obs1 <- c(4,1,2,1,1,19,5,1,1,1,6,19,1,1,1,10,1,1,4...)

syntax error
Execution halted

  • 仕様では見かけていない(見落としているだけかもしれませんが)のでメモリとアーキテクチャ次第でしょうかね
    syntax errorなのだから文法が間違っているんじゃないでしょうか?少なくとも16000くらいだったら平気でしょう~ 手元の環境では
    eval(parse(text=paste("c(",paste(1:16000,collapse=","),")")))
    は問題なく通ります
    ソースコードを載せれば誰かがどこが間違っているか指摘してくれると思います(16000個のデータをそのまま載せたりはしないで下さいね)-- takahashi? 2005-09-06 (火) 13:39:56
  • 私のところはMacintosh OS 10.4 で,R 2.1.1 ですが,10倍の 16000 個の引数をとっても,平気でしたよ。あなたの環境 OS, R のバージョンはどのようになっているのですか? -- 青木繁伸 2005-09-06 (火) 13:44:23
  • バッファの問題は異なるエラーになります(見落としが無ければ).syntaxなら、括弧やカンマが全角文字の可能性が高いですね。 -- なかま 2005-09-06 (火) 14:03:00
  • OSは TurboLinux? WorkStation8でRのバージョンは2.0.1 メモリは2Gです.R -no-save < program.Rとして読み込んでいます.program.Rのないようですが,obs1 <- c(4,1,2,1,1,19,5,1,1,1,6,19,1,1,1,10,1,1,4...(省略))と16000ぐらいのデータが並んでいます.eval(parse(text=paste("c(",paste(1:16000,collapse=","),")")))は通ったので,バッファの問題ではないようです.括弧やカンマかもしれないので -- 初心者ち? 2005-09-06 (火) 17:20:40
  • 同じデータでWindows XPで実験したところSyntaxは出ませんでした.どうやら括弧やカンマではないようです. -- 初心者ち? 2005-09-06 (火) 17:41:13
  • 改行コードかな?そうであればprogram.Rをobs1<-c()としてもエラーになるはずなので試してみたらいかがでしょう?データが多すぎることが問題なのかそれ以外のことが問題なのかわかるとおもいますので -- takahashik? 2005-09-06 (火) 18:10:19
  • 何より、16000 個もの数字を書き連ねたソースファイルというのが混乱の源のような気がしますね。分割して読み込み、必要なら R 内部で結合するか、それがだめなら、ソースファイル自体を見やすい(つまり間違いが見付けやすい)形に清書しておくか。天下りデータでそうはいきませんか? -- [[ MKR]] 2005-09-06 (火) 23:40:16
  • obs1<-c()はうまくいきます.私自身も改行コードかなっと思ったのですがnkfで変換してもエラーが出るのでそうではないようです.根本の解決にはなっていませんが,分割読み込みでやってみます.(実はJavaでobs1 ->c(...)を自動的に生成してまして半角文字などは入っていないのです)皆様いろいろありがとうございましたm(= =)mまたなにか原因が分かりましたら報告させていただきます. -- 初心者ち? 2005-09-07 (水) 12:38:10
  • 追記:Rを起動したあとにsource("foo.R")とするとうまくいきますが,R --no-save < foo.R だとうまくいきませんでした.JavaからRを呼び出すのにsource()のほうを使って解決できそうです.ありがとうございました&ご迷惑をおかけしました -- 初心者ち? 2005-09-07 (水) 13:03:16
  • 一行に長く書きすぎたのではないでしょうか.http://cran.r-project.org/doc/manuals/R-intro.html#R-commands_003b-case-sensitivity-etc (Command lines entered at the console are limited to about 1024 bytes) -- numata? 2007-03-09 (金) 15:22:11

連番の変数の参照

Katsura? (2005-09-05 (月) 13:22:09)

連番の変数(m1〜m30)のそれぞれに対数線形モデル(loglin)の結果を代入したのですが、その成分m1$lrt,m2$lrt,...,m30$lrtという値を一括して取り出したいのです。いままでは、いちいちm1$lrt、m2$lrtと30行続けて入力していたのですが、もっと賢い方法がありそうな気がします。ご教唆いただけないでしょうか。

  • sapply(1:30,function(x){eval(parse(paste("m",x,"$lrt",sep="")))}) #まあ普通はmをリストにしたほうがいいと思いますが -- takahashi? 2005-09-05 (月) 14:18:26
  • sapply(1:30, function(i) eval(parse(text=paste('m', i, '$lrt', sep="")))) とか? -- 青木繁伸 2005-09-05 (月) 14:18:43
  • もろかぶりしましたね。
    ちなみに,教唆とは「1 他人をそそのかすこと。けしかけること。扇動すること。哄唆。2 犯意のない他人をそそのかして犯罪を行なわせること。」です -- 青木繁伸 2005-09-05 (月) 14:23:03
  • あーparseのtext=忘れてました(別に「教唆」したわけではないです) 青木先生の方が正しいです  -- takahashi? 2005-09-05 (月) 14:29:36
  • m<-as.list(NULL);m[[ 1]]<-loglin(なんちゃら);...;lapply(m,function(x){x$lrt}) 最初のはリストの器の用意で,あとは器の中に入れ込んで,最後にapplyで取り出すのが多分正統派「前のは教唆された亜流...:-)」だと思います. -- なかま 2005-09-05 (月) 14:32:49
  • ありがとうございます。これでだいぶ楽になります。ちなみに「教唆」という意味、いままでずっと間違えてメールなどで使っていました。恥ずかしい。 -- Katsura? 2005-09-05 (月) 16:23:43

凡例に+-記号を

mitsu5? (2005-09-04 (日) 11:34:42)

凡例に mean +/- SD と書きたいのですが、このプラス マイナス を1つの記号としてどうして入力するのでしょうか? 使用しているのは英語版のRです。faqかも知れませんが、うまく探せませんでしたのでよろしくお願いします。

  • plot(1:10,1:10) ;legend(8,8,expression(mean %+-% SD)) どんな数式が使えるかは,demo(plotmath)等としてみてください -- なかま 2005-09-04 (日) 11:57:27
  • どうもありがとうございます。help("plotmath")ではぴんと来ないことが、demo(plotmath)でよく判りました。 -- mitsu5? 2005-09-04 (日) 12:46:38

横軸の下とか縦軸の左に文字を描き込みたい

青木繁伸 (2005-09-04 (日) 11:25:25)

横軸の下に,何か描こうとしても,描けません。

x <- 0:10
y <- x^2
plot(x, y)
text(5, -2, "a")

のようにしても,ウインドウサイズにもよりますが,描けません(上の方だけ見えたり,全く見えなかったり)。~
どうしたらいいでしょうか。

  • par(xpd=T) として,textしてみてください. -- なかま 2005-09-04 (日) 11:48:03
  • ありがとうございました。できました。安易に聞いてしまいました。反省。 -- 青木繁伸 2005-09-04 (日) 12:33:56

重回帰分析におけるマルチコについて

TAK? (2005-08-29 (月) 20:46:25)

こんばんは。申し訳ありませんが、質問させてください。

Y X1 X2 X3
18.235 13.256 0.632 0.011
30.444 20.668 0.592 0.017
10.529 7.167 0.654 0.018
19.579 13.307 0.668 0.025
27.143 20.413 0.534 0.006
15.000 12.362 0.515 0.010
17.433 13.061 0.448 0.002
30.250 22.299 0.692 0.008
30.000 22.000 0.500 0.007
18.571 13.582 1.311 0.832
31.026 20.697 1.569 0.585
11.028 7.558 1.526 0.702
19.785 14.200 1.138 0.436
27.452 20.470 1.187 0.525
15.598 12.604 0.784 0.525
17.692 13.388 0.684 0.863
30.928 22.752 1.320 0.266
30.086 22.558 0.688 0.434
19.507 14.512 2.239 1.784
31.082 21.684 2.031 0.671
11.577 7.855 1.751 0.888
19.952 14.945 1.900 0.909
27.970 20.639 1.310 0.671
16.577 13.166 1.094 1.473
18.017 13.711 1.623 1.791

上記のデータに対して、目的変数をY、説明変数をX1、X2、X3とした時に、全説明変数で回帰させたところ、X2の偏回帰係数の符号はプラス、Yに対するX2の単相関係数はマイナスとなりました。そこで、マルチコが起こっているものだと考え、http://aoki2.si.gunma-u.ac.jp/R/find_multico.html にあるfind.multico()関数を用いて、X1、X2、X3に対してマルチコのチェックを行ったところ、結果は"no multico"となりました。気になって、VIFも計算したところ、どの変数も皆3より小さい値となっていました。このような場合に、「マルチコという現象は起こっていない」とみなして良いのでしょうか?また、このデータから得られたY=aX1+bX2+cX3+dという重回帰式は、安定した、将来の予測に有効な式となりうるのでしょうか?どなたかお考えを教えてください。よろしくお願い致します。

  • 参照された find.multico() は,固有値・固有ベクトルの方からマルチコをチェックするものです。若干基準が強すぎる。つまり,ややマルチコのうたがいありの場合でもチェックから逃れてしまう気配?
    一方,http://aoki2.si.gunma-u.ac.jp/R/tolerance.html の方によれば,
       tolerance    VIF
    X1   0.92483 1.0813
    X2   0.44101 2.2675
    X3   0.42018 2.3799
    となり,TAK さんがおっしゃるように X2もX3 もVIF は3以下ですよね。しかし,3 以下なら良いという基準もそんなに確かなことでしょうか?
    単相関係数と回帰係数の符号が違うのがもっとも明確な理由になると思います。 もっとも,X2 と Y の単相関は -0.0255 と非常に小さいので,このような独立変数を重回帰式に入れる必然性は全くない。もっと言えば,X3 も必要かどうか。。。
    ということで,ステップワイズ変数選択(Pin=Pout=0.05)をやってみると X1 しか採用されない。そうですよね Y と X1 の単相関係数が 0.988 もあるのですからね。自由度調整済みの決定係数を見てみると,
    X1 のみ  0.97442
    X1, X3  0.97401
    X1, X2  0.97466
    X1,X2,X3 0.98074
    ということです。
    感じとしては,X1 のみで十分かな。 -- 青木繁伸 2005-08-29 (月) 21:51:44
  • 予測に徹するなら、計算機が計算を放棄しない限りマルチコ完全無視で変数を投入していいのでは?因果分析等がしたい訳じゃなさそうですし。 -- 2005-08-30 (火) 08:56:34
  • 青木先生およびもう一方、有意義なわかりやすいコメントをしていただきまして大変感謝しております。お二方のおっしゃることで頭がすっきりしました。本当にありがとうございました。 -- TAK? 2005-08-30 (火) 09:09:10
  • そのような方針もありかもしれませんね。しかし,マルチコ状態では,たとえシンギュラーマトリックスということで計算がストップしなくても,得られる結果は不安定です。相関係数のいくつかがほんの少し変化しても結果がかなり大きく違うと言うこともあります。そういうことで,マルチコを避けるのは,安定した結果が得られることを担保することにもなるのです。 -- 青木繁伸 2005-08-30 (火) 09:10:53
  • 「結果が不安定でも予測制度は良い」というのが最近の潮流のような気もしますね。 -- 2005-08-30 (火) 23:58:10

プロットの色分けの方法

阿部? (2005-08-29 (月) 14:34:46)

質問です。

  代金   緯度   経度
 48000	139.77379 35.70892
103000	139.77705 35.69388
 68000	139.7803  35.69616
103500	139.7728  35.696

のようなデータをRのplotを用いて緯度経度をプロットし、代金によって色分けしたいのですがどのような方法がありますでしょうか。
たとえば代金が、50000未満の点を青、50000以上,80000未満の点を赤、80000以上の点を黒という風に表示したいと思っております。よろしくお願いいたします。

  • この例題は理解できますか?
    x <- 1:20
    color <- ifelse(x < 6, "blue", ifelse(x < 15, "red", "black"))
    plot(x, x^2, col=color)
    2行目で color という変数に付値してから使っていますが,直接使ってもよいですよ。
    わかりにくければ,for 文と if 文を使ってプログラムすればいいです。 -- 青木繁伸 2005-08-29 (月) 15:36:15
  • ご回答ありがとうございました。
    先生の回答を参考にさせていただき、無事色分けすることができました。
    ご指導ありがとうございました。 -- 阿部? 2005-08-29 (月) 15:59:00

t.testでの母平均算出

Akira? (2005-08-17 (水) 20:06:29)

matrixデータに対して、row毎にt.testを適用し、estimateを取り出そうとしています。
すると、【以下にエラーt.test.default(x) : データは本質的に定数です】とエラーになります。
おそらく、rowの値が全て等しいからだと思うのですが、t.testのヘルプにはそれを回避するオプションがあるか分かりませんでした。ご教授いただけないでしょうか?
今は以下のようにしていますが、上手く行きません。

> # このようなデータです。
> x
        x1       x2       x3
a 7.864811 7.402886 7.633848
b 8.210283 7.838808 8.024545
c 5.624960 5.624960 5.624960
d 3.451917 3.491341 3.471629
e 9.493148 8.908723 9.200936
f 7.759156 8.401847 8.080501

> #関数を定義しました
> fun <- function(x, y){
+   if(all(x==mean(x, na.rm=TRUE), na.rm=TRUE)){
+     x <- mean(as.numeric(x), na.rm=TRUE)
+    }else{
+     x <- t.test(x)[[y]]
+    }
+   return(x)
+ }
>
> #そしてapplyでrow毎にfunを実行しました
> apply(x,1,fun,5)
>
> #そうすると、以下のエラーが出てしまいます
以下にエラーt.test.default(x) : データは本質的に定数です

3行目の値が全て同じなのが、原因と分かりました。

> x.row3 <- x[3,]
> x.row3
       x1      x2      x3
c 5.62496 5.62496 5.62496
> apply(x.row3,1,fun,5)
以下にエラーt.test.default(x) : データは本質的に定数です

しかし、人為的に作成したデータでは上手く動きます

> test <- matrix(c(1:12, rep(3,3)), ncol=3, byrow=T)
> test
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9
[4,]   10   11   12
[5,]    3    3    3
> apply(test,1,fun,5)
[1]  2  5  8 11  3

all(x==mean(x, na.rm=TRUE), na.rm=TRUE)の条件分岐が上手く動いていないと思いますが、Rではx.row3は全て同じ値に見えます。小数点以下の表示していない桁でおかしくなるのでしょうか?

> str(x.row3)
`data.frame':   1 obs. of  3 variables:
 $ x1: num 5.62
 $ x2: num 5.62
 $ x3: num 5.62
> x.row3[1]
       x1
c 5.62496
> x.row3[2]
       x2
c 5.62496
> x.row3[3]
       x3
c 5.62496

よろしくお願いします。

  • 誰かが, 行を潰してた。。ブツブツ。。
    > mean(x, na.rm=TRUE)
          x1       x2       x3 
    7.067379 6.944761 7.006070 
    となるのは分かっているのでしょうか。
    そもそも,あなたがやろうとしていることがよく分からないのですが「row毎にt.testを適用し、estimateを取り出そうとしています」というのはどういうことでしょう。estimate ってなんですか。うまく動く例の結果を見ると,行ごとの平均値を求めているのですか??
    だとしたら,エラーチェック・エラー対処は別として,apply(x, 1, mean) ではなぜいけないんでしょうか? -- 青木繁伸 2005-08-17 (水) 20:56:12
  • ありがとうございます。「estimateを取り出す」という部分ですが、「t.test(x)$estimateを取り出す」ということをお伝えしたかったのです。青木先生のご指摘通り、実作業ではapply(x, 1, mean)で解決します。難しく考えすぎていました。
    上のエラーはt.test(x)$p.valueを取り出す関数を使いまわしたら発生したので、疑問に思って質問させていただきました。
    > mean(x, na.rm=TRUE)
          x1       x2       x3 
    7.067379 6.944761 7.006070 
    については、理解しております。x <- mean(as.numeric(x), na.rm=TRUE)としておりますが、それならばx==mean(as.numeric(x), na.rm=TRUE)にするべきでした。
    ここが【以下にエラーt.test.default(x) : データは本質的に定数です】の原因なのでしょうか? -- Akira? 2005-08-17 (水) 23:17:12
  • 確かに意味不明ですが、もし質問の趣旨が(本質的に同値な)行を弾き出すという趣旨なら-- QDU? 2005-08-17 (水) 23:22:41
    if(all(x==mean(x, na.rm=TRUE), na.rm=TRUE))
    の部分を
    if(identical(all.equal(var(x[2,]), 0), TRUE))
    にでもかえれば良いはずです。投稿が前後してしまいましたが、本質的に等しい(つまり浮動小数演算の意味で本質的に等しい)データの分散はゼロですから、t 検定など出来るはずがありませんから、エラーになるのは当然です。
  • QDU先生ありがとうございます。エラーなく動作しました。
    if(identical(all.equal(var(x[2,]), 0), TRUE))
    のような使い方を知りませんでした。
    「xに含まれる要素が全て等しいかどうか」を判定したかったのですが、identical(a, b)という風な2つの組合せを作れずにいました。a=identical(all.equal(var(x[2,])でb=0とすれば良かったのですね。事前に検索をしっかりします。
    どうでも良いと思いますが、x==mean(as.numeric(x), na.rm=TRUE)ではエラーになって動きませんでした。これは青木先生のご指摘にありましたmean(x, na.rm=TRUE)の意味が分かっていなかったからです。すみません。 -- Akira? 2005-08-18 (木) 08:37:35

度数分布の求め方

A380? (2005-08-16 (火) 23:40:38)

 度分布を計算するRの関数はないでしょうか?

  • table とか hist とか。
    > x <- rbinom(10000, 10, 0.5)
    > table(x)
    x
       0    1    2    3    4    5    6    7    8    9   10 
       9   87  439 1142 2140 2483 1993 1162  422  112   11 
    > hist(x, breaks=0:11-0.5)$counts
     [1]    9   87  439 1142 2140 2483 1993 1162  422  112   11
    オンラインヘルプを良く読んで使ってください。 -- 青木繁伸 2005-08-17 (水) 08:10:42
  • 御回答ありがとうございました。 -- A380? 2005-08-17 (水) 08:50:46

計算結果の画面表示の保存

take? (2005-08-16 (火) 20:21:32)

お世話になります。
例えば下記のようなリストの表示そのものを保存するコマンドはどなたかご存知ですか?
コピー&ペーストがもっとも手っ取り早いですが、大量にある場合、a b 別々にwrite.table()と保存してDOS窓でcopy/bによって結合するという方法もありますが、R上で行ってしまいたいです。
capture.output()などではベクター形式で表示されてしまいますし、sinkではRで使うには便利ですが、テキストで別のソフトに移すときは不便です。
なお、当方の環境は2.01 WinXPsp1です。どうぞよろしくお願いします

> a <-matrix(1:4,ncol=2); b <-matrix(4:1,ncol=2)
> list <-list(a,b)
> list
[[1]]
     [,1] [,2]
[1,]    1    3
[2,]    2    4

[[2]]
     [,1] [,2]
[1,]    4    2
[2,]    3    1
  • どのような形でファイルに出力したいのか不明なので,なんともいえませんし,あなたが候補に上げたが棄却した理由もはっきりしない。あなたの望むようにはき出す関数を作れば良いだけでしょう。その関数を print メソッドで準備すれば,いちいち意識して print しなくても良いし,cat を使えばどのような形式ででもファイルに書き出すことができますよ。
    > a <-matrix(1:4,ncol=2); b <-matrix(4:1,ncol=2)
    > list <-list(a,b)
    > class(list) <- "your.class"
    > print.your.class <- function(x)
    + {
    + 	cat(list[[1]][1,1], list[[1]][1,2], list[[1]][2,1], list[[1]][2,2],"?n")
    + 	cat(list[[2]][1,1], list[[2]][1,2], list[[2]][2,1], list[[2]][2,2],"?n")
    + }
    > list
    1 3 2 4 
    4 2 3 1 
    cat の出力先をファイルにする,ファイルに追加出力するとか,なんでもかんでも後は自分で色々改造。自分でなんでもかんでも好き勝手にできると言うところが R のよいところ。 -- 青木繁伸 2005-08-16 (火) 21:14:16
  • 青木先生早速ありがとうございます。わかりにくい書き方ですみませんでしたが、質問の出力形式は上記質問のRの画面表示自体でした(>list以下の9行)。cat()によって解決しますしこの方法は応用がきき大変助かります。また、Rでのオブジェクト思考的(?)な概念の例は初めてみましたので大変勉強になりました。 -- take? 2005-08-16 (火) 22:06:23

パッケージが見つからない。

econome007? (2005-08-15 (月) 17:19:25)

初心者です。
MASS,nnet,ts,nlsといったパッケージが見あたらないのですが、どこからダウンロードできるのでしょうか。

  • MASSとnnetについては、VRというパッケージの中にバンドルされていることが分かりました。 -- econome007? 2005-08-15 (月) 17:38:47
  • ts も nls もパッケージではなく,関数ですが。いずれも stats に含まれていますよ。 -- 青木繁伸 2005-08-15 (月) 17:50:39
  • 質問にお答え頂き、ありがとうございました。 -- econome007? 2005-08-16 (火) 10:49:15

ポアソン方程式を扱える関数

ラプラス? (2005-08-10 (水) 14:29:03)

 Rでポアソン方程式を扱える関数はあるのでしょうか?

 いろいろ探してみましたが見つからないのですが?

  • R は統計システムですから、これは中華料理店でフランス料理を注文するようなものでは。それでも、常微分方程式を数値的に解くパッケージ odesolve というのがありますが、ポアソン方程式は偏微分方程式ですね。他のソフトを検討された方が賢そうです。 -- QDU? 2005-08-10 (水) 20:04:00
  • ないものは,自分で書けば良いだけです。それも,何も R で書かなくても良いなら,C で書かれたプログラムがあるようですよ。その方面には全く暗いので,そのプログラムであなたの望むことができるのかどうか,それより前にそもそもあなたが C を解するのかどうかさえ分からないのですが。 -- 青木繁伸 2005-08-10 (水) 20:18:09
  • ROctaveでOctaveを呼んで,ん中でポアソンなりラプラスなり処理するとかは?.(R向けの処理で無いことはたしか) -- なかま 2005-08-11 (木) 09:10:11

R のスクリプトで再帰アルゴリズムを書く方法

Lost? (2005-08-05 (金) 17:17:45)

 R のスクリプトで、再帰アルゴリズムを使ったサンプルはないでしょうか?

  • どの程度のサンプルといっておられるのか分かりませんが。よくある,階乗を求める再帰関数の例は,以下のようになるでしょう。
    > fact2 <- function(n)
    + {
    +     if (n == 0) 1
    +     else n*fact2(n-1)
    + }
    > fact2(8)
    [1] 40320
    > fact2(3)
    [1] 6
    base に本来ある factorial 関数でチェックしてみましょうか・・
    > factorial(8)
    [1] 40320
    > factorial(3)
    [1] 6
    いかがでしょうか(レポートの手助けだったらやだな)。
    else n*fact2(n-1)
    のところは,
    else n*Recall(n-1)
    のほうがいいらしい。 -- 青木繁伸 2005-08-05 (金) 17:44:25
  • 青木先生ご返答ありがとうございます。
     上の関数を以下のように加算に変更してみました。
    fact3 <- function(n)
    {
    if (n == 0) 1
    else n+fact3(n-1) #加算に変更
    }
     R2.1.0 での実行結果は以下のようになりました。
    > fact3(998)
    [1] 498502
    > fact3(999)
    エラー:protect():プロテクションスタックが溢れました
     結果がオーバーフローしているようではないですが。関数内の"fact2"の替わりに"Recall"を使用すると、引数がより小さな値で、上のエラーが出てしまいます。
     上記のエラーへの対処方法はないのでしょうか?-- >Lost? &new {2005-08-05 (金) 20:04:44};
  • 加算に変更って。。。これは何を計算することになるんでしょうね。
    あなたは,分かってますよね。だとすると,この例の関数名は不適切ですね。そして,関数定義も不適切です。
    解が求まるまでに,何回再帰する必要があるかを計算することは可能ですね。あなたが変更した関数は,階乗を計算する場合に比べて遙かに大きいということになるでしょう。というか,fact2 はそんなに大きな引数に対応できないだけですが。
    再帰をするということについて,コンピュータサイエンスを分かっておられれば,スタック(作業領域)が必要ということはおわかりいただけると思います。要するに,用意していたスタック領域が足りなくなったということでしょう。そしてたぶん Recall はもっと余分なメモリーを要求するんでしょうね。
    解決法?
    再帰関数はスマートな定義ですが,計算上はちっともスマートではないので,再帰呼び出しをしない(繰り返し計算)に転換すると言うことだと思いますよ。
    スタックの変更方法は私は知りません。誰かが書いてくれるかも。 -- 青木繁伸 2005-08-05 (金) 21:35:07
  • 一般(?)に再帰とは鼠講に例えられて、階層が浅いうちは上手くいくんですが,資源には限りがあるので深くなると破綻します.src/include/Defn.hのR_PPSSIZEを増やせば良いですが,上記のような物ではあっと言う間に他の制限に引っかかってすぐに破綻するでしょう.言語の仕様や今日の計算機資源からしても、今のサイズでも十分です.上のようなものはたんなるループ処理が適切です. -- なかま 2005-08-06 (土) 00:19:24
  • 前にパーコレーション問題を解くプログラムをRを使って書けという問題を出したら、それまでRを一度も使ったことの無い学生が、翌週には立派なプログラムを提出して驚かされました。計算機言語の基礎を学んだ学生にはRはすぐに使えるようになるようです。ただし、全員が再帰的アルゴリズムを使っていて、どうしたことかとおもったら、アルゴリズム論の講義でしっかり刷り込まれていたようです。当然ごくごく小さなサイズで無いと破綻します。再帰的アルゴリズムを使わない(私なりに工夫した)模範例を紹介したら、すぐに数十倍も早いプログラムを書いてきた学生がいて、これまた驚かされました。話が脱線気味ですが、再帰的アルゴリズムをつかった(Rの)実用的な関数、再帰的アルゴリズムでないとどうしようも無い関数ってはたしてあるのでしょうか? ちなみに、R の階乗関数 factorial(x) は再帰的に計算したりせず、単に gamma(x+1) を計算するだけです。 -- 間瀬? 2005-08-06 (土) 02:00:23
    > factorial
    function (x) 
    gamma(x + 1)
    <environment: namespace:base>
  • 末尾再帰になってないからかと思ったらどうやらRは末尾最適化をサポートしてないんですね 純粋なfunctional language だと思ってたんですがそうでもないのか 通常の関数型言語だと
    function(a,b=1){ifelse(a==0,b,Recall(a-1,b+a))}
    とかにすればスタック消費しなくなるはずですけどね 例えばschemeとか -- takahashi? 2005-08-06 (土) 03:10:06
  • リコンパイル無しでも,R --max-ppsize=100000 とかの設定が可能のようです. Recallは,べ_な_ぶ_る_ずとり_ぷ_り_ー(spamのようだ)の仕組んだ黒○術?だそうで(訳には自信がもてないが...伏せとこう).実際R内の統計系のソース内では殆どRecallの使用は認められず(kernel,ts.smspline,mosaicplotぐらい)システムの維持用(○魔術の方が主に書かれている...)に多く使われています. -- なかま 2005-08-06 (土) 23:30:51

標準パッケージ

issei? (2005-08-04 (木) 14:35:12)

標準パッケージのリファランスはどこで得ることができるのでしょうか?
CRAN拡張パッケージについては、ダウンロードするところにあるのはわかるのですが。
spatialを使いたいのですが困っています。
よろしくお願いします。

  • spatial パッケージはパッケージバンドル VR の一部です。VR がインストールされていれば、library(spatial) で使えるようになります。内容一覧は命令 library(help=spatial) で得られます。但し、spatial パッケージよりも現在ではもっと良いパッケージ(例えば spatstat)等が複数ありますから、そちらを検討されることをお勧めします。 -- 間瀬? 2005-08-04 (木) 18:23:50
  • すばやい返答ありがとうございます。そうですか、他のパッケージを検討してみます。 -- issei? 2005-08-05 (金) 10:38:39

OSX R起動できない

いしだ? (2005-08-01 (月) 16:56:10)

あまりに初歩的なので、「何でも掲示板」に質問しましたが、指摘を受けて、こちらに移動しました。また、しばらく前に、青木先生の掲示板で質問したのですが、解決しませんでした。
MacOSX10.3.9でR2.1.1が起動しません。具体的には、起動させても、ドックの中でバウンドしたまま、予期せず終了してしまいます。何度インストールしてもだめです。
 新しいアカウントからログインすると、普通に使えます。また、自分の通常のアカウントでも、ターミナルで"r"と打つと、ターミナルの中で使えます。x11でも同様に使えます。
 なぜでしょうか。ライブラリの中のRとついてるものをけしたり、正常に動くパソコンからコピーしてきたりしたのですが、いっこうに、改善しません。何か、解決方法がわかれば、教えてください。

  • 以下をターミナルから実行するとどうなんでしょう, -- なかま 2005-08-01 (月) 18:35:12
    $ /Applications/R.app/Contents/MacOS/R
    で何かメッセージが出ないか.
    $ otool -L /Applications/R.app/Contents/MacOS/R
    ついでにライブラリのリンクがどうなるか. にわかMacユーザですんで,あれですが...
    $ gdb /Applications/R.app/Contents/MacOS/R
    デバッガのめっせーじ
    (gdb) run   # デバッガで実行
    なんか出てくる実行可能なら,Rの方が動くが、もし落ちてるようなら
    (gdb) bt    # とりあえず落ちるならバックトレース
    なにもかえって来ないようなら
    (gdb) [control+C] # シグナルを送って
    (gdb) bt    # とりあえずバックトレース
    よくわからなかったらすみません.
  • ":"は点々,"(hoge)"はこめんとのつもりだったんですが,疎通できませんようなので,もう一度上の部分を直しましたので見てもらえますか,"(gdb)はプロンプトですので入力しなくてよいです."#"以降はコメントになるはず. -- なかま 2005-08-01 (月) 21:13:35
  • お手数おかけして、申し訳ありません。以下のように試しました。お願いします。
    一文目(起動しようとして、予期せず終了)
    Trace/BPT trap
    二文目
    /Applications/R.app/Contents/MacOS/R:
    /System/Library/Frameworks/Cocoa.framework/Versions/A/Cocoa (compatibility version 1.0.0, 
    current version 9.0.0)
    /Library/Frameworks/R.framework/Versions/2.1.1/Resources/lib/libR.dylib (compatibility version 
    2.1.0, current version 2.1.1)
    /System/Library/Frameworks/WebKit.framework/Versions/A/WebKit (compatibility version 1.0.0, 
    current version 1.0.0)
    /System/Library/Frameworks/Security.framework/Versions/A/Security (compatibility version 1.0.0,
    current version 177.0.0)
    /System/Library/Frameworks/ExceptionHandling.framework/Versions/A/ExceptionHandling 
    (compatibility version 1.0.0, current version 4.9.0)
    /usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 71.1.3)
    第三文
    GNU gdb 5.3-20030128 (Apple version gdb-309) (Thu Dec  4 15:41:30 GMT 2003)
    Copyright 2003 Free Software Foundation, Inc.
    GDB is free software, covered by the GNU General Public License, and you are
    welcome to change it and/or distribute copies of it under certain conditions.
    Type "show copying" to see the conditions.
    There is absolutely no warranty for GDB.  Type "show warranty" for details.
    This GDB was configured as "powerpc-apple-darwin".
    Reading symbols for shared libraries ....... done
    runすると
    Starting program: /Applications/R.app/Contents/MacOS/R 
    Reading symbols for shared libraries +.................................................................... done
    Reading symbols for shared libraries . done
    Reading symbols for shared libraries . done
    Program received signal SIGTRAP, Trace/breakpoint trap.
    0x90a8d318 in _NSRaiseError ()
    (Rは長時間かけて起動したものの、全く反応無し。 ウィンドウも表示されず、メニューバーも出ない)
    btすると
    #0  0x90a8d318 in _NSRaiseError () 
    #1  0x90a8d1fc in +[NSException raise:format:] ()
    #2  0x90a2d35c in -[NSString stringByAppendingString:] ()
    #3  0x00006ce0 in -[RController doLoadHistory:] ()
    #4  0x00003ce0 in -[RController awakeFromNib] ()
    #5  0x90a5f0e8 in -[NSSet makeObjectsPerformSelector:] ()
    #6  0x92ea2150 in -[NSIBObjectData nibInstantiateWithOwner:topLevelObjects:] ()
    #7  0x92f93c2c in loadNib ()
    #8  0x92eeae24 in +[NSBundle(NSNibLoading) _loadNibFile:nameTable:withZone:ownerBundle:] ()
    #9  0x92f69d28 in +[NSBundle(NSNibLoading) loadNibFile:externalNameTable:withZone:] ()
    #10 0x92f7b51c in +[NSBundle(NSNibLoading) loadNibNamed:owner:] ()
    #11 0x92f69b90 in NSApplicationMain ()
    #12 0x00002dd0 in _start (argc=1, argv=0xbffffe60, envp=0xbffffe68) at /SourceCache/Csu/Csu
    -47/crt.c:267
    #13 0x8fe1a278 in __dyld__dyld_start ()
     -- いしだ? 2005-08-02 (月) 09:40:35
  • どこかでR.appがオーバーランでもしているような気がするんですが,ターミナルから
    export LANG=C
    export LC_ALL=C
    /Applications/R.app/Contents/MacOS/R
    では動かないでしょうか,動くなら,"C"を"ja_JP.UTF-8"に置き換えるといかがでしょうか -- なかま 2005-08-02 (火) 10:50:40
  • だめですね。"ja_JP.UTF-8"に変えても、~"Trace/BPT trap"と出て、R.appはすぐに予期せず終了です。とにかくいろいろと、ありがとうございます。 -- いしだ? 2005-08-02 (火) 11:54:28
  • なんとなくですが、これはリンク時の shared library と実行時の shared library のバージョンが違うからでは?動作する場合の環境変数と、動作しない場合の環境変数で違うやつ(特に PATH 等)を比較してみるとわかるかも。 -- 後藤 2005-08-02 (火) 13:08:31
  • うーん。残念ながら、それを確かめる方法がわかりません。 -- いしだ? 2005-08-02 (火) 15:17:21
  • ところで、同じR.appファイルとR.frameworkフォルダを使っているのに、なぜ、R.appが起動するアカウントと起動しないアカウントがあるのでしょうか。最悪、新しいアカウントに自分の環境を移行しても、同じようにRが起動しなくなるのではと不安です。 -- いしだ? 2005-08-02 (火) 15:21:14
  • 起動するアカウントとしないアカウントがあるんですか。立ち上がらないアカウントのホームディレクトリに .R で始まるファイルがいくつかありますか。もしかしたら,そられのファイルが壊れているとかいうことが関係するのかなぁ。いきなり消去は恐いので,どっかべつの所へ mv して,起動してみればいかがでしょう。。 -- 青木繁伸 2005-08-02 (火) 20:00:33
  • RController doLoadHistory? でひっかかってる...R(Mac GUI)側のバグのようですが、回避策としてはとりあえずターミナルで mv ~/.Rhistory ~/.Rhistory.bak でどうでしょう。 -- 岡田 2005-08-03 (水) 01:29:07
  • .Rで検索されるファイルは、すべてHD直下のライブラリ/Frameworks/R.Framework内にあるようでした。
    岡田さまのご指摘の方法で、起動しました。
    これで、解決ですね?皆さん、ありがとうございました。-- いしだ? 2005-08-03 (火) 9:41:14
  • すみません。まだ一部未解決でした。日本語が表示されません。
    ターミナルから
    export LANG=ja_JP.UTF-8
    export LC_ALL=ja_JP.UTF-8
    /Applications/R.app/Contents/MacOS/R
    としないかぎり、コンソールに日本語が表示されません。
    R.appのダブルクリックでは、コンソール内が英語表記です。-- いしだ? 2005-08-03 (火) 10:32:14
  • OS X自体は日本語環境(システム環境設定>言語環境で日本語が一番上)になっているのですよね。/Applications/R.app/Contents/Resources/ja.lproj に、9個の.nibファイルとLocalizable.strings がありますか?もしなければR.appが壊れているので再インストールして下さい。 -- 岡田 2005-08-03 (水) 10:43:12
  • OSX自体は、日本語環境で、Rもコンソール内以外は、日本語です。ja.lproj内には、おっしゃる通りファイルがあります。一応再インストールもしました。別アカウントでログインすると、コンソール内も日本誤表示されます。 -- いしだ? 2005-08-03 (水) 11:19:25
  • > .Rで検索されるファイルは、すべてHD直下のライブラリ/Frameworks/R.Framework内にあるようでした
    っていうのは,普通の状態ですかね?R関連のフォルダ,ファイル( .ファイルを含む)を全部消去してから,再インストールはいかが。 -- 青木繁伸 2005-08-03 (水) 12:56:35
  • ふむふむ, そうするとMac GUIの問題ではなく、R本体が日本語環境であることを認識していませんね。普通OS Xでは環境変数はいじりませんが、何か他のアプリなりツールによってFinder起動時の環境変数LANGがデフォルトから変更されてしまっている可能性が高いです。R.appをダブルクリックして起動したR内から Sys.getenv("LANG") とか Sys.getenv("LC_CTYPE") とか Sys.getenv("LC_ALL") を実行するとどうなりますか? -- 岡田 2005-08-03 (水) 12:59:04
  • ("LANG")でLANG (改行)"en_GB.UTF-8", ("LC_CTYPE")でLC_CTYPE(改行)"", ("LC_ALL")もLC_ALL(改行)""です。 -- いしだ? 2005-08-03 (水) 14:47:00
  • 以前、R関連のフォルダは、すべて捨てて再インストールしたことがありますが、改善しませんでした。改めて確かめましたが、.Rファイルは、HD直下のライブラリフォルダ何しかありません。 -- いしだ? 2005-08-03 (水) 14:54:09
  • 起動項目すべて削除し、再起動、最初にRを立ち上げてもコンソール内は英語のままです。コンパネも標準以外使ってません。 -- いしだ? 2005-08-03 (水) 15:06:44
  • おそらく解決しました。私は、言語環境では日本語を最上位にしていたのですが、日付の書式をカスタム設定にしていました(曜日の表記だけ英語にしたかったので)。何度か、日本とカスタムを切り替えて試したところ、Rコンソール内の言語は、それに連動していました。皆様、ほんとうにありがとうございました。 -- いしだ? 2005-08-03 (水) 15:39:50
  • en_GB!本当ですね。MacOS Xは、言語環境の「言語」メニューではなく、「書式」メニューでロケールを設定していたとは.....新しい発見でした。Rではないですが別のアプリの日本語化がうまくいかない原因がこれで判明しそうです。 -- 岡田 2005-08-03 (水) 23:41:40

要素に変数を持つ行列(関数)の扱いについて

生物系院生? (2005-07-27 (水) 13:00:57)

格子データの回帰分析に使用するパラメータ行列を作る際、function機能の中でfor文を使って変数を入れる位置を指定するという作業がうまくいかなかったので質問させて頂きました。
状況を取りうる最小サイズの格子データで説明すると

123
456
789

といったデータ(数字は変数名)の場合、例えば5のデータはfirst neighborのみを考慮するマルコフ確率場を仮定した場合2,8のデータと4,6のデータによって推測される。このときそれぞれが縦のデータの影響の度合いと横のデータにの影響の度合い関係するパラメータをもっており、それは以下のように表される。(データは周期条件を仮定して、要素内の数字の1が横要素、2が縦要素に関係するパラメータ、行の数字が予測対象の変数名、列の数字が説明変数の変数名に相当する行列)

    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,]    0    1    1    2    0    0    2    0    0
[2,]    1    0    1    0    2    0    0    2    0
[3,]    1    1    0    0    0    2    0    0    2
[4,]    2    0    0    0    1    1    2    0    0
[5,]    0    2    0    1    0    1    0    2    0
[6,]    0    0    2    1    1    0    0    0    2
[7,]    2    0    0    2    0    0    0    1    1
[8,]    0    2    0    0    2    0    1    0    1
[9,]    0    0    2    0    0    2    1    1    0

この行列自身は以下のコマンド

el <- 9  #格子データの総要素数 パラメータ行列では行数、列数に相当
col <- 3  #格子データの一行の要素数、すなわち列数
mat <- matrix(0, nrow = el, ncol = el)
xiA <- 1
xiB <- 2  #文字列操作では不都合が多いため、パラメータを数値で仮置き
for(i in 1:6){
    W <- i - 1
        if(ceiling(W/col) != ceiling(i/col)){W <- W + col}
    E <- i + 1
        if(ceiling(E/col) != ceiling(i/col)){E <- E - col}
    N <- i - col
        if(N < 1){N <- N + el}
    S <- i + col
    mat[i,W] <- xiA
    mat[i,E] <- xiA
    mat[i,N] <- xiB
    mat[i,S] <- xiB
}
for(i in 7:9){
    W <- i - 1
        if(ceiling(W/col) != ceiling(i/col)){W <- W + col}
    E <- i + 1
        if(ceiling(E/col) != ceiling(i/col)){E <- E - col}
    N <- i - col
        if(N < 1){N <- N + el}
    S <- i + col
        if(S > col){S <- S - el}
    mat[i,W] <- xiA
    mat[i,E] <- xiA
    mat[i,N] <- xiB
    mat[i,S] <- xiB
}

で作ることができるのですが、このコマンドにおいてxiA,xiBの部分を変数のままにした関数を作ってから1,2を代入しようとしてもうまくいかないのです(ちなみに実際のデータを使用する際には400×400の行列を使用する予定になっています)。そのうまくいっていないコマンドは以下のような感じで

el <- 9  #格子データの総要素数 パラメータ行列では行数、列数に相当
col <- 3  #格子データの一行の要素数、すなわち列数
matB <- function(a,xiA,xiB){
    for(i in 1:6){
        W <- i - 1  #W,Eは横方向に関するパラメータの位置
            if(ceiling(W/col) != ceiling(i/col)){W <- W + col}
        E <- i + 1
            if(ceiling(E/col) != ceiling(i/col)){E <- E - col}
        N <- i - col  #N,Sは縦方向に関するパラメータの位置
            if(N < 1){N <- N + el}
        S <- i + col
        a[i,W] <- xiA  #パラメータを代入する位置の指定
        a[i,E] <- xiA
        a[i,N] <- xiB
        a[i,S] <- xiB
    }
    for(i in 7:9){
        W <- i - 1
            if(ceiling(W/col) != ceiling(i/col)){W <- W + col}
        E <- i + 1
            if(ceiling(E/col) != ceiling(i/col)){E <- E - col}
        N <- i - col
            if(N < 1){N <- N + el}
        S <- i + col
            if(S > col){S <- S - el}
        a[i,W] <- xiA
        a[i,E] <- xiA
        a[i,N] <- xiB
        a[i,S] <- xiB
    }
}
base <- matrix(0, nrow = el, ncol = el)
matB(Base,1,2)

for文で要素の位置指定をすることができていないようでした(エラーメッセージは出ず、全角スペースが混じっていないことは検索済)。

手持ちの資料やネット上の解説では、関数定義で行列を扱っている例は少なく、他の言語の知識にも乏しいため行き詰って質問させていただきました。宜しくお願いします。

  • Note that any ordinary assignments done within the function are local and temporary and are lost after exit from the function. Thus the assignment X <- qr(X) does not affect the value of the argument in the calling program. です. 関数の末尾でreturn(a)として
    base<-matB(base,1,2)
    でよいのでは?forのインデクス云々というのは関係ないと思いますよ -- takahashi? 2005-07-27 (水) 13:23:08
  • returnの使いどころがちゃんと理解できていなかったので助かりました(結果は問題なく作動)、ありがとうございます。また、結果的に非常に単純な質問となってしまい失礼しました。 -- 生物系院生? 2005-07-27 (水) 15:48:26
  • 不審な繰り返しは無駄の徴候・バグの温床
    el <- 9  #格子データの総要素数 パラメータ行列では行数、列数に相当
    col <- 3  #格子データの一行の要素数、すなわち列数
    mat <- matrix(0, nrow = el, ncol = el)
    xiA <- 1
    xiB <- 2  #文字列操作では不都合が多いため、パラメータを数値で仮置き
    for (i in 1:9) {
        W <- i - 1
            if (ceiling(W/col) != ceiling(i/col)) {W <- W + col}
        E <- i + 1
            if (ceiling(E/col) != ceiling(i/col)) {E <- E - col}
        N <- i - col
            if (N < 1) {N <- N + el}
        S <- i + col
            if (i > el-col && S > col) {S <- S - el}
        mat[i,W] <- xiA
        mat[i,E] <- xiA
        mat[i,N] <- xiB
        mat[i,S] <- xiB
    }
    まだ無駄があるような。 -- 2005-07-27 (水) 16:55:39
  • 生物系院生が R の使用法で悩んだときは北大の久保先生に相談するのではなかったですか (^^) 恐らくギブスサンプリング等への利用を考えているような気がしますが、違いますか。各格子の添字 (i,j) に対し本質的に必要な情報は、上下左右の格子の座標とそれに付随する重み値だけですから、お尋ねの行列を作るのは冗長で無駄なのでは。参考例をあげておきます。(「変数を要素に持つ行列」といういい方は意味不明ですね) -- QDU? 2005-07-27 (水) 18:52:52
    > no <- 9
    > X <- matrix(1:81, c(no,no))   # テスト用格子データ
    > N <- W <- E <- S <- vector("list", no*no) # リスト行列の用意
    > dim(N) <- dim(W) <- dim(E) <- dim(S) <- c(no, no)
    > wN <- wW <- wE <- wS <- matrix(0, no, no) # 重み用行列
    >
    > for (i in 1:no)
    +    for (j in 1:no) {
    +       N[[i,j]] <- matrix(c(ifelse(i==1, no, i-1), j), c(1,2)) # 北(上)側格子の座標(周期条件付き)
    +       wN[i,j]  <- 2                                           # その重み
    +       S[[i,j]] <- matrix(c(ifelse(i==no, 1, i+1), j), c(1,2)) # 南(下)側格子の座標(周期条件付き)
    +       wS[i,j]  <- 2                                           # その重み
    +       E[[i,j]] <- matrix(c(i, ifelse(j==no, 1, j+1)), c(1,2)) # 西(右)側格子の座標(周期条件付き)
    +       wE[i,j]  <- 1                                           # その重み
    +       W[[i,j]] <- matrix(c(i, ifelse(j==1, no, j-1)), c(1,2)) # 東(左)側格子の座標(周期条件付き)
    +       wW[i,j]  <- 1
    +    }
    >
    > test1 <- function() {     # テスト用関数
    +             for (k in 1:1000)
    +                for (i in 1:no)
    +                   for (j in 1:no) { # (i,j) の上下左右の格子の値と重み
    +                      c(X[N[[i,j]]], X[E[[i,j]]], X[S[[i,j]]], X[W[[i,j]]])
    +                      c(wN[i,j], wE[i,j], wS[i,j], wW[i,j])
    +                   }
    +          }
    > gc(); system.time(test1()) # 実行時間
             used (Mb) gc trigger (Mb) max used (Mb)
    Ncells 168482  4.5     350000  9.4   350000  9.4
    Vcells  63009  0.5     786432  6.0   280690  2.2
    [1] 2.83 0.03 2.89 0.00 0.00
    
    ## 関数形で表現すると
    > N <- function(i,j) matrix(c(ifelse(i==1, no, i-1), j), c(1,2))
    > S <- function(i,j) matrix(c(ifelse(i==no, 1, i+1), j), c(1,2))
    > E <- function(i,j) matrix(c(i, ifelse(j==no, 1, j+1)), c(1,2))
    > W <- function(i,j) matrix(c(i, ifelse(j==1, no, j-1)), c(1,2))
    > wN <- wS <- function(i,j) 2
    > wE <- wW <- function(i,j) 1
    
    > test2 <- function() {
    +             for (k in 1:1000)
    +                for (i in 1:no)
    +                   for (j in 1:no) { # (i,j) の上下左右の格子の値と重み
    +                      c(X[N(i,j)], X[E(i,j)], X[S(i,j)], X[W(i,j)])
    +                      c(wN(i,j), wE(i,j), wS(i,j), wW(i,j))
    +                   }
    +          }
    > gc(); system.time(test2()) # 但し時間は45倍かかる
             used (Mb) gc trigger (Mb) max used (Mb)
    Ncells 167600  4.5     350000  9.4   350000  9.4
    Vcells  61585  0.5     786432  6.0   280690  2.2
    [1] 125.99   0.07 131.18   0.00   0.00
    

whileによる繰り返しについて

小山? (2005-07-25 (月) 11:04:15)

uniroot()を使用したく、Rでプログラムを書いています。元データとして約5000行,4列のデータがあります。その4列目のデータを1行目から1つずつ関数に代入させ、uniroot()を用いて関数の最適解を求め、それをまた別の関数に代入し、最終的な解を得たいのですが、繰り返し関数がうまく走りません。
プログラムとしては、まず4列目のデータの個数を調べ(x)、算出した解を格納する配列を定義し、xがデータの個数以下の場合は繰り返し配列に答えを格納していく様にしようと、whileを使用しました。
作ったプログラムは以下の様になっております。

no <- length(data[[4]])
ha <- array(data=NA, dim = c(nd, 7))
x <- 0
while (x <= nd){
   x <- x+1
   a  <- data[x, 4]
   b  <- 7
   c  <- 42
   d  <- 10
   e  <- 50
   f  <- 40
   obj <- function(phi,a)
     1/cos((phi*pi)/180)-cosh(a*tan((phi*pi)/180))
   res <- uniroot(obj,c(0.1,89.9),a=a)
   phi <- res$root
   func <- res$f.root
   depth <- rep(NA,b)
   i <- 1:b
   depth[i] <- d+c+L*(((1/tan(phi*pi/180))^2+1)^0.5
               -(1/(tan(phi*pi/180))^2+(1-(f+e*(i-1))/L)^2)^0.5)
   ha[x, ] <- c(depth, a)
 }

しかし配列は綺麗に出てくるのですが、データはすべてNAのままで、計算結果が格納されていません。Rの参考書やネット等で調べましたが、どの部分が間違っているのかわからず、質問させていただくことにしました。つたない説明で申し訳ありませんが、アドバイスよろしくお願いします。

  • インデンテーションは、tab キーで行うべきです。全角空白は使ってはいけません。追試できなくなります。今回は直しておきました。なお、nd が不定なのと(L もこのプログラムリスト中では定義されていない)、data もおおむねどんな感じのデータなのかわからないとこのプログラムを動かしてみることができないのです。
    たぶん,nd は,一行目の no なんでしょうね。
    関係ないところでは,
    x <- 0
    while (x <= nd) {
       x <- x+1
    というのはおかしいですね。x が nd のときに,新たなループが生じますが,その中で x は nd より1大きな値になりますよ。そして配列要素のアクセスに失敗するはずですが。
    また,while ループ中で,b,c,d,e,f,i という定数(ベクトル)を定義していますが,プログラム全体を通じて不変なものですからループの外で定義する方がよいでしょう。obj という関数の定義もループの外の方がいいでしょう(実際に複数回関数定義が行われるのかどうかは知らないが)。
    また,
       i <- 1:b
       depth[i] <- d+c+L*(((1/tan(phi*pi/180))^2+1)^0.5
                   -(1/(tan(phi*pi/180))^2+(1-(f+e*(i-1))/L)^2)^0.5)
    というところですが,左辺は depth だけでよいのではないでしょうか( [i] はいらないということ。またdepth <- rep(NA, b) も不要になる)
    ha[x, ] <- c(depth, a) としているが,ha の列数は 7 で,depth が要素数7ならaを加えて8個の要素をhaの1行に入れようとするとエラーになる。
    示されたプログラムのままだと気づきにくいかもしれないが,プログラムを整理していくと見通しがよくなって,矛盾に気づきやすくなると思いますよ。 -- 2005-07-25 (月) 11:30:57
  • コメントありがとうございます!申し訳ありません、ndはnoの誤りです。また、Lは簡潔に書こうとするあまりに省いてしまいました。
       L <- ((en-1)*e+2*f)/2
    ご指摘いただいた点を参考にして、以下の様に直してみました(長くなるので行間を詰めました)が、やはり出てくるのはNAだけでした。プログラムを実行ごすると 「ERROR: 構文解析エラーです」と出ます。質問なのですが、このプログラム実行後正確に計算されたかを見るには「ha]と入力し実行すれば良いのですよね?
       no <- length(data[[4]])
       ha <- array(data=NA, dim = c(no, 8))
       b  <- 7
       c  <- 42
       d  <- 10
       e  <- 50
       f  <- 40
       obj <- function(phi,a)
         1/cos((phi*pi)/180)-cosh(a*tan((phi*pi)/180))
       x <- 0
       while (x < no){
         x <- x+1
         a  <- data[x, 4]
         res <- uniroot(obj,c(0.1,89.9),a=a)
         phi <- res$root
         func <- res$f.root
         L <- ((en-1)*e+2*f)/2
         i <- 1:b
        depth[i] <- d+c+L*(((1/tan(phi*pi/180))^2+1)^0.5
                   -(1/(tan(phi*pi/180))^2+(1-(f+e*(i-1))/L)^2)^0.5)
         ha[x, ] <- c(depth, a)
        }
    また、データは以下の様なデータになっております。1列目から3列目までは時間データで、本プログラムでは使用しません。
       8	40	10	0.927 
       8	40	20	0.965 
       8	40	30	0.962 
       8	40	40	0.964 
       8	40	50	0.967 
    • 小山? 2005-07-25 (月) 12:38:11
  • えーと,
    depth[i] <- d+c+L*(((1/tan(phi*pi/180))^2+1)^0.5
    の前に全角空白が1つ残っていますね。左辺は depth だけ。
    L の定義の中の en も未定義です。
    no <- length(data[[4]])
    no <- nrow(data)
    で。 en を適当に定めたら,ha に,一応数値は入ったようですよ。 -- 2005-07-25 (月) 13:51:20
  • そうそう。付け加えとして,while を使うよりは for を使う方が構文的にははっきり・すっきりします。つまり,
    x <- 0
    while (x < no) {
         x <- x+1
    を,
    for (x in 1:no) {
    とするのです。そうすれば,x を1,2,...,no として順次くりかえすのだなとすぐ分かりますね。no は対象になるのかならないのかとか,余計な心配をしなくてもよい(エラーの種が減る)。
    なお,制御変数としては(数学のときと同じような慣例に従って),i,j,k,l あたりを使うのがよいでしょうね。 -- 2005-07-25 (月) 17:26:30
  • またしてもうっかりしていましたが、en=bでした。ご指摘の通り訂正したところ、haに値が入りました!ありがとうございます!しかし、ここでまたしても問題が発生してしまいまして、重ね重ね申し訳ありませんが質問させてください。aの値が1.0以上の場合に、unirootがおかしくなるようで、 「以下にエラーuniroot(obj, c(0.1, 89.9), a = a) : f() の端点での値が異なった符号を持ちません」と、出ます。 そこで、a > 1.0の場合にはa = 1.0(もしくは1.0では走らない場合0.99)にしようと、以下の様に付け足しました。しかし、また同じエラー文が出てしまいます。
    	for (j in 1:no){
    	a  <- data[j, 4]
    	if(a > 1.0){
    		a == 1.0
    		}else{
    		a == a
    		}
    	res <- uniroot(obj,c(0.1,89.9),a=a)
    原因はaの値なのでしょうか。a > 1.0の場合は次のようなデータになっています。
    	14	15	0	1.015 
    	14	15	10	1.018 
    	14	15	20	1.017 
    	14	15	30	1.014 
    	14	15	40	1.011 
    	14	15	50	1.013 
    • 小山? 2005-07-25 (月) 18:37:17
  • uniroot は二分法による求解なので,初期状態で二つの端点において目的関数が正・負の状態でないと解を求められません。a = 1.015 あたりのとき,0.1〜89.9 のobj関数の値を見てみたらいかがでしょう。縦軸にobj関数値,横軸に0.1〜89.9をとって,図を描くなりしてみましょう。横軸と交わるところがありますか?ないように思いますが。どうでしょう。 -- 2005-07-25 (月) 19:56:33
  • 以下のような例を考えましょう。関数は二次関数で,二つの解,x=0.1234, x= 0.5678 を持ちます。

    #ref(): File not found: "quad.png" at page "初級Q&A アーカイブ(3)"

    # 関数定義
    > quad <- function(x) (x-0.1234)*(x-0.5678)
    # 0,1 の範囲から解を求めようとしたとき
    > uniroot(quad, c(0, 1))
    以下にエラーuniroot(quad, c(0, 1)) : f() の端点での値が異なった符号を持ちません
    # つまり,quad(0)= 0.07006652, quad(1)=0.3788665
    # 異なった符号を持たない(両方とも正)ということ
    # 関数値は正,負,正と変化し区間内に二つの解を持つ
    
    > uniroot(quad, c(0, 0.1))
    以下にエラーuniroot(quad, c(0, 0.1)) : f() の端点での値が異なった符号を持ちません
    # つまり,quad(0)= 0.07006652, quad(0.1)=0.01094652
    # この場合は,関数値は,いつも正,つまり,この区間内に解はない!
    
    > uniroot(quad, c(0, 0.5))
    $root
    [1] 0.1233996
    
    $f.root
    [1] 1.843689e-07
    
    $iter
    [1] 7
    
    $estim.prec
    [1] 6.103516e-05
    
    # つまり,quad(0)= 0.07006652, quad(0.5)=-0.02553348
    # 異なった符号を持つということ
    # 関数値は正,負と変化し区間内に一つの解を持つ
    #(もしかしたら3個とか5個とか奇数個の解を持つのかもしれないが)
    よく吟味してみてください。 -- 2005-07-25 (月) 21:50:39
  • よく解りました。確かにその通りでした。条件文を書いてエラー値をはじき出したらプログラムはちゃんと走りました。エクセルのソルバーでやると気の遠くなるような時間がかかるであろう計算が数秒で出来、感激しております。アドバイスありがとうございました。 -- 小山? 2005-07-26 (火) 14:03:36
  • 一つ気になることが。pi/180 というのは円周率が弧度法と思い込みしているのでは? 今回は数秒でできたということで必要ないかも知れませんが、シミュレーション等時間のかかる作業をくり返し行う際は、途中でエラーが起きると泣くことに。try 関数はエラーが起きても、そこで実行中断せずに、次以降も実行してくれますから精神衞生に良いです。参考例をあげておきます。 -- QDU? 2005-07-26 (火) 15:15:30
    data <- rbind(c(8,40,10,0.927),
                  c(8,40,20,1.20 ),  # エラーになるようにした
                  c(8,40,30,0.962),
                  c(8,40,40,0.964),
                  c(8,40,50,1.20 ) ) # エラーになるようにした
    
    A <- data[,4]; no <- length(A)
    b  <- 7; c  <- 42; d  <- 10; e  <- 50; f  <- 40
    
    obj <- function(phi, a) 1/cos((phi*pi)/180)-cosh(a*tan((phi*pi)/180))
    
    doit <- function(x) {  # 一回の試行用の関数
              res <- uniroot(obj, c(0.1,89.9), a=A[x])
              phi <- res$root
              # func <- res$f.root  # 不要ですね
              L <- ((e-1)*e+2*f)/2
              y <- (1:b)-1
              depth <- d+c+L*(sqrt((1/tan(phi*pi/180))^2+1)
                        -sqrt(1/(tan(phi*pi/180))^2+(1-(f+e*y)/L)^2))
              return(c(depth, A[x]))
        }
    res <- vector("list", no) # 結果を納めるリストを予めサイズ指定で用意
    for(i in 1:no) res[[i]] <- try(doit(i), TRUE) # 全試行
     
    > res
    [[1]]
    [1]  75.40177 103.94505 131.67010 158.54344 184.53087 209.59754 233.70809
    [8]   0.92700 
    
    [[2]]
    [1] "以下にエラーuniroot(obj, c(0.1, 89.9), a = A[x]) : ?n?tf() の端点での値が異なった符号を持ちません?n"
    attr(,"class")
    [1] "try-error"
    
    [[3]]
    [1]  69.74911  91.28434 112.08025 132.11997 151.38677 169.86406 187.53553
    [8]   0.96200
    
    [[4]]
    [1]  69.32652  90.34202 110.62861 130.17047 148.95193 166.95752 184.17200
    [8]   0.96400
    
    [[5]]
    [1] "以下にエラーuniroot(obj, c(0.1, 89.9), a = A[x]) : ?n?tf() の端点での値が異なった符号を持ちません?n"
    attr(,"class")
    [1] "try-error"
  • ↑のエラーメッセージ中のエスケープシーケンス /n,/t が変ですね。 -- QDU? 2005-07-26 (火) 19:58:25
  • /n, /t ではなくバックスラッシュ n, バックスラッシュ t ですが
    表示関数がエスケープシークェンスを解釈していないようですね。 -- 2005-07-26 (火) 21:10:42

whichについて

田中? (2005-07-21 (木) 22:41:35)

n×n行列をすでに持っており、各行ごとに最大値を与える添え字を返させ、それをn×2の行列に格納してやりたい(つまりは、ある行について最大値をもつのはここですよというのを、全ての行について記録させたいということです。)のですが、いろいろな方法を試しても上手く出来ませんでした。
解決法を教えていただければと思います。よろしくお願いします。 一応、作ってみたのは以下のようなものです

Q <- matrix(0,1,n)
	for(b in 1:n){
		for(c in 1:n){
				Q[1,c] <- x[b,c]
				P[b] <- which(max(Q), arr.ind=TRUE)
				}
			}

xが調べたいn×n行列で、Pに格納するように作りました。

  • 最大値の位置でいいのでしょうか? -- なかま 2005-07-22 (金) 01:16:25
    col=4
    MM<-matrix(runif(col^2),ncol=col)
    SQM<-function(matrix){
            ml<-function(n){ncol=length(n);((1:ncol)[n==max(n)])[1]}
            rbind(apply(matrix,1,ml),apply(matrix,2,ml))
    }
    SQM(MM)
  • 質問中のプログラムが不完全なのは仕方ないですが,質問本文の意図を読み取ると以下のようなことかな?
    (x <- matrix(rnorm(12), 3, 4))
    which(x == apply(x, 1, max), arr.ind=TRUE)
    本質は二行目。
    なお,which の説明を読めば分かりますが,質問者の場合では一行中の最大値が同値の場合には,結果として返ってくるのは元の x の行数より多いですからご注意。
    > (x <- matrix(c(1,2,3,2,1,2,3,3,2,1,2,3),3,4))
         [,1] [,2] [,3] [,4]
    [1,]    1    2    3    1
    [2,]    2    1    3    2
    [3,]    3    2    2    3
    >  which(x == apply(x, 1, max), arr.ind=TRUE)
         row col
    [1,]   3   1
    [2,]   1   3
    [3,]   2   3
    [4,]   3   4
    ところで,arr.ind=TRUE のとき,結果の1列目でソートされない結果が返ってくるのは仕様かな。
    と思って,ソースを見たら。あら。結構泥臭いことやってるんだ。 -- 青木繁伸 2005-07-22 (金) 10:45:04
  • ご回答ありがとうございます。質問の意図はおっしゃるとおりです。
    x <<- mix2_16 #n×n行列
    P <- which(x == apply(x, 1, max), arr.ind=TRUE)
    P <<- P
    gn <- sprintf("result_mix%g_16.txt",a)
    write.table(P,file=gn, sep="?t", row.names=FALSE, quote=FALSE)
    print(P)
    のようなプログラムををつくり実行したところ、
          row col
    [1,]   4   4
    [2,]   5   4
    [3,]   8   4
    [4,]   9   4
    [5,]  13   4
    [6,]  15   4
    [7,]   7   7
    [8,]  19  25
    [9,]   1  29
    [10,]   2  29
    のような結果が返ってきました。(n=38のときで11以下省略)行ごとに最大値を与える添え字を出力させたつもりなので、rowが1,2,3・・・,10・・・と並ぶはずなのですがそうではありません。 なぜなのか分かりますでしょうか? よろしくお願いいたします。-- 田中? 2005-07-22 (金) 18:19:24
  • 一列目が順に並んでいないのは,私がコメントしたように,そのようにプログラムされているからです(^_^;)
    TRUE のある場所を列単位に探索しているからなんですね。which と入力すれば,ソースを見ることができますよ。
    並べ替えた結果が欲しければ,結果が P に入っているなら P[order(P[,1]),] とすればよいでしょう。
    ちなみに,<<- は使う必然性があるんですか? P <<- P は特に意味不明です。 -- 青木繁伸 2005-07-22 (金) 18:39:41
  • 早速の回答ありがとうございます。 P <<- P は先ほど書いたプログラムを自分で作成した関数内のプログラムとして組み込んでいるので、あとからPを取り出しやすくするために(いちいち、またそのプログラムを実行しないでよいように)P <<- Pとしてやっています^^; 教えて頂いた方法をさっそく試してみます。 -- 田中? 2005-07-22 (金) 22:33:47
  • 好みの問題かもしれませんが,解析が難しいバグの温床になる可能性もあるので,return(list(P=P,x=x)) 等として値を返す方が良いと思います。 -- 青木繁伸 2005-07-22 (金) 22:49:40
  • 日の下に新しいものはない(旧約聖書)のたとえの通りすでに max.col (列ごとの最大値の位置の添字を返す)が用意されています。但しタイがあるとランダムに一つが選ばれるようです。which -> which.min -> max.col とヘルプの See Also を注意するとたどり着けます。 -- QDU? 2005-07-23 (土) 00:27:36
> set.seed(31415); (x <- matrix(runif(20), c(4,5)))
          [,1]      [,2]      [,3]       [,4]      [,5]
[1,] 0.9502223 0.0607455 0.9365462 0.77713002 0.1302115
[2,] 0.3357378 0.3579473 0.8376269 0.08262812 0.9671198
[3,] 0.1330718 0.2271160 0.1275758 0.81371308 0.8880396
[4,] 0.4901114 0.4000740 0.7251304 0.30108535 0.5821669
> cbind(1:dim(x)[1], max.col(x)) # 列ごとの最大値位置の添字
     [,1] [,2]
[1,]    1    1
[2,]    2    5
[3,]    3    5
[4,]    4    3
> rbind(1:dim(x)[2], max.col(t(x))) # 行ごとの最大位置の添字
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    1    4    1    3    2
> cbind(1:dim(x)[1], max.col(-x)) # 最小位置が欲しければ -x とするだけ
     [,1] [,2]
[1,]    1    2
[2,]    2    4
[3,]    3    3
[4,]    4    4

lowess変換前と変換後のデータの対応方法

sugimoto? (2005-07-14 (木) 18:08:02)

初めて質問させていただきます。
R tipsを購入し、昨日から使い始めました。
今、lowess関数を使っています。
この関数を使って変換されたデータと、変換前のデータをつなぐにはどのような方法がありますでしょうか。もしくは、lowess関数に与えるx,yにIDのような値を持たせることは可能なのでしょうか。
大変つたない質問で申し訳ありませんが、よろしくお願いします。

  • 「つなぐ」というのは,変換前の値の隣の列に変換後の値をくっつけるということですか。cbind(cars[,2], lowess(cars)$y) などでよいと思いますが。 -- 青木繁伸 2005-07-14 (木) 18:58:28
  • 意味が分かりにくいですね。lowess の結果はあくまで平滑化で、そのものとしては、元の散布図中の各点と一対一に対応はしていませんから、変換後のこの点が変換前のこの点に対応と言うことは無いはずですね。 -- QDU? 2005-07-14 (木) 21:34:22
  • cars というデータセットに,lowess(cars) をして結果は $x, $y となり $x は cars[,1] になっていて,plot(lowess(cars))で平滑化後の曲線が描画できるので,一応は $y は $x のときの平滑後の予測値・推定値というべきものになっていると理解しましたが。間違っているんでしょうか? -- 青木繁伸 2005-07-14 (木) 23:49:26
  • すみません…。諸悪の根源は、The R-Tips で分かりにくい説明をしているせいだと思われます…。データハンドリングの章で、もっと分かりやすい説明をすべきでした。以下ではまず、変換前のデータと変換後のデータを、ID で紐付けしてデータをマージ(併合)しております。次に、マージ後に、不要な列を削除しております。(2つの操作は一度に出来るでしょうが、解説上、2つに分けております)-- 舟尾? 2005-07-15 (金) 01:00:01
# 変換前のデータをデータフレームにする(マージ用にIDを付加)
( before <- data.frame(ID=1:nrow(cars), cars) ) 
  ID speed dist
1   1     4    2
2   2     4   10
3   3     7    4
................
# 変換後のデータをデータフレームにする(マージ用にIDを付加)
( after  <- data.frame(ID=1:nrow(cars), lowess(cars)) ) 
   ID  x         y
1   1  4  4.965459
2   2  4  4.965459
3   3  7 13.124495
..................
# ID で紐付けしてデータをマージ
result <- merge(before, after, by="ID")
# x と dist は重複しているので x を削除
result$x <- NULL
# 結果を表示
result
   ID speed dist         y
1   1     4    2  4.965459
2   2     4   10  4.965459
3   3     7    4 13.124495
..........................
  • あと、青木先生のおっしゃることはご尤もでして、$y は $x のときの平滑後の予測値・推定値というべきものになっていると思います。sugimotoさんのジレンマは、「変換前のデータには存在したdistが、変換後のデータにない!」ということではないでしょうか??(間違っていたらすみません・・・) -- 舟尾? 2005-07-15 (金) 01:07:28
  • data.frame(ID=1:nrow(cars), cars, y=lowess(cars)$y) でよいということですね?(IDがいるかどうは別として) -- 青木繁伸 2005-07-15 (金) 07:29:32
  • 丁寧な解説、ありがとうございました。私が実行したかったのは、まさに舟尾さんが解説していただいた事そのものです。本当に大変参考になりました。ありがとうございました。 -- sugimoto? 2005-07-15 (金) 09:57:06
  • > data.frame(ID=1:nrow(cars), cars, y=lowess(cars)$y) でよいということですね?
    あ・・・・・。そうです!何を面倒なことをやっていたのでしょうか!!青木先生、ご指摘ありがとうございます。 -- 舟尾? 2005-07-15 (金) 14:06:27

ファイルの読み込みについて

田中? (2005-07-08 (金) 18:50:56)

以前、他の件で質問させて頂きました。
今回は、また別の事柄になるのですが、ファイル読み込みの際のファイル名の指定についての質問です。
具体的には、C:/sp01.txt,C:/sp02.txt,・・・C:/sp30.txtのようなファイルに対してそれぞれ同じプログラムを実行したいのですが、ファイル数が多くなることも考えて、

x01 <-read.table("C:/sp01.txt")

をfor文で回せないかと考えております。つまり、以下のようにしてfor文で変数?(ただ、分かりやすいように?にしただけです。)を変えてやり、ファイル名の変更を自動化したいわけです。(?の部分を01、02、03・・・と変更するのをfor文で回したい。また、program(a,b)は各ファイルに対し実行したいプログラムのことです。)

kurikaesi <- function(a,b) {
    for(? in 1:30) {
        x <- read.table("C:/sp?.txt")
        y? <- program(a,b)
     }
}

自分で調べては見たものの、結論を得ることが出来なかったので何か方法があればおしえてください。 よろしくお願いします。

  • x <- read.table("C:/sp?.txt") は,自分でファイル名を作ればいいだけでしょう
    for (i in 1:2) {
        fn <- sprintf("test%i.data", i))
        x <- read.table(fn)
    }
    のような感じ。
    y? <- program(a,b) は結果をy[i] に入れればいいと思うが,program() が返す形にもよる。質問をもう少し詳しく書いた方がよい。 -- 青木繁伸 2005-07-08 (金) 19:08:43
  • 早速のご回答ありがとうございます。 質問が分かりにくくてすみませんでした。質問の意図はご回答の通りで、 無事解決いたしました。一度、fnに渡して、それから代入する方法を思いつきませんでした。 ありがとうございました。-- 田中? 2005-07-08 (金) 20:25:20
  • 補注:もし test1.data でなく test01.data の様にしたければ -- QDU? 2005-07-09 (土) 08:26:19
    > i=1; sprintf("test%02i.data", i)  # 「整数二桁、一桁なら最初を0で埋める」の意味
    [1] "test01.data"
    > i=10; sprintf("test%02i.data", i)
    [1] "test10.data"
    こんなのもあり(R のループ範囲はベクトル、リストでもよい) -> 青木さんには怒られるかも :-)
    > for(file in lapply(1:20, function(i) sprintf("file%02i.data", i)))
      cat(file, "?n")
    file01.data 
    file02.data 
    ...........
    file19.data 
    file20.data
  • こういうのは好きなんですよ。書式にこだわらないなら,
    for (file in paste("file", 1:20, ".data", sep=""))
    のほうが短い。 -- 青木繁伸 2005-07-09 (土) 11:15:09
  • 様々なご意見ありがとうございます。非常に参考になりました。 -- 田中? 2005-07-11 (月) 23:44:01
  • 既にあるファイル名をリストに得れば事足りるらしいので,list.files("C://",pattern="sp[0-9]+.txt")でもいいのでは? -- 中澤? 2005-07-12 (火) 12:44:07

RODBCを使ったExcelの読み込み

Riemann? (2005-07-04 (月) 20:59:13)

R Bookに載っているExcelとの連携のところを読んで、「これは便利!」と思いやってみたのですが、1つのセルに長いテキストが入っていると途中できれてしまいます。ncharで測るとちょうど255文字になっています。どこかのオプションで制限サイズを変えればよさそうに思うのですが、どこをどういじったよいのかわかりません。ご教示ください。

  • うちでは300文字かつ300バイト以上でも問題なく読めますが?Windowsですよね? -- なかま 2005-07-05 (火) 00:25:15
  • はい、Windows XPで、Rのバージョンは1.91です。それから、Excelのセルのなかにあるのは日本語の2バイト文字です -- Riemann? 2005-07-05 (火) 01:01:56
  • 特に理由が無ければ2.1.1にバージョンアップしてください. -- なかま 2005-07-05 (火) 01:47:40
  • 2.1.1でやってみましたが、症状は同じでした。テキストが同じ箇所で切れてしまいます。 -- Riemann? 2005-07-05 (火) 14:46:49
  • Excel のバージョンは関係ないの?最新版ですか? -- 2005-07-05 (火) 15:48:55
  • 使っているのはExcel 2002 SP3でした。最新ではないですが、そんなに古いということもないと思います。 -- Riemann? 2005-07-05 (火) 16:57:16
  • こちらでは何ら問題は(改行も読めるし)無いんですが?もちろん多バイト文字も問題ありません. -- なかま 2005-07-05 (火) 17:11:33
  • Cドライブ直下にフォルダを作ってそこにデータを移動してみたら、途切れずに全部読めました。問題がおきた時の置き場所haには2バイト文字や空白が含まれているのですが、そのせいでしょうか?でも、途中で切れるのはちょっと不思議です。255という意味ありげな数字も??? -- Riemann? 2005-07-05 (火) 17:51:11
  • Cドライブ直下にフォルダを作ってそこにデータを移動してみたら、途切れずに全部読めました。問題がおきた時の置き場所はかなり深い場所で途中のフォルダ名には2バイト文字や空白が含まれているのですが、そのせいでしょうか?でも、途中で切れるのはちょっと不思議です。255という意味ありげな数字も??? -- Riemann? 2005-07-05 (火) 17:52:26
  • だぶってますね。パスに二バイト文字を含めないのは常識だと思います。いずれにせよ,二バイト文字を含めないようなパスにしたらうまくいくなら,今後すべからくそうすべきですね。 -- 2005-07-05 (火) 18:56:17
  • パスが2byteでも大丈夫で〜す.(T_T) 何が違うかもし判ったら教えてください. -- なかま 2005-07-05 (火) 19:31:36

デジタル信号処理

K? (2005-07-04 (月) 00:58:04)

RにはMatlabのSignal Processing Toolboxに相当するような
信号処理に適したパッケージはありますでしょうか?
fftやウェーブレット関係のパッケージはいくつかあるようですが、
ChebychevやButterWorth?など古典的なフィルタ群がないように思います。

「dsp」や「digital filter」、「butterworth」などのキーワードで
検索してみましたが、レスのついていないMLの投稿にヒットするようです。

既存の関数(filterなど)で十分実現可能だから必要ないということなのか
それとも誰も開発していない(需要がない)ということなのか・・。

どなたか非公式パッケージでもご存じでしたら是非教えてください。

  • Octave にはないのでしょうか? -- FlyMeToTheMoon? 2005-07-07 (木) 15:18:41

ベジェ曲線による補間

Krokodile? (2005-06-29 (水) 16:48:38)

 R では akima パッケージで、スプライン補間ができますが、
ベジェ曲線による補間はできるのでしょうか?

 このサイトを検索しましたが、「ベジェ曲線」はヒットしませんでした。

  • トップページからたどれる Rsite search でキーワード "Bezier" で検索すると Hmisc パッケージの bezier という関数がヒットしました。お試しあれ。なお Krokodile ではなく crocodile では。検索するときはローマ字は使えませんので、老婆心ですが一言。 -- 2005-06-29 (水) 19:37:41
  • ありがとうございます。"Krokodile"はドイツ語です。英語では"crocodiles"です。個人で複数形を使うのは問題でしたが。 -- Krokodil? 2005-06-29 (水) 20:19:13
  • ↑ 失礼しました -- 2005-06-29 (水) 22:01:35

行列について

田中? (2005-06-27 (月) 17:30:29)

以下の質問のご回答よろしくお願いいたします。
k×(m*i)行列Aを作成するときに、1×mの行列Bkm(k、mは添え字で、B行列の要素はk,m,iの関数として与えられる)を合成することで作成したい(つまり、行列Aの要素として行列Bを利用したい)のですが、

centroid <- function(i,K,m) {
    C <- matrix(0,K,m*i)
    U <- matrix(0,1,m)

    for (a in 1:K) {
        for (b in 1:i) {
            for (w in 1:m) {
                U[a,b][1,w] <- (a,b,wの3変数による関数)
            }
        }
    }
    for (a in 1:K) {
        for (b in 1:i*m) {
            t <- ceiling(b/m)
            C[a,b] <- U[a,t][1,(b-(t-1)*m)]
        }
    }
}

としたのですが、エラーがでます。Uの定義の仕方に問題があるのであろうとは思うのですが、調べても(行列に2変数の添え字を付け、さらにその行列の要素を先の2変数+1変数で指定してやる方法について)分からなかったので、助けていただければと思います。
(分かりにくい気がするので、まとめると、C[i,j]=Uij Uij[1,k]=f(i,j,k)という場合の、Uの記述方法について教えていただきたいということです)
よろしくお願いいたします。

  • 行列でなければいけないのでしょうか?一番分かりやすい方法は3次元配列を使うことだと思うのですが。 次のようにすれば x[i,j,k] はご質問の C の (i,j) ブロック行列 U_{ij} の (1,k) 成分になります。もう一つの方法は C を U 行列からなるリストとして表現し、
    C[[i,j]][1,k] 
    のようにアクセスすることでしょうか。-- MKR? 2005-06-27 (月) 18:35:22
    x <- array(0,c(3,3,3))
    for (i in 1:3) for (j in 1:3) for (k in 1:3) x[i,j,k] <- f(i,j,k)
    x[i,j,k]
  • C のポインタみたいな機能が欲しいということでしょうかねえ。
    例えば PTR というベクトルがあって,その要素は別の行列を指している。
    PTR【n】 <- matrix(0, n, n)
     【】で表される要素はメモリ領域の先頭を指すポインタを表す(^_^;)
    PTR【2】 は 2×2行列を意味し,PTR【10】は 10×10行列。。
    こういうのがあったら便利だなと思うこともありますね(あるよ!という答えが返ってきたりして)。 -- 青木繁伸 2005-06-27 (月) 19:13:23
  • 早速のご返答ありがとうございます。
    私の目的としては、C行列をk×i個の部分行列Uijで考えることで、あとでUijに対していくつかの計算を施して、新たにk×(m*i)のD行列を作成するということを考えています。そのため、C行列の中に部分行列を考える方針にしました。(具体的に言うと、C行列の中身と言うのはm個のデータ列←U (これで1つのまとまり)のk×i観測分データ(さらにはkは要素、iは次元)となっており、のちに、k×i個のデータ列に対してそれぞれ同様な操作を行い1×mの行列Tとし、これを合成してk×(m*i)のD行列を作りたいと思っています。)この目的に適した方法があれば教えていただければと思います。言葉足らずですみませんでした。よろしくお願い致します。 -- 田中? 2005-06-27 (月) 22:16:24
  • 特に高速性を要求しなければ、一番柔軟な方法はリストの行列を考えることのような気がします。こうすればベクトルでも、行列でも、任意の R オブジェクトを行列風にアクセスで来ますし、後から変更も自由自在です。 -- MKR? 2005-06-28 (火) 01:13:00
    C <- as.list(rep(NA, 3*3)) 
    dim(C) <- c(3,3)                # リストを行列化 
    for (i in 1:3) for (j in 1:3)   # リスト行列の各成分に(なんでも)ほうりこむ
       C[[i,j]] <- runif(i*j) 
    C    
              [,1]      [,2]      [,3] 
    [1,] 0.6133985 Numeric,2 Numeric,3 
    [2,] Numeric,2 Numeric,4 Numeric,6 
    [3,] Numeric,3 Numeric,6 Numeric,9 
    C[[2,3]][4]   # (2,3) 番目のブロック(ベクトル)の4番目の成分 
    [1] 0.7488104
  • 青木さん、MKRさん、ありがとうございました。何とか問題は解決に近づきそうです^^ ただ、1つつまづいていることがあります。それはリストCの宣言の仕方です。
    centroid <- function(i,K,m){
    
    C <- ????
    for(a in 1:K){
    	for(b in 1:i){
    		for(w in 1:m){
    			C[[a,b]][w] <- (a,b,wの関数)
    				}
    			}
    		}
    	
    print(t(C))
    
    }
    とプログラムを作りました。K×i個のベクトル(要素数m)からK×iのリスト行列(?)を出力したいのですが、Cの宣言の仕方がよく分かりませんでした。(いろいろ調べて試してみたのですが、エラーばかりでてしまいました。)これについて教えていただけませんでしょうか。 たびたびになり申し訳ないのですが、よろしくお願いいたします。-- 田中? 2005-06-29 (水) 01:19:36
  • C <- lapply(rep(m, K*i), numeric); dim(C) <- c(K,i) とするか、先の例のように C <- as.list(rep(NA, K*i)); dim(C) <- c(K,i) とし、各成分にサイズ m の「ベクトル」を付値します。後者の場合、最初まだ各成分はベクトルになっていませんから、
    C[[a,b]][w] 
    という構文は使えません。前者なら使えるはずです。ついでに始めて気づきましたが、リストは行列(風)にアクセスできるだけでなく、dim(C) <- c(3,3,3) 等とすれば配列(風)にアクセスすることもできるようです。なお、リストへのアクセスは純粋な行列・配列へのアクセスに比べれば遅めですから、特に高速性が要求されるときは問題有りかもしれません。-- MKR? 2005-06-29 (水) 06:51:01
  • 回答ありがとうございます。とりあえず、問題は解決しました。 あまりにも、計算に時間がかかるようでしたらまた別の方法も模索してみようと思います。 回答して頂けた方、ありがとうございました。 -- 田中? 2005-06-29 (水) 17:59:35

for文を使わず行列の要素毎に関数を適用する方法

13m? (2005-06-26 (日) 06:26:05)

ある行列もしくはベクトルの要素それぞれに対して、1から任意の整数aまでのx/nの総和を計算するf(x,a)=Σ[n=1:a](x/n)のような関数を適用したいと思い、 以下のようなものを書いてみました。

f <- function(x, a)
{
    res <- 0
    for(n in 1:a) res <- res + x/n
    return( res )
}

実際に使ってみた様子はこんな感じです

> m=matrix(c(1,2,3,4),ncol=2)
> m
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> f(m,10)
         [,1]      [,2]
[1,] 2.928968  8.786905
[2,] 5.857937 11.715873

これと同じ結果を返す関数を、for文(または要素毎のapply)を使わず、もっと簡潔に、R的に、スマートに書けますでしょうか?すぐにでも書けそうな気がするのですが、私の力不足のせいかなかなかいい方法が思いつきません。なんだかお題みたいな質問になってしまってすみません。

  • 御質問を文字通りに理解すると単に f(x,a)=x*(1/1+1/2+...+1/a) を計算することになると思うのですが、本当にそうなのですか。 -- MKR? 2005-06-26 (日) 10:04:58
temp <- function(x,a) x*sum(1/(1:a))
temp(m,10)
         [,1]      [,2]
[1,] 2.928968  8.786905
[2,] 5.857937 11.715873
  • あ!ごめんなさい。ついつい問題を単純化しすぎました。私が知りたいのはxが外にくくりだせないような場合です・・いまからこの下に新たに定義した問題を付け加えます。申し訳ないです。ご指摘感謝します。
    ○「改訂版」
    f(x,a)=Σ[n=1:a](1/(x+n))のような関数を、行列またはベクトルの要素毎に適用する方法を知りたい。
    そこで以下のようなものを書いてみた。
    g <- function(x, a)
    {
        res <- 0
        for(n in 1:a) res <- res + 1/(x+n)
        return( res )
    }
    実行結果
    > m=matrix(1:4,ncol=2)
    > g(m,5)
             [,1]      [,2]
    [1,] 1.450000 0.8845238
    [2,] 1.092857 0.7456349
    これと同じ結果を返すものを、いかにもfor文を使わずに書けそうな気がするので、もしあれば教えていただきたいと思います。 -- 13m? 2005-06-26 (日) 13:13:57
  • 以下のプログラムは R らしいでしょうか。注釈なしで一目でわかるでしょうか。
    > g <- function(x, n)
    + {
    + 	ncol <- ncol(x)
    + 	nrow <- nrow(x)
    + 	k <- array(rep(1:n, each=ncol*nrow), dim=c(nrow, ncol, n))
    + 	x <- array(x, dim=c(nrow, ncol, n))
    + 	apply(1/(x+k), c(1,2), sum)
    + }
    > m <- matrix(1:4, 2)
    > g(m, 5)
             [,1]      [,2]
    [1,] 1.450000 0.8845238
    [2,] 1.092857 0.7456349
    我ながらわかりにくいと思います。for を使ったプログラムの方が遙かに分かりやすい。 -- 青木繁伸 2005-06-26 (日) 20:20:33
  • 元の記事を書き直してしまったのはMKRさんに失礼でした。お詫びします。 -- 13m? 2005-06-26 (日) 21:40:35
  • 青木さん>丁寧な回答いただきありがとうございます。私は、「いかにもありそう」な気がしていたので、もしかすると案外もっと素直な書き方やそれ専用の関数があるのでは…と考えていただけです。特段、アクロバティックなコードを求めていたわけではなかったので、私が例示した書き方でもうすでに十分だということをご指摘いただいても私にとっては全く同様でありがたいことです。 -- 13m? 2005-06-26 (日) 22:04:19
  • 青木さん>とは言え、アクロバティックなコードでも、もしあれば是非一度見てみたいという気持ちも(大いに)あります…。rep(1:n, each=ncol*nrow)の使い方はとても参考になりました。私も、array()で3次元方向にa個の行列を生成して、最後に和をとるという方法でやろうとしたのですが、どうしても要素毎(c(1,2)のような)applyやfor文が必要になりそうなので断念してました。 -- 13m? 2005-06-26 (日) 22:20:21
  • 記事を追加せずに,投稿済みの記事の途中に挿入・改変するから,話の流れがおかしくなる。そのようにするのなら,私の発言も幾つか削って,あなたの発言の中から,私の削除分についても削除させて頂く。話の流れはねじ曲がってはいないと思うので,あしからず。
    要するに,コメントの付いた記事は誤字の訂正程度以外はいじるのは良くないと言うことですよね。 -- 青木繁伸 2005-06-26 (日) 23:31:16
  • 良いコードかどうかは別にして、例えば -- MKR? 2005-06-27 (月) 00:06:11
    > matrix(colSums(1/mapply(seq, m+1, m+5)), dim(m))
             [,1]      [,2]
    [1,] 1.450000 0.8845238
    [2,] 1.092857 0.7456349
  • おおお、一行コード。mapplyとseqの組み合わせで要素ごとの数列が縦に作れるんですね。おもしろいなあ…。mapplyや匿名関数使えば色々凝った一行コードができそうですね。 -- 13m? 2005-06-27 (月) 01:45:41
  • MKRさんの一行コードを参考に(パクって)sapplyバージョンで作ってみました。
    > h <- function(x,n) matrix(colSums(sapply( x, function(x) 1/( x+seq(n) ))), dim(x))
    > h(m,5)
             [,1]      [,2]
    [1,] 1.450000 0.8845238
    [2,] 1.092857 0.7456349
    上のやり方だとかなりアレなので、きちんと別途関数定義して、x,nに関するある程度どんな数列でも扱えるよう工夫してみました。ついでに引数はベクトルにも行列にも対応。
    eess <- function(FUN,x,n) 
    { # 引数xのそれぞれの要素に対して一般項FUNに関する1からn項までの総和を求める
       res <- colSums(sapply(x, FUN, seq(n)))
       if(is.matrix(x)) res <- matrix( res, dim(x) )
       return( res )
    }
    
    feess <- function(FUN,x,n)
    { # 上記と全く同様の機能をfor文で書き直したもの
       FUN <- match.fun(FUN)
       res <- 0
       for(i in seq(n)) res <- res + FUN(x,i)
       return(res)
    }
    使用例とsystem.time()
    > gn <- function(x,k) 1/(x+k) # 一般項を適当な関数で定義
    > eess(gn,m,5)
             [,1]      [,2]
    [1,] 1.450000 0.8845238
    [2,] 1.092857 0.7456349
    > eess(gn,1:5,5)
    [1] 1.4500000 1.0928571 0.8845238 0.7456349 0.6456349
    > eess("^",m,5) # 二項演算子を渡して奇抜な使い方も(Σ[k=1:n]x^k)
         [,1] [,2]
    [1,]    5  363
    [2,]   62 1364
    > system.time(eess(gn,rnorm(10000),100))
    [1] 0.51 0.04 0.59   NA   NA
    > system.time(feess(gn,rnorm(10000),100))
    [1] 0.06 0.02 0.08   NA   NA
    sapply()の演算は思ったより高速なんですね。驚きました。速度的にはそれほどでなくとも(なんだかんだで最初のfor文が一番速い)、いつか何か別の用途に使えるかもしれません。可読性があるかと言われればあまりないかもしれませんが。 -- 13m? 2005-06-27 (月) 02:25:07

画像の貼り付け

坂田? (2005-06-24 (金) 22:05:29)

主成分得点の平面に固体を名前で表示するよりも顔写真を貼り付けたい
と思います(顔に限りませんけど)。このための方法をご教授ください。
巷の本ではプロ野球選手の成績に対して主成分分析を行っているものが
見受けられますが、名前ではなく顔を張りつければインパクトも大きいと思
うのですが。

  • 適当な画像ファイルに落して、あとは適当な画像編集ソフトを使い手作業で重ねるしか方法はないのでは。 -- 2005-06-25 (土) 00:51:38
  • 早速のレス有難うございます。 むー 出来ればRの世界で閉じていてほしいのですが、そういうのって想定外なのかもですね。  -- 坂田? 2005-06-25 (土) 13:28:02
  • 好みの問題ですが,画像をアスキーアートに変換して小さなフォントでテキスト出力するとか.:-) -- なかま 2005-06-25 (土) 20:00:06
  • 回答に出遅れてしまった(^o^; 任意の位置に画像を貼る方法はあります。このwikiのどこかに以前書き込んだ気がしますが、pixmapのaddlogoなんていかがでしょう。詳しくはaddlogoのヘルプを見てください。-- 谷村 2005-06-27 (月) 15:05:28
    > library(pixmap)
    > r <- read.pnm(system.file("pictures/logo.ppm", package = "pixmap")[1])
    > x <- runif(10)
    > y <- runif(10)
    > plot(x,y)
    > for ( i in 1:length(x)) addlogo(r,px=c(x[i]-0.05,x[i]+0.05),py=c(y[i]-0.05,y[i]+0.05))

    #ref(): File not found: "addlogo.png" at page "初級Q&A アーカイブ(3)"

  • なかま様、谷村様 ご教授有難うございました。Rでも可能なことが分かりました。ここまで表現できればすばらしいですね。主成分分析がずいぶん楽しくなりました。 -- 坂田? 2005-06-28 (火) 02:29:56

plot.pointって???

? (2005-06-22 (水) 21:22:30)

はじめまして。
あるデータの平均と、そのなかの特定の1つのデータの図を同じ図にplotするのにpointを使用するときいたのですが、どうすればいいのでしょうか???

  • そう言った人に聞くのが一番だと思いませんか。
    points 関数のことではないかと思いますが?
    そうであるとしても,やりたいことをちゃんと書かないと,他人には分かりませんよ。
    > x <- rnorm(100) # 100個の正規乱数を発生させて
    > plot(1:2, c(min(x), max(x)), type="n") # プロット領域を定義して
    > points(1.5, mean(x)) # 平均値を描いて
    > points(1.5, x[3], col="red") # 3番目のデータを赤でプロットして
    なんて馬鹿馬鹿しいことをやるわけじゃないんでしょう?
    投稿法もちゃんと確認してから投稿しましょう。 -- 青木繁伸 2005-06-22 (水) 22:35:11

ダミー変数の表示について

? (2005-06-22 (水) 18:35:56)

重回帰分析をする際にダミー変数を用いました。
しかし、実測値と推定値をplotする際にどうしてもエラーになっていしまいます。どのようにして表現するのでしょうか?
ちなみにダミー変数は数字からアルファベットに変えてやりました。

  • どのようにプログラムしたのか,再現できるプログラムと必要なら最小限のデータを提示すべきでしょう。
    実測値と推定値のプロットにダミー変数の扱いが原因のエラーが出るとは思えませんよ。
    ハンドルネームはなんでも良いとはいうものの,「?」というのはおかしいでしょう。
    > x1 <- factor(c("foo", "bar", "bar", "foo", "baz", "baz"))
    > x2 <- factor(c("hoge", "hoge", "wao", "wao", "hoge", "hoge"))
    > y <- c( 2, 1, 2, 4, 3, 4)
    > ans <- lm(y ~ x1+x2)
    > cbind(y, predict(ans))
      y     
    1 2 2.25
    2 1 0.75
    3 2 2.25
    4 4 3.75
    5 3 3.50
    6 4 3.50
    > plot(y, predict(ans))

    #ref(): File not found: "regplot.png" at page "初級Q&A アーカイブ(3)"

    というようなことになったが? -- 青木繁伸 2005-06-22 (水) 20:15:21
  • すみません。名前まで考えてなかったもので…
    データはCSVで読み込んで、ある従属変数を5つの説明変数で重回帰しようとして、その前の段階で6つのダミー変数により全体を説明しようとしました。説明とPCが下手ですみません。 -- question? 2005-06-22 (水) 21:12:31
  • で,結局,何がやりたくて,どんなプログラムを書いたら,どんなエラーが出たというのでしょうか。回答が不要なら,質問は無用です。 -- 青木繁伸 2005-06-22 (水) 22:47:56

RMySQLのインストールについて

森野? (2005-06-17 (金) 19:40:38)

初歩的な質問で,申し訳ありません。
CRANからRMySQLをダウンロードして,RGuiのパッケージから(ローカルにあるZIPファイルからのパッケージのインストール)を実行しました。
library()で確認したところ,

** No title available (pre-2.0.0 install?) **

と表示されてしまい,使用する事ができません。
WindowsXPサービスパック1,R 2.1.0を使用しているのですが,どなたかご教授頂けないでしょうか。

  • RODBCを使いましょう. -- なかま 2005-06-17 (金) 20:32:25
  • RODBCをダウンロードして頑張ってみます. -- 森野? 2005-06-17 (金) 21:33:03
  • ご教授,本当に有難うございました. -- 森野? 2005-06-17 (金) 21:35:50
  • 投稿法も、良く読んでくださいね。最初の投稿と今では、見た目が全然違うでしょ?どんなシステムでも、ちゃんと使うための方法を学んでから使うべきなんですね。ちゃんと読んで使わないから、ちゃんと使えないわけです(あたりまえですね)。みんな同じ原理なんですね。そのあたりを理解しないと、新しいことやろうとするたんびに、同じような(根元的に同じ)トラブルに遭遇するんです(トラブルに遭遇しないとすれば、それは偶然の幸運なんですね) -- 2005-06-17 (金) 22:28:25
  • 以後気をつけます.ご迷惑お掛けしました. -- 森野? 2005-06-18 (土) 00:45:20
  • RODBCインストール成功しました!有難うございます. -- 森野? 2005-06-20 (月) 20:08:02

VimでRを使用したい

山下? (2005-06-15 (水) 17:52:32)

WindowsのVimでRを使用したいのですが、
皆様使用されている方はいませんでしょうか?

エディタ版でVimファイルは見つけたのですが、Syntaxも上手くいかず、
デフォルトの :set Syntax=rを使用しています。

xyzzy並みに使っておられる方がいましたらご教授下さい。

  • 私は Vim でスクリプトを編集し,Emacs(ESS) 上の R で source() で読み込ませる,という使い方をしています (OSはLinuxです).自動でsyntaxを判別しないようならば,ファイルの末尾に # vim: set syn=r: と付け加えてみてはいかがでしょう? -- ivan? 2005-06-22 (水) 22:50:51
  • ご返答ありがとうございました。ご提示頂いた方法を使用させて頂きます。 -- 山下? 2005-06-24 (金) 19:59:18
  • 途中で切れてしまったので、続きを記載いたします。やはりVim上からRを実行させるというようなことは出来ないのでしょうか?どなたそのような環境を作った!!という方がおられましたらご教授下さい。Vi使いは以外に多いと思いますので有用な情報かと・・・私が調べた限りはわからなかったので・・・どうぞ宜しく御願いします。 -- 山下? 2005-06-24 (金) 20:01:46

ATLASを使いたい

どうむ (2005-06-15 (水) 12:13:12)

LinuxでRを使用しております。
ATLASを使用したいのですが、どのようにしたらいいのでしょうか、、、
以下「うまくいかない」現象です

  • おこなったこと(ATLAS)
    • http://math-atlas.sourceforge.net/からソースをダウンロード、解凍
    • 一般ユーザーで
      make
    • rootユーザーで
      make install arch=Linux_HAMMER64SSE2_2
      • CPUは、に上からもわかるとおり、AMD64bitです(汗
    • 作成された「.a」をパスの通っているディレクトリへ、コピー
  • おこなったこと(R)
    • バージョン > R-2.0.1
    • ソースをダウンロード、解凍
    • 一般ユーザーで
      ./configure --with-blas=f77blas
      • この時点でログには
        checking for sgemm_ in -lf77blas... no
        checking for sgemm_... no
        checking for ATL_xerbla in -latlas... no
        checking for sgemm_ in -lblas... yes
      • f77blasを使用することにはなっていません、、、
    • 一般ユーザーで
      make
      make check
    • rootユーザーで
      make install
      • 案の定lddで実行ファイルを見ても普通のlibblas.soがリンクされている

これを、解決するには

  • ATLASがデフォルトで作成するlibf77blas.aというスタティックライブラリをRのmakeのときにリンクさせるか
  • libf77blas.soというシェアードライブラリを作成する
    ということになると思うのですが、どうなのでしょうか?_?
    blas or atlas で、当サイトを検索して、調べてみましたが、みなさん「.so」ファイルは当たり前のようにできていたりするのでしょうか
    ご教授よろしくお願いいたします
    (次はgotoにチャレンジだ!)

p.s.
Win版ATLASはdllをダウンロードしてきて、既存のものと置き換えるだけでしたので簡単でした(汗
行列演算が、計算の中身にもよりますが、4分の1になったりして、驚いています
Linuxでもこのパフォーマンスを体験したい、、、

  • パスの通ったディレクトリとありますが(ライブラリの検索はパスとは無関係です),野良ビルドモノなので/usr/local/libあたりが適切かと思います(このとき必要なのはlibatlas.a,libf77blas.aおまけでlibclbas.aです) ./configure --with-blas='-L/usr/local/lib -lf77blas -latlas' (-L:ここのライブラリ見てね, f77blas:ATLASのFortranシンボル窓口, atlas:そう、本体が無い!)としてください. 動的共有ライブラリを作るなら,Makefile内にatlas.soの作り方が書いてあります. 動的共有ライブラリを作る場合はビルド時に-fPICが必要(デフォルトはそうなってたかも)かもしれません. 動的共有ライブラリを使う場合は,/etc/ld.so.confに格納場所を追加して/sbin/ldconfig するか, Rビルド時に -rpath を加えるかするといいかと思います. -- なかま 2005-06-15 (水) 13:34:05
  • ご回答ありがとうございます。現在頂いた情報を元に、いろいろ試しています。(Rのmakeがなぜか、うまくいきません、、、涙) -- どうむ 2005-06-15 (水) 20:26:29
  • ATLAS関係以外なら,たぶんヘッダー関係(hoge-devel or hoge-dev)が無いのでは? -- なかま 2005-06-15 (水) 22:41:58
  • R の configure スクリプトは単に sgemm_ をリンクできるかどうかをチェックしているだけなので、libf77blas.a に sgemm_ が含まれているかを調べた方がいいです。商用コンパイラで生成した ATLAS の場合、sgemm ("_" がないのがミソ) となる場合があります。チェック方法は nm libf77blas.a | grep sgemm_ で T のシンボルが出れば大丈夫です。 -- 後藤 2005-06-16 (木) 11:44:47
  • ついでに config.log の該当部分を書き写してくれると一発で幸せになれます(それとも R のバージョンが古いのが原因かな?2.1.0 を使ってみましょう)。 -- 後藤 2005-06-16 (木) 12:03:23
  • なかま様・後藤様、情報ありがとうございます!結果から申し上げますと何とか(若干の問題も残しながら)できました! -- どうむ 2005-06-16 (木) 15:43:16
    • 環境
      • AMD64では、64ビット系?のところでうまくいかないようでしたので、オーソドックスなCPU?に変更しています
      • Xeon3.2GHz×2
      • メモリ2Gbyte
      • Red Hat Enterprise Linux WS release 3 (Taroon) Kernel 2.4.21-4.ELsmp on an i686
    • ATLASのコンパイル(CPU一個しか使わない設定で、、)
      make install arch=Linux_UNKNOWNSSE2
      • できた「libatlas.a」「libf77blas.a」を「/usr/lib」へコピー
    • ./configure --with-blas="-L/usr/lib -lf77blas -latlas"
      make
      make check
      make install
      • コンフィグのログ該当部分?
        checking for sgemm_ in -L/usr/lib -lf77blas -latlas... yes
      • 静的リンクに成功!!
    • 結果
      • test2関数
        test2 <- function ( n=500 )
        {
        	A<-array(rnorm(n^2), dim=c(n,n))
        	B<-array(rnorm(n^2), dim=c(n,n))
        	C<-array(rnorm(n^2), dim=c(n,n))
        	D<-array(rnorm(n^2), dim=c(n,n))
        	BA <- B%*%A 
        	return (system.time(A%*%solve(t(BA)%*%BA+C)%*%BA%*%D ))
        }
      • 導入前
        > test2()
        [1] 2.24 0.08 2.33 0.00 0.00
        > test2()
        [1] 2.26 0.06 2.33 0.00 0.00
        > test2()
        [1] 2.27 0.05 2.32 0.00 0.00
      • 導入後
        > test2()
        [1] 0.46 0.07 0.54 0.00 0.00
        > test2()
        [1] 0.45 0.06 0.51 0.00 0.00
        > test2()
        [1] 0.43 0.05 0.48 0.00 0.00
      • やった!!
      • が、、、当マシン(Xeon×2)では、ATLASのコンパイル時にマシンに負荷をかけないと、途中でエラーが出てコンパイルできません、、、たまたまなのでしょうか???
      • 3回失敗、負荷(CPU使用率50%-ひとつのCPUが100%の状態)をかけてコンパイル2回成功、、、(いえ、一回目はたまたま負荷がかかっている状態だっただけなのですが、、、)
      • ↓これが途中でエラーが出てとまります。どういうエラーだったかは、確認中です(涙
        make install arch=Linux_UNKNOWNSSE2
      • AMD64では、こんなことはなかったのですが、、、
  • HAMMER64SSE2_2の2を見落としていました.Rのビルド時に --disable-R-profiling --with-blas='-L/usr/local/lib -lptf77blas -latlas -lpthread' を付与してみてください. -- なかま 2005-06-16 (木) 16:19:33
    # ATLAS
    make xconfig
    ./xconfig -F f '-fomit-frame-pointer -O -m64 -fPIC' ?
              -F c '-fomit-frame-pointer -O -mfpmath=387 -m64 -fPIC' ?
              -F m '-fomit-frame-pointer -O -mfpmath=387 -m64 -fPIC'
    make install arch=Linux_HAMMER64SSE2_2 # rootである必要はありません
    
    # R (ATLASのスレッドとR-profilingは不仲!)
    ./configure --disable-R-profiling --with-blas='-L/usr/local/lib -lptf77blas -latlas -lpthread'
  • AMD64も上記のとおりすると、できました。ありがとうございました! -- どうむ 2005-06-16 (木) 19:24:06
    • 導入前は、大体4-5秒かかっていたのが
      > test2()
      [1] 0.76 0.05 0.57 0.00 0.00
      > test2()
      [1] 0.75 0.07 0.58 0.00 0.00
      > test2()
      [1] 0.75 0.06 0.56 0.00 0.00
    • すばらしい!
  • XEONなら3.7.10(developer)あたりが良いかとおもいます. 普通にxconfigした後,-fPICを加えるか上のように自分で指定してもかまいません. -- なかま 2005-06-16 (木) 21:05:15
  • 再々、貴重な情報ありがとうございます。当方業務にてRを使っておりまして、RのパフォーマンスUPは社内外に影響が大きいのです(もちろん良いほうに)。 -- どうむ 2005-06-17 (金) 08:31:10
  • ちょっと大変ですが,出来たらATLAS+Rのビルドページ作りませんか.>どうむさん. 私もそんなに沢山のCPUは持っていませんので.(笑) -- なかま 2005-06-17 (金) 11:35:37
  • どのような構成にすればよいのか、まったく見当もつきませんが、サマリーページ、チャレンジしてみます、、、 -- どうむ 2005-06-17 (金) 13:34:14

データをソートしたい

学生? (2005-06-14 (火) 18:45:06)

はじめまして、とっても初歩的な質問ですが、
R環境で、入力データ(数値)をソートしたいです、
できれば、出力は順番づけられてほしい
ご存知の方、教えていただければありがたいです。
よろしくお願いいたします

  • 「出力は順番づけられてほしい」というのは,どういうことなんでしょう。対面して対話しているのでないので,ある程度きちんと書かないと通じないです。
    sort 関数というのはご存じなんでしょうか。 -- 青木繁伸 2005-06-14 (火) 19:09:29
  • 早速返事して頂きありがとうございました!今までRを使ったことがないから、sortという関数は知りません...その関数の使い方を教えていただけませんでしょうか、後、順番づけの件ですが、私の方が勘違いしましたので、なかった事にしといて頂ければと思います!よろしくお願いいたします! -- 学生? 2005-06-15 (水) 11:23:25
  • やりたいことをやってくれる関数名だけが分かっていて,詳しい使い方が分からないときには,help という関数を使うことができます。
    help(sort) # または簡略化された ? sortでもよい
    今後は自分で調べてもらうことにして,sort の使い方は以下の通り。
    > x <- c(3,6,5,2,8)
    > sort(x)
    [1] 2 3 5 6 8
    Net で得られる R の使い方についての資料は,舟尾さんの R-tipsが優れていますが,本を買ってあげてください。
    書名(著者・編者等) The R Tips データ解析環境Rの基本技・グラフィック活用集(舟尾暢男著)
    発行元        株式会社 九天社
    発行年        2005
    価格         3500円(税別)
    ISBN         4-86167-039-x
    です。 -- 青木繁伸 2005-06-15 (水) 11:27:17
  • ○○したことがないから○○出来ない、というフレーズは感心しません。世の中の問題、解法は自分から編み出さなきゃ。みんなが他力本願では日本はつぶれるぜ。sortを知らなければ、x[order(x)]やfor2重ループを使えばよろしい。要は、自分の知っている範囲をまず使いこなすことです。その次にあちこち探してみる。それでへたばったら初めて泣きを入れることです。forループすら知らなければ、まずプログラミングとアルゴリズムの勉強から始めるべきです。Rは言語です。ワープロみたいなアプリではありません。それなりの訓練が必要です。 -- 2005-06-15 (水) 14:53:53
  • 丁寧な回答をどうもありがとうございました!後は、自分で手を動かしてみたいと思います -- 学生? 2005-06-16 (木) 18:27:22

heatmapの縦横比変更 及び 複数の図を描画する方法について

小島? (2005-06-14 (火) 17:40:44)

はじめまして。
heatmap関数を使用して作図をしているのですが、行き詰ってしまいました。
以下の事柄についてご存知の方がいましたら解決方法を教え下さい。

質問1:heatmapの縦横比の変更について

heatmap関数を用いて解析を行うと、その縦横比(X軸とY軸の長さ)が必ず1:1の正方形になるのですが、この比率は変更できないのでしょうか?
扱っているデータがX軸の数の方が圧倒的に多いので、X軸を長くした長方形の図を得たいと思っています。
通常、出力デバイス(png, pdf等)の大きさを変更するとそれに合わせて図も伸縮してくれるのですが、heatmap関数の場合この比率が変更されずに余白が出来てしまいます。
例を挙げると以下のようになります。

#サンプルデータ(データをX軸の方が多くなるように加工)
data(mtcars)
sampledata<-rbind(c(mtcars[,2]), c(mtcars[,3]), c(mtcars[,4]), c(mtcars[,5]))
#デバイスの大きさを横長の長方形に指定する
pdf("test.pdf", width=12, height=5)
heatmap(sampledata)
dev.off()

marginsをいじれば形は変化するのですが、余白が出来る上、xlabの値が残ってしまいます。

heatmap(sampledata, xlab="XLAB", margins=c(20,5))

質問2:heatmapを使った複数の図の描画について

複数の図を描画する時、

#layoutを用いた画面分割(下図を大きく)
mat <- matrix(c(1,2,2), 3, 1, byrow=TRUE)
layout(mat)
plot(sin)
plot(cos)

など、plot関数を用いた場合は描画できるのですが、
heatmap関数を用いると、heatmapの図だけ独立して新規に出来てしまいます。

mat <- matrix(c(1,2,2), 3, 1, byrow = TRUE)
layout(mat)
#データは上記と同じものを使用
plot(sampledata)
heatmap(sampledata)

何か回避する方法はあるのでしょうか?

ヘルプを見たり検索をしてみたのですが、私自身の力不足とheatmapに関する記述が全体的に少なかったため方法が見つかりませんでした。
ご存知の方がいらっしゃいましたらよろしくお願いいたします。

  • まず質問1について
    heatmap を例えば heatmap2 に付値して,heatmap2 の中にある
        layout(lmat, widths = lwid, heights = lhei, respect = TRUE)
    という行の,respect を FALSE にします。つまり,
        layout(lmat, widths = lwid, heights = lhei, respect = FALSE)
    とすれば,いいようですが。そして当然使うのは heatmap2(sampledata) のように。以下のようになります(デバイスの大きさを示すため,余分な背景も切り取りました)。

    #ref(): File not found: "heatmap.png" at page "初級Q&A アーカイブ(3)"

    自己責任で試してみてください。
    質問2は,あ,自明かな。heatmap が上のように layout を呼ぶからでしょう。
    望むとおりのグラフを描くには,相当の手入れが必要かと。R でグラフを合成するより,個々のグラフを別のソフト(たとえば Word ?)に取り込んで合成(レイアウト)する方が簡単でしょう。 -- 青木繁伸 2005-06-14 (火) 19:56:08
  • 早速の返信ありがとうございました。
    こちらの返信が遅くなってしまいすみません。heatmap.2のインストールに思いのほか時間がかかってしまいました・・・。
    #上にはheatmap2()となっていますが、gregmiscのheatmap.2()でいいのですか?

    再度質問なのですが、heatmap.2にしたところ、時間軸(X軸)を固定していた部分がエラーで弾かれてしまいました。
    # heatmapでは時間軸の固定が可能
    heatmap(sampledata, Colv=NA)
    
    # heatmap.2ではColvに欠損値(NA)は入れられないとエラーが出る
    heatmap.2(sampledata, Colv=NA, trace="none", key=FALSE)
    以下のように dendrogram を指定すると似たようなグラフは出てくるのですが順序がバラバラになってしまいます。
    heatmap.2(sampledata, dendrogram="row", trace="none", key=FALSE)
    ヘルプの Note: 欄に書いてあった時間軸固定に関する記述がなくなっているので、heatmap.2 ではこのような設定はできないのでしょうか?-- 小島? 2005-06-15 (水) 15:35:47
  • ?? gregmiscのheatmap.2()でいいのですか?
    全然よくありません。gregmisc にheatmap.2 という関数があるんですか。。。。
    私が上に書いたように,あくまでも heatmap を使うのだが,関数の一部を書き換えないといけない。
    直接書き換えるのはお勧めできないので,適当な名前に付値して(だから heatmap2 なんてのでなくてどんな名前でもいいのだが),その付値された名前で定義される関数の書き換えをしましょう。と。
    書き換える場所は,これこれですよと。
    このような説明が必要と言うことは,自己責任でその過程を踏むことができなさそうなので,今回は,あきらめた方がよろしいのではないかと思います。
    教えられればなんでもできるというものでもありませんから,しかたないです。 -- 青木繁伸 2005-06-15 (水) 19:58:51
  • と,言ってしまうのも何ですから。
    heatmap2 <- heatmap
    fix(heatmap2)
     ここで,前記の通りの行を修正後,保存し編集ウインドウを閉じる
    heatmap2(sampledata)
    これで,うまく行くと思うのだが。。。 -- 青木繁伸 2005-06-15 (水) 20:35:42
  • 丁寧な回答をありがとうございます。
    もともとheatmap.2という関数名を知っていたので勘違いをしてしまい、お手数をおかけしました。
    heatmap.2()でも、
    rowInd <- 1:nr
    colInd <- 1:nc
    という記述を適当な場所に加えればいいようです。
    もっとRを使いこなせるように今後も勉強していきたいと思います。
    青木先生、ありがとうございました。-- 小島? 2005-06-16 (木) 11:12:35

重回帰分析でのカテゴリー変数:内部表現の値

ショーゾー? (2005-06-06 (月) 13:20:11)

一つのカテゴリー変数(FとM)と5つの数値変数を説明変数としてyの予測式を重回帰分析(ステップワイズ)で求めようとしています。性別のカテゴリー変数(sex)はfactorとして読み込み,str()でsex : Factor w/ 2 levels "F","M": 1 1 1 1 2 2 ... と確認できました。重回帰分析のsummary()ではこのsexの係数(estimate)が,sexM 6.3007という表示(有意でした)がされました。
ここで,予測値を計算する時はsex=Mのときは6.3007 * 2としてもちいるのでしょうか? sex=Fのときは6.3007 * 1でいいのでしょうか?
チェックのために,Fのときは0,Mのときは1のダミー変数(sexcode)をつくって実行したところ,おなじ結果がでました。このときの予測値の計算の時はsex=Mのときは6.3007 * 1,sex=Fのときは6.3007 * 0となるように思います。すると,factorの内部表現はどうなっているのでしょうか?

  • factor は,アルファベット順で数値を割り当てられるのだったと思います。すなわち,M,F なら(常識に反して)M=2,F=1だと思いますが。
    > x <- c("A", "B", "C")
    > as.integer(as.factor(x))
    [1] 1 2 3
    > x <- c("C", "B", "A")
    > as.integer(as.factor(x))
    [1] 3 2 1
    という結果を精査すべし。 -- 青木繁伸 2005-06-06 (月) 13:37:25
  • 関係するので書いておきます(どこかに書いた気がするが。。。)
    複数のレベルを持つ factor を lm で使うと,ちゃんとダミー変数を作って分析してくれる。
    > x <- as.factor(sample(3, 10, replace=T))
    > x
     [1] 1 3 2 1 1 3 1 3 1 1
    Levels: 1 2 3
    > y <- rnorm(10)
    > lm(y ~ x)
    
    Call:
    lm(formula = y ~ x)
    
    Coefficients:
    (Intercept)           x2           x3        # x2,x3 は何かといえば,,,
       -0.04008      0.80732     -0.18625  
    
    > x2 <- ifelse(x ==2, 1, 0) # こういう風に解釈してくれている
    > x3 <- ifelse(x ==3, 1, 0)
    > lm(y ~ x2+x3)
    
    Call:
    lm(formula = y ~ x2 + x3)
    
    Coefficients:
    (Intercept)           x2           x3  
       -0.04008      0.80732     -0.18625  # で,同じ答えになる
    こんな風になる。賢いなぁ R は。 -- 青木繁伸 2005-06-07 (火) 00:32:00
  • 青木さんありがとうございます。factorをlmで使った時に青木さんのコメントで,x2 <- ifelse(x == 2, 1, 0)のようにRが解釈してくれることが分かりました。factorを使った時と,代りに1,0を使ったときの結果が一致することから,factorの内部表現の数字ではなくて論理値(でいいのかな?)の0, 1が使われているような感じをもっていました。ヘルプでlm, formula, termを当たってみたのですが,このような解説はみつけることができませんでした。-- ショーゾー? 2005-06-08 (水) 13:50:01
  • カテゴリー変数をfactor化するときは、factor()にlevelsオプションをつけないとアルファベット順の順序化がされてしまいます。hist()などでもこの順序が使われるので、factor化するときはできるだけlevelsオプションをつけておいたほうがよいです。-- 2005-06-08 (水) 15:56:01
    > c <- c("Low", "Medium", "High")
    > f <- factor(c)
    > f
    [1] Low    Medium High  
    Levels: High Low Medium
    > as.integer(f)
    [1] 2 3 1
    > f2 <- factor(c, levels=c("Low", "Medium", "High"))
    > f2
    [1] Low    Medium High  
    Levels: Low Medium High
    > as.integer(f2)
    [1] 1 2 3
    > 
  • 青木さん「どこかに書いた気がするが...」ありました。なんでも掲示板の2005-04-18の記事です。検索時のキーワードが不適切で見逃しておりましたが,今回みつけました。 -- ショーゾー? 2005-06-09 (木) 12:51:17

R commanderを使いたい

g? (2005-06-06 (月) 10:48:48)

R commanderを使いたくて,パッケージRcmdrを読み込んだら
以下のメッセージが出てしまいました.
同様な状況を解決した人がいましたら解決方法を教えてください.
ちなみにRのバージョンはR2.1.0で
R commandeのバージョンはRcmdr1.0-2です.
また,以下のパッケージも読み込んでいます.
(バージョンもチェック済み)
abind, car (>= 1.0-15), effects (>= 1.0-7), foreign,
grid, lattice, lmtest, MASS, mgcv, multcomp,
mvtnorm, nlme, nnet, relimp, sandwich, strucchange, zoo

*****エラーメッセージの始まり*****

local({pkg <- select.list(sort(.packages(all.available = TRUE)))

  1. if(nchar(pkg)) library(pkg, character.only=TRUE)})
    要求されたパッケージ tcltk をロード中です
    要求されたパッケージ rgl をロード中です
    要求されたパッケージ zoo をロード中です
    要求されたパッケージ strucchange をロード中です
    要求されたパッケージ sandwich をロード中です
    要求されたパッケージ relimp をロード中です
    要求されたパッケージ nnet をロード中です
    要求されたパッケージ nlme をロード中です
    要求されたパッケージ mvtnorm をロード中です
    要求されたパッケージ multcomp をロード中です
    要求されたパッケージ mgcv をロード中です
    This is mgcv 1.3-0
    要求されたパッケージ MASS をロード中です
    要求されたパッケージ lmtest をロード中です
    要求されたパッケージ lattice をロード中です
    要求されたパッケージ grid をロード中です
    要求されたパッケージ foreign をロード中です
    要求されたパッケージ effects をロード中です
    要求されたパッケージ car をロード中です
    要求されたパッケージ abind をロード中です
    以下にエラーparse(file, n, text, prompt) : 構文解析エラーです
    エラー:.onAttach は 'attachNamespace' で失敗しました
    エラー:'Rcmdr' に対するパッケージもしくは名前空間のロードが失敗しました
  • それは、なんと言うOS上でのお話しでしょうか? -- なかま 2005-06-06 (月) 14:05:06
  • 失礼しました.WinXP(pro)のsp2です.ちなみにこのエラーは,R commanderのフレームとメニューが表示され,ワークスペースが表示されません.また,Rcmdr09-11では,エラーメッセージは表示されませんが,フレームしか表示されません(メニューの表示なし).Rの再インストールとWinの再起動も試みましたが同様の結果でした. -- g? 2005-06-06 (月) 14:33:38
  • Win2000ですが、R2.1.0およびR2.1.0 patched(6月4日)ではこの現象が出ていません。以前のpatchedバージョンで、同じ現象があった記憶があります。 -- a? 2005-06-06 (月) 14:38:15
  • Win版ならlocaleToCharset?の影響かも.日替わりR(パッチ適用版)ではどうでしょうか.(日本語不要ならLC_ALL=Cとかでもいいかも) -- なかま 2005-06-06 (月) 14:42:39
  • 質問当初は,Rterm上でも同様のエラーが起こったのですが,JGRのインストールの後Rterm上で確かめたら, -- g? 2005-06-06 (月) 19:06:32
  • (続きです.)R commanderが起動しました.使用目的からすると,これで十分です.(GUIの方は,相変わらずダメです.JGRが解決してくれたのかどうかも定かではありません.)ご回答してくださった皆様,ありがとうございました. -- g? 2005-06-06 (月) 19:09:30
  • あれからさらにいろいろ試してみたところ,インストール時に作成されるショートカットでGUIを起動するとエラーが起き,Rgui.exeそのもの,またはRgui.exeから作ったショートカットから起動するとエラーは起こりませんでした.偶然なのかバグなのかわかりませんがとりあえず報告まで. -- g? 2005-06-07 (火) 10:53:19

散布図のプロットについて

ken? (2005-06-03 (金) 23:03:02)

はじめまして
plot関数を使って散布図を作成しているのですが、x軸、y軸に対応した項目を直接図の中に表示するにはどうすればよいでしょうか? 例えば、

項目 X Y
A    3  4
B    7  5  # 投稿法を良く読んで投稿しなさい

のような表の散布図を作成した場合、A、Bを直接図の中に表したいと思っています。一通り検索してみましたが、いい方法が見つからなかったので、質問させていただきました。 よろしくお願いします。

  • こんなのはどうでしょうか?
    > dat <- matrix(c(3,4,7,5), nrow = 2, byrow=TRUE, dimnames = list(c("A", "B"), c("x", "y")))
    > plot(dat, pch=rownames(dat))
    松村俊和 (2005-06-04 (土) 00:33:00)
  • 早速の回答、ありがとうございます。 この方法を探していました。 プログラミング自体が初めてなので、このような場があると本当に助かります。ありがとうございました。 -- ken? 2005-06-04 (土) 14:33:46
  • 上記のプログラムでは、A,B等、1文字だけは出力されますが、2文字になるとグラフに表現されなくなります。 -- ken? 2005-06-04 (土) 16:50:39
  • plot(dat, pch=" ");text(dat, rownames(dat)) -- 2005-06-04 (土) 18:21:37
  • データファイルがちゃんと作られていれば,データを読んだ後で行に名前を付与する必要もなく,以下のようにするだけでよいでしょう。
    d <- read.table("test.dat")
    plot(d, pch=" ")
    text(d, rownames(d))
    プロット位置にマークを書いたり,text で書かれる文字列がマークと重ならないようにするにはどうしたらよいかは,オンラインヘルプを参照すればよいでしょう。 -- 青木繁伸 2005-06-04 (土) 21:33:35
  • 無事、出来上がりました。 松村様、青木様、いろいろとご指導いただきありがとうございました。 -- ken? 2005-06-05 (日) 00:26:02

標本平均のSD

Mari? (2005-06-02 (木) 13:14:43)

以下のように乱数を発生させ、その平均値の分布のSDがσ/sqrt(n) になることを確かめたかったのですが、なぜかσ/sqrt(n)にはなりません・・。
途中で計算のさせかたを間違っているのかと思いましたが、どうしても分かりません。

> dat <- numeric(10000)
> for ( i in 1:10000) dat[i] <- mean(rnorm(100,1,3))
> sd(dat)
[1] 0.30414
> 3/sqrt(10000)
[1] 0.03

どなたかご教示いただけませんでしょうか・・・

  • 標本の大きさは100でしょう?( mean(rnorm(100,1,3)) って書いてあるもの)だから,平均値の標準誤差は 3/sqrt(10000) ではなく,3/sqrt(100)=0.3 で,あってますね。
    > x <- matrix(rnorm(1000000,1,3), 10000,100)
    > sd(apply(x, 1, mean))
    [1] 0.29930224896266783
    ですね。 -- 青木繁伸 2005-06-02 (木) 13:34:23
  • ありがとうございました。nを勘違いしておりました。もっと勉強します! -- Mari? 2005-06-02 (木) 13:45:07

ヒストグラムの階級の区切りと重なった値の扱い

<ふ>? (2005-06-02 (木) 00:04:40)

R2.0.1を使ってます。

x <- c(150,160,170,180,190)

というデータをつくり、

> hist(x)

でグラフをかかせますと、150と160が同じ区間にカウントされます。
help(hist)を読みましたら、
'right = TRUE' (default), では、'(a, b]'で値をカウントするけど、左端の区間は、別、と書いてありました。
それなら、と、

> hist(x,xlim=c(140,200))

で描画させてみたのですが、左端とは、データの左端なんですね。
次に、

> hist(x,right = FALSE)

としますと、180と190が同じ区間にカウントされます。
手でグラフを描いたら、こんなふうにはせずに、平らな図にします。このように、最小値なり最大値の扱いを他と変える理由はどこにあるのでしょうか。
よろしくご教示ください。

> stem(x)
 The decimal point is 1 digit(s) to the right of the |

  15 | 0
  16 | 0
  17 | 0
  18 | 0
  19 | 0

> hist(x,xlim=c(140,200))
> hist(x,xlim=c(140,200),right = FALSE)
  • breaks を使って
    hist(x,breaks=15:20*10,right=FALSE)
    とすれば良いのでは?
    help にちゃんと書いてありますよ。
    Note that xlim is not used to define the histogram (breaks),
    but only for plotting (when plot = TRUE).
    しかしまぁ,エラーまでペーストする必要はないでしょう。
    先頭1文字空白行の使い方もまずかったし。-- 青木繁伸 2005-06-02 (木) 00:08:57
  • ありがとうございます -- <ふ>? 2005-06-03 (金) 00:18:18
  • 青木先生、すみません、ここの使い方、まだ、よくわかってません。先のエラー部分の送信も含めて、ご迷惑おかけしました。breaksの使い方、ありがとうございます。やってみます。 -- <ふ>? 2005-06-03 (金) 00:20:30

kruskal.test()の使い方間違っているのでしょうか?

RIRISU? (2005-05-31 (火) 01:52:02)

クルスカル・ワリスの検定についてお聞きします。

   1回群  2回群  4回群  8回群
no.1 9    13  19   24
no.2 11   15  19   23
no.3 11   14  21   20
no.4 14   16  22   24
no.5 20   17  19   19
no.6 13   21  20   24

このデータから群によって評定の差があるかどうかを検定したいのですが、

> x<-c(9,11,11,14,20,13)
> y<-c(13,15,14,16,17,21)
> z<-c(19,19,21,22,19,20)
> e<-c(24,23,20,24,19,24)
> x<-c(x,y,z,e)
> g<-factor(rep(1:4,c(6,6,6,6)),
+ labels=c("1group","2group","4group","8group"))
> kruskal.test(x,g)

で合っているでしょうか?手計算(表計算ソフトで公式に数値を当てはめるやり方)でやると、カイ二乗値は1494.129となり、自由度は24-1で23になります。
手計算でやると、0.1%の危険率で帰無仮説を棄却するのですが、うえの実行結果ではカイ二乗値=14.4297で自由度=3となり、1%の危険率で帰無仮説を棄却する。となり結果がちがってしまいます。

  • そもそも,自由度についてさえ,なぜ23になるのか。計算が間違えている,あるいは計算方法を誤解しているだけでしょう。R での結果は正しいと思います。
    なお,データがどのように取られているか書かれていないので,本当にクラスカル・ウォリス検定が適切なのかどうかは保留しておきます。
    投稿記事は,(コメントが付く前なら,訂正・編集して良いと思いますよ。今回は私が校正しておきます(^_^;)) -- 青木繁伸 2005-05-31 (火) 08:33:52

post-hoc Fisher's LSD test

Sybock? (2005-05-26 (木) 18:50:12)

はじめまして、質問させていただきます。

anova を3群のデータに適用し、有意差が出たものについて post-hoc analysis を行いたいと思っております。
中澤港様の「Rによる統計解析の基礎」P.109 を参考に、

# 例
firingRatio <- c(52,30,35,41,29,…, 33)
types <- c("A", "C", "A", "B", …, "C")
# 両ベクトルとも同じ length
pairwise.t.test(firingRatio, types, p.adjust.method = "holm")

といった形で実行しましたが、この pairwise.t.test という関数の"pairwise"は、t.test(paired = TRUE)とは意味が違うものなのでしょうか。

また、古典的な Fisher の制約つき LSD 法 も行って見たいのですが、RjpWiki 内で Fisher, LSD 等で単語検索を行ったところ、それらしいものを発見できず、yahoo 等で検索しても R における実現方法を探り当てることができませんでした。
知識不足のため、検索キーワードが不適切であった可能性もあり、申し訳ありませんが、上記の件につきましてアドバイスをよろしくお願い致します。

  • pairwise.*.test は,K 群の全ての二群の組み合わせについて*.test を行う(対比較を行う)という意味です。 -- 青木繁伸 2005-05-26 (木) 19:10:22
  • google でキーワード "filetype:R LSD" で検索すると12例見つかります。R site search でも探してみて下さい。-- 2005-05-26 (木) 19:44:55
  • 青木先生>ご回答いただき、ありがとうございます。pairの意味が全く違ったのですね。 -- Sybock? 2005-05-27 (金) 11:53:05
  • Fisher's LSD test についても検索して勉強してみます。 -- Sybock? 2005-05-27 (金) 12:06:42

日付値の取り扱い

さかな? (2005-05-25 (水) 21:35:34)

質問させて頂きます。

例えば、開始日(2002/05/03)と終了日(2004/03/21)を指定した時に、
期間中の一連の日付を生成する良い方法はないでしょうか?

2002/05/02
2002/05/03
....
2004/03/20
2004/03/21

アドバイスよろしくお願いします。

  • help(seq)とか,seq(as.POSIXct("2002/05/03"),as.POSIXct("2004/03/21"),by = "day")とか,日付のクラスはちとややこしいですね. -- なかま 2005-05-25 (水) 21:58:16
  • 月/日/年の順で良ければ以下が一番簡単か?
    library(chron)
    seq.dates("05/03/2002", "03/21/2004")
    。。。-- 青木繁伸 2005-05-25 (水) 23:12:54
  • base パッケージだけで済ますなら(ちなみに始点と終点に日付オブジェクトを与えると seq 関数はメソッド関数 seq.POSIXt を使います) -- 2005-05-25 (水) 23:20:23
    x <- seq(as.POSIXct("2002/05/03"),as.POSIXct("2004/03/21"),by = "day")
    gsub("-","/",as.character(x))
  • gsub 使うなら, format(x,"%Y/%m/%d") の方が素直かも.%Aとかもちゃんと(OSに依存するけど)国際化されています.:-) -- なかま 2005-05-25 (水) 23:46:14
  • seqも総称関数だったんですね。アドバイスありがとうございました。 -- さかな? 2005-05-26 (木) 21:12:37

表示の問題

こじろう? (2005-05-12 (木) 01:21:32)

Rのカーソルが|ではなく■の形になっていてとても使いにくいです。|の形にしたいのですが、どなたかご教示いただけませんか??

  • これは R の問題ではなく、OS の問題ではと思います。だとすれば、ここで質問する内容ではありません。そうでなくても、使っている OS がなんであるぐらい書くものですよ。 -- 2005-05-12 (木) 06:37:07
  • コメントありがとうございます。そして、お気に触る質問の仕方、失敬しました。Windows MEです。OSの問題ならあきらめるしかないですね。 -- こじろう? 2005-05-12 (木) 14:02:59
  • カーソルはたしか, src/gnuwin32/console.cあたりで描画しています. CTRL+Oで挿入/上書きの切替が出来ますから,それに合わせてカーソルの描画を変えれると,うれしいかもしれませんね. -- なかま 2005-05-12 (木) 18:21:30
  • アドバイス、ありがとうございました。でも、やはり無理でした。残念です。 -- こじろう? 2005-05-17 (火) 21:54:05

表と図を一枚のPDFに出力したいです

Akira? (2005-05-04 (水) 0:01:12)

LaTeXを勉強しろと言われそうですが、グラフと表を一枚のPDFファイルで出力する方法はあるのでしょうか?
LaTeXRによるポストスクリプト画像のLaTeXでの利用mimetexを確認しました。理解不足もありますが、少し違うように思います。
非効率的とはわかりながらも、今は座標を探しながらmtextで表を作成しています。

pdf(file="test.pdf")
layout(c(1:2))
plot(x=0, y=0, xlim=range(0,4), ylim=range(0,5), pch=&quot;.", main="test", xaxs="i", 
               yaxs="i", axes=FALSE, xlab="", ylab="", col="white")
for(i in seq(0, 6)){
 segments(0.5, i, 14.5, i)
}
mtext(text="column1", side=3, line = -1, adj=0, at=2, font=2)
mtext(text="column2", side=3, line = -1, adj=0, at=3, font=2)
mtext(text="row1", side=3, line = -2.2, adj=0, at=1, font=3)
mtext(text="row2", side=3, line = -3.4, adj=0, at=1, font=3)
mtext(text="row3", side=3, line = -4.5, adj=0, at=1, font=3)
mtext(text="row4", side=3, line = -5.7, adj=0, at=1, font=3)
mtext(text="data1.1", side=3, line = -2.2, adj=0, at=2, font=1)
mtext(text="data2.1", side=3, line = -3.4, adj=0, at=2, font=1)
mtext(text="data3.1", side=3, line = -4.5, adj=0, at=2, font=1)
mtext(text="data4.1", side=3, line = -5.7, adj=0, at=2, font=1)
mtext(text="data1.2", side=3, line = -2.2, adj=0, at=3, font=1)
mtext(text="data2.2", side=3, line = -3.4, adj=0, at=3, font=1)
mtext(text="data3.2", side=3, line = -4.5, adj=0, at=3, font=1)
mtext(text="data4.2", side=3, line = -5.7, adj=0, at=3, font=1)
plot(1:10)
dev.off()
  • 任意のワープロソフトで表組みして,Rで作成した画像(pdf でも png でもなんでも)を貼り込んで,print 時に pdf ファイルとして書き出す。。。これは,Macintoshでないとできませんか?
    LaTeX で作表,フロートで画像を取り込んで,pdf 出力する。。それなら,Windows でもできる? -- 青木繁伸 &new{2005-05-04 (水) 00:29:43};
  • ご要望の趣旨とは違いますが、Sweave と呼ばれるツールを使うと Latex 風の原稿に R のコードを埋め込み、コードと出力(eps, pdf画像とともに)を Latex ソースに変換してくれます。これをコンパイルした dvi ファイルを例えば dvipdfmx で pdf ファイルに変換するというやりかたは考えられるかも知れません。つまり画像と解説テキストを一ページの dvi ファイルに収まるようにあんばいし、pdf に変換するわけです。どんな風になるかは example(Sweave) を実行してみて下さい。以下は Sweave ソースファイルsweavetest.Rnw の中身と、R で Sweave("sweavetest.Rnw") を実行してできた sweavetest.tex (および関連画像ファイル)を latex でコンパイルしてできた dvi ファイルを dvipdfmx で pdf ファイルに変換したものです(ここに表示するため更に imagemagic で png ファイルに変換しています。) -- 2005-05-04 (水) 15:48:04
    ?documentclass[a4paper]{article}
    
    ?SweaveOpts{echo=FALSE}
    ?usepackage{a4wide}
    ?usepackage{Sweave} 
    
    ?begin{document} 
    
    <<>>=
    x <- 1:20
    y <- x/10+rnorm(20)
    summary(z <- lm(y ~ x))
    @
    
    ?SweaveOpts{echo=true}
    
    ?begin{figure}[htbp]
      ?begin{center}
    <<fig=TRUE>>=
    par(mfrow=c(2,2))
    plot(z)
    @
        ?caption{linear regression}
      ?end{center}
    ?end{figure}
    
    ?end{document}

#ref(): File not found: "sweavetest.png" at page "初級Q&A アーカイブ(3)"

  • 表と図ですから,xtableとSweaveで. -- なかま 2005-05-04 (水) 16:44:38
    ?documentclass[a4paper]{article}
    ?usepackage{graphicx}
    ?begin{document}
    ?SweaveOpts{echo=false,pdf=FALSE}
    
    <<results=hide>>=
    library(xtable)
    hoge<-data.frame(Tokyo=c(100,200),
    	Osaka=c(80,69),
    	Sapporo=c(80,69),
    	Fukuoka=c(80,69),
    	Naha=c(80,69),
    	Kyoto=c(80,69),
    	Hiroshima=c(80,69),
    	Nagoya=c(90,250)
    	)
    @ 
    
    ?section{table and graph}
    
    &lt;<results=tex>>=
    xtable(hoge, caption="table dayo")
    @ 
    
    ?setkeys{Gin}{width=1.0?textwidth}
    <<fig=TRUE,eps=TRUE>>=
    boxplot(hoge)
    @
    
    ?end{document}

#ref(): File not found: "test1.png" at page "初級Q&A アーカイブ(3)"

  • Sweave は最終的なまとめだけでなく、結果の一時的メモをとるにも便利かも知れませんね。RjpWIki にも使い方の紹介があってしかるべきですね。それにしても、元記事投稿者は Latex は使いたくないということでしょうか。もしそうなら単独の PDF ファイルをどのように使う予定なのでしょう。 -- 2005-05-04 (水) 17:57:12
  • 発言意図がなかなかくみ取れない質問ではありましたが,作表の部分は,本当に必要ならばせめて関数にするなり,もっと汎用性のある関数にしておけばよろしいのではないかと思います。
    temp.pdf というファイルもアップロードしておきましたが表示はできないので,表示用には temp.png に変換したものをアップして表示します。 -- 青木繁伸 2005-05-04 (水) 21:35:58

    #ref(): File not found: "temp.png" at page "初級Q&A アーカイブ(3)"

  • 元記事の参考例を実行してみましたが、やはり表組みの質がいまひとつですね。それともっと複雑な図だと画面を分割すると見にくいかも。やはり R だけでやろうというよりも、Latex 利用や青木さんの例を真似たほうが生産的な気がします。 -- 2005-05-04 (水) 22:10:34
  • 青木先生、なかま様、その他の皆様、ありがとうございます。Sweavesを存じ上げず、LaTeXの最後の記載を理解しておりませんでした。
    青木先生のご指摘どおり、関数の作成とWordでの出力(AcrobatがないためPDF出力はできていませんが)をこれまで行っておりました。
    現在、「測定データから測定対象の品質レポートを作成する」という作業をを別の人にお願いしています。今のWordの作業だとレポート作成に時間がかかってしまい、別の仕事をお願いできない状況です。そこで測定データからPDFの最終レポートまでのルーチンを関数一発で完了させたいと思った次第です。
    私がLaTeXを使いこなせないことも原因ですが、スクリプトの知識のない人でも簡単に表+グラフのレポートを出力できるような関数を作りたいと思った次第です。しかし、今の関数はmtextを使っており、汎用的とは言いがたい状態でした。Sweavesをこれから勉強させていただきます。 -- Akira? 2005-05-06 (金) 09:57:51
  • 私も噂だけで使ってみたことのない Sweave を今回勉強できて、参考になりました。LATEX 環境を整備するのはすこしハードかも知れませんが、マスターすればお話のような定型的作業にはうってつけかと思います。 -- 2005-05-06 (金) 11:38:47
  • PDFLATEX というのを使えば dvi ファイルでなく、PDF ファイルで出力してくれるらしい。ただし、日本語はつかえない? -- 2005-05-10 (火) 05:50:46
  • 揚げ足取りのようなコメントですが、Acrobatがないとpdf出力できないというわけではなく、pdfを作成するたくさんの廉価・無料ソフトがあります。PDFの作り方を参照してください。それからpdflatexは日本が使えません。しかし、platex+dvipdfmxでpsを経由させずにpdfが作成できますので、pdflatexに魅力を感じません。 -- 2005-05-10 (火) 07:01:31
  • 素朴に疑問に思ったことだったのですが、関連するページをご準備していただけて感激しております。今は目下、Texを使えるようになることが目標です。Sweavesまでたどり着けるように頑張ります。
    ちなみに私もplatex+dvipdfmxを使っております。 -- Akira? 2005-05-13 (金) 19:58:08
  • Sweaveを使って結果をまとめることができるようになりました。ありがとうございます。 -- Akira? 2005-10-07 (金) 09:16:51

緯度・経度について

? (2005-05-03 (火) 11:52:12)

 東京の地図をmaptoolsで表示させた上に、地価のポイントを表示させたいのですが、やり方がわかりません。
 東京の地図は表示できたのですが、地価ポイントの緯度・経度のエクセルファイルのみ(シェープファイルなし)で、東京の上に表示させることはできないのでしょうか。

  • できます。その東京の地図と同じ座標系のシェープファイルにその緯度経度を変換し、オーバレイさせます。その際は、緯度経度の測地系に注意してください。私の演習講義ではgen2shpを使ってシェープファイルを作成させていますが、他にもいろいろなツールがあります。 -- 谷村 2005-05-10 (火) 07:15:20
  • 補足です。もし、その東京の地図がUTMや平面直行座標系などではなく経緯度で表示され、なおかつ地価の座標の測地系が同じであるならば、そのエクセルファイルを読み込んで、points()でオーバレイできます。ご参考になれば幸いです。 -- 谷村 2005-05-10 (火) 11:45:57

Mac 版 R 2.1.0 のインストール

初心者? (2005-05-02 (月) 22:14:17)

Windows 版 R 2.1.0 ならば,East Asian Languages でインストールした場合に Rconsole と Rdevga を修正する必要がありますが,Mac 版 R 2.1.0 はそのような作業を行う必要は無いのでしょうか?

  • どういう意味でしょう。何もせずに日本語化されたRが起動されましたけど。 -- 青木繁伸 2005-05-02 (月) 22:46:29
  • 何もせずに日本語化されたRが起動したので,「あれ?何もしなくても日本語化された R が起動したけど,本当に『Rconsole と Rdevga を修正するような作業』はいらないのかな?」と不安になったので質問させていただきました.ご回答ありがとうございます. -- 初心者? 2005-05-02 (月) 23:17:19
  • MacOS X版のGUIはフォントはロケールに応じて使い分けられるので追加作業は全く必要ありません:-) -- 岡田 2005-05-06 (金) 11:31:31

debian(woody)におけるR-2.1.0のapt-get用sources.list

さとう? (2005-05-02 (月) 20:01:21)

debian(woody)においてR-2.0.1をapt-getするためのsources.listファイルの記述はどのようにすればいいのでしょうか。CRANやそのmirrorのサーバーを指定すると,

Ign http://******* Release

のようになりapt-getできません。

  • 今、手元にDebianが無いのでなんですが、 "deb http://cran.r-project.org/bin/linux/debian/stable ./" ではどうでしょうか? -- なかま 2005-05-02 (月) 22:26:05
  • CRAN にある R の deb パッケージは確か保守されていないはず。以下を sources.list に加えて試して下さい。最新 R 2.1.0 日本語版はもしかすると unstable にしかないかも(未確認) -- 間瀬? 2005-05-02 (月) 23:08:23
    # Stable
    deb http://ftp.de.debian.org/pub/debian stable main contrib non-free
    deb http://ftp.de.debian.org/pub/debian-non-US stable/non-US main contrib non-free
    
    # Testing
    deb http://ftp.de.debian.org/pub/debian testing main contrib non-free
    deb http://ftp.de.debian.org/pub/debian-non-US testing/non-US main contrib non-free
    
    # Unstable
    deb http://ftp.de.debian.org/debian unstable main contrib non-free
    deb http://ftp.de.debian.org/debian-non-US unstable/non-US main contrib non-free
  • 念のため確認したら CRAN にも最新版がありました。しかし Debian 本家にもありますからそちらを使う方が便利かも。それに多数のパッケージの deb ファイルもありますしね。 -- 間瀬? 2005-05-02 (月) 23:16:54
  • コメントありがとうございます。Debian.orgのパッケージ一覧をみると、Unstable,TestingはR-2.1.0ですがStableはR-1.5.1になってます。また、CRANにR-2.1.0がありますが、それをapt-getしたくてもうまくいきません。 -- さとう? 2005-05-03 (火) 05:15:18
  • souces.list に以下の一行を加えたら、(少なくとも) at-get update ではエラーが起きませんでした。 -- 間瀬? 2005-05-03 (火) 07:58:26
    deb http://cran.md.tsukuba.ac.jp/bin/linux/debian stable main
  • 再度ありがとうございます。ところで、上述のもので'apt-cache show r-base'するとversionは1.5.1となりはしないでしょうか? -- 2005-05-03 (火) 09:24:08
  • つくば大学の CRAN mirror は本家の完全コピーかと思いましたが、そうでもないようですね。以下の一行を sources.list に加えて apt-get update し apt-cache search r-base としたら 2.1.0 の woody 版が表示されました。最後の stable/ (スラッシュ)が大事のようです。-- 間瀬? 2005-05-03 (火) 23:31:51
    http://cran.au.r-project.org/bin/linux/debian stable/ 
    これでも駄目なら deb ファイルをダウンロードし dpkg -i r-base 等と手作業でインストールするか、ソースからコンパイルしたほうがてっとり早いですね。ところで woody 版に固執する理由があるのでしょうか。
  • 結局、sources.listをRに関連するCRANミラーのもの一行だけにし、'apt-get install'することでwoodyの最新版のものをインストールできました。多分CRANのパッケージとdebian本家のパッケージが衝突(重複)していて、最新版がうまくインストールできなかったのではと思ってます(sources.listの列挙順序の問題でしょうか。)。リリースファイルが無いという警告メッセージ'Ign'はべつにインストール上問題ないことがわかりました。 -- さとう? 2005-05-04 (水) 07:58:09
  • 筑波ミラーにも2.1.0-woodyはきていました。.ac.jp/bin/linux/debian stable/ でも2.1.0-woodyがインストールできるはずです。一方、.ac.jp/bin/linux/debian (sarge|sid|woody|stable|unstable) main でインストールされるバージョンは1.8で古いので本家からは削除されましたが、ミラーには残ってしまっていたようです。 -- 岡田 &new{2005-05-06 (金) 11:16:05}

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