COLOR(green){SIZE(20){初心者のための R および RjpWiki に関する質問コーナー}}~
新規投稿はできません
----
-[[初級Ｑ＆Ａ アーカイブ(4)]] (元記事が 2005-11-09 より 2006-01-25 まで)~
-[[初級Ｑ＆Ａ アーカイブ(3)]] (元記事が 2005-05-02 より 2005-11-07 まで)~
-[[初級Ｑ＆Ａ アーカイブ(2)]] (元記事が 2004-12-13 より 2005-04-26 まで)~
-[[初級Ｑ＆Ａ アーカイブ(1)]] (元記事が 2004-08-03 より 2004-11-30 まで)~
----
#contents
----
**自作関数データのファイルへの取り込み [#xc14f3e7]
>[[maechan]] (2006-06-26 (月) 14:04:44)~
~
ある行列データがありそれを一つのヒストグラムで見ようと思い行列を一列にするプログラムを書きました。それを画面上に表示することはできるのですがファイルに取り込む方法がわかりません。~
関数の値をファイルに書き込む方法が知りたいです。~
//
-関数は -- [[maechan]] &new{2006-06-26 (月) 14:06:36};
-?sink または ?capture.output -- &new{2006-06-26 (月) 14:12:18};
-ありがとうございました。sinkでファイルに書き込むことができました -- [[maechan]] &new{2006-06-26 (月) 14:27:15};
-> 行列を一列にするプログラムを書きました~
プログラムを書くって，as.vector(行列オブジェクト）だけでことが済みますよね。~
それと，ファイルに取り込むの？書き込むの？「関数の値」って何でしょうか。何を書き込んだり取り込んだりするのか。よくわかりませんでした。 -- &new{2006-06-26 (月) 14:59:11};
-適当なエディタを開き，画面（コンソール）に表示されたものをコピーしてエディタのウインドウにペーストして，ファイルセーブするのが，一番簡単では？それに，プログラムだろうと，データだろうと，計算結果だろうと，何でもファイルにできますよ。エディタといえば，R だって自前のエディタを持っているんですよ。 -- &new{2006-06-26 (月) 15:07:58};
**R-2.3.1のインストール手順10のファイルの上書き [#m9ece0ae]
>[[初心者]] (2006-06-24 (土) 12:03:31)~
~
インストールの手順の最後にRconsole，Rdevga，Rprofile.siteをダウンロードし、同名ファイルに上書きとありますが、これらのファイルはどこからダウンロードしたらよいのでしょうか？~
//
-ファイル名のところにリンクが張ってありません？？ -- &new{2006-06-24 (土) 12:05:15};
-インターネットでは（そのほかでも），半角カタカナは使わないようにしましょう。（お節介をして，直しておきました） -- &new{2006-06-24 (土) 12:31:10};
-ファイル名の所を押しても文字が出てくるだけなのですが？ -- &new{2006-06-24 (土) 16:09:24};
-MIMEタイプの設定がうまくできていないサーバーなんでしょうかね。~
このような場合の対処法は，覚えておくと，今後得しますよ。~
そのような場合は，うんとね，ういんどーずのようですから，右ボタンクリックで，リンク先の内容を保存する（ブラウザの種類・バージョンで言葉は違うだろうが）を選べばよいのではないかなと。（それにしても，ユーザに追加の処理を要求しないように，「私でも，インストールできます！」というような，ワンタッチ方式にはできないのだろうか？CRANとしてはむりでも，それをミラーリングする日本のサイトで独自のバイナリファイルをおいてもいいんじゃないか？特に，ういんどーず） -- &new{2006-06-24 (土) 21:00:59};
**関数c()内でのargs[]の使用 [#e44553c1]
>[[チョコボール]] (2006-06-23 (金) 07:21:00)~
~
現在、perlでのスクリプトをバッチモードで使用しています。~
変数を用いているのですが、c()の中では使用できません。~
(こんなかんじ↓)
Y <- X[c(args[5])]
実行は
# R --vanilla --quiet --args test.pdf 3 5 3 < test.R
とやっています。~
他の変数は得られています。~
どなたかおしえてください。
//
-あのね。追試できるように，必要なファイルを全て用意して，それらをどのように操作するかが分かるように記述した方が良いですよ。~
なにも，あなたが実際にやったときのそのままのファイルを見せて！というんじゃなくて（全部見せられたら余計に困る），問題が再現できる短いファイルということです。~
回答してみようかと思う人でも，あなたがどういう風にやったか分からないので，手も足も出せないんじゃないかなと思います。決して，回答は分かっているけど答えないというような，意地悪じゃないと思いますよ。~
全体が見通せれば，「そういうことをやりたいのなら，こんな風にすれば良いんですよ」という別解（しかも，エレガントな別解）が提案されるかもしれません。 -- &new{2006-06-24 (土) 21:26:44};
-すみませんでした。質問のしかたから勉強します。 -- &new{2006-06-24 (土) 22:09:43};
**関数名を文字列に [#l2aaf4df]
>[[ショーン]] (2006-06-22 (木) 19:50:20)~
~
以下のように関数オブジェクトを変数に入れて、関数として呼び出すことができると思いますが、~
> a = cos
> a
.Primitive("cos")
> a(0)
[1] 1
後でaの元の関数名である、"cos"という文字列を得る方法はありますでしょうか？つまり、~
> f2s(a)
[1] "cos"
のf2sようなことはできますでしょうか？~
//
-出来ないことは無いですが、まずなぜそんなことが必要なのか一言説明してもらわないと、答える気にはなりませんね。それに、例のようなことは内部関数にたいしてだけ出来ることですが。 -- &new{2006-06-22 (木) 20:08:53};
-理由がないと教えられないと言うことでもないようには思いますが。~
もっとも，私にはわかりませんでした。~
変数に関数定義を付値したら，その変数の class は function になるので，元の名前と同じように使えますが，「そのようなことは内部関数に対してだけできること」というのもちょっと私にはわかりませんでした。 -- &new{2006-06-22 (木) 21:16:11};
> cubic.root <- function(x) x^(1/3)
> cubic.root(27)
[1] 3
> a <- cubic.root
> a(27)
[1] 3
> a
function(x) x^(1/3)
> class(a)
[1] "function"
-答えるにはそれなりの手間がかかりますから、必要性が私には感じられない質問には（私は）たとえ答えを知っていても答える気にはなりませんよ、という意味です。そんなことに拘らず答えてあげようという親切な方はご自由に。また、直前の例では関数 a には元の関数の情報（名前） "cubic.root”が保存されませんから、a から文字列 "cubic.root" を復元することは原理的に不可能ですね。 -- &new{2006-06-22 (木) 21:40:58};
-勉強の機会を与えてくれて，ありがとうございます。後者の場合に，名前を得る方法（に近いものを）さがしてみました。けど，是では，名無しさんさんが持っている解の方つまり，内部関数に対しては無力です。 -- &new{2006-06-22 (木) 22:25:03};
> f2s <- function(arg1, arg2)
+ {
+ for (i in arg2) {
+ cat(sprintf("ref <- %s?n", i), file="temp")
+ source("temp")
+ if (identical(ref, arg1)) print(i)
+ }
+ }
> cubic.root <- function(x) return(x^(1/3))
> a <- cubic.root
> f2s(a, ls())
[1] "a"
[1] "cubic.root"
> b <- sin
> f2s(b, ls())
[1] "b"
勉強になりました。内部関数の場合でも，たとえば以下のようなものでも，できないことはないですね。~
来るか来ないか分からない答えを待つより，何とか考えてみるのもよいもんですね。
> f2s <- function(a)
+ {
+ sink("temp")
+ print(a)
+ sink()
+ res <- readLines(con="temp")
+ res <- sub("?.Primitive???(???"", "", res)
+ res <- sub("???"?)", "", res)
+ return(res)
+ }
> a <- sin
> f2s(a)
[1] "sin"
> b <- cos
> f2s(b)
[1] "cos"
-うむ、あなたは初級者などでは決してない。何につかうのか秘密のようですが、第二の例は次のようにも出来ます。 -- &new{2006-06-22 (木) 23:33:51};
f2a <- function(x) {capture.output(x ,file=(filename <- tempfile())); unlist(strsplit(readLines(filename),"?""))[2]}
-私は，原質問者じゃないよ。~
だから，原質問者は用途が秘密だとは言っていない。~
質問者をもてあそぶヒマはあっても，答えるヒマはないというのが，鼻についただけ。 -- &new{2006-06-22 (木) 23:41:11};
-暇のある無しの問題ではありません。質問するからには、やはりなぜそんなことを知る必要があるのか（多くの自明な場合は除き）述べる方が、回答を得る可能性が高まりますよ、といいたかっただけ。 -- [[某常連回答者]] &new{2006-06-22 (木) 23:55:28};
-回答ありがとうございます。おかげさまでうまくいきそうです。f2sの使い道ですが、ユーザ次第でいろいろあるとは思いますが、私の場合は、
funcs = c(cos,sin,tan)
のようにして、ループでいろいろな関数のグラフを描こうと思っています。複数のグラフファイルができますが、そのときのファイル名を関数名にするためにf2sを使います。 -- [[ショーン]] &new{2006-06-23 (金) 10:47:09};
-2006-06-22 (木) 23:33:51 は，もっと簡単になるでしょう？ -- &new{2006-06-23 (金) 10:47:31};
f2a <- function(x)
{
unlist(strsplit(capture.output(x),'"'))[2]
}
**引数...の長さ [#u9a72ee2]
>[[ショーン]] (2006-06-22 (木) 12:37:36)~
~
可変長引数...の長さを知りたいのですが、以下のようにうまくいきません。~
> func = function(...){print(length(...))}
> func(1,2,3)
以下にエラーprint(length(...)) : 'length' に対する引数の個数が正しくありません
「[[Rの関数定義の基本]]」を見るとlength(...)という書き方もあるように思えるのですが、何が間違っているのでしょうか？~
//
- 次のようにする。nargs() は関数中でつかわれ、実引数の数を与える。 なぜ length(...) が失敗するかというと、今の場合 length(c(1,2,3)) でなく length(1,2,3) を実行しようとしたから。 -- &new{2006-06-22 (木) 13:09:43};
func <- function(...) print(nargs())
-なるほど。~
純粋に ... の個数を知りたい場合には，以下の方がよいのかな？ -- &new{2006-06-22 (木) 16:36:19};
> func = function(a, ...){print(length(c(...)))}
> func(1,2,3)
[1] 2
-【「Rの関数定義の基本」を見るとlength(...)という書き方もあるように思えるのですが】は記述が間違えているというか抽象的に書かれているので，実際にプログラムの中に書く場合には，可変長引数の一番最後の引数の参照は ...[length(c(...))] と書く必要がある。 -- &new{2006-06-22 (木) 16:46:28};
-ショーンさんはlength(1:3)とlength(c(1,2,3))、length(1,2,3)の区別がついていなかったということですね。 -- &new{2006-06-22 (木) 17:10:19};
-COLOR(blue){【「Rの関数定義の基本」を見るとlength(...)という書き方もあるように思えるのですが】は記述が間違えているというか抽象的に書かれている}......済みません。ちなみに、... 引数の各項目を取り出すには次のように list(...)[[.]] とすれば良いようです。 -- [[参照記事投稿者]] &new{2006-06-22 (木) 17:20:22};
> test <- function(i,...) print(list(...)[[i]])
> test(1,1,"abc",list(1:4),matrix(1:4,2,2))
[1] 1
> test(2,1,"abc",list(1:4),matrix(1:4,2,2))
[1] "abc"
> test(3,1,"abc",list(1:4),matrix(1:4,2,2))
[[1]]
[1] 1 2 3 4
> test(4,1,"abc",list(1:4),matrix(1:4,2,2))
[,1] [,2]
[1,] 1 3
[2,] 2 4
つまり、質問に対する汎用的な答えは length(list(...)) でしょうか。
-ありがとうございました。よくわかりました。 -- [[ショーン]] &new{2006-06-22 (木) 19:27:52};
**データ行列の行をまたぐレコード参照 [#b26a6989]
>[[ginga]] (2006-06-21 (水) 15:19:35)~
~
以下のような時系列データがあるとします
time ID1 ID2 val1 val2
1 1077525908 1 2 35.256 90.188
2 1077529016 1 3 33.664 90.509
3 1077532131 1 3 30.543 68.935
4 1077535240 1 3 36.345 123.351
5 1077538354 1 5 36.408 93.581
6 1077541476 1 5 37.916 106.016
7 1077544589 1 2 35.809 93.094
... . . ... ...(数十万行)
ここで各行について，以下の処理をしたいというのは，R ではどのように書くのでしょうか?~
for を使ってやるしかないでしょうか?~
~
result[i] = (((他の行のtime - i行めのtime)< 10000) であり，
かつ (val1 < 35.0) となるレコードについて
ID2 のラベルが何種類あるか(ID2 が同じものは無視して可))
例えば 4行めだと ID2,3 が該当するので答えは 2 という感じです~
(time では一応ソート済なので，処理の際に前後100行のみ，という感じにある程度スコープを絞ることはできますが，それがまたどう書くのか分からないです)~
~
大量のデータ列の行間の関係を洗うような処理について，なにか良い資料というか，例などがあれば有り難いです．~
~
配列の添字に条件をいろいろ書けばできるような話なのか，添字部分に書く書き方を探すというのがそもそも間違いなのか...~
~
R向きの考え方でないのであれば改めて別言語で検討します．~
~
よろしくお願いします~
//
-lapplyと比較演算子を組み合わせればできそうに思いますが -- &new{2006-06-21 (水) 15:30:07};
-例に挙げたデータについての例解が間違えているし。条件の記述も間違えている。(他の行のtime - i行めのtime)< 10000)は引き算の順序が逆？val1<35なのに val1>35 のものを対象にしているし。 -- &new{2006-06-21 (水) 17:36:37};
-まずは，val1 < 35.0 という条件なら，それを満たさないものは事前に除外した方が探索時間は短くなる。データフレーム名を d とすれば，d <- d[d[,4]<35,] でとりあえず。 -- &new{2006-06-21 (水) 17:38:54};
-しかしまあ，val1 < 35.0 は，第 i 行に関する条件なのか，比較される他の行に関する条件なのか，それとも両方なのか。曖昧だねぇ。曖昧だとプログラムが書けません。とりあえず，上では，両方ということにしたけど，解釈を間違えていそうだねぇ。 -- &new{2006-06-21 (水) 17:42:30};
-いろいろ面倒だが，骨格を作ったから直してみて？val1 についてのチェックは，定義が曖昧だから，プログラムに含めていない。定義がどうであろうと，簡単に組み込めるだろう。~
たぶん遅いから，R じゃなきゃいけないとか R の機能をふんだんに使っているというのではないのだから，Cなんかで書く方が良いでしょう。 -- &new{2006-06-21 (水) 17:53:05};
nr <- nrow(d) # 行数
sapply(1:nr, function(i)
{
ti <- d[i] # i 行目の time
lo <- max(1, i-100) # i 行目の前100行（前か後だけでよいのだろうが）
hi <- min(nr, i+100) # i 行目の後100行（前か後だけでよいのだろうが）
temp <- d[lo:hi,] # とりあえず対象可能性のある行列を取り出す
ok <- abs(temp[,1]-ti) < 10000 # 時間差をチェック（abs は要らないのだろうが，そのときは不等号の向きを確認）
temp <- d[ok, 3] # ID2 の列を取り出して
length(table(temp)) # 度数分布を求めてその長さを見れば，カテゴリーが何種類あったかわかる
})
-ありがとうございます．lapply() 調べてみます&ご提示頂いた sapply() を参考に練ってみます．~
その他の皆様も，アドバイスありがとうございます(たしかに例が正確ではなかったかも) -- [[ginga]] &new{2006-06-21 (水) 18:52:15};
-上のコードと本質的に同じコードで50万行の人工例を（一部実行で）試してみたら、所要時間19時間余りと予測されました（笑）。本当に全部を一度に処理する必要があるかどうか熟考した方がよさそうです。普通そんなことは無いものですが。 -- &new{2006-06-21 (水) 22:46:41};
-一昼夜コンプータを動かしておけば解が出るなら，安いモンでしょう(^_^) -- &new{2006-06-21 (水) 23:08:33};
**Williamsの多重比較のパーセント点 [#p8461227]
>[[tm]] (2006-06-20 (火) 22:23:00)~
~
Rとは直接関係なく、大変申し訳ありませんが、Williamsの多重比較のパーセント点の表を、web経由で入手できるところはありませんでしょうか？~
身近で入手できる本をあさったのですが、5%点くらいしか手に入りませんでした。~
できれば、1%点の表が欲しいのですが、どなたかご存知ありませんか？~
//
-その本に，生成アルゴリズムなどは書いてなかったですか。。。あるいは，その本の底本となっている文献にも？Williams の検定はRでもできるようなので，そのソースをハッキングしてみては？ -- &new{2006-06-20 (火) 22:34:48};
**.Internalの中身 [#o72ed84b]
>[[のぞき]] (2006-06-20 (火) 14:35:37)~
~
.Internalが呼ばれる関数dnormなどの
dnorm
function (x, mean = 0, sd = 1, log = FALSE)
.Internal(dnorm(x, mean, sd, log))
<environment: namespace:stats>
の実際のコードを覗いてみることはできないのでしょうか。~
//
-過去に同じ類の質問が沢山あります。ちゃんと調べましょう。 -- &new{2006-06-20 (火) 14:48:11};
-.Internal については，あったかな？ -- &new{2006-06-20 (火) 16:33:37};
R-2.3.x/src/nmath ディレクトリに
dnorm.c pnorm.c qnorm.c rnorm.c という C ソースがある
**主成分分析の主成分と元の変数の関係 [#sb27b706]
>[[ショーン]] (2006-06-19 (月) 09:40:45)~
~
prcompを使って求まる主成分PC1,PC2...と、元の変数xの関係を知りたいです。以下のサイトを見ると、行列の掛け算であるようなので、実際に試してみました。~
http://aoki2.si.gunma-u.ac.jp/lecture/PCA/pca1.html~
しかし、以下のようにうまくいきません。~
~
print(d)
x1 x2 x3
1 0.1991 -15.0277 0.2790
2 0.2003 -16.8090 0.3462
3 0.2169 -16.6008 0.2990
4 0.1429 -17.3694 0.2586
5 0.1840 -16.7424 0.2597
6 0.2231 -15.0885 0.3075
7 0.1834 -15.7661 0.1550
8 0.1516 -15.1187 0.2940
9 0.1478 -14.4152 0.3411
10 0.2071 -17.1567 0.2789
11 0.1613 -13.5147 0.2592
12 0.2288 -15.0654 0.2165
res = prcomp(d,scale = TRUE)
prims = res[["x"]];
l = as.matrix(res[["rotation"]]);
x = as.vector(as.matrix(d[1,]))
pr = as.vector(as.matrix(prims[1,]))
pr2 = l%*%x;
print(x)
[1] 0.1991 -15.0277 0.2790
print(pr)
[1] -0.1190277 0.1568528 -0.6743579
print(pr2)
[,1]
x1 -0.793312
x2 -5.952815
x3 13.779837
prとpr2の値が同じになると思っていたのですが、違います。どこかで計算ミスをしているのでしょうか？~
//
-主成分負荷量を求めるのでしょうか？だとしたら，貴方は，説明を読み違えているのでしょう。 -- &new{2006-06-19 (月) 13:14:23};
> d <- structure(c(0.1991, 0.2003, 0.2169, 0.1429, 0.184, 0.2231, 0.1834,
+ 0.1516, 0.1478, 0.2071, 0.1613, 0.2288, -15.0277, -16.809, -16.6008,
+ -17.3694, -16.7424, -15.0885, -15.7661, -15.1187, -14.4152, -17.1567,
+ -13.5147, -15.0654, 0.279, 0.3462, 0.299, 0.2586, 0.2597, 0.3075,
+ 0.155, 0.294, 0.3411, 0.2789, 0.2592, 0.2165), .Dim = c(12, 3
+ ), .Dimnames = list(c("1", "2", "3", "4", "5", "6", "7", "8",
+ "9", "10", "11", "12"), c("x1", "x2", "x3")))
> # eigen 関数から求める
> res <- eigen(cor(d))
> print(t(sqrt(res$values)*t(res$vectors)))
[,1] [,2] [,3]
[1,] 0.7725894 -0.04950733 -0.6329728
[2,] -0.7142710 -0.37685279 -0.5897448
[3,] -0.2484622 0.92942168 -0.2728404
> # prcomp 関数から求める
> res2 <- prcomp(d, scale.=TRUE)
> print(t(res2$sdev*t(res2$rotation)))
PC1 PC2 PC3
x1 0.7725894 0.04950733 -0.6329728
x2 -0.7142710 0.37685279 -0.5897448
x3 -0.2484622 -0.92942168 -0.2728404
-ご回答ありがとうございます。すみません。基本的なことを確認させてください。教えていただいた、
t(res2$sdev*t(res2$rotation))
は、ここのサイト(http://aoki2.si.gunma-u.ac.jp/lecture/PCA/pca1.html )でいう相関係数(因子負荷量)aですよね。それから、
res2$rotation
が行列Lですよね。それから、
res2$x
が主成分zですよね。投稿した質問のように、元の変数xから主成分zを求めるには、xにaではなくLをかけるのではないのでしょうか？ですので、res2$xとL%*%xが一致するはずだと思っているのですが、どこで勘違いをしているのでしょうか？ -- [[ショーン]] &new{2006-06-19 (月) 17:39:26};
-だから，貴方が何を求めたいのかを聞いたのに。貴方が求めたいのは，主成分得点ですね。「主成分PC1,PC2...と、元の変数xの関係」なんて言うから主成分負荷量かと勘違いしましたよ。~
主成分得点を求めるときに L に掛けるのは標準化したデータですよ。~
主成分得点の計算については，~
http://aoki2.si.gunma-u.ac.jp/lecture/PCA/pca3.html ~
http://aoki2.si.gunma-u.ac.jp/lecture/PCA/pca4.html ~
を見ないと。 -- &new{2006-06-19 (月) 18:34:20};
> res2 <- prcomp(d, scale.=TRUE)
> res2$x
PC1 PC2 PC3
1 -0.11902772 0.15685279 -0.6743579
2 0.58870153 -1.58534961 -0.1294253
3 1.07405220 -0.65752984 -0.3523574
4 -0.07378289 -0.30135908 1.9991581
5 0.54804949 -0.05898904 0.7081565
6 0.35523899 -0.32780390 -1.3586354
7 0.46001105 2.09814940 0.7988055
8 -1.25606299 -0.21446936 0.3827542
9 -1.93793258 -0.83621938 -0.1796872
10 1.23546343 -0.49028619 0.2885350
11 -1.75203822 0.91644637 -0.5043241
12 0.87732771 1.30055785 -0.9786220
> scale(d)%*%res2$rotation
PC1 PC2 PC3
1 -0.11902772 0.15685279 -0.6743579
2 0.58870153 -1.58534961 -0.1294253
3 1.07405220 -0.65752984 -0.3523574
4 -0.07378289 -0.30135908 1.9991581
5 0.54804949 -0.05898904 0.7081565
6 0.35523899 -0.32780390 -1.3586354
7 0.46001105 2.09814940 0.7988055
8 -1.25606299 -0.21446936 0.3827542
9 -1.93793258 -0.83621938 -0.1796872
10 1.23546343 -0.49028619 0.2885350
11 -1.75203822 0.91644637 -0.5043241
12 0.87732771 1.30055785 -0.9786220
-ありがとうございました。よくわかりました。 -- [[ショーン]] &new{2006-06-22 (木) 12:35:00};
**Sweave環境でsink()は使えないのでしょうか？ [#fda5d2d5]
>[[akira]] (2006-06-16 (金) 14:11:29)~
~
Sweave環境でxtableを書き換えたい場合、~
(1)xtableの出力先をtexファイルからsink()で別のファイル（例えば"temp"）へ切り替えて、~
(2)readLinesでtempを読み込んで修正して、~
(3)texファイルへ書き出す~
ということを試してみました。しかし、sinkは受け付けてくれませんでした。~
何かオプションがあるのでしょうか？~
?documentclass{jarticle}
%
?begin{document}
?section{sinkでtestが書き出せない}
<<dat1, echo=F>>=
sink("test")
cat("ここはsinkでtestに出るはず?n")
cat("testに出てますよ")
cat("ここまで?n")
sink()
cat("test2?n")
x <- readLines("test")
cat("読み出したtestを書き出すよ?n")
cat(x,"?n")
cat("ここまで?n")
@
?section{table出力もできない}
<<dat2, echo=F, results=tex>>=
library(xtable)
x <- data.frame(a=1:3, b=5:7)
sink("test")
cat("ここはsinkでtestに出るはず?n")
print(xtable(x, caption="変更前のtable"))
cat("????", "?n")
cat("ここまで?n")
sink()
try(x <- readLines("test")) # ここがエラーでとまる
cat("????", "?n")
x[grep("caption", x)] <- sub("変更前", "変更後", x[grep("caption", x)])
x <- paste(x, "?n")
cat("読み出したtestを書き出すよ?n")
cat(x,"?n")
cat("????", "?n")
cat("ここまで?n")
@
?end{document}
そこで、xtableの出力をxに入れなおして、readLinesの部分を書き換えると上手く動きます。
# xtableをxに入れる
x <- xtable(x, caption="変更前のtable")
# readLinesの部分はこうする
x <- unlist(strsplit(x, split="?n")) # testがないからxtableの出力を変更する
'RweaveLatex', 'Rtangle'当たりを見てもそれらしい記述が見当たりません。
//
-[[print.xtable2:http://sapporo.cool.ne.jp/matsut/r02.txt]]を発見。しかし、これを修正するのは大変そうだ。 -- [[akira]] &new{2006-06-19 (月) 16:36:45};
**multcomp パッケージでの Williams 検定について [#pa7f07ff]
>[[tm]] (2006-06-15 (木) 21:10:37)~
~
multcomp パッケージ中の simtest の Williams 検定についてですが、うまく動作していない気がします。~
Rのバージョンは 2.0.1、multcomp のバージョンは 0.4-8、mvtnorm のバージョンは 0.7-2 です。~
> library(mvtnorm)
> library(multcomp)
> data(angina)
> summary(simtest(response ~ dose, angina, type="Williams", alternative="greater"))
Simultaneous tests: Williams contrasts
Call:
simtest.formula(formula = response ~ dose, data = angina,type = "Williams",
alternative = "greater")
Williams contrasts for factor dose
Contrast matrix:
dose0 dose1 dose2 dose3 dose4
C1 0 -1 0.00 0.0000000 0.0000000 1.0000000
C2 0 -1 0.00 0.0000000 0.5000000 0.5000000
C3 0 -1 0.00 0.3333333 0.3333333 0.3333333
C4 0 -1 0.25 0.2500000 0.2500000 0.2500000
Absolute Error Tolerance: 0.001
Coefficients:
Estimate t value Std.Err. p raw p Bonf p adj
C1 10.499 6.778 1.549 0 0 0
C2 7.747 5.775 1.341 0 0 0
C3 6.297 4.979 1.265 0 0 0
C4 5.247 4.284 1.225 0 0 0
といった結果です。また、Interval の計算結果も~
> summary(simint(response ~ dose, angina, type="Williams", alternative="greater"))
Simultaneous 95% confidence intervals: Williams contrasts
Call:
simint.formula(formula = response ~ dose, data = angina, type = "Williams",
alternative = "greater")
Williams contrasts for factor dose
Contrast matrix:
dose0 dose1 dose2 dose3 dose4
C1 0 -1 0.00 0.0000000 0.0000000 1.0000000
C2 0 -1 0.00 0.0000000 0.5000000 0.5000000
C3 0 -1 0.00 0.3333333 0.3333333 0.3333333
C4 0 -1 0.25 0.2500000 0.2500000 0.2500000
Absolute Error Tolerance: 0.001
95 % quantile: 1.978
Coefficients:
Estimate 5 % -- t value Std.Err. p raw p Bonf p adj
C1 10.499 7.435 Inf 6.778 1.549 0 0 0
C2 7.747 5.093 Inf 5.775 1.341 0 0 0
C3 6.297 3.795 Inf 4.979 1.265 0 0 0
C4 5.247 2.824 Inf 4.284 1.225 0 0 0
となってしまいます。p値はすべて 0 ですし、Interval に Inf が含まれていますし、結果がおかしい気がするのですが…~
WEBで検索しても、multcomp の Williams 検定については情報が見つけられませんでした。~
使用方法や結果の見方などご教授いただけないでしょうか？~
よろしくお願いいたします。~
//
-あるプログラムにバグがあるというとき，確証がないといけません。そのデータは本来どのような結果になるべきものなのか分からないと，プログラムにバグがある，ちゃんと計算されていないということはできないと思います。「p値はすべて 0 ですし、Interval に Inf が含まれていますし」というのは，根拠にならないと思います。以下に示すものも，結果を検証したものではありませんが，少なくとも，pはすべて0でもないし，Interval に Inf が含まれるのは，片側推定(greater)だからですよね（less の場合には -Inf になるし）。「使用方法や結果の見方などご教授いただけないでしょうか」ということでは，「バグがあるよ〜」とはいえないのではないでしょうか。~
プログラムの使い方などについては，教科書などに載っている答えの分かっているデータを分析してみれば，どこのどの数値がどれに対応しているんだなあとか，使い方を間違えていないんだという確信が得られるとか。そう言う風にやっていくのがよろしいのではないでしょうか。 -- &new{2006-06-15 (木) 21:39:54};
> library(mvtnorm)
> library(multcomp)
> set.seed(12345678)
> df <- data.frame(g = factor(rep(1:6, each=20)), y = rnorm(120)+rep(1:6/6, each=20))
> summary(simtest(y ~ g, df, type="Williams", alternative="greater"))
Simultaneous tests: Williams contrasts
Call:
simtest.formula(formula = y ~ g, data = df, type = "Williams",
alternative = "greater")
Williams contrasts for factor g
Contrast matrix:
g1 g2 g3 g4 g5 g6
C1 0 -1 0.0 0.00 0.0000000 0.0000000 1.0000000
C2 0 -1 0.0 0.00 0.0000000 0.5000000 0.5000000
C3 0 -1 0.0 0.00 0.3333333 0.3333333 0.3333333
C4 0 -1 0.0 0.25 0.2500000 0.2500000 0.2500000
C5 0 -1 0.2 0.20 0.2000000 0.2000000 0.2000000
Absolute Error Tolerance: 0.001
Coefficients:
Estimate t value Std.Err. p raw p Bonf p adj
C1 0.771 2.751 0.280 0.003 0.017 0.008
C2 0.464 1.911 0.243 0.029 0.117 0.047
C3 0.299 1.309 0.229 0.097 0.290 0.124
C4 0.242 1.094 0.221 0.138 0.290 0.156
C5 0.092 0.426 0.217 0.335 0.335 0.335
Warning message:
full precision was not achieved in 'pnt'
> summary(simint(y ~ g, df, type="Williams", alternative="greater"))
Simultaneous 95% confidence intervals: Williams
contrasts
Call:
simint.formula(formula = y ~ g, data = df, type = "Williams",
alternative = "greater")
Williams contrasts for factor g
Contrast matrix:
g1 g2 g3 g4 g5 g6
C1 0 -1 0.0 0.00 0.0000000 0.0000000 1.0000000
C2 0 -1 0.0 0.00 0.0000000 0.5000000 0.5000000
C3 0 -1 0.0 0.00 0.3333333 0.3333333 0.3333333
C4 0 -1 0.0 0.25 0.2500000 0.2500000 0.2500000
C5 0 -1 0.2 0.20 0.2000000 0.2000000 0.2000000
Absolute Error Tolerance: 0.001
95 % quantile: 1.976
Coefficients:
Estimate 5 % -- t value Std.Err. p raw p Bonf p adj
C1 0.771 0.217 Inf 2.751 0.280 0.003 0.017 0.008
C2 0.464 -0.016 Inf 1.911 0.243 0.029 0.146 0.057
C3 0.299 -0.153 Inf 1.309 0.229 0.097 0.483 0.166
C4 0.242 -0.195 Inf 1.094 0.221 0.138 0.690 0.227
C5 0.092 -0.336 Inf 0.426 0.217 0.335 1.000 0.476
-pが『0』に見えるのは出力の有効桁数だけの問題では？全ての要素は確認済みですか？ -- [[akira]] &new{2006-06-16 (金) 13:00:34};
-ご返信、アドバイスありがとうございます。 ~
まず、いいたかったのは「バグがあるよ。」ということではなく、「私の環境ではうまく動作していないようです。」あるいは「使い方がわかりません。」ということで、使用方法等インターネットで検索してみましたが、見つかりませんでしたので、皆様に助言いただこうとしたものでした。~
おっしゃられるように、正解例がないとよくわかりませんので、群馬大の青木先生のページでも参考にしておられた「統計的多重比較法の基礎」（サイエンティスト社）より例題を引っ張ってきました。
> library(mvtnorm)
> library(multcomp)
> set.seed(12345)
> df <- data.frame(
+ Dose = ordered(rep(c("Cont", "Low", "MidLow", "MidHigh", "High"), each=7), c("Cont", "Low", "MidLow", "MidHigh", "High")),
+ Response = c(
+ 415, 380, 391, 413, 372, 359, 401, # Cont
+ 387, 378, 359, 391, 362, 351, 348, # Low
+ 357, 379, 401, 412, 392, 356, 366, # MidLow
+ 361, 351, 378, 332, 318, 344, 315, # MidHigh
+ 299, 308, 323, 351, 311, 285, 297 # High
+ ))
> summary(simtest(formula=Response~Dose, data=df, type="Williams", alternative="less"))
Simultaneous tests: Williams contrasts
Call:
simtest.formula(formula = Response ~ Dose, data = df, type = "Williams",
alternative = "less")
Williams contrasts for factor Dose
Contrast matrix:
DoseCont DoseLow DoseMidLow DoseMidHigh DoseHigh
C1 0 -1 0.00 0.0000000 0.0000000 1.0000000
C2 0 -1 0.00 0.0000000 0.5000000 0.5000000
C3 0 -1 0.00 0.3333333 0.3333333 0.3333333
C4 0 -1 0.25 0.2500000 0.2500000 0.2500000
Absolute Error Tolerance: 0.001
Coefficients:
Estimate t value Std.Err. p raw p Bonf p adj
C1 -79.571 -7.076 11.245 0 0 0
C2 -63.500 -6.520 9.739 0 0 0
C3 -45.571 -4.963 9.182 0 0 0
C4 -39.714 -4.467 8.890 0 0 0
となります。前述の本によると、~
t1 = 7.076~
t2 = 4.218~
t3 = 1.417~
となるらしく、Cont と High の比較ですでに t value の値が正負逆ですし、その後の比較では値がまったく異なり、さらに、本では t3 のところで帰無仮説が保留されるのに対し、multcomp の Williams 検定によると、すべての比較において棄却されています。~
この違いは、私の R 環境の問題でしょうか？それとも、実装されている手法が本のものと違うのでしょうか？~
ご教授いただけますよう、よろしくお願いいたします。
-- [[tm]] &new{2006-06-20 (火) 10:32:31};
-Contrast マトリックスと件の本を比べると違いがわかるのでは？件の本は，コントロール以外の群の複数の組み合わせで線形結合を作ってその最小値を与える群とコントロールを比較している。simtestでは単に線形結合とコントロールの比較でしょう。tの符号が違うのは，反応が小さい順に並んでいるか大きい順に並んでいるかとか検定の方向が greater なのか less なのかによる。 -- &new{2006-06-20 (火) 12:02:36};
-contrast マトリックスの意味がなんとなく理解できました。やりたい検定に対して、入力の方法が間違っているのか、それとも simtest の contrast マトリックスが間違っているのかまではわかりませんが。。。~
ただ、おっしゃられるように反応の並びや greater、 less あたりをいじっても正負は変わりませんが、greater、less の変更による p 値の変化は納得できるので、こうゆうものなのでしょうか？~
実装が、本と異なるのは理解できましたが、contrast マトリックスを見る限り、行いたい検定をするには、
> df <- data.frame(
+ Dose = ordered(rep(c("Cont", "Low", "MidLow", "MidHigh", "High"), each=7), c("Cont", "High", "MidHigh", "MidLow", "Low")),
+ Response = c(
+ 415, 380, 391, 413, 372, 359, 401, # Cont
+ 387, 378, 359, 391, 362, 351, 348, # Low
+ 357, 379, 401, 412, 392, 356, 366, # MidLow
+ 361, 351, 378, 332, 318, 344, 315, # MidHigh
+ 299, 308, 323, 351, 311, 285, 297 # High
+ ))
> summary(simtest(formula=Response~Dose, data=df, type="Williams", alternative="less"))
Simultaneous tests: Williams contrasts
Call:
simtest.formula(formula = Response ~ Dose, data = df, type = "Williams",
alternative = "less")
Williams contrasts for factor Dose
Contrast matrix:
DoseCont DoseHigh DoseMidHigh DoseMidLow DoseLow
C1 0 -1 0.00 0.0000000 0.0000000 1.0000000
C2 0 -1 0.00 0.0000000 0.5000000 0.5000000
C3 0 -1 0.00 0.3333333 0.3333333 0.3333333
C4 0 -1 0.25 0.2500000 0.2500000 0.2500000
Absolute Error Tolerance: 0.001
Coefficients:
Estimate t value Std.Err. p raw p Bonf p adj
C4 -39.714 -4.467 11.245 0.000 0.000 0.000
C3 -26.429 -2.878 9.739 0.004 0.011 0.007
C1 -22.143 -1.969 9.182 0.029 0.058 0.042
C2 -15.929 -1.636 8.890 0.056 0.058 0.056
とするのが近いのかも。~
本と同じ検定をするためには、simtest の中を調べる必要がありそうで、ちょっとがっかりです。。。 -- [[tm]] &new{2006-06-20 (火) 13:31:22};
**RからWebへのアクセス [#t4afdaac]
>[[Goeppingen]] (2006-06-14 (水) 10:42:34)~
~
Ｒ内より、ＷｅｂにアクセスしてHTMLやＸＭＬを取得できるのでしょうか？~
//
-?url -- &new{2006-06-14 (水) 11:00:37};
-? download.file -- &new{2006-06-14 (水) 11:06:49};
**"fitted probabilities numerically 0 or 1 occurred" [#ica65f0c]
>[[kaz]] (2006-06-12 (月) 16:00:47)~
~
glmを用いておよそ100万件のデータをロジスティック回帰しています。説明変数は５個です。たまに首記のエラーが発生してしまいます。日本語版のエラーは「数値的に 0 か 1 である確率が生じました」です。~
~
インターネットで検索したところ、perfect separation、つまりある説明変数XによってYが完全に説明されてしまう（Y=1 if and only if X<X0 for an X0）場合に当該エラーが出るようですが、100万件もありますし、またいくつかデータをチェックしましたが、そのようなことは起こっていないようです。（perfect separation参考：https://stat.ethz.ch/pipermail/r-help/2005-July/074693.html）~
~
他にこのエラーの生じる原因をご存知であれば教えていただけますでしょうか。なお、100万件もデータを用いることの是非については別問題としてください（データが多すぎることがこのエラーを引き起こしているのならば別ですが、同じ数のデータで回帰してもエラーが起こらないこともあります）~
//
-また、5つの説明変数は全て連続量です。Rと回帰分析の初心者で申し訳ありませんが、どんなコメントでも大変うれしく思います。 -- [[kaz]] &new{2006-06-12 (月) 16:04:19};
-Rのバージョンは2.2.1（日本語化版）です。 -- [[kaz]] &new{2006-06-12 (月) 16:21:33};
-Quasi-complete separation とか Complete separation というのを調べましょう。~
簡単に説明すれば，ある変数（または複数の変数の線形結合）の分布を見たとき，二つのグループが完全に分かれるとき（分布が重ならないとき）問題のエラーが生じるのです。それは，連続変数であろうと無かろうと起こります。以下に簡単な例を挙げましょう。3変数で，従属変数が0を取るケースが25，1を取るケースが25としましょう（多くても少なくても良い）。乱数を発生させてロジスティック回帰を行って見ましょう。うまくいきます。分析の前に第一列の変数の最大値と最小値を群ごとに示してあります。分布は重なっていますね。一方の群にだけこの変数に定数を加えて分布が重ならないようにしてやると，件のエラーメッセージがめでたく出現します。~
データが何百万，何千万あろうとも，ある変数の分布が重ならない（つまり平均値が相当違う等）ときにはよくあること。~
共用パソコンに入っているなら仕方ないけど，最新版を使うのが吉。 -- &new{2006-06-12 (月) 16:31:10};
> require(MASS)
[1] TRUE
> set.seed(12345)
> y <- rep(0:1, each=25)
> x <- matrix(rnorm(150), ncol=3)
> cat(min(x[1:25,1]), max(x[1:25,1]),"?n")
-1.817956 1.817312 ここと
> cat(min(x[26:50,1]), max(x[26:50,1]),"?n")
-2.380358 2.196834 ここの数値を見ると分布は重なっていない
> summary(res <- glm(y ~ x, family=binomial))
Call:
glm(formula = y ~ x, family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.65228 -1.06530 -0.06326 1.03977 1.68695
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.1365 0.3115 -0.438 0.661
x1 0.3894 0.2823 1.379 0.168
x2 0.2135 0.2689 0.794 0.427
x3 -0.4449 0.2768 -1.607 0.108
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 69.315 on 49 degrees of freedom
Residual deviance: 64.623 on 46 degrees of freedom
AIC: 72.623
Number of Fisher Scoring iterations: 4
> x[1:25,1] <- x[1:25,1]+4.2 定数を加えて分布が重ならないようにしてやる
言い換えると，x[,1]の分布を描くと，二群が完全に分離しているのだ
> cat(min(x[1:25,1]), max(x[1:25,1]),"?n")
2.382044 6.017312 ここと
> cat(min(x[26:50,1]), max(x[26:50,1]),"?n")
-2.380358 2.196834 ここの数値を見れば，分布が全く重なっていないことがわかる
> summary(res <- glm(y ~ x, family=binomial))
Call:
glm(formula = y ~ x, family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-5.647e-05 -2.107e-08 0.000e+00 2.107e-08 5.560e-05
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 229.574 100367.677 0.002 0.998
x1 -99.764 42439.569 -0.002 0.998
x2 -6.802 19068.857 -0.000357 1.000
x3 19.210 15815.467 0.001 0.999
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 6.9315e+01 on 49 degrees of freedom
Residual deviance: 9.7282e-09 on 46 degrees of freedom
AIC: 8
Number of Fisher Scoring iterations: 25
Warning messages:
1: アルゴリズムは収束しませんでした in: glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,
2: 数値的に 0 か 1 である確率が生じました in: glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,
-大変わかり易い具体例を有難うございます。頭がすっきりしました。確認させていただけないでしょうか。Y~X1+X2+X3の回帰で、当該エラーが（perfect separationを理由として）発生した場合、Y~X1, Y~X2, Y~X3の回帰のうち少なくとも一つで同様のエラーが発生しなければならない。これは正しいでしょうか？ -- [[kaz]] &new{2006-06-12 (月) 17:28:43};
-私のやり方が間違っているのかもしれませんが、上記のように単回帰に分解してみたものの、当該エラーは発生しませんでした。 -- [[kaz]] &new{2006-06-12 (月) 17:29:50};
-先ほど追加した部分と関連しますが，複数の変数の線形結合値の分布が重ならないという場合には，個々の分析がうまくいったのにというのは理由になりません。個々の分析はうまくいったのに，全部使うとダメになるということもあるでしょう。 -- &new{2006-06-12 (月) 17:32:28};
-set.seed も含めて。上の通りやってみましたか？set.seed がないなら，データが重ならないように加える数値を調整しないと再現できませんよ。x[,1]だけ使った結果を以下に載せておきましょうか？ -- &new{2006-06-12 (月) 17:35:58};
> set.seed(12345)
> y <- rep(0:1, each=25)
> x <- matrix(rnorm(150), ncol=3)
> x[1:25,1] <- x[1:25,1]+4.2
> summary(res <- glm(y ~ x[,1], family=binomial))
Call:
glm(formula = y ~ x[, 1], family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.256e-04 -2.107e-08 0.000e+00 2.107e-08 1.282e-04
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 460.8 118212.0 0.004 0.997
x[, 1] -201.3 51633.8 -0.004 0.997
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 6.9315e+01 on 49 degrees of freedom
Residual deviance: 3.2194e-08 on 48 degrees of freedom
AIC: 4
Number of Fisher Scoring iterations: 25
Warning messages:
1: アルゴリズムは収束しませんでした in: glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,
2: 数値的に 0 か 1 である確率が生じました in: glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,
-また、（これは単にRのバージョンの問題かもしれませんが）warning messageはご提示していただいた具体例の２．のみで１．の「アルゴリズムは収束しませんでした」は私の回帰では表示されませんでした。 -- [[kaz]] &new{2006-06-12 (月) 17:41:48};
-なるほど、確かに単回帰のエラーは重回帰のエラーの十分条件ではあるが必要条件で無いと理解しました。タイムリーなご回答有難うございます。 -- [[kaz]] &new{2006-06-12 (月) 17:49:59};
-set.seed で生成される乱数列は，バージョンによっても違うのかな？だからこそ，よほどの事情がない限り最新バージョンを使うべし。それと，コメントはよいのだけど，たとえば先の私のコメント set.seed は入れましたかとか，追試結果に関係するコメントに対する答えもちゃんと書いた方がいいよ。後で，他の人が見るときに，消化不良になるよ。 -- &new{2006-06-12 (月) 17:53:52};
-個々の変数の分布は重なるが，線形結合の分布が重ならない例 -- &new{2006-06-12 (月) 18:07:32};
> set.seed(12345)
> y <- rep(0:1, each=25)
> x <- round(matrix(rnorm(150)*10+50, ncol=3), 1)
> z <- apply(c(1,2,3)*t(x), 2, sum) # 線形結合を作って
> oz <- order(z)
> x <- x[oz,] # その大きい順に並べ替えた
> z <- z[oz] # Z は確認のために表示するだけ
> cbind(x, z) # 確認
z
[1,] 26.2 41.4 37.1 220.3
[2,] 35.9 44.9 34.5 229.2
[3,] 68.2 36.6 36.8 251.8
[4,] 56.1 54.8 28.8 252.1
[5,] 40.8 26.5 53.2 253.4
中略
[47,] 55.2 65.9 63.2 376.6
[48,] 34.5 73.3 65.3 377.0
[49,] 61.2 55.2 71.8 387.0
[50,] 52.5 59.7 76.6 401.7
> # それぞれの変数を使って分析するとうまくいくが
> summary(res <- glm(y ~ x[,1], family=binomial))
Call:
glm(formula = y ~ x[, 1], family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.66515 -1.06756 0.05372 1.06560 1.69950
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.50275 1.60011 -2.189 0.0286 *
x[, 1] 0.06747 0.03016 2.237 0.0253 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 69.315 on 49 degrees of freedom
Residual deviance: 63.530 on 48 degrees of freedom
AIC: 67.53
Number of Fisher Scoring iterations: 4
> summary(res <- glm(y ~ x[,2], family=binomial))
Call:
glm(formula = y ~ x[, 2], family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.79526 -0.77141 0.04637 0.90898 1.81508
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -6.78401 2.06861 -3.280 0.001040 **
x[, 2] 0.12770 0.03841 3.324 0.000886 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 69.315 on 49 degrees of freedom
Residual deviance: 52.668 on 48 degrees of freedom
AIC: 56.668
Number of Fisher Scoring iterations: 4
> summary(res <- glm(y ~ x[,3], family=binomial))
Call:
glm(formula = y ~ x[, 3], family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.65849 -0.81491 -0.05309 0.67276 1.95259
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -7.58755 2.15363 -3.523 0.000426 ***
x[, 3] 0.15294 0.04298 3.558 0.000373 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 69.315 on 49 degrees of freedom
Residual deviance: 48.393 on 48 degrees of freedom
AIC: 52.393
Number of Fisher Scoring iterations: 4
> # 全部を使ったらこける
> summary(res <- glm(y ~ x, family=binomial))
Call:
glm(formula = y ~ x, family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-9.418e-05 -2.107e-08 0.000e+00 2.107e-08 6.798e-05
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2164.803 719087.033 -0.003 0.998
x1 5.672 2238.150 0.003 0.998
x2 14.675 5396.431 0.003 0.998
x3 21.969 7225.664 0.003 0.998
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 6.9315e+01 on 49 degrees of freedom
Residual deviance: 1.7709e-08 on 46 degrees of freedom
AIC: 8
Number of Fisher Scoring iterations: 25
Warning messages:
1: アルゴリズムは収束しませんでした in: glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,
2: 数値的に 0 か 1 である確率が生じました in: glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,
-詳しい例示ありがとうございます。「個々の変数の分布は重なるが，線形結合の分布が重ならない例 」ということですが、逆では無いでしょうか。」私は依然、「単回帰のエラーは重回帰のエラーの十分条件ではあるが必要条件で無い」と理解しています。もしこの理解が正しければ非常にすっきりと理解できたかと思います。有難うございました。今は、「では、(quasi)complete separation」となった場合の処方箋を調べています。追ってご報告いたします。 -- [[kaz]] &new{2006-06-13 (火) 10:37:37};
-> 「個々の変数の分布は重なるが，線形結合の分布が重ならない例 」ということですが、逆では無いでしょうか。~
上の例を見ればおわかり頂けると思ったのですが。分布が重ならないから完全に分割できるのであって（complete separation）その場合に，上述のようなエラーが出るんですよ。~
その方が理解しやすいということならかまいませんが，必要条件とか十分条件とかいう言葉を持ち出す必要もないと思います。 -- &new{2006-06-13 (火) 11:01:00};
上の x から合成される z を，y 別に描いてみる
重なっていないことがわかる
boxplot(z ~ y, horizontal=T)
#ref(boxplot.png)
- ずいぶん時間が経ってから見させていただきましたが、参考になりました。ありがとうございます。ところで、結局このような場合はどうしたらいいのでしょうか？記載のように3つの変数全て使うとWarningが出る場合は、３つの変数を使った解析はしてはいけないのですか、それともよく説明されていてすばらしい、ということなのでしょうか？素人で済みません -- [[unta]] &new{2008-07-25 (金) 14:48:54};
-Firth bias-correctionという方法があるようです -- nu &new{2015-02-17 (火) 04:19:46};
--参考
---http://cran.r-project.org/web/packages/brglm/index.html
---http://www.ats.ucla.edu/stat/mult_pkg/faq/general/complete_separation_logit_models.htm
**Sweaveでxtableをminipage環境に適用したいです [#x624f617]
>[[Akira]] (2006-06-09 (金) 10:40:42)~
~
latexコマンドをRスクリプトに埋め込んで、Sweaveで結果をまとめています。~
xtableで作成したtableをminipage環境で2つ横に並べたいのですが、Sweave()が出力するtexファイルの該当箇所を探して修正するのは面倒なので、minipage環境でxtableを出力したいと思っています。
> library(xtable)
> x <- data.frame(matrix(1:12, 2))
> xtable(x)
% latex table generated in R 2.3.0 by xtable 1.3-2 package
% Fri Jun 09 10:27:57 2006
?begin{table}[ht]
?begin{center}
?begin{tabular}{rrrrrrr}
?hline
& X1 & X2 & X3 & X4 & X5 & X6 ??
?hline
1 & 1.00 & 3.00 & 5.00 & 7.00 & 9.00 & 11.00 ??
2 & 2.00 & 4.00 & 6.00 & 8.00 & 10.00 & 12.00 ??
?hline
?end{tabular}
?end{center}
?end{table}
この出力にminipage環境のコマンドを埋め込みたいです。
% latex table generated in R 2.3.0 by xtable 1.3-2 package
% Fri Jun 09 10:27:57 2006
?begin{table}[ht]
?begin{minipage}{0.5?hsize}
?begin{center}
?begin{tabular}{rrrrrrr}
?hline
& X1 & X2 & X3 & X4 & X5 & X6 ??
?hline
1 & 1.00 & 3.00 & 5.00 & 7.00 & 9.00 & 11.00 ??
2 & 2.00 & 4.00 & 6.00 & 8.00 & 10.00 & 12.00 ??
?hline
?end{tabular}
?end{center}
?end{minipage}
#2列の場合は?begin{minipage}〜?end{minipage}が繰り返される
?end{table}
こんな感じです。~
汎用性を持たせる技術がないので、今は
> cat("??begin{table}[ht]?n?n")
> y <- print(xtable(x)) #(1)
> y <- sub("????begin??{table??}??[ht??]?n", "", y)
> y <- sub("????end??{table??}?n", "", y)
> print(cat(y))
> cat("??end{table}?n?n")
としていますが、この場合(1)のところの出力がtexファイルに書き込まれてしまうので、結局修正しなればなりません。~
printのオプションを確認しましたが、出力しない設定が分かりませんでした。RSiteSearchでもminipage環境の情報が少なく、困っています。~
何か他に良い方法がありますでしょうか？
R-2.3.0、WinXPSP1を使用しています。
//
-print.xtableをみましょう. カラクリ屋敷(ぼそっ)... -- [[なかま]] &new{2006-06-09 (金) 14:00:32};
?documentclass[a4j]{jsarticle}
?usepackage[dvipdfm]{graphicx}
?title{Happy xtable}
?date{}
?begin{document}
<<minipage,echo=T>>=
library(xtable)
set.seed(123)
X<-array(runif(10*5),dim=c(10,5))
@
?begin{table}[ht]
?begin{minipage}{0.5?hsize}
<<fig=F,echo=F,results=tex>>=
print.xtable(xtable(X)[1:5,],floating=F)
@
?end{minipage}
?begin{minipage}{0.5?hsize}
<<fig=F,echo=F,results=tex>>=
print.xtable(xtable(X)[6:10,],floating=F)
@
?end{minipage}
?end{table}
?end{document}
-理想通りの結果です。ありがとうございます。print.xtalbeはsee alsoにありました。すみません。欲しい機能は大概あるものですね。 -- [[akira]] &new{2006-06-09 (金) 15:25:14};
-なかまさまの前に記入されていたsinkを利用する方法がなくなっています。こちらも別の目的で使えそうなので、戻してほしいのですが。 -- [[akira]] &new{2006-06-11 (日) 21:35:38};
-バックアップをたどればありますよ. (良く見ろとか怒られるかも...) で, たぶん戻して欲しいと言うのはややアレ(消した理由もあるんでしょう)なので, Sweave関連のTIPS頁にサマリーとしてakiraさんの方でまとめると言うのはいかがでしょうか. -- [[なかま]] &new{2006-06-12 (月) 01:24:29};
-勝手にしていいものかと・・・でも、Tipsに乗せてしまいました。sinkを教えてくれた方、ご了承ください。 -- [[akira]] &new{2006-06-12 (月) 12:17:02};
**dates関数でout.formatをyyyy/mm/ddにするには [#aa1c7ffd]
>[[mame]] (2006-06-08 (木) 20:13:58)~
~
dates関数を使い数値データをyyyy/mm/ddに整形するにはどうすればいいでしょうか？~
下記のように実行してみたのですが、月だけmon形式になってしまいます。
> dates(as.numeric(13149),out="yyyy/mm/dd")
[1] 2006/Jan/01
因みに"yyyy/m/dd"でも結果は同じでした。~
R var2.3.0~
//
-パッケージ survival には数値を様々な形式で出力する関数が用意されていますが、ご希望の順番に並べるものは無いようです。でも date.mmddyyyy という関数がありますからそのコードをすこしいじくればすみます。面倒なら date.mmddyyyy で我慢しましょう。ところで私の Linux 版 R 2.3.0 には dates という関数は見当たらないのですが、MSW 版固有の関数でしょうか？ -- &new{2006-06-08 (木) 20:46:56};
> date.yyyymmdd <- function(sdate, sep="/") {
require(survival) # もちろんシステムにあらかじめインストールしておきましょう
temp <- date.mdy(sdate) # date.mdy は survival の関数
ifelse(is.na(sdate), "NA", paste(temp$year, temp$month, temp$day, sep=sep) )}
> date.yyyymmdd(100:105)
[1] "1960/4/10" "1960/4/11" "1960/4/12" "1960/4/13" "1960/4/14" "1960/4/15"
ふむ、何か変かな。次のようになるのだけれど。
> date.yyyymmdd(13149)
[1] "1996/1/1"
-パッケージ data にも survival にも全く同じ関数がありますね。 -- &new{2006-06-08 (木) 21:17:11};
-dates関数って，どのライブラリにあるのか。。。。~
date.yyyymmdd は，以下のように微修正すればよいでしょう。 -- &new{2006-06-08 (木) 21:53:32};
> date.yyyymmdd <- function(sdate) {
+ require(survival) # もちろんシステムにあらかじめインストールしておきましょう
+ temp <- date.mdy(sdate) # date.mdy は survival の関数
+ ifelse(is.na(sdate), "NA", sprintf("%4i/%02i/%02i", temp$year, temp$month, temp$day ))
+ }
> date.yyyymmdd(100:105)
[1] "1960/04/10" "1960/04/11" "1960/04/12" "1960/04/13" "1960/04/14" "1960/04/15"
> date.yyyymmdd(13149)
[1] "1996/01/01"
-13149 て Julian day (1960/01/01) からの経過日数ですよね。最初の例の結果が何か変なような気が？ -- &new{2006-06-08 (木) 22:40:31};
-変も何も，そういう前提で経過日数後の日付を yyyy/mm/dd の形式で表示したいと言うことなんでしょう？意図まで疑うときりがない。ちなみに，月の名前が英字ではいるのが気に入らないと言うのは関数仕様ではどうしようもないので，自分で関数を書くか，他の書式を制御できる関数を使うか，不満足な関数の出力を修正するか（つまり，英字の月名部分を切り出してそれを数値01〜12の二文字で置き換えて結果とする関数を書くか。。まあ，どれでも大差ない。あるいは，ユリウス日の関数から自分で書くか。 -- &new{2006-06-08 (木) 22:48:20};
> dates2 <- function(x)
+ {
+ s <- as.character(dates(as.numeric(13149),out="yyyy/mm/dd"))
+ sprintf("%s%02i%s",substring(s, 1, 5),which(month.abb==substr(s,6,8)),substr(s,9,11))
+ }
> dates2(13149)
[1] "2006/01/01"
-うん？
> dates(0,out="yyyy/mm/dd")
[1] 1970/Jan/01
なんだから，dates 関数では 1970/01/01 からの日付では？~
結局は，関数の定義によるようで，別の関数で代替処理するときにはチェックが必要。自分で書くのが一番確か。 -- &new{2006-06-08 (木) 22:52:20};
-datesって S-Plusの関数では? -- &new{2006-06-09 (金) 07:31:39};
-google で少し調べたら Julian date とは本来紀元前4713/1/1の正午を起点とする主に天文学で使われる日付単位とか(そして Julian とはユリウス暦の Julius とは無関係）。Unix 界では1970/1/1/0:00:00 からの経過秒で表す慣習が。さらに S では 1960/1/1 からの経過日数を使うようで、R はこれを継承しているのか。ややこしい。 -- &new{2006-06-09 (金) 07:55:32};
-元の発言にも記載がないので，いったいどのdates関数だと，議論する人も混乱するのでは？標準でないライブラリの関数については，ちゃんとどのライブラリであるか書くべきですね。私がみたdates関数はchronライブラリのもの。これは，1970/01/01からカウントしているようですよ。ああ，本当にややこしい!! -- &new{2006-06-09 (金) 08:20:13};
> library(chron)
> dates(0)
[1] 01/01/70
> dates(0,out="yyyy/mm/dd")
[1] 1970/Jan/01
> dates(13149,out="yyyy/mm/dd")
[1] 2006/Jan/01
-皆さんどうもありがとう御座いました。~
dates関数はchronライブラリのものです。現在SからRへの移植の仕事をしていますが、どちらもど素人で苦労しています。始めはdates関数はRにない！！と思い、別の関数を探していましたが、いろいろ調べているうちにchronライブラリが使えるのかも･･･と考えた次第です。それで筑波大学さんから[chron_2.3-3.zip]をダウンロードさせていただきました。~
きちんとdates関数についての記述をきちんとせず、皆さんを混乱させてしまって済みませんでした。~
今回は自分で関数を書いてこの問題を解決しようと思います。 -- [[mame]] &new{2006-06-09 (金) 10:24:29};
**分散の違い [#hfeda6c6]
>[[kanako]] (2006-06-07 (水) 14:21:24)~
~
異なる平均と標準偏差を持つ分布を3種類用意して、それらを順番に足したものを3列作成する際、別々に分けて作成したものと、一度に作成したものでは分散が違うようなのですが、何が違うのでしょうか？Rは2.3を使っていますが、2.1でも同じでした。いろいろ調べてみたのですが、何が悪いのか良く分かりませんでした。またどちらが正しいのでしょうか？
> a <- rnorm(100,0,1)
> b <- rnorm(100,3,3)
> c <- rnorm(100,5,5)
> abc <- cbind(a,a+b,a+b+c)
> var(abc)
a
a 1.0695965 1.092034 0.7527298
1.0920344 9.323825 7.0911109
0.7527298 7.091111 30.6758695
> test <- cbind(rnorm(100,0,1),
rnorm(100,0,1)+rnorm(100,3,3),
rnorm(100,0,1)+rnorm(100,3,3)+rnorm(100,5,5))
> var(test)
[,1] [,2] [,3]
[1,] 0.9833900 0.10720269 0.56137488
[2,] 0.1072027 10.47080201 -0.03146996
[3,] 0.5613749 -0.03146996 34.98877571
~
//
-乱数は文字通りランダムですから、毎回同じ物が作られるとは限りません。敢えて同じものを発生したければ乱数種を同じにしないといけません。例えば set.seed(1234) をさせると同じ乱数を発生させられますが、それでも二度目の例では、総計６組の乱数ベクトルを発生させていますから、やはり同じにはなりませんよ。つまり、三種類の rnorm(100,0,1) はちがうものです。ためしに次の例を実行してみてください。（これってもしかして学校の宿題？） -- &new{2006-06-07 (水) 14:45:14};
> rnorm(3,0,1)
> rnorm(3,0,1)
> set.seed(1234)
> rnorm(3,0,1)
> set.seed(1234)
> rnorm(3,0,1)
-コメントありがとうございました。説明が不足していました。乱数なので、毎回違い値になることは想像していたのですが、2つ目の方法で作った時にほとんどいつもマイナスがどこかに含まれるので、（最初の方はだいたいいつも正の数）何か意味が違うのかなと思い、こちらへ投稿しました。（ちなみに宿題ではないです。）コメントありがとうございました。 -- [[kanako]] &new{2006-06-07 (水) 14:52:29};
-ちなみに、二つの例で分散（と言うよりは共分散行列ですね）が大幅にことなるのは、第一の例では相関を持つ（同じ a,b が含まれる）のに、第二の例では三つのベクトルは６種類の独立なベクトルから作られる和だから互いに独立になるからです（老婆心）。var の代わりに cor を使うとこの点がもっとはっきり分かりますね。 -- &new{2006-06-07 (水) 15:03:44};
**データから円グラフを作成 [#y864b83a]
>[[ヤキソバパンマン]] (2006-06-06 (火) 21:58:03)~
~
例えば、
2,1,3,2,2,2,1,3,1,2,2,2,1,2
という.csvファイルがあるとして、そのデータを比率に直し、~
円グラフを書くスクリプトがわかりません。~
2が8個、1が4個、3が2個なので、2の割合が大きいですが…~
この計算などを行って、円グラフで出力したいです。~
どなたか教えてください。
~
//
-このサイトのグラフィック参考実例集の pie chart を見てください。データは度数で与えても自動的に比率（角度）に変換されて描画されます。cvs ファイルを読み込み、度数データに集計するには table 関数を使います。たとえば x がデータなら pie(table(x)) とするだけ. 注釈はどうする、などと関連質問がありそうですが、まず自分で考えましょう。-- [[MKR]] &new{2006-06-06 (火) 23:22:04};
-ありがとうございます。頑張ってみます。 -- [[やきそばぱんまん]] &new{2006-06-08 (木) 13:43:50};
**プロビットモデルでの疑似決定係数の計算 [#o080db1a]
>[[イマイ]] (2006-06-05 (月) 23:18:44)~
~
プロビットモデルの回帰分析を行っていますが、疑似決定係数、あるいは尤度比インデックスを計算させるにはどうすればよいのか、教えていただけませんか。~
超初心者なりに考え付く限り検索をかけてみたのですが、みつけられませんでした。よろしくお願い致します。~
//
-前にも同じ質問をどこかでしてませんでしたっけ？というのは，どうでもよいが。~
「プロビットモデル 疑似決定係数」でググってみると，20件くらいしか出てこないけど，[[http://www.jil.go.jp/jil/zassi/200110abe.html:http://www.jil.go.jp/jil/zassi/200110abe.html]]には，
「12）プロビットモデルで疑似決定係数を計算する方法はいくつか存在するが，ここでは１−（対数尤度−制限つき対数尤度）で計算した。」なんて書いてありますね。それでわからなければ，この著者「阿部正浩」さんに連絡してみればいかがですか？ -- &new{2006-06-06 (火) 13:47:16};
-尤度比インデックスに関しては存じ上げませんが，疑似決定係数に関しては以下の操作で解決すると思います。第1に，google で，「probit pseudo filetype:r」と検索すると7件表示されます。Robert Andersen の講義HP (http://socserv.mcmaster.ca/andersen/esrc.html ) からSolutions to Logit and Probit Exercises にアクセスします。疑似決定係数を算出する関数を定義してくれていますので，どのタイプの疑似決定係数か，教科書でも見ながら確認すると，解決すると思います。 -- [[mint]] &new{2006-06-06 (火) 17:15:37};
-検索がまだまだ不十分だったようです。アドバイス、本当にありがとうございました。 -- [[イマイ]] &new{2006-06-06 (火) 21:15:39};
**Mac OSXでのR 2.3.1でパッケージアップデートができない [#d91636d7]
>[[高井]] (2006-06-05 (月) 17:21:12)~
~
PowerPC G5（デュアル 2.3GHz）, OSX 10.4.6でR 2.3.1のGUI版の「パッケージとデータ」-->「パッケージインストーラ」を選んでバイナリーを「すべてアップデート」しようとすると，「パッケージのアップデートに失敗しました。詳細についてはコンソールを参照してください」となります。コンソールには~
以下にエラーupdate.packages(lib = ""/Library/Frameworks/R.framework/Resources/library/"", :
オブジェクト "Library" は存在しません
とあります。ところが，ターミナルから起動してupdate.packages()するとエラーはでません。バグでしょうか？~
//
-補足です：エラーメッセージが出た後に，Rパッケージインストーラのウィンドウにはバージョンが表示されます。現在の私の状態は最新の導入バージョンですので，エラーメッセージだけが出ているだけかもしれません。（エラーメッセージが出てもアップデートされるのかどうかは不明です） -- [[高井]] &new{2006-06-05 (月) 17:25:03};
-R コンソールから，update.packages() すると，バージョンの新しいライブラリを検索してくれて，インストールするかどうか答えてくれますね。ask=FALSE を指定すると黙って更新してくれます。 -- &new{2006-06-05 (月) 17:37:11};
-はい，そうですね。Rコンソールの結果もターミナルからの結果も同じです。でも，メニューバーからupdateする時にエラーメッセージがでるのです。 -- [[高井]] &new{2006-06-05 (月) 20:04:22};
-Libraryの前のダブルクォートが２つになっているのが怪しいですね...バグかもしれません。あとからソースみてみますが、最近きた翻訳アップデート依頼でパッケージマネージャまわりがあったので、最新版ではfixされているかもです -- [[岡田]] &new{2006-06-07 (水) 23:05:29};
-翻訳はアップデートされてたけどアップデートは直ってないですね...デバッグしなくちゃ(^_^;; -- [[岡田]] &new{2006-06-09 (金) 13:17:35};
-間違いなくバグでした。PackageInstaller.mの中でupdate.packages()を呼ぶときにクォートが多すぎです。レポートしておきます。 -- [[岡田]] &new{2006-06-09 (金) 16:15:56};
-私の環境のせいではなくて一安心です。 -- [[高井]] &new{2006-06-09 (金) 21:14:07};
**Step関数の実行 [#n2eba3cc]
>[[R初心者]] (2006-06-01 (木) 19:56:59)~
~
Windows XPでR2.2.1を用いて回帰分析を行っています。~
step関数でAICによる交互作用も含めた変数選択を行おうとしたのですが、
> slm2 <- step(lm(y~(A+B+C+D+E)^2))
エラー：サイズ 1027089 Kb のベクトルを割り当てることができません
追加情報: Warning messages:
1: Reached total allocation of 734Mb: see help(memory.size)
2: Reached total allocation of 734Mb: see help(memory.size)
> summary(slm2)
以下にエラーsummary(slm2) : オブジェクト "slm2" は存在しません
となって計算が実行できません。~
ちなみにメモリーは
> memory.size(T)
[1] 61349888
です。~
Rではハードディスクを用いた仮想メモリを利用できる機能等はないのでしょうか?~
//
-help(memory.size) はしてみましたか？ この wiki 内を「メモリ」で検索するだけでも有用な情報は得られると思うんですがね。 -- &new{2006-06-02 (金) 02:36:35};
-help(memory.size)も実行済みです。memory.limit(size=1024)も行ってみましたが、NULLが返ってきます。 -- [[R初心者]] &new{2006-06-02 (金) 14:14:31};
-NULLが返ってきた後に、memory.limit()をすると上限が増えてませんか？ -- &new{2006-06-02 (金) 15:31:49};
-(A+B+C+D+E)^2を独立変数にした回帰で，しかも変数選択をしたいというのが適切な状況があまり思い浮かばないので，モデル自体の吟味をすることを，まずお薦めしますが，Windows版では，起動コマンドラインで，--max-mem-size=2500000000などと2GBより大きな値を指定しておくと，最大限の2GBメモリが使えるようになります（物理メモリとして実装していなくても）。それでもメモリが不足する場合，[[Windows版のFAQの2.9:http://cran.md.tsukuba.ac.jp/bin/windows/base/rw-FAQ.html#There-seems-to-be-a-limit-on-the-memory-it-uses_0021]]によれば，WindowsのBoot.iniに/3GBスイッチをつけると仮想メモリで3 GBまで使えるようになるらしいですが，未確認です。それでもダメならLinux版を使うのも一案でしょう。 -- [[中澤]] &new{2006-06-02 (金) 15:38:26};
-memory.limit()で確かに上限が増えています。が、計算が上手くいきません。~
--モデル自体の吟味をすることを，まずお薦めします~
そうかもしれません。実は回帰と言うより、実際は数量化I類で、Aが10水準、Bが2水準、Cが20水準、Dが4水準、Eが20水準のカテゴリーに分かれていて、総データ数が138,920件です。y~A+B+C+D+Eで回帰を行ったところ、全ての説明変数でのAICが一番低く、そう言った意味では変数削除は出来ないのです。かつ、自由度調整済み寄与率が0.9922と大変高い値を示しているので、ひょっとしたらこのままで構わないのかもしれませんが、念の為に交互作用を検出できるかどうか調べてみたかったのです。~
ところで｢起動コマンドライン｣とは何でしょうか?単語検索でも引っ掛かりませんでした。--[[R初心者]] &new{2006-06-02 (金) 16:40:56};
-Rを起動するアイコンの上でマウスを右クリックしてプロパティを選ぶと，「リンク先」のところに"C:?Program Files?R?R-2.3.1?bin?Rgui.exe"となっていると思います。その後にスペースを空けて，オプションとして--max-mem-size=2500000000と打てばいいです。 -- [[中澤]] &new{2006-06-02 (金) 17:29:36};
-それを｢起動コマンドライン｣と言うのですね。実は--sdiにしているんですが、重複させて --max-mem-size=2500000000も打ち込んで構わないのでしょうか? -- [[R初心者]] &new{2006-06-02 (金) 18:18:41};
-起動アイコンのコマンドラインパラメータなので，そう呼んでいるだけです。重複というか，２つ以上つけても大丈夫です。 -- [[中澤]] &new{2006-06-02 (金) 20:26:04};
-なるほど。Tipsありがとうございます。数量化I類は結局どのみちメモリ不足で計算できませんでしたが、メモリ容量を確保できる方法論が分かったので良かったです。 -- [[R初心者]] &new{2006-06-02 (金) 20:44:53};
-カテゴリー変数をダミーデータに変換して重回帰分析を行う場合，変数選択によって落とされたカテゴリー（ダミー変数）がある場合，その元締めであるカテゴリー変数自体を落とさないと，引き続く変数選択過程に齟齬を来すのではなかったでしょうか。あるいは，ダミー変数が一つ落とされたら，落とされたもの以外と落とされたものをプールしたダミー変数というようにデータ行列自体を再構成しないいけないのでは。駒沢先生が数量化における変数選択について本を書いていたような記憶が。。。。それにしても，カテゴリーが20個もあるとか，データが14万ケースとか，自由度調整済みの決定係数が0.9922 とか，すごいですね。このレベルになると，普通の統計解析レベルの話で収まるのでしょうか。 -- &new{2006-06-03 (土) 00:45:17};
-Rのfactor関数を使ってダミー化したものですから、元締めの変数が落とせる物なのかどうかちょっと分かってません。ただ、総カテゴリカル変数の内、係数では有意じゃなかったのは6つしかありませんでした。ですから、この6つを削除してまた回帰をやりなおせばいいのでしょうか?駒沢先生の本は未読なもので・・・・・・。~
#ref(QQplot.jpg)
あれ、画像添付こんなのでいいのかな･･････?~
ご覧のように、正規Q-Qプロットを一応試みたのですが、大部分は直線に一致していますが、右に行くほど外れ値が出てきてしまいます。~
そんなわけで、交互作用の存在を疑ってみたのです。-- [[R初心者]] &new{2006-06-03 (土) 03:19:55};
-step関数は大変便利な方法だと思いますが、強引な面も持っています。~
まずは、主効果のみモデルのAICを求めて、~
その後、2次の交互作用項を一つ入れたモデルを作って~
AICが改善されるかを見ていくのもいいと思います。-- [[MK]] &new{2006-06-03 (土) 10:33:52};
-なるほど。力技ですね。試してみます。 -- [[R初心者]] &new{2006-06-03 (土) 15:37:30};
-バグってRが止まってしまいました。。。。。。○|￣|＿~
ただ、交互作用A:Bを入れたモデルだけは計算できたんですが、やはり予想通りAICは減少しています。と言うことは、何らかの形で交互作用が存在する可能性がある、って事ですね。何とか調べられればいいのですが･･････。~
ところで、AIC関数を用いるとstep関数で言うRSSが表示されるみたいですが、この2つは何が違うのでしょう? -- [[R初心者]] &new{2006-06-03 (土) 16:43:19};
-その交互作用項のt値はどうでしたか？有意であれば、交互作用項は影響を与えているといえます。 -- [[MK]] &new{2006-06-05 (月) 16:13:30};
-交互作用項A:BのP値は全て 2×10^-16以下を示しています。
#ref(QQplot2.jpg)
交互作用A:Bを考慮した場合の正規Q-Qプロットです。あんまり見た目が変わりません。と言うことは、依然とまだどっかに交互作用がある、と言う事でしょうね。 -- [[R初心者]] &new{2006-06-05 (月) 16:25:19};
-交互作用を疑うよりは，81510番のデータを良く精査して，問題があるようなら除いて分析するのが先決かと（14万件もの内の一例なんで，影響はないのかも知れませんが，非常に気持ち悪い）。14万件もあると，係数の有意性検定もほとんど意味がないのではないだろうか。 -- &new{2006-06-05 (月) 16:44:46};
-81510番は確かにイレギュラーに見えるので、次のようにして削除してみました。~
df <- df[-c(81510),]
それで元の交互作用項を含まないモデル式に戻して、正規Q-Qプロットを試みてみました。
#ref(QQplot4.jpg)
今度は乖離が派手に見えます。指数関数的に増加していってますね。-- [[R初心者]] &new{2006-06-05 (月) 18:30:46};
-これは，交互作用と言うよりは，独立変数と従属変数が直線相関で内政ではないでしょうか？横軸を Theoretical Quantiles ではなくて，予測値にして，回帰分析なさるのがよいのでは？非常に得意な標準化残差の様相を呈していると思います。予測値の大小にかかわらず，標準化残差が0を中心として上下にばらつくのが重回帰モデルで期待される結論でしょう。（図がちょっと大きすぎませんか） -- &new{2006-06-05 (月) 18:51:43};
-大きすぎると思います。wiki上で図を小さくする方法が分からないのです。すみません。~
正規Q-Qプロットより残差プロットかS-Lプロットの方が宜しいのでしょうか? -- [[R初心者]] &new{2006-06-05 (月) 18:54:17};
-最初から小さな図を描けばいいのです。ま，ちょっと描いてみればいかがですか？データは貴方しか持っていないのだから，我々はいろいろ分析したり図を描いたりすることができないのですから。（なんか，R と無関係なまま）~
残差は以下のような分布になることが期待されています。同じデータ（結果に対する）q-q plot は2つめのグラフのようになるんです。 -- &new{2006-06-05 (月) 19:14:48};
#ref(temp.png)
#ref(temp2.png)
-大量のデータを全部使って分析しなければ得られない結果というのはそうはないように思います。あるモデルで表現できるシステムならば，データ数は少なくても同じように表現できるはず。データ数が少ないと困るのは，たとえば回帰係数の有意性検定で有意になりにくいというようなこと。逆に，データ数が多すぎれば，なんでもかんでも有意にはなるはろくなことはないように思いますよ。近代統計学はサンプルを用いて母集団を推測するという点を切り開いたわけで，データ・マイニングというのは，まあ，母集団を分析しようとしているようなものではないのでしょうかね。たとえば，統計調査というのは，ほぼ2000〜3000人に調査すればよい，それ以上である必要はないということになっているようですね。しかし，一方で，たとえば因子分析で安定した有効な因子構造を見るには変数の数十倍（だったかな）のデータが必要とか。それにしても，項目100個で数万でしょ？ -- &new{2006-06-05 (月) 20:22:44};
-確かにRとあんまり関係無くなってますね。・・・・・・どっかで大容量メモリーのコンピュータがあればいいのですが。
#ref(回帰診断.jpg)
Rで画像を小さくする、と言うのはこんな感じで宜しいのでしょうか?~
間瀬先生の本に従って、見よう見まねで回帰診断プロットを行ってみました。標準化残差は0を中心にして上下にバラついてるとは思いますが････。~
--あるモデルで表現できるシステムならば，データ数は少なくても同じように表現できるはず。~
--近代統計学はサンプルを用いて母集団を推測するという点を切り開いたわけで~
仰る事は分かります。その点も気になって、逆に、｢データが多ければ多いほど推定精度が高くなる｣と言ったアイディアである、｢ベイズ推定｣での数量化I類も実際試してみたのです。最尤推定じゃなくっても、理論的にはベイズ推定での数量化I類も可能なのではないか?と。そこでMCMCpackでのMCMCregressを試してみました。結果は･･････計算に30分以上かかり、しかもバグって止まってしまったようです。~
場合に拠ってはWinBugsを試してみようかな、とは思っていますが。 -- [[R初心者]] &new{2006-06-05 (月) 20:16:44};
-0を中心にしてばらついているというのではなく，予測値が70〜150ぐらいのときには非対称（残差ですか？が10を超えるようなものがいくつかある。Q-Qプロットの右の方の急に立ち上がっているものが是に由来しているのでしょう。これらは，要するに外れ値です。その外れ値がどのような場合に現れているのかは把握しましたか？その外れ値が出てくる要因を独立変数で定義できればQ-Qプロットも直線上に載るのでしょう。sqrt(Standardized residuals)のグラフをみると，若干線形性からの乖離が見られるようではありますがこの図は外れ値を含んでいるので薄められている（実際はもっと線形からはずれている）のかもしれません。 -- &new{2006-06-05 (月) 20:36:04};
-外れ値の影響が大きいんですね。という事は・・・・ロバスト/抵抗回帰で数量化I類を行ってみるのも一つの手でしょうか?｢S-PLUSによる統計解析｣によると、lqs関数が推奨されていましたが、これはRのライブラリーにも存在するのでしょうか? -- [[R初心者]] &new{2006-06-05 (月) 20:54:30};
-とか言ってましたがダメでしたね。
エラー：lqs failed: all the samples were singular
ロバスト/抵抗回帰と数量化I類は相性が悪い?連続量じゃないとダメなんですね。一つ賢くなりました。 -- [[R初心者]] &new{2006-06-05 (月) 21:01:57};
**windows版Rのインストール [#sae95e3a]
>[[ありんこ]] (2006-06-01 (木) 16:51:01)~
~
すみません、中級者の所にしゃしゃり出てしまって…。インストールの説明のところに、2.2.1、2,2,0、2,1,0、2,1,1の四つがあったのですが…。初心者にはどれが良いのでしょう。~
//
-最新版の[[R-2.3.0:http://cran.md.tsukuba.ac.jp/]]を使いましょう。 -- &new{2006-06-01 (木) 17:11:41};
-と書いた矢先にR-2.3.1が。上のリンク先（筑波大ミラーサイト）も2.3.1になってますね。 -- &new{2006-06-02 (金) 09:57:51};
-確かにインストールの説明のところに新しい2.3.1が追加されていませんね。 -- &new{2006-06-05 (月) 23:25:06};
-「インストールの説明のところに新しい2.3.1が追加されていませんね」って，すぐには対応できないのでしょう。そのページを書いた人だって，ボランティアでやっているんだから，責任取る必要もないんだから。気づいた人がやれば？ -- &new{2006-06-06 (火) 13:17:53};
-実はみんな違う人が書いたんだったりして:-) -- [１つは書いた覚えがある者] &new{2006-06-09 (金) 17:31:12};
**ノンパラメトリックなモデルでの変数選択 [#i7bc81e2]
>[[金井]] (2006-06-01 (木) 01:35:15)~
~
【やろうとしていること】~
20個程度の説明変数があり，その説明変数を利用して，目的変数（0/1）を予測しようとしています．~
~
【現状】~
現在，R2.2.1を，Windows XPで利用中です．~
線形判別分析，ロジスティック回帰分析の変数選択に関してはAICに基づく変数選択を過去ログから発見できたのですが，分類木（treeパッケージ），ニューラルネット（nnetパッケージ）に関してはここの過去ログを含め，Webを調べても情報がありませんでした．~
~
【お聞きしたい点】~
お聞きしたのは，２点です．~
~
１．ニューラルネットや分類木で，線形判別分析で提供されている関数 step() のような簡単な仕組みで変数選択を行う方法はありますか？~
~
２．そもそもニューラルネットや分類木で0/1の２値を予測する際に，変数選択を行うことが間違いでしょうか．~
//
-質問２は R プロパーの話題ではありませんね。統計プロパーの話題については、たまたま知識のあるかたが親切で回答する場合をのぞき、このサイトでは期待しないでください。質問１については既に自分で調査済み、ということですね。 -- &new{2006-06-01 (木) 17:31:20};
-回答しないのに、わざわざへこませるようなこと書かなくてもいいんじゃない？ -- &new{2006-06-02 (金) 14:59:56};
-ほっとくよりいいんじゃない？回答者の意気をくじくようなことこそ書かなくても良いんじゃない？~
聞く人は多いけど，答える人は少ないよ。答える人がいなくなると，このWikiはなりたたないよ。~
2. については，変数選択というか予測に重要な変数を選択するのが悪いわけがないでしょう。~
既存の仕組みがないと言うのは障害かもしれないが，やり方はいろいろある。~
そもそも，重回帰の変数選択の仕組み自体を評価しない人もいるわけで（理論的根拠と一致しない結果が得られても困るということ），結果が良ければ理論根拠はなくても良いという立場はあるだろうが。~
総当たり法も，変数が20もあると無理だけど，有望そうな変数セットというのは，そう多くはないのではないか？~
理論根拠に基づき，試行錯誤的に変数選択するのも良いのではないか。それが，そのデータセットで最適なものである保証はないが，それでも十分な結果が得られるのならそれで満足しようよ。 -- &new{2006-06-02 (金) 23:37:12};
**rglがOSXでコンパイルできない [#g1afb318]
>[[高井]] (2006-05-30 (火) 12:30:55)~
~
R 2.3.0のフルセットをMac OSX 10.4.6, Xcode 2.3，PowerMac G5で使っていますが，rgl 0.66のコンパイルに失敗します。エラーメッセージは以下のようなものです。R 2.2.1のときはwarningはでるものの，rgl 0.66はコンパイルできて動きました。~
~
api.cpp: In function ‘void rgl_user2window(int*, int*, double*, double*, double*, double*, int*)’:
api.cpp:608: error: invalid conversion from ‘int*’ to ‘const GLint*’
api.cpp:608: error: initializing argument 6 of ‘GLint gluProject(GLdouble, GLdouble, GLdouble,
const GLdouble*, const GLdouble*, const GLint*, GLdouble*, GLdouble*, GLdouble*)’
api.cpp: In function ‘void rgl_window2user(int*, int*, double*, double*, double*, double*, int*)’:
api.cpp:633: error: invalid conversion from ‘int*’ to ‘const GLint*’
api.cpp:633: error: initializing argument 6 of ‘GLint gluUnProject(GLdouble, GLdouble, GLdouble,
const GLdouble*, const GLdouble*, const GLint*, GLdouble*,
GLdouble*, GLdouble*)’
make: *** [api.o] Error 1
chmod: /Users/takai/Library/R/library/rgl/libs/ppc/*: No such file or directory
ERROR: compilation failed for package 'rgl'
//
-608行，633行を修正すればよいだけでしょう（誰でもソースをいじれるのが，フリーソフトウエアのよいところ） -- &new{2006-05-30 (火) 13:20:03};
-CRANからソースをダウンロードして，ダブルクリックで解凍し，api.cppというファイルをダブルクリックしました。Xcodeというアプリケーションが開いたようですが，そこから先は初心者にはどうにもお手上げです。const GLintなる個所も無いようだし... rglのバイナリーが出るまで待つしかないようです。「誰でもソースをいじれる...」の誰にもには私は含まれないようです。 -- [[高井]] &new{2006-05-30 (火) 16:28:29};
-[[R-develメーリングリストのアーカイブの中:http://tolstoy.newcastle.edu.au/R/devel/06/05/5238.html]]にDuncan Murdockのパッチがあります。 -- [[中澤]] &new{2006-05-30 (火) 17:02:49};
-ちゃんとやれば，誰にでもできるはずなんですけどね。~
レスポンスがないようなのですが，以下は，ターミナルで操作します（大丈夫？）~
中澤さんのコメントを参考に，まず解凍された rgl ディレクトリのあるディレクトリへ移動し，rgl/src/api.cpp を修正します（直接 src ディレクトリへ行ってもいいけど，その場合にはインストール時に上位ディレクトリに戻ること）。
601a602
> GLint viewport[4]; 追加
605a607
> for (int i=0; i < 4; i++) viewport[i] = view[i]; 追加
607c609
< gluProject(point[0],point[1],point[2],model,proj,view, 次行のように変更
---
> gluProject(point[0],point[1],point[2],model,proj,viewport,
624a627
> GLint viewport[4]; 追加
628a632
> for (int i=0; i < 4; i++) viewport[i] = view[i]; 追加
632c636
< gluUnProject(pixel[0],pixel[1],pixel[2],model,proj,view, 次行のように変更
---
> gluUnProject(pixel[0],pixel[1],pixel[2],model,proj,viewport,
その後
$ R CMD INSTALL rgl
でしょうか？？ -- &new{2006-05-30 (火) 18:24:22};
-上記のお二方のコメントに従って再トライしました。api.cppのソースの修正をしてから，R CMD INSTALL rglを試しました。が，WARNING: invalid package 'rgl'とERROR: no packages specifiedが出て失敗しました。Mac OSXではこうはいかないのではと思い，GUIのRを起動-->「パッケージとデータ」-->「パッケージインストーラ」-->「このコンピュータ上のパッケージディレクトリ」から修正したrglを指定しました。こうしたらコンパイルは実行しましたが，結果は
collect2: ld シグナル 6 [Abort trap] で終了させられました
/usr/bin/ld: warning internal error: output_flush(offset = 462632, size = 2352) overlaps with
flushed block(offset = 464972, size = 3068)
calling abort()
make: *** [rgl.so] Error 1
chmod: /Library/Frameworks/R.framework/Versions/2.3/Resources/library/rgl/libs/ppc/*:
No such file or directory
ERROR: compilation failed for package 'rgl'
となってコンパイル失敗です。一歩前進しましたがエラーメッセージの意味が分からず（リンクローダーでのエラーでしょうか？）立ち往生です。-- [[高井]] &new{2006-05-31 (水) 11:32:56};
-R CMD INSTALL rgl を実行したときの pwd が，rgl ディレクトリを含むディレクトリでしたか？ Mac OSX 10.4.6 で確認済（インストール済）なんですけど変ですね。環境が微妙に違うのかな？ -- &new{2006-05-31 (水) 12:01:59};
-ああ，「上位のディレクトリでR CMD INSTALL rgl」とあったので，rglディレクトリの中でR CMD INSTALLやっていました。rglディレクトリの上位のディレクトリでコンパイルしたら走りました。しかし，GUIでRを起動してメニューバーからの指令でコンパイルした時と同様のエラーになりました。環境設定が違う故の差異でしょうかね〜？ -- [[高井]] &new{2006-05-31 (水) 21:32:24};
-/usr/bin/ld の internal errorが出るのは, RかDeveloper Toolsのインストールに問題がありそうですね。R-2.3.1+PMG5 OS10.4.6では中澤先生ご紹介のパッチだけで、GUIから「このコンピュータ上のパッケージディレクトリ」でインストールできました。[[バイナリ:http://www.okada.jp.org/~mokada/R/Mac/rgl_0.66_R_powerpc-apple-darwin8.6.0.tar.gz]]一応おいときますね。 -- [[岡田]] &new{2006-06-09 (金) 18:06:38};
-岡田さん，rglのバイナリーありがとうございます。無事にインストールできました。Rのほうは何度もインストールし直したのですが... 今度はDeveloper Toolsを再インストールしてみます。謝謝 -- [[高井]] &new{2006-06-09 (金) 21:16:02};
-3時間くらいかかって，本家からDevelopper Tools 2.3をダウンロードして再インストールしてみました。念のために，uninstall.plを実行しました。が，やはり同じエラーが出ます。 -- [[高井]] &new{2006-06-11 (日) 08:43:25};
-環境依存の問題であることはほぼ明らかになりましたから、あとはOSの構成ファイル自体がどこかおかしくなっているか、メモリ等のハードウェアエラーぐらいでしょうか。外付けHDDにTigerをクリーンインストールしてみて、それでも再現するようならハードウェアでしょうね。 -- [[岡田]] &new{2006-06-11 (日) 18:22:42};
-岡田さんに励まされ（？），ネジ捲かれ（？）ハードウェアチェックとクリーンインストールをやってみました。Appleのハードウェアチェックディスクでチェックしたところ，問題なし。で，外付けHDDにOSX10.4をクリーンインストールして，アップデートしました。Developper Tools 2.3とR 2.3.1をインストールしてrgl 0.66修正版をインストールしたところ，みごとに（？！）コンパイルが通りました。ということは，今の環境のどこかに問題があるということです。クリーンインストールしたファイルシステムと現ファイルシステムを比べて見ながら検討してみます。皆様，ありがとうございました。 -- [[高井]] &new{2006-06-13 (火) 12:14:21};
-どうやらgccがおかしいんじゃないかと当たりをつけました。/usr/local/bin/gccがインストールされてまして，PATHの先頭でこれを探すようになっていました。このgccは[[HPC:http://hpc.sourceforge.net/]]のgccでした。で，このgcc関係（g77, gfortran）のファイルを消してみました。これでうまくいきました。 -- [[高井]] &new{2006-06-13 (火) 15:28:36};
**文字をx軸に持つデータのプロット [#jd511ae0]
>[[Kamo]] (2006-05-26 (金) 16:27:22)~
~
x軸が文字列,y軸が数値のデータを点と線でプロット(type="b")したいです。~
次のようにすると、強制的に棒グラフ？になってしまいます。~
plot(factor(x),y,type="b")
次のようにするとプロットはできますが、今度は軸の値が数値になってしまいます。~
plot(as.numeric(factor(x)),y,type="b")
軸の値は上のケース、プロットは下のケースのようにするにはどうしたらいいのでしょうか？~
//
-plot(1:3,1:3,xaxt="n");axis(1,1:3,c("a","b","c")) -- [[takahashi]] &new{2006-05-26 (金) 18:12:22};
-x が実際にどうなっているのかわからないのだけど。他の人が追試するなりプログラムするなりが可能なように，端的に表現すべし。で，適当に仮定して，以下のようにすればよい。
x <- letters
y <- rnorm(length(x))
plot(as.numeric(factor(x)), y, type="b", xaxt="n")
axis(1, at=factor(x), labels=x)
老婆心ながら，このような場合には折れ線グラフは不適切と，統計学の教科書に書いてないか？ -- &new{2006-05-26 (金) 18:18:26};
-ああ、これ苦労した記憶があります。結局出来なかったんですよね。一般的には不適切か否かはともかくとして、時系列データで特に間隔が一定じゃなくってイベントを強調したグラフを作りたいとき(例えばバレンタインデーとかクリスマスとかのイベント)、こう言ったグラフを確かに作りたい。でもRだと余計な操作してくれて、棒グラフになった覚えがあります。結局面倒臭くなっちゃってExcelでグラフを作りました。 -- &new{2006-06-06 (火) 11:38:04};
-これだとlabels内文字列でX軸の順番がソートされてしまいますね。次の方がいいかもしれない -- [[nishi]] &new{2006-07-04 (火) 10:34:18};
-y<-c(1,2,3,4);plot(1:length(y),y,xaxt="n");axis(1,1:length(y),c("a","b","c","d")) -- &new{2006-07-04 (火) 10:47:13};
**スクリプトでread.delimを使うとエラーになってしまいます。 [#hadb5af0]
>[[KT]] (2006-05-25 (木) 17:52:50)~
~
コンソールから~
data <- read.delim("data.txt", header=T)
と入力・実行すると正しく実行されるのに、スクリプトファイル("test.R")に同じ内容を記述して、~
source("test.R")
と実行すると、~
以下にエラーeval.with.vis(expr, envir, enclos) :
関数 "raed.delim" を見つけることができませんでした
とエラーになってしまいます。対策をご教授いただければ幸いです。よろしくお願いします。~
//
-x raed.delim / o read.delim -- [[takahashi]] &new{2006-05-25 (木) 17:55:35};
-takahashi様、コメントありがとうございます。ですが、超初心者の私にはどう直せばよいのか分かりません。もう少し詳しく教えていただけましたら幸いです。 -- [[KT]] &new{2006-05-25 (木) 18:21:17};
-あのですね。「つづりが間違えてますよ」と指摘されたわけです。 -- &new{2006-05-25 (木) 18:26:36};
-scriptの方がtypoなんじゃないですか？ x raed / o read -- [[takahashi]] &new{2006-05-25 (木) 18:27:02};
-まさにスペルミスでした。何たるケアレスミス。どうもありがとうございました。 -- [[KT]] &new{2006-05-25 (木) 18:31:42};
**英語版R1.3の入手法は？ [#a767b2b7]
>[[Yabo]] (2006-05-25 (木) 08:32:16)~
~
R1.3を本家からダウンロードしても、インストールすると日本語表示になってしまいます。すでに"RとJava"のところでR1.2について報告があるのと同じく、JGRを動かすと文字化けします。残念ですが日本語が使えなくてもかまわないことにして、JGRを動かしたいです。英語版Rの導入法をどなたかご教示下さい。~
//
-R1.3って(@.@)何？ -- [[okinawa]] &new{2006-05-25 (木) 09:56:23};
-えっと
set LANG=C
set LC_ALL=C
"*****?jgr.exe" #JGRのパス
というバッチファイルを作って，そこから起動すればO.K.です． -- [[g]] &new{2006-05-25 (木) 10:29:28};
-R3.1.0 のことでしょう。未来のお話です(^_^;) というのは冗談で，JGR 1.3 のことでしょうが？？？ -- &new{2006-05-25 (木) 10:30:52};
-32bitかつUTF-8のロカールを持つOS上なら日本語でも動きますが. -- [[なかま]] &new{2006-05-25 (木) 12:16:11};
-ｇさん、ご教示の方法でうまくいきました。ひとりぼっちでやっているので、すぐ答えていただき感激です。R1.3はR2.3no -- [[yabo]] &new{2006-05-25 (木) 17:11:22};
-R1.3はR2.3の誤りでした。すんません。 -- [[yabo]] &new{2006-05-25 (木) 17:16:32};
-gさんのバッチファイルの set LC_ALL=C のCをUTF-8にしても動きました。 -- [[yabo]] &new{2006-05-25 (木) 17:27:07};
**条件が長さが2以上なので [#wb1cb369]
>[[レッド]] (2006-05-23 (火) 04:30:21)~
~
次ののメッセージがでて条件式が効きません．
TSi <- amp$Si
> if (TSi < 8) {
+ TAl <- 8 - TSi
+ } else {
+ TAl <- 0
+ }
Warning message:
条件が長さが2以上なので，最初の一つだけが使われます in: if (TSi < 8) {
超初心者で，Win．Version 2.2.1 ~
//
-if()を調べましょう。~
cond: A length-one logical vector that is not 'NA'. Conditions of
length greater than one are accepted with a warning, but only
the first element is used.
if(x)ではis.logical(x)がTRUE、かつis.na(x)がFALSEかつ、length(x)が1（ただし2以上は一つ目が使われるようです）である必要があります。レッドさんのTSiはこのどれかを満たしていないのです。~
ただ、この条件を満たしているのでしたら、amp$Siがどういったものか提示した方が良いと思います。~
エラーが出る場合は、いきなり応用問題はやめて基礎からはじめた方が近道ですよ -- [[akira]] &new{2006-05-23 (火) 08:10:59};
-TSi がベクトルなんでは？ c(1,3) < 2 は論理ベクトル c(TRUE,FALSE) になります。 -- &new{2006-05-23 (火) 10:07:04};
-このページの上の方の注意事項を読んでから質問・投稿しましょうね。~
amp$Si がベクトルなんでしょうね。だとしたら，TAl <- ifelse(TSi < 8, 8-TSi, 0) だけで，貴方がやりたいことができるでしょう。~
ifelse 関数を調べましょう。 -- &new{2006-05-23 (火) 10:41:59};
-amp$Siはベクトルです． ex.amp$Si <- c(8.123, 7.897, 7.977・・・）.ifelseで上手くいきました．ベクトルについて勉強必要ですね. 親切で的確なコメントありがとうございます． -- [[レッド]] &new{2006-05-23 (火) 14:26:33};
**平均と標準偏差を使ったboxplotは描けませんか？ [#r2b5ff6e]
>[[くに]] (2006-05-22 (月) 18:35:13)~
~
箱ひげ図には平均値と標準偏差を用いたものもありますが、これをRで作図することは可能でしょうか？
//
-平均値・標準偏差を求める関数があり，plot, lines, points, rect 関数 などがありますので，ちゃんと描けますよ。 -- &new{2006-05-22 (月) 18:42:04};
-ありがとうございます。できたら、具体例か参考になるウェブページを示していただけると大変ありがたいのですが。よろしくお願いします。 -- [[くに]] &new{2006-05-22 (月) 21:30:56};
-可能ですが，平均値と標準偏差を使うなら，ストリップチャートに平均と標準偏差の線を重ねがきする方がいいと思います。http://phi.med.gunma-u.ac.jp/swtips/sumfunc.htmlをご覧ください。少なくとも平均と標準偏差の線を打つ方法は参考になると思います。 -- [[中澤]] &new{2006-05-22 (月) 21:46:05};
-可能かどうか聞かれたので「可能ですと」お答えしましたが。。。~
何だってできますから，どこまでやらなきゃならないかを決めることが大事（ダイジ，あるいは，オオゴトと読んでも良いが）なんですね。~
元のデータを記号で書き加えたり，そのほかなんでもあなたのやりたいことを付け加えればよい。 -- &new{2006-05-22 (月) 22:24:13};
bp <- function(x, px, wx=0.25)
{
m <- mean(x)
s <- sd(x)
mx <- max(x)
mn <- min(x)
x1 <- px-wx
y1 <- m-s
x2 <- px+wx
y2 <- m+s
rect(x1, y1, x2, y2)
lines(c(px-wx, px+wx), c(m, m))
arrows(px, m+s, px, mx, angle=90, length=2*wx)
arrows(px, m-s, px, mn, angle=90, length=2*wx)
}
x <- rnorm(100)
plot(c(0,1), c(min(x), max(x)), type="n", xaxt="n", xlab="", bty="n", ylim=c(min(x), max(x)))
bp(x, 0.5, 0.4)
z <- matrix(rnorm(300), nr=100, nc=3)
plot(c(0,3), c(min(z), max(z)), , type="n", xaxt="n", xlab="", bty="n", xlim=c(0,4), ylim=c(min(z),
max(z)), ylab="foo bar baz")
for (i in 1:3) {
bp(z[,i], i, wx=0.12)
}
#ref(fig1.png)
fig 1. example1
#ref(fig2.png)
fig 2. example2
-R の boxplot 関数にご希望のオプションが無い意味も考えましょう。EDA (探索的データ解析) の精神に反する邪道だからです。 -- &new{2006-05-22 (月) 23:20:33};
-ボックスプロット類似のグラフは，探索的データ解析のためだけに描かれるものでもないと思います。確かに，ボックスプロットを提唱したチューキーは文字通り，探索的データ解析のためにボックスプロットを提唱しましたが，類似のデータ表現（というか，それより詳細なグラフ）つまり，データ点と平均値と標準偏差，最大・最小値を群毎に表現するグラフは医学分野では相当昔から使われていました。もっとも，それは "boxplot" ではないわけですが。ストリップチャートよりは洗練されていると思います。デー点の重なりを避ける工夫（丸めるわけですがね）とか，データ点を中心線に対して左右対称に配置するとか。~
棒グラフに標準偏差を示すひげを付けた，「ダメダメ」なグラフもあります。このグラフの発祥は，ボックスプロットより前なのか後なのか。今でもよく見られるようではありますが，だめさ加減は無限大。~
ちなみに，「EDA (探索的データ解析) の精神」とはどのようなことでしょうか。お教え頂けるとありがたく思います。 -- &new{2006-05-22 (月) 23:47:26};
-「平均と標準偏差を使ったboxplotを描くにはどうすればよいですか？」とお尋ねするべきだったかもしれませんね。皆さん、どうもありがとうございました。 -- [[くに]] &new{2006-05-23 (火) 00:23:02};
-「EDA (探索的データ解析) の精神に反する」、というのはデータの分布を見ようとしているのに、外れ値の影響が大きい平均値を使うのはナンセンスだ、という風に解釈してよろしいのでしょうか？ -- [[くに]] &new{2006-05-23 (火) 01:27:09};
**ksvmでのstring kernelの処理 [#o5939ea9]
>[[北陸人？]] (2006-05-18 (木) 22:42:37)~
~
ksvmでのstring kernelを使いたいと考えています。~
データセットを次のようにセットしました。(sktest.txt)~
AA 1
AB -1
BB 1
コマンドは次のようにしました。~
> library(kernlab)
> sktest<-read.table("C:ファイルの場所")
> sktest
V1 V2
1 AA 1
2 AB -1
3 BB 1
> set.seed(50)
> sktest.num<-sample(3,2)
> sktest.train<-sktest[sktest.num,]
> sktest.test<-sktest[-sktest.num,]
> sktest.svm<-ksvm(V2~.,data=sktest,kernel="stringdot",kpar=list (lambda=0.5))
ここで、エラーメッセージが次のようにでます。~
以下にエラーmatch.arg(kernel, c("rbfdot", "polydot", "tanhdot", "vanilladot", :
'arg' は以下の一つでなければなりません： rbfdot, polydot, tanhdot, vanilladot, laplacedot,
besseldot, anovadot, splinedot, matrix
stringdotが選択肢にないようです。
パッケージのマニュアルには、stringdotは選択肢にあり、ソース(stringk.c)という
ものがあることは確認しました。それともデータセットの与え方がまずいのでしょうか？ご教示いただければ幸いです。~
//
-よく知らないのだけど，ヘルプを見る限り，
## S4 method for signature 'list':
ksvm(x, y = NULL, type = NULL, kernel = "stringdot",
:
x a symbolic description of the model to be fit. When not using a formula x
can be a matrix or vector containg the training data
or a kernel matrix of class kernelMatrix of the training data or a list of
character vectors (for use with the string kernel).
なんだから，x はリストでないといけないのでは？ -- &new{2006-05-18 (木) 23:16:13};
-下記のように試してみました。 -- [[北陸人?]] &new{2006-05-19 (金) 04:32:29};
> x<-c("aaa","abb","bbb")
> xl<-list(x)
> class(xl)
[1] "list"
> xl
[[1]]
[1] "aaa" "abb" "bbb"
> s.svm<-ksvm(xl,kernel="stringdot",kpar=list(lambda=0.5))
ここで、エラーメッセージが次のようにでます。
以下にエラーas.double.default(t(x)) : (list)オブジェクトは 'double' に変換できません
ラベルの与え方がまだ理解できていませんが、ラベルyとすると、y=NULLというのは
指定しないということでしょうか？
-なかなか，コメントがつきませんね。「わたしは，その方面は興味がないし，ひつようなら自分で調べてね」ということでしょうか。 -- &new{2006-05-20 (土) 22:12:09};
-そうですね。string kernelが正確に使えるSVMは少ないですね。 -- [[北陸人]] &new{2006-05-21 (日) 15:24:34};
**ksvmでの入力データ(listから行列の変換) [#q4e1ede1]
>[[北陸人]] (2006-05-17 (水) 13:40:12)~
~
ksvmを使いたいと思っていますが、~
入力するデータの変換についてお伺いします。
testdata <- read.table("C:ファイルの場所") #データ読み込み
testdatam <- data.matrix(testdata) #list(data.frame)から行列に変換
set.seed(50) #乱数の設定
testdatam.num <- sample(14965, 9000) #訓練データの数の設定
testdatam.train <- testdatam[testdatam.num,] #訓練データ
testdatam.test <- testdatam[-testdatam.num,] #テストデータの設定
testdatam.svm <- ksvm(type~., data=testdatam.train, kernel="rbfdot", kpar=list(sigma=0.001))
#SVMでの訓練データの学習
ここで
以下にエラーmodel.frame.default(data = list(V1 = c(33,19,13,24, 21,23,: オブジェクトが行列ではありません
とエラーメッセージがでます。~
確認してみますと
> class(testdatam.train)
[1] "matrix"
> is.list(testdatam.train)
[1] FALSE
> is.matrix(testdatam.train)
[1] TRUE
とでます。~
データは
7 11 4 1 2 3 0 2 （略） 9 4 3 8 5 3 1 2 0 3 1 -1
というスペース区切りの数値データです。~
どういう操作が足らないのかご教示いただけませんでしょうか？~
//
-問題を再現できる必要最小限のデータと共に，再度質問される方がよいでしょう。データも，どういう風にファイルになっているかわかりませんよね。スペース区切りで長い一行として用意されているのでしょうか？ -- &new{2006-05-17 (水) 17:32:33};
-そのデータは１行しかなくて、行列型に旨く変換されなかったのでしょう。ncolsやnrowsでtestdatam.trainがどうなったか覗いてみましょう。もちろん再現しやすい情報を載せるのが吉ですが。 -- [[なかま]] &new{2006-05-17 (水) 19:26:20};
-説明不足で申し訳ありません。 -- [[北陸人]] &new{2006-05-17 (水) 20:26:51};
データセット(test1.txt)データ数を少なくしました。
0 1 -1
2 3 -1
1 5 +1
コマンドは次のようにしました。
> test1<-read.table("C:ファイルの場所",header=F)
> test1
V1 V2 V3
1 0 1 -1
2 2 3 -1
3 1 5 1
> test1m<-data.matrix(test1)
> class(test1m)
[1] "matrix"
> set.seed(50)
> test1m.num<-sample(3,2)
> test1m.train<-test1m[test1m.num,]
> test1m.test<-test1m[-test1m.num,]
> test1m.svm<-ksvm(type~.,data=test1m.train,kernel="rbfdot",kpar=list (sigma=0.001))
以下にエラーmodel.frame.default(data = list(V1 = c(1, 0), V2 = c(5, 1), V3 = c(1,
:オブジェクトが行列ではありません
> dim(test1m)
[1] 3 3
> class(test1m.train)
[1] "matrix"
> dim(test1m.train)
[1] 2 3
上記のようになりました。ksvmにかけたときのオブジェクト(test1m.train?)が行列ではないというメッセージでした。dim(test1m.train)で2*3の行列だというメッセージがでていますので、行列型になっていると考えました。初歩的なミスだと思いますがどの点が間違っているでしょうか？長文になりご迷惑をお掛けして申し訳ありません。
-第一引数のtypeがデータの列名に見つからないからじゃないですか？ V3~.（reponse vectorはV3ですよね？）に変更してみては。 -- [[takahashi]] &new{2006-05-17 (水) 21:03:56};
-library(kernlab)であることを書いて頂いた方が早く返事が来たと思います。エラーメッセージが変なのかもしれません。data=で指定するのはデータフレームでなくてはいけません。example(ksvm)で出てくるspamの例でも，data=で指定されるのはデータフレームです。typeという変数はspamというデータフレームの中の変数なので，あなたの例では使えません。ちなみに，ksvm(V3~.,data=test1,kernel="rbfdot",kpar=list(sigma=0.001))とするとエラーは出ません。 -- &new{2006-05-17 (水) 21:15:04};
-どんぴしゃり解決しました。アホでした。どうもありがとうございます。m(_ _)m -- [[北陸人]] &new{2006-05-17 (水) 21:53:04};
**MacでのRコマンダーの使用について [#y76f2c61]
>[[どいま]] (2006-05-16 (火) 14:42:13)~
~
MacOS XでのRコマンダーのインストールが出来ず困っています。~
Rコマンダーをロードしようとすると、以下のようなメッセージが表示されます。~
Rのバージョンは2.3.0です。
要求されたパッケージ tcltk をロード中です
Loading Tcl/Tk interface ... 以下にエラーdyn.load(x, as.logical(local), as.logical(now)) :
共有ライブラリ '/Library/Frameworks/R.framework/Resources/library/tcltk/libs/i386/tcltk.so'
を読み込めません
dlopen(/Library/Frameworks/R.framework/Resources/library/tcltk/libs/i386/tcltk.so, 10):
Library not loaded: /usr/X11R6/lib/libX11.6.dylib
Referenced from: /Library/Frameworks/R.framework/Resources/library/tcltk/libs/i386/tcltk.so
Reason: image not found
追加情報: Warning message:
ライブラリ ‘/Users/doima/Library/R/library’ はパッケージを含んでいません in: library()
エラー：.onLoad は 'tcltk' のための 'loadNamespace' に失敗しました
エラー：パッケージ 'tcltk' をロードできませんでした
//
-Macを使っていないので予想ですが、Tcl/Tkがinstallされていないのでは？R Commanderには"For Mac OS/X use, it's necessary to run R under X-Windows, and to install Tcl/Tk."と書いてありました。 -- [[Akira]] &new{2006-05-16 (火) 14:52:08};
-R Commanderには"X11"が必須ですが、これは標準ではインストールされていません。OSのディスクからX11だけを追加インストールしてみてください。 -- &new{2006-05-17 (水) 15:35:00};
**Rで使えるデータ一覧 [#y1977666]
>[[究極超具あーる]] (2006-05-12 (金) 12:57:54)~
~
Rで例えば、~
library(stats)
data(morley)
とすると、データmorleyが使えるようになりますよね。このような既にRで利用可能なデータの一覧と、その内容を知りたいのですが、どこかにまとめられたページはありますか？もしくは調べ方がわかればいいのですが、教えてください。~
//
-データセットmorleyはstatsパッケージではなくbaseパッケージの中にあります。[[base(基本)パッケージ中のオブジェクト一覧 :http://www.okada.jp.org/RWiki/index.php?base%28%B4%F0%CB%DC%29%A5%D1%A5%C3%A5%B1%A1%BC%A5%B8%C3%E6%A4%CE%A5%AA%A5%D6%A5%B8%A5%A7%A5%AF%A5%C8%B0%EC%CD%F7]]を見ましょう。RjpWikiの中だと『オブジェクト一覧』で検索すると、一次情報には行き着くと思います。 -- &new{2006-05-12 (金) 13:15:12};
-morley などは，data(morley) などとしなくても，使えます。base パッケージや stats パッケージなどは library 関数でロードしなくても最初からロードされている（あるいは，必要になった時点で即ロードされる）ので，それらに含まれるデータも，data 関数でロードしなくても使えます。 -- &new{2006-05-12 (金) 13:21:56};
-インストール済みの全パッケージの全データ一覧を知りたいです。いろいろ自分でも調べた結果、
data(package="パッケージ名")
とするとパッケージ中のデータ一覧が得られることがわかりました。
ただ、この中のItem列を抜き出そうとすると、以下のようなエラーが出ます。
datas = data(package="accuracy")[["results"]];
datas[["Item"]]
以下にエラーdatas[["Item"]] : 添え字が許される範囲外です
datasのデータ構造が良くわかりません。Item列を抜き出すにはどうしたらいいでしょうか？ -- [[究極超具あーる]] &new{2006-05-12 (金) 14:03:22};
-datas[,"Item"]とするか、data(package="パッケージ名")[["results"]][,"Item"]なら一発でしょう。"パッケージ名"の部分にインストールしている全てのパッケージを入れるのは、パッケージの読み込みの際に使われる
local({pkg <- select.list(sort(.packages(all.available = TRUE)))
+ if(nchar(pkg)) library(pkg, character.only=TRUE)})
の部分を書き換えると対応できるのではないかと。 -- &new{2006-05-12 (金) 14:20:37};
-単純ミスでした。ありがとうございました。 -- [[究極超具あーる]] &new{2006-05-12 (金) 14:34:20};
-次のようにすると、究極超具あーるさんの望みは叶うと思います。
x<-sort(.packages(all.available = TRUE))
lspkg <- function(pkgname){
data(package=pkgname)[["results"]][,"Item"]}
sapply(x,lspkg)
データを含んでないパッケージまで一々character(0)と返すのはご愛嬌。 -- &new{2006-05-12 (金) 14:51:45};
**filterの使い方 [#w410852a]
>[[ワース]] (2006-05-09 (火) 19:59:58)~
~
下のフィルタをfilterで実装したいと思っています。~
h = 0.01/(1-z^-1)
以下のようにスクリプトを書きました。~
x = ts(rep(1, 1000), freq=100, start=0)
x = lag(x,k=-100)
y = filter(filter(x, 0.01, method="con", sides = 1), -1, method="rec")
plot(y)
フィルタhがうまく動けば、yは(t,y) = (0,1), (10,9)を通る直線となるはずなのですが、上記の例ではそうなりません。何が間違っているのでしょうか？~
アドバイスをいただけると幸いです。よろしくお願いいたします。~
//
-それぞれの行ごとに，x, y がどのようになっているか，確認したのでしょうか？？~
x は全部1なんですけどいいのかなぁ。 x = ts(1:1000, freq=100, start=0) じゃないのかなぁ。もっとも，そうであっても，直線にはなるけど(0,1), (10,9)は通らないしなぁ。 -- &new{2006-05-10 (水) 17:58:27};
-典型的なまずい質問の例ですね。自分のしたいことを人にもわかるようにまず説明できなければ、良い回答は得られません。 自分にはわかっていても人にはほとんどわからない（例えば filter h って具体的に何をするもの？）ことが多いのです。それでも、何をしたいのか理解させられれば、ヒントが得られる可能性は高くなります。わたしも、少し考えましたが、そもそも何をしたいのかわからないので、匙を投げました。そもそも、内側の filter は単に時系列全体を 0.01 倍しているだけのように見えますが？ -- &new{2006-05-10 (水) 23:11:41};
-まずいですか？「フィルタをfilterでどう実装するか」という問題なのですから、フィルタhの式が明確になっていれば十分かと思ったのですが…。hは積分器1/sをサンプリング周期0.01で離散化したものです。filterによる実装対象として、簡単な例を挙げたつもりです。積分器の使い方に関してはよろしいかと思います。xが全部1なのは単なるステップ入力だからです。問題はfilterの使ったフィルタの実装方法なので、入力xは別に何でもいいのです。入出力に関係なく、フィルタの式をfilterでどう記述するかという話ですから。が、やはりなるべく簡単な例を挙げたつもりです。この情報が役に立つことを願います。 -- [[ワース]] &new{2006-05-11 (木) 01:17:05};
-あなたの説明はフィルターやら積分器とやらをすでにある程度知っているひとにしかわからない説明です。このサイトの閲覧者（特に常連回答者）にはそうした分野の知識はおそらく皆無でしょう。何をしたいか、とは ?filter で得られる説明中の y[i]=f[1]*x{i]+... という類の元の時系列を、どういう時系列に変換するのかという形式的関係式です。R (というより時系列理論)ではフィルターとはまさにこうした入力・出力時系列間の関係式にすぎません。それが何を表しているのかなどはある意味でどうでも良いことなのです。 -- &new{2006-05-11 (木) 05:23:26};
-お勉強の結果（1）~
filter の引数 method が "convolution" のときは，「（重みつき）移動平均（移動和）」を求める。~
二番目の引数は項数ごとの重み，つまり，n 項移動平均なら rep(1/n, n)~
ただし，本当は移動和なので，重みは等しくなくても良いし和が1にならなくてもかまわない。~
使用例 a3 と a4 の定義と結果を比較のこと。重みが「逆順」であることもわかる。~
sides=1 と sides=2 の違いは使用例 a1 と a2 の比較で明白。
> x <- ts(c(1,3,4,2,1,4,6,2), freq=2, start=0)
> a1 <- filter(x, rep(1/3, 3), method="convolution", sides=1)
> a2 <- filter(x, rep(1/3, 3), method="convolution", sides=2)
> a3 <- filter(x, c(2,0.5,1), method="convolution", sides=2)
> a4 <- rep(NA, 8)
> for (i in 2:7) a4[i] <- 1*x[i-1]+0.5*x[i]+2*x[i+1]
> cbind(x, a1, a2, a3, a4)
Time Series:
Start = c(0, 1)
End = c(3, 2)
Frequency = 2
x a1 a2 a3 a4
0.0 1 NA NA NA NA
0.5 3 NA 2.666667 10.5 10.5
1.0 4 2.666667 3.000000 9.0 9.0
1.5 2 3.000000 2.333333 7.0 7.0
2.0 1 2.333333 2.333333 10.5 10.5
2.5 4 2.333333 3.666667 15.0 15.0
3.0 6 3.666667 4.000000 11.0 11.0
3.5 2 4.000000 NA NA NA
ということですか。 -- &new{2006-05-11 (木) 11:15:05};
-お勉強の結果（2）~
filter の引数 method が "recursive" のときは，ある時点の値はそれ以前の時点のフィルター出力に重みを掛けたものとの和（重みは逆順で）。~
重みの長さが1で値も1なら，cumsum と同じ結果になる（当然だが）。a1 と a2 を比較。~
重みの長さが1で値が2なら，a4[1] = x[1], a4[2] = x[2]+2*a4[1], a4[3] = x[3]+2*a4[2] この意味で「再帰的」~
重みの長さが2以上で，値が違っても同様。a5[4] = x[4]+1*a5[3]+2*a5[2]
> x <- ts(c(1,3,4,2,1,4,6,2), freq=2, start=0)
> a1 <- filter(x, 1, method="recursive")
> a2 <- cumsum(x)
> a3 <- filter(x, -1, method="recursive")
> a4 <- filter(x, 2, method="recursive")
> a5 <- filter(x, 1:2, method="recursive")
> cbind(x, a1, a2, a3, a4, a5)
Time Series:
Start = c(0, 1)
End = c(3, 2)
Frequency = 2
x a1 a2 a3 a4 a5
0.0 1 1 1 1 1 1
0.5 3 4 4 2 5 4
1.0 4 8 8 2 14 10
1.5 2 10 10 0 30 20
2.0 1 11 11 1 61 41
2.5 4 15 15 3 126 85
3.0 6 21 21 3 258 173
3.5 2 23 23 -1 518 345
ということですか。~
どうも，書かれたプログラムはあなたがやろうとしたこととは別なことをやっている。どこがどう違うかは，あなたにならわかるはず。 -- &new{2006-05-11 (木) 11:53:28};
-お勉強の結果（3）
> x1 <- ts(rep(1, 1000), freq=100, start=0)
> x2 <- lag(x1,k=-100)
> y1 <- filter(x2, 0.01, method="con", sides = 1)
> y2 <- filter(y1, -1, method="rec")
> cbind(x1, x2, y1, y2)
Time Series:
Start = c(0, 1)
End = c(10, 100)
Frequency = 100
x1 x2 y1 y2
0.00 1 NA NA NA
すっぱり省略
0.99 1 NA NA NA
1.00 1 1 0.01 0.01
1.01 1 1 0.01 0.00
1.02 1 1 0.01 0.01
1.03 1 1 0.01 0.00
1.04 1 1 0.01 0.01
以下略
プログラムに書かれたとおりの動きをしているようだ。 -- &new{2006-05-11 (木) 12:01:24};
-お勉強の結果（4）~
あの〜。もとの発言にあった z^-1 というのは，「遅延素子」なんでしょうか。で，-1, method="rec" でそれを表そうとしたのでしょうか？つまり，y[i] = y[i-1] ??~
そうだとすると，関数の定義上それはあなたの意図とは違ったことをやってくれますね。-- &new{2006-05-11 (木) 12:08:43};
-お勉強の結果（5）~
遅延素子かどうかはともかくとして，あなたの最初の質問で，「yは(t,y) = (0,1), (10,9)を通る直線となるはず」と書いてありましたが，(0,1) ではなくて，(1,0) ではないのですか？~
もしそうだとすると，
> y <- filter(filter(x, 0.01, method="con"), 1, method="rec")
でいいのでは。つまり，外側の filter の第2引数が -1 ではなくて 1 って，あっけない。 -- &new{2006-05-11 (木) 12:22:22};
#ref(filter.png)
-いろいろありがとうございました。おっしゃるとおり、ｚ^-1は遅延素子で、(0,1)ではなく(1,0)の間違いです。おかげさまで良くわかりました。 -- [[ワース]] &new{2006-05-11 (木) 15:44:04};
**ksvmに関連してのcsvファイルの読み書き [#c846105c]
>[[こねこ]] (2006-04-29 (土) 17:36:31)~
~
ksvmに関連してcsvファイルの読み書を行おうとしているのですが、詰ってしまいました。~
環境はWindowsXP + R(2.3.0)です。~
アドバイスをいただけますと幸いです。~
~
まずlibrary(kernlab)の後、下記コマンドを実行すると正常に処理が行われます。~
~
> x<-as.matrix(iris[51:150,3:4])
> y<-as.matrix(iris[51:150,5])
> iris1 <- ksvm(x,y,kernel="anovadot",kpar=list(sigma=0.5))
> table(y,predict(iris1,x))
~
ところが下記のように一度データをcsvファイルに保存してから、読み込みなおして実行するとエラーが発生してしまいます。~
~
> x<-as.matrix(iris[51:150,3:4])
> y<-as.matrix(iris[51:150,5])
> write.table(x, file = "svmsample_x.csv", sep = ",", col.names = NA)
> write.table(y, file = "svmsample_y.csv", sep = ",", col.names = NA)
> rm(list=ls(all=TRUE))
> x<-read.table("svmsample_x.csv", header = TRUE, sep = ",", row.names=1)
> y<-read.table("svmsample_y.csv", header = TRUE, sep = ",", row.names=1)
> iris1 <- ksvm(x,y,kernel="anovadot",kpar=list(sigma=0.5))
以下にエラーksvm(x, y, kernel = "anovadot", kpar = list(sigma = 0.5)) :
no direct or inherited method for function 'ksvm' for this call
> table(y,predict(iris1,x))
エラー：オブジェクト "iris1" は存在しません
以下にエラーpredict(iris1, x) : unable to find the argument 'object' in
selecting a method for function 'predict'
~
読み書きでデータ型が変化したわけでもないようです。~
申し訳ありませんがご教示宜しくお願いします。~
//
-自己解決しました。失礼しました。 >x<-as.matrix(read.table("svmsample_x.csv", header = TRUE, sep = ",", row.names=1)) -- [[こねこ]] &new{2006-04-29 (土) 20:45:44};
-read.table 関数はデータフレームを返す一方、ksvm は第一引数が formula でなければベクトルを必要とするということでしょうか。 -- &new{2006-04-29 (土) 20:58:04};
**成長曲線の当てはめ [#ja28eabb]
>[[たつ]] (2006-04-29 (土) 15:09:09)~
~
y(売上高）とx(キャンペーン額)の関係を成長曲線に当てはめたいと思っています。~
以上の目的に適当なパッケージはありますでしょうか？~
その際に売上高の上限と下限を指定できたりすると好都合なのですが････~
~
アドバイスよろしくお願いします。~
//
-glm関数で出来るかもと思い色々と試行錯誤もしましたが、自分の知識では無理でした。 -- [[たつ]] &new{2006-04-29 (土) 15:10:45};
-まず RSiteSearch("growth curve") を実行し、ヒットしたたくさんの記事を（うんざりしなが）らチェックしてください。 -- &new{2006-04-29 (土) 17:50:11};
-ヒットしました･･･覚束ない英語力でチェックしていたところ、上限と下限を指定したい場合は比較的メジャーなパッケージであるnlsで目的を果たせる事がわかりました。
お陰さまで解決できたのですが最後に質問をさせてください。
上限と下限を指定しない場合はglmでも可能なのでしょうか？-- [[たつ]] &new{2006-04-29 (土) 19:55:27};
-glm は非正規誤差にたいする線形回帰手法ですから、成長モデル等の非線形モデルは扱えないのでは。成長モデルにもいろいろあるようですが、R の基本パッケージにはそうした特殊非線型回帰専用の関数（nls を用いるのですが、適切な初期パラメータを自動推定してくれる SS(SelfStart) モデルと呼ばれるものです）が用意されています。この初期推定値はしばしば最終推定値として使えるほど優れています。非線型回帰分析でいちばん困るのは適切な初期推定値を得ることです（しばしば試行錯誤が必要となります）。例えば以下のヘルプを見てください。なお、nls はパッケージでは無く、基本パッケージ stats中の一関数です。 -- &new{2006-04-29 (土) 20:45:46};
~
SSlogis
SSgompertz(stats) Gompertz Growth Model
SSweibull(stats) Weibull growth curve model
-素人の私にはSSlogisが使い勝手が良さそうです。これから詳しい使用法を学習しようと思います。ありがとうございました。 -- [[たつ]] &new{2006-04-29 (土) 22:30:54};
**最小絶対偏差回帰 [#r1f44837]
>[[たろう]] (2006-04-26 (水) 20:19:04)~
~
R初心者ですが、Splusをしばし使用しておりました。~
Splusでは、最小絶対偏差回帰は~
l1fitですが、Ｒで対応する関数を見出すことできずにおります。~
ご存知のかた、ご教示宜しくお願いします。 ~
//
- r-help ML の過去記事を検索したらつぎのような(二件の別個の)やり取りが。&new{2006-04-26 (水) 20:49:36};
> is there any package or method which enables to compute a linear
> regression with the leas absolute value fit-criterion?
Try rq() from the quantreg package.
~
You can also look at package nlrq, if you intend to do non linear LAD regression.
-より頑健な回帰手法（おそらく単なる l1fit よりベター？）が望みならパッケージ MASS （パッケージバンドル VR の一部）にいくつもあります。 -- &new{2006-04-26 (水) 21:18:13};
-ご教示ありがとうございます。おかげさまで無事できました。 -- [[たろう]] &new{2006-04-27 (木) 10:39:52};
-ちなみに二件の記事は Rsitesearch でそれぞれキーワード "least abs", "l1fit" で見つかったものです（上級者も何でも知っているわけではありません、調べ方が上手なだけです）。組み込み命令 RSiteSearch("l1fit") を試してみてください。非線形版があることは私も今回始めて知りました。 -- &new{2006-04-27 (木) 16:06:59};
-なるほど。RSiteSearch、初めてしり早速、確認いたしました。今後つかいたいと思います。 -- [[たろう]] &new{2006-04-27 (木) 16:33:44};
-「より頑健な回帰手法」の実行方法についても、お教え頂きありがとうございました。身の回りに問題関心の近いRユーザが少ないこともあり、レスを頂いたことがとても励みにもなりました。 -- [[たろう]] &new{2006-04-27 (木) 16:37:59};
**部分一致検索によるtableからの行の抽出 [#i01f05f8]
>[[hogehoge]] (2006-04-25 (火) 23:21:19)~
~
R初心者です。~
以下のようなtableからNameに関する部分一致検索で行を抽出したいのですが、どのように書けば良いのでしょうか？
>in1
Name Score
AAA 25
AAB 50
ABB 75
BBB 100
完全一致の場合は以下のように書けるのを見てなるほどと思い
>in2 <- subset(in1, Name=="AAB")
>in2
Name Score
2 AAB 50
次のように書いたのですが
>in2 <- subset(in1, match("AA", in4$Name))
>in2 <- sunbset(in1, match("AA", Name))
etc...
'subset' は論理値として評価されなければなりません。
と怒られてしまいます。またhelp(match)の内容も良く分かりませんでした。RはWin用2.3.0を使用しています。~
すみませんが、ご教示宜しくお願いします。~
//
-以下に示すように，
in1[grep('AA', as.character(in1$Name)),]
でできますが、これだと、"BAA" があれば、それも摘出してしまいます。もし"AA" で始まる物だけ、ということなら、
in1[grep('AA', substr(as.character(in1$Name),1,2)), ]
でできます。 -- [[さら]] &new{2006-04-25 (火) 23:50:27};
-"AA" で始まる物だけということなら，行頭を表す正規表現 ^ を使って
in1[grep('^AA', as.character(in1$Name)),]
で良いと思います。 -- [[青木繁伸]] &new{2006-04-26 (水) 00:19:45};
-ちなみに「'subset' は論理値として評価されなければなりません。」とは、subsetの引数がTRUE or FALSEでないとだめということ。match()の返り値はnumericなので、errorが出たのだと思います。 -- &new{2006-04-26 (水) 08:25:06};
-さらさん、青木先生、どうもありがとうございました。as.characterを作用せるのがポイントですね。 -- [[hogehoge]] &new{2006-04-26 (水) 10:44:56};
-行頭を表す正規表現 ^ 。 そんな便利な物があったんですね。いろいろ使えそうです。 -- [[さら]] &new{2006-04-27 (木) 01:22:50};
-Rjpwikiにも正規表現の詳しい説明がある -- &new{2006-04-27 (木) 09:40:24};
**常微分方程式系のバラメータ推定をするには？ [#le048e5a]
>[[使用1ヶ月目]] (2006-04-19 (水) 11:44:48)~
~
常微分方程式系のバラメータ推定を考えています。常微分方程式系のバラメータ推定をするには、おそらくodesolveの中の関数rk4とnlmを組み合わせる必要があると思います。まだ、Rの操作に慣れていないので、うまく記述することができません。どこかにコーディング例がありましたらご紹介ください。よろしくお願いします。~
//
-それぞれの関数の help に example がありますでしょ？そのなかに参考文献も書かれてますでしょ？ -- &new{2006-04-19 (水) 23:18:19};
-いくら教えてる立場だとは言え、全く関係ない人が見ても不愉快になるような不用意な言葉遣いは控えた方が… -- &new{2006-05-10 (水) 09:23:41};
**データファイルを読み込んでcor.testをする時のエラー [#f633e39a]
>[[tensho]] (2006-04-17 (月) 23:26:15)~
~
基本的な質問で申し訳ありません。~
データファイルを読みこんでcor.testを行うと、「長いオブジェクトの長さが短いオブジェクトの長さの倍数になっていません」というエラーが出てしまいます。対処法を宜しくお願いします。~
//
-文字とおりの意味では？これ以上の回答を期待するなら、読み込んだ（はずの）データをこの場に紹介しない限りコメントのしようがないと思います。 -- &new{2006-04-17 (月) 23:30:09};
-データ入力は以下の通りです。
> data <- read.delim("test.txt")
> data
A B C D E
1 3.0 4.0 8.0 1.0 7.0
2 0.3 0.6 0.1 0.7 0.2
> cor.test(data[1,],data[2,],method="p")
これだと上手くいくのですが。
> x <- c(3,4,8,1,7)
> y <- c(0.3,0.6,0.1,0.7,0.2)
> cor.test(x,y,mehtod="p")
-データフレームの意味を誤解していませんか。行がケース、列は変数です。data[1,] は一つのケース（例えば個人）にたいする A,B,C,D,E ５つの変数の値（つまりこれもデータフレームで、ベクトルでは無い）という意味になります。期待する結果は x <- t(data); cor.test(x[,1],x[,2], method="p") で得られます。データフレームは転置すると行列に強制変換されます。 -- &new{2006-04-18 (火) 07:16:56};
-ご親切に教えていただき、どうもありがとうございました。 -- [[tensho]] &new{2006-04-18 (火) 13:27:02};
**GLMでの決定係数 [#n25e8784]
>[[ごう]] (2006-04-11 (火) 06:30:08)~
~
GLMの基本的な部分かもしれないのですが、情報をいただけませんでしょうか。~
正規分布を前提とした普通の重回帰モデルでは決定係数が算出できますが、GLMでのポアソン回帰などでは普通の決定係数は算出できませんよね。~
近似した値は「誰々（人名）」の決定係数といった雰囲気で算出できるのかもしれないのですが、これはロジスティック回帰を見る限りではどうも低く算出されているような気がします。~
正規分布重回帰モデルでもサンプルサイズが大きければ有意なモデルとして評価されますが、R^2=0.1ではあまり意味がありません。GLMではこのあたりをどのように解釈すれば良いのでしょうか。Rで step(model)とすればAICで評価されたモデルが手に入りますが。特にGLMMではもっとめんどうで、モデルが成り立つかどうかで判断するというような記載を読んだ気がします。GLM、GLMMではサンプルサイズ至上主義なのでしょうか？~
重複投稿ご容赦願います。他のサイトでも質問したのですが返答頂けなかったもので、こちらで再度質問させて頂きました。~
//
-＞ 近似した値は「誰々（人名）」の決定係数といった雰囲気で算出できるのかもしれないのですが、~
できるのかもしれないのですがとは，できるんですか，できないんですか？~
誰々の決定係数というのは，あるんですか，ないんですか？~
あるのなら，具体的に何なんですか？~
＞ これはロジスティック回帰を見る限りではどうも低く算出されているような気がします。~
気がするだけなんですか，実際に低いんですか？~
その判断は，あなたがしたものですか，そのように書いてあったんですか（どこに）？~
＞ モデルが成り立つかどうかで判断するというような記載を読んだ気がします~
気がするだけなのですか？ -- &new{2006-04-11 (火) 08:55:27};
-不明瞭な質問でした。決定係数の例：Nagelkerke R^2, Cox&Snell R^2 決定係数が低く見積もられているとの記載 ロジスティック回帰分析におけるモデルの適合度指標に関する考察と提案 2004 東京情報大学研究論集 8(1) 9〜14 GLMMでのR^2については、うろ覚えで勘違いかもしれません。URLのような記載を読んで、難しいものだと思ったところまでは思い出せたのですが。https://stat.ethz.ch/pipermail/r-help/2006-March/090658.html http://finzi.psych.upenn.edu/R/Rhelp02a/archive/48023.html -- [[ごう]] &new{2006-04-11 (火) 13:24:40};
**bartlett.test関数の中身を見るには？ [#j47456dd]
>[[使用1ヶ月目]] (2006-04-10 (月) 17:17:00)~
~
通常、関数の中身を見るには、関数名をそのまま入力すればよいのですが、bartlett.test関数は、この方法で中身を見ることができません。S-plusには、bartlett検定の関数がないので、Rのbartlett.test関数を移植しようと思っています。bartlett.test関数の中身を見る方法をご存知の方がいらっしゃいましたら、よろしくお願いします。~
//
-このスレッドの下の方に同じような質問がありますよ。 -- &new{2006-04-10 (月) 18:07:21};
**尤度のプロットについて [#l09ecac0]
>[[にゃあ]] (2006-04-08 (土) 14:26:34)~
~
例えば二項分布binomの確率分布はdbinomとして計算できますが、同様な手法で簡単に尤度を計算させる命令というのはあるのでしょうか?~
また尤度のプロットというのも可能なんでしょうか?~
//
-尤度の定義通りに計算すればよろしいのでは。ただし、普通は対数尤度を計算する方が賢明です -- [[QDU]] &new{2006-04-08 (土) 19:41:47};
N <- 20
X <- rbinom(100, N, 0.55) # データ
loglikelihood <- function(p) {
log.prob <- dbinom(0:N, N, prob=p, log=TRUE) # 対数二項確率
LL <- 0
for (x in X) LL <- LL + log.prob[x+1] # x=n に対する確率は log.prob[n+1] であることに注意
# もしくは
# for (i in 0:20)
# # LL <- LL + sum(X[X=i])*log.prob[i+1] #ではなくて(冷汗)
# LL <- LL + sum(X==i)*log.prob[i+1] # ですね（失礼）
return(LL)
}
-後で消去しますが，というか，もし指摘が正しいなら，記事修正後，この記事もあなたが消去してくださる方が良いのですが，「もしくは」の前と後の定義では結果が違ってませんか？ -- &new{2006-04-08 (土) 21:32:18};
-ベルヌーイ尤度にベータ分布を掛ける作業を視覚化してみたかったんです。~
関数を作ってみたんですが、以下の様になりました。~
mybayes1 <- function(n, prob){
theta <- seq(0, 1, length=100)
distribution <- dbeta(theta, 1, 1)
plot(theta, distribution, ylim=c(0, 1), type="l")
count <- 0
for(i in 1:n){
.x <- runif(1)
if (.x <= prob) cointoss <- 1
else cointoss <- 0
if (cointoss == 1) count <- count + 1
likelihood <- dbinom(count, n, theta)
posterior <- likelihood*distribution/sum(likelihood*distribution)
plot(theta, posterior, ylim=c(0, 1), type="l")
distribution <- posterior
}
}
かなり雑なんですが、見よう見まねでやってみました。~
しかし、実行速度にかなり問題があります。~
どなたかいいアイディアがありますでしょうか?-- [[にゃあ]] &new{2006-04-08 (土) 22:01:55};
-関数を書いて，それに対する意見を求めて，回答があったとしてそのアイデアをインストールしてということをやっているよりも，今のプログラムでやった方が実行時間は遙かに短いのですから。実行時間が1年間が1日になるのはうれしいだろうが，365秒が1秒になってもあんまりうれしくはないでしょう。歯でも磨いていれば，あっという間に365秒は過ぎていきますよ。 -- &new{2006-04-08 (土) 22:35:12};
-なるほど。言われてみればその通りかもしれません。 -- [[にゃあ]] &new{2006-04-08 (土) 22:58:44};
-毎回プロットしていれば時間はかかるし、途中の結果も見られない（par(ask=TRUE) としていれば毎回 CR おさねばならぬ)し。 -- &new{2006-04-09 (日) 18:13:36};
-そうですね。ただし、根本的な考え方の間違いに気付いたので、このプログラムはオクラとします。これじゃあ伝統的な標本理論である事に気付きました。 -- [[にゃあ]] &new{2006-04-09 (日) 20:11:32};
-関連して気になったこと。非負整数値データの度数分布ベクトルを作るスマートな方法は何でしょう。table 関数は欠損した値をパスしてしまいますが、欠損した値に対応する度数はゼロとして欲しい場合です。 -- [[QDU]] &new{2006-04-10 (月) 21:58:06};
> x <- sample(0:10, 20, replace=TRUE)
> x
[1] 0 2 0 10 0 6 0 8 10 3 10 3 5 6 2 4 3 1 4 6
> table(x) # 欠損した値 7,9 は無視される
x
0 1 2 3 4 5 6 8 10
4 1 2 3 2 1 3 1 3
> sapply(0:max(x), function(i) sum(x==i)) # 例えばこうすればよいがもっと直接的な方法があってもよさそうな？
[1] 4 1 2 3 2 1 3 0 1 0 3
-中級Ｑ＆Ａに同じような質問がありますよ -- &new{2006-04-10 (月) 22:36:57};
-参照先を示してあげるか，再度例示してあげると，おんぶにだっこですね。~
> x <- sample(0:10, 20, replace=TRUE)
> table(x)
x
0 2 3 4 5 6 7 8 9 10
1 3 2 1 3 3 3 1 2 1
> x2 <- factor(x, levels=0:10)
> table(x2)
x2
0 1 2 3 4 5 6 7 8 9 10
1 0 3 2 1 3 3 3 1 2 1
かな？ -- &new{2006-04-10 (月) 22:45:32};
-なるほど table(factor(x, levels=0:max(x)) ですか。ありがとうございます。しかし（負け惜しみではけっしてないですが）一発で計算できる基本関数があってもよさそうな( たとえば table 関数のオプションとして)。 -- [[QDU]] &new{2006-04-10 (月) 23:02:06};
-ですよね(^_^;)要望してみればいかがでしょう。 -- &new{2006-04-10 (月) 23:33:23};
**ライブラリ [#i1aeaca4]
>[[超初心者]] (2006-04-03 (月) 15:18:03)~
~
Package source: evd_2.1-7.tar.gz ~
Windows binary: evd_2.1-7.zip ~
ライブラリには、上記のように、拡張子がgzのものと、zipのものがあります。当方はWINDOWS環境でRを使用することを考えていますが、2種類の拡張子は使い分けが必要なのでしょうか？あるいは、zipだけで十分なのでしょうか???~
//
-拡張子より，その前の情報を見ましょう。「Package source」はソースプログラム，「Windows binary」はウインドウズ用のバイナリプログラムです。インストールして，何にもせずにすぐ使うにはバイナリプログラムの方がよいでしょう。 -- [[青木繁伸]] &new{2006-04-03 (月) 15:28:02};
-R(最新2.2.1)では、CRANサイトにあるパッケーッジは動作するのでしょうか？パッケージにおいて、対応するRのバージョンは確認できないのでしょうか？~
何卒よろしくお願いします。 -- [[超初心者]] &new{2006-04-03 (月) 14:54:13};
-パッケージの添付ファイルに，その情報が書かれているものが入っていると思います。~
あなたも,
[[パッケージの説明:http://cran.md.tsukuba.ac.jp/src/contrib/PACKAGES.html]]はご覧になったのではないかと思いますが，そこに掲載されているパッケージ名をクリックすると，「Depends: R (>= 2.0.0)」のように，書いてありませんでしたか？ -- [[青木繁伸]] &new{2006-04-03 (月) 15:29:23};
-度々申し訳ありません。~
ライブラリの拡張子gz「Package source」はソースプログラムについてですが、これは、Ｒと同じディレクトリに保存すればよいのでしょうか？インストール方法を何卒ご教示下さい。~
（zipの場合は、Rメニューから呼び出し後インストール機能があるようなのですが??) -- [[超初心者]] &new{2006-04-07 (金) 10:40:05};
-ソースパッケージをビルドするにはRをビルドするのとほぼ同じ環境が必要になります. [[rw-FAQ:http://cran.r-project.org/bin/windows/base/rw-FAQ.html#Can-I-install-packages-into-libraries-in-this-version_003f]]に従って下さい. かなり面倒な作業であることは保証します. -- [[なかま]] &new{2006-04-07 (金) 11:31:40};
-Windows版のRをお使いなら、パッケージのインストールはメニューから選択することで実行できます（インターネットに接続していることが必要です）．具体的な手順は[[舟尾さんのR-tips:http://cse.naro.affrc.go.jp/takezawa/r-tips.pdf]]（29ページ〜）に分かりやすい説明があります．-- &new{2006-04-07 (金) 11:39:24};
-早速ご指導深謝します。舟尾様のページを見ました。 -- [[超初心者]];
-なかま様早速ご指導深謝します。~
状況としては、ライブラリの拡張子gz「Package source」がＰＣのＣドライブに保存されてます。~
ＲはＰＣにインストールされています。~
この場合、ソースをＲからインストールする場合
>install.packages("c:/***.gz",CRAN=NULL)
では、パッケージ***がインストールできないのでしょうか？~
誠に申し訳ありません。 -- [[超初心者]] &new{2006-04-07 (金) 12:57:16};
-Rtools.zipと言うのと,Mingw(C,Fortranコンパイラ)を入れてですね,[[FAQのリンク先:http://www.murdoch-sutherland.com/Rtools/]]に従って開発環境構築を行ってください. install.packagesでは無く,Rcmd install hoge.tar.gz をしましょう. 構築した環境が正しく機能するかは, R本体のビルドを行えれば間違いはありません. (だから面倒な作業なんですよ) -- [[なかま]] &new{2006-04-07 (金) 13:16:00};
-[[超初心者]]さんがインストールしたいのは[[evdパッケージ:http://cran.md.tsukuba.ac.jp/src/contrib/Descriptions/evd.html]]ではないのですか？普通にバイナリをインストールすれば、わざわざ開発環境構築を行う必要はないのでは？ -- &new{2006-04-07 (金) 14:13:00};
-超初心者さんなら，バイナリをインストールすればいいでしょう。あなたの環境（ビルドに必要なものがインストールされているのかどうかを含め），あなたの発言だけではわからないのですから。示された URL を読んでね。 -- &new{2006-04-07 (金) 14:58:47};
-礼儀作法がわからずに誠に申し訳ありません。具体的に申し上げますとインストール致したいのはrobustbase_0.1-4.tar.tarです。これをRで使用していのですが、どのようにインストールできるか具体的にご指導下さい。アドレスは、http://cran.md.tsukuba.ac.jp/src/contrib/robustbase_0.1-4.tar.gzです。何卒ご指導をお願いします。この度は申し訳ありません。 -- [[超初心者]] &new{2006-04-07 (金) 15:32:55};
-上のなかまさんのコメントの「[[FAQのリンク先:http://www.murdoch-sutherland.com/Rtools/]]」は見ました？そして，開発環境の構築はできました？その上で，やってみたけどできないというなら，どのようにやってみたらどのような状況・エラーメッセージが出たと，具体的に書きましょう。「具体的にご指導ください」なんていわれても，どこまで詳しく書かないといけないのかわからないし，もしコンパイラのインストールから説明しないといけないと言うことになったら，たいていの人は二の足踏みますよ。 -- &new{2006-04-07 (金) 16:09:29};
-[[ビルドのログ:http://cran.us.r-project.org/bin/windows/contrib/2.2/check/robustbase-check.log]]を見たんですが,ソースを開くと
#include <sys/types.h>
/* --> int64_t ; if people don't have the above, they can forget about it.. */
/* #include "int64.h" */
とか言ってるので, <sys/types.h>では無く<inttypes.h>にすれば, あなたも私も幸せになれるとお伝え下さい. -- [[なかま]] &new{2006-04-07 (金) 16:40:02};
-#include "int64.h" では？で，sys/types.h はある場合の方が多いのではないかなぁ？もうずいぶん C から離れているので．．．不確かだが。私もちょっとやってみたけど，上とは別のところで文句を言っている *.o がないって。あるんだけどな？ -- &new{2006-04-07 (金) 17:20:00};
-[[データサイズ:http://www.unix.org/version2/whatsnew/datasize.html]]です.mingwにも含まれます. int64.hはその場凌ぎでしょう. ビルドした環境と若干でも違うとバイナリ提供の環境では辛いかも. -- [[なかま]] &new{2006-04-07 (金) 19:04:22};
-カレントディレクトリにある *.h は " でくくるんじゃなかった？ということでしたが， -- &new{2006-04-08 (土) 21:38:14};
-それはそうです. ただ int64.hの中で何を定義しているかと言えばコメントに書いてあるとおりで, int64_tは何をインクルードするべきかと言う話です. stdint.h(C99)もありますが, ポータビリティ(otherOS)の点においてinttypes.hがよかろうと言う話です. ざっと斜めに読んでも他に問題は無いようです. sys/types.hもint64_t欲しさ(特定の環境では定義してしまいますが)に読んでいるようです. -- [[なかま]] &new{2006-04-08 (土) 23:16:53};
-ライブラリ(zip)でインストール時にエラーが発生します。
何卒よろしくご指導をお願いします。
一例をあげますと、copula_0[1].3-5~
http://cran.r-project.org/bin/windows/contrib/r-release/copula_0.3-6.zip
エラーメッセージは
> utils:::menuInstallLocal()
以下にエラーgzfile(file, "r") : コネクションを開くことができません
追加情報: Warning message:
圧縮されたファイル 'copula_0[1].3-5/DESCRIPTION' を開くことができません
その他、e1071_1[1].5-13、やkernlab_0[1].6-2もインストールしたいのですが、同様の現象です。~
これらzip形式ライブラリのインストール方法をご教示下さい。-- [[超初心者]] &new{2006-04-11 (火) 09:32:33};
-まだ質問の仕方がおわかりになっていないようですね。関係ないと思っているのでしょうが，今回の件では，あなたのインターネット環境でプロキシーサーバーが必要なのかどうかが重要な情報なんです。「単語検索」で「プロキシーサーバー」を検索すると，まさに，「Rをインストールしたのですが、CRANにアクセスできません。」という質問とそれに対する回答が書いてあります。 -- &new{2006-04-11 (火) 11:05:48};
-返す言葉もありません。正直なところ該当箇所を見ましたが、〜御恥ずかしい話ですが、それが答えになっているか〜どうかすらよくわかりません。システム担当者に教えを請います。〜この度はご迷惑をお掛けしました。 -- [[超初心者]] &new{2006-04-11 (火) 16:34:39};
-ご指導頂いたように、アーカイブの「プロキシーサーバー」を熟読し、.Rprofileにはoptions(CRAN="http://cran.md.tsukuba.ac.jp") Sys.putenv("http_proxy"="http://172.21.12.109:8080/")と入力しますとRのプロンプトで「エラー："options(CRAN="http://cran.md.tsukuba.ac.jp") Sys.putenv" に構文エラーがありました」と表示されました。何度も繰り返し試行していますが同じ結果です。なにとぞ、ご指導ください。 -- [[超初心者]] &new{2006-04-13 (木) 15:09:47};
-options 関数と Sys.putenv 関数は同じ行に入力したのですか？同じ行に入力したのならSys.putenv の前にセミコロンがないといけないのでは？2行に分けたのなら，別の原因でしょうが。 -- &new{2006-04-13 (木) 16:49:24};
-2段に別々に入力しました。本件に関連しているかどうかわかりかねますが、インターネットを閲覧するときに、当方の環境では、IDとパスワードを求められます。これは何か影響しているのでしょうか？誠に申し訳ありません。 -- [[超初心者]] &new{2006-04-14 (金) 14:33:16};
-何の ID とパスワードなんでしょう。どのページを見ても要求されるんですか？よくわからないことは，身近の人に聞くのがいいでしょう。あなたは何らかの組織の中にいるようであるし，「システム担当者に聞いてみる」とも書いていたのだから，その人に聞いたらいいでしょう。しかしまあ，よくまあ，こんなまだるっこしいやりとりが我慢できるものだと思いますよ。 -- &new{2006-04-14 (金) 14:34:42};
**数値の置き換え [#g6c16a4c]
>[[青木繁伸]] (2006-03-31 (金) 23:03:18)~
~
整数ベクトル x を，別の数値に置き換える関数を作る必要がありました。~
条件はn行3列の行列で与えようと思います。xの要素が1列目以上で2列目以下であるときに3列目の数値に変換しようと言うことになります。~
> y <- matrix(c(1,3,1, 4,7,5, 8,10,9), byrow=TRUE, nc=3)
> y
[,1] [,2] [,3]
[1,] 1 3 1
[2,] 4 7 5
[3,] 8 10 9
xの要素が，2のときはy[1,1] <= x <= y[1,2] ですので y[1,3] つまり1に置き換える，8のときは9になるなどです。~
xが次のようなとき~
> x <- c(7, 4, 3, 1, 8, 6, 5, 5, 8, 1, 1, 10, 7, 6, 3, 5, 10, 8, 3, 3)
答えが~
5, 5, 1, 1, 9, 5, 5, 5, 9, 1, 1, 9, 5, 5, 1, 5, 9, 9, 1, 1
になることを期待されます。~
そこで作ったのが，~
recode <- function(x, y)
{
sapply(x, function(z) y[,3][y[,1] <= z & z <= y[,2]])
}
ですが，もっとスマートな定義があるでしょうか...~
//
-library(car)の関数recode()の定義の長さを考えると（あちらは変換ルールを文字列として与える点で汎用性が高いので直接比べるのは不公平かもしれませんが），1行でできてしまうというのは驚異です。ちなみにcarライブラリのrecodeでは，recode(x,"1:3=1;4:7=5;8:10=9")となります（回答としてはズレているので匿名にしておきます）。 -- &new{2006-04-03 (月) 20:08:24};
-コメント，ありがとうございます。私の悪い癖で，探すより作る方に走ってしまいました。ご紹介頂いた recode も候補としたいと思います。 -- [[青木繁伸]] &new{2006-04-03 (月) 21:08:00};
-青木さんの例のように置き換える数値が単調増加正整数であれば、次のような（少し前処理が必要なのでスマートではないですが）やりかたも可能ですね。（自分で関数書くのが嫌いなたなぼた期待不精者。） base パッケージ中の関数 findInterval は第一引数ベクトル中の数字が第二引数で定義される区間列中の何番目に入るかを計算する関数です。区間は半開区間 [a.b) とされますからわざと [a,b+0.5) となるようにし、結果が1,5,9番目の区間になるようにわざと途中にダミーの区間をいれています。-- [[間瀬茂]] &new{2006-04-04 (火) 09:31:44};
> x <- c(7, 4, 3, 1, 8, 6, 5, 5, 8, 1, 1, 10, 7, 6, 3, 5, 10, 8, 3, 3)
# 1,2,3 は第一の区間 [1,3.5)、4,5,6,7 は第5の区間 [4,7.5)、8,9,10 は第9の区間 [8,10.5) に入る
> y <- c(1,3.5, 3.6, 3.7, 4,7.5, 7.6,7.7, 8,10.5)
> findInterval(x,y)
[1] 5 5 1 1 9 5 5 5 9 1 1 9 5 5 1 5 9 9 1 1
**OR？ [#z1565397]
>[[内藤]] (2006-03-24 (金) 10:57:13)~
~
X中のF2=1，かつF3=1であるF1を抽出するには，~
下記のとおりでいいと思うのですが，~
X$F1[X$F2=1 & X$F3=1]~
X中のF2=1，または，F3=1であるF1を抽出するには，~
どのようにすればよいでしょうか？~
//
-思うぐらいならやってみる. -- [[なかま]] &new{2006-03-24 (金) 12:32:12};
> X<-data.frame(F1=sample(3,10,replace = TRUE),
+ F2=sample(3,10,replace = TRUE),
+ F3=sample(3,10,replace = TRUE))
> X$F1[X$F2==1 & X$F3==2] # AND
[1] 1 3
> X$F1[X$F2==1 | X$F3==2] # OR
[1] 1 3 1 2 2
>?"&"
>?"=" # これと
>?"==" # これの違いも
> X$F2==1 | X$F3==2 # これが何を出力するかとかも
-どうもありがとうございます。助かりました。 -- [[内藤]] &new{2006-03-24 (金) 20:00:20};
**図をpdfファイルに保存するには [#j0387b16]
>[[初心者ST]] (2006-03-20 (月) 16:50:24)~
~
皆様．初めて質問させていただきます．初心者STです．宜しくお願いいたします．
x <- c(0:100)
y <- 10000 / (1+0.05*x)
z <- 10000 / (1+0.10*x)
plot (x,y, ylim=c(0, 10000), main="hypothetical function", sub="Fig 1",pch=2)
par(new=T)
plot (x,z, ylim=c(0, 10000), ann=F, pch=4)
上記の図を描いた後，ファイル−別名で保存より，pdfを選択すると，次のようなエラーが生じて，pdfファイルとして正常に保存できません．
エラー：invalid character sent to 'PostScriptCIDMetricInfo' in a single-byte locale
追加情報: Warning messages:
1: invalid string in 'PostScriptStringWidth'
2: invalid string in 'PostScriptStringWidth'
3: invalid string in 'PostScriptStringWidth'
何が原因でしょうか？基本的な質問で大変申し訳ありませんが，もしお分かりの方がいらっしゃれば，ご教示頂ければ幸いです．なお，使用環境は，R version 2.2.1, 2005-12-20, i386-pc-mingw32です．以上，宜しくお願い申し上げます．~
//
-?pdfをしてください。
pdf()
x <- c(0:100)
y <- 10000 / (1+0.05*x)
z <- 10000 / (1+0.10*x)
plot (x,y, ylim=c(0, 10000), main="hypothetical function", sub="Fig 1", pch=2)
par(new=T)
plot (x,z, ylim=c(0, 10000), ann=F, pch=4)
dev.off()
初心者ならどんな基本的なことを質問しても良いということではありません。まずは検索とヘルプで調べましょう。-- [[Akira]] &new{2006-03-20 (月) 17:28:36};
-Akira様．初心者STです．準備不足でいろいろとご迷惑をおかけして申し訳ありませんでした．それにもかかわらず，早速に上記のアドバイスを頂き，有難うございました．ご指摘のデバイス関数を用いたpdf(file="*.pdf")とdev.off()の組み合わせによる方法はThe R tipsという著書でも紹介されていましたので，先に試してみたのですが，やはり，同じエラーが生じるようです．取り急ぎ，御礼とご報告まで． -- [[初心者ST]] &new{2006-03-20 (月) 18:01:18};
-今回初めて pdf ファイルとして書き出そうとしたのでしょうか。他のフォーマットでは書き出せますか？インストールをやり直してみるというのも，基本的ではありますがやってみる価値はあるかもしれません。 -- [[青木繁伸]] &new{2006-03-20 (月) 18:39:00};
-R version 2.2.1, 2005-12-20, i386-pc-mingw32 512M WXPSP2で初心者STさんと同じことをやってみましたが、正常にPDF出力されました。青木先生の言われたとおりインストールをやりなおしてみたらいかがでしょうか？ -- [[okinawa]] &new{2006-03-20 (月) 20:57:42};
-R version 2.2.1, 2005-12-20, i386-pc-mingw32 256MB WXPSP1でもOKでした。出力できないのはpdfだけでしょうか？やはり再インストールが良いと思います。 -- [[Akira]] &new{2006-03-21 (火) 00:02:11};
-多分インストール時に東アジアを選択しなかったのでは？(この選択は2.3.0以降は無くなるはず) -- [[なかま]] &new{2006-03-21 (火) 00:11:07};
-青木様，okinawa様，Akira様，なかま様．初心者STです．再インストールし直したところ，メニューバーを使う方法，デバイス関数を使う方法のいずれにおいても，pdfファイルに正常に書き出せることが確認されました．また，これまでうまく書き出せなかった原因は，なかま様のご推察のように，インストールの途中で，Version for East Asian languageを選択しなかったことにあるようです（再インストールの結果，メニューの文字化けが解消されたことも，この説の正当性を裏付けました）．皆様の有益なコメントのおかげで，問題解決に至りました．お忙しいところ，本当に，有難うございました． -- [[初心者ST]] &new{2006-03-21 (火) 16:34:29};
**大量のグラフを描くためには [#a6a1ddb0]
>[[antt]] (2006-03-14 (火) 00:08:58)~
~
はじめまして。大量のグラフを描きたいのですが，以下に記したようにエラーが出て解決できません。どなたか解決方法をご存知の方，お教え頂けると幸いです。~
【環境】R version 2.2.1, i386-pc-mingw32, メモリ 1GB~
【問題】大量のグラフを描き続けると，オブジェクトは別に増えていないのに，メモリを消費したままで，その内，Rがフリーズする~
【知りたいこと】大量のグラフを描きたい場合に，どのようなプログラムを書けば，メモリの消費を抑えられるか知りたい~
【プログラム例】
alternative <- 1:5 #選択肢の数
nitem <- 100 #項目の数
nschool <- 10 #施設の数
n <- 10000 #標本サイズ
##仮想データを作成
data <- data.frame(matrix(alternative, ncol=nitem, nrow=n))
data$school <- rep(1:nschool, each=n/nschool) #1施設 1000人
##画像作成用の関数を定義する
plotitem <- function(i, j) {
dir.create(paste("school", i, sep=""), showWarnings=F) #施設ごとにディレクトリを作成
filename <- paste("./school", i,"/item", j, ".ps", sep="") #施設ごとに，各項目のファイルを出す
postscript(filename) #保存する
subdata <- which(data$school==i) #施設ごとのデータを抽出する
hist(data[subdata,i], breaks=c(.5 + 0:5))
dev.off()
}
##まずは，1施設，10項目を描画してみる
for (i in 1) {
for (j in 1:10){
plotitem(i,j)
}
}
##次に，1施設，100項目を描画してみる
##メモリサイズにより，この時点でRがフリーズする
##postscript(); dev.off() が終了した時点でも，メモリが消費しているようである
for (i in 1) {
for (j in 1:100) {
plotitem(i,j)
}
}
//
-がんばって整形してくれましたが、無意味に長すぎるのも読む気を無くすものです。少し短くしてみました。 ついでにコメントすると Linux, 1Gb の私の環境では j in 1:1000 としても問題なくファイルが作れます。 MSW 固有の問題？-- [[MKR]] &new{2006-03-14 (火) 07:51:26};
-R version 2.2.1, 2005-12-20, i386-pc-mingw32 メモリ512M WindowsXPSP2で実行しましたが正常に出力されましたよ。 -- [[okinawa]] &new{2006-03-14 (火) 08:16:34};
-闇版ならFMの処理のリークです.setHook(packageEvent("grDevices", "onLoad"),function(...) grDevices::ps.options(family="Helvetica")) 等にすればCIDのメトリックをロードしないので大丈夫だと思います.R-develでは正月開けにこの処理はwcwidthに置き換えられたのでこの問題自体無くなりました.(^_^; -- [[なかま]] &new{2006-03-14 (火) 10:09:39};
-日本語をどうしても使いたい場合は, Postscriptを吐くならTeXの環境をお持ちだと思うので, 一端複数頁のPostScriptを作成してpsselect -p1 複数P.ps 1頁.ps とかして頂ければ回避出来るかと思います. -- [[なかま]] &new{2006-03-14 (火) 10:20:22};
-ご回答ありがとうございました。仰るとおり闇版を使っております。「FMの処理のリーク」という意味を理解しかねるのですが，ようするに，日本語を含めるグラフを作成したい場合は，上記の問題が発生するという理解で良いのでしょうか。複数ページのpsやpdfを作成して回避することを検討してみます。 -- [[antt]] &new{2006-03-14 (火) 16:22:02};
-フォントのメトリック情報の処理によってリークが発生するということです.(あんまり変わってないような...滝汗) Type1に比べてCIDのサイズが大きいから顕著になるだけで(リークの原因は他にもあるんで),普通に使う分には問題ないかと.(いいかげんですみません) -- [[なかま]] &new{2006-03-15 (水) 23:44:58};
-リーク => メモリーリーク => 使われたまま回収再利用不可能なメモリが生じる、てな理解で良いのかしら（？）。 -- [[半可通]] &new{2006-03-16 (木) 01:07:40};
-はい。 -- [[なかま]] &new{2006-03-17 (金) 00:08:31};
**imageの色を絶対指定するには [#i140c8f9]
>[[子牛]] (2006-03-13 (月) 19:39:19)~
~
はじめまして。imageで得た画像にcolをつかって色を指定したいのですが、たとえば次の例にあるような二つのグラフの違いを明らかにしたいのです。この例ではグラフaはグラフa-1よりもいつも1だけ大きくなります。~
数値にたいして絶対的な方法で色を指定するにはどうすればよろしいのでしょうか？ ~
~
> a <- rnorm(100)
> dim(a) <- c(10,10)
> b<- a-1
>
> x <- 1:10; y <- 1:10
>
> op <- par(mfrow=c(1,2))
> image(x,y,a, col=cm.colors(100), main=paste("a"))
> image(x,y,b, col=cm.colors(100), main=paste("a-1"))
>
> par(op)
~
以下は使用環境のメモです。~
> sessionInfo()
R version 2.1.1, 2005-06-20, i386-pc-mingw32
attached base packages:
[1] "methods" "stats" "graphics" "grDevices" "utils" "datasets"
[7] "base"
//
-100段階のパレットで一段階の違いを表現するという例題は，効果がわかりにくいので，10段階の違いを110段階のパレットで表現する例に変換してみました。要するに colors(n) を多めに用意して，二つの行列に共通して適用する，ただし「オフセット付きで」と言うことですが，~
a <- rnorm(100)
dim(a) <- c(10,10)
b<- a-10 # ***
x <- 1:10; y <- 1:10
op <- par(mfrow=c(1,2))
image(x,y,a, col=cm.colors(110)[11:110], main=paste("a")) # ***
image(x,y,b, col=cm.colors(110)[1:100], main=paste("a-10")) # ***
par(op)
これでいかが？ -- [[青木繁伸]] &new{2006-03-13 (月) 20:45:28};
#ref(pict2.png)
-ありがとうございました、色の使う範囲をこれで指定できるのですね。~
ただ、教えていただいた方法ですと、オフセットを手で調整せねばなりません。
image(x,y,a, zlim=c(-2,2), col=cm.colors(110), main=paste("a"))
この方法で絶対的な値を固定できるようです。（すみません、昨夜の段階ではHelpの説明が理解できてませんでした）。ちなみに上下の指定をはみだした値は、透明のイメージになるようです。はみ出さないようにしておかないと、落としてしまいそうですね。-- [[子牛]] &new{2006-03-14 (火) 14:47:55};
**SVMのマージン距離出力方法について [#o5f5dbe6]
>[[ちょめ夫]] (2006-03-13 (月) 12:10:44)~
~
現在SVMに関するパッケージの検討を行っております。~
検討対象としているパッケージは以下の２つです。~
１）e1071~
２）kernlab~
なお、判別したいクラスは2クラスの分類となっております。上記二つのパッケージのどちらでも良いのですが、SVMの予測結果において、与えた各サンプルと各クラスのマージンまでの距離を知りたいと考えておりますが、e1071,kernlabのpredict関数ともに、予測されたクラスやProbability、vote結果は出力可能ですがマージンまでの距離はどのように求めればよいか分かりません。イメージとしては、以下ウェブページの~
http://www.msi.co.jp/vmstudio/materials/tech/classification.html#1-5~
「図 Support Vector Machine 出力画面」のSVM.出力.1、２のような
数値情報です。もしご存じの方いらっしゃったらアドバイスいただければ幸いです。
どうぞよろしくお願いいたします。
> sessionInfo()
R version 2.1.1, 2005-06-20, i386-pc-mingw32
attached base packages:
[1] "methods" "stats" "graphics" "grDevices" "utils" "datasets" [7] "base"
other attached packages:
colorspace kernlab e1071 class
"0.9" "0.6-2" "1.5-11" "7.2-16"~
//
-投稿法についてご指摘いただきまして、誠に申し訳ありませんm(_ _)m -- [[ちょめ夫]] &new{2006-03-13 (月) 13:40:16};
-当てずっぽうですが、もし話題の量が SVM アルゴリズムにとり重要な中間量ならば、関数の出力に暗に含まれている可能性が高いと思います。 summary 関数で表示されないならば、str 関
![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
