初級Ｑ＆Ａ アーカイブ(5)
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
COLOR(green){SIZE(20){初心者のための R および RjpWiki に...
新規投稿はできません
----
-[[初級Ｑ＆Ａ アーカイブ(4)]] (元記事が 2005-11-09 より 2...
-[[初級Ｑ＆Ａ アーカイブ(3)]] (元記事が 2005-05-02 より 2...
-[[初級Ｑ＆Ａ アーカイブ(2)]] (元記事が 2004-12-13 より 2...
-[[初級Ｑ＆Ａ アーカイブ(1)]] (元記事が 2004-08-03 より 2...
----
#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...
-ありがとうございました。sinkでファイルに書き込むことがで...
-> 行列を一列にするプログラムを書きました~
プログラムを書くって，as.vector(行列オブジェクト）だけで...
それと，ファイルに取り込むの？書き込むの？「関数の値」っ...
-適当なエディタを開き，画面（コンソール）に表示されたもの...
**R-2.3.1のインストール手順10のファイルの上書き [#m9ece0ae]
>[[初心者]] (2006-06-24 (土) 12:03:31)~
~
インストールの手順の最後にRconsole，Rdevga，Rprofile.site...
//
-ファイル名のところにリンクが張ってありません？？ -- &ne...
-インターネットでは（そのほかでも），半角カタカナは使わな...
-ファイル名の所を押しても文字が出てくるだけなのですが？ -...
-MIMEタイプの設定がうまくできていないサーバーなんでしょう...
このような場合の対処法は，覚えておくと，今後得しますよ。~
そのような場合は，うんとね，ういんどーずのようですから，...
**関数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{...
**関数名を文字列に [#l2aaf4df]
>[[ショーン]] (2006-06-22 (木) 19:50:20)~
~
以下のように関数オブジェクトを変数に入れて、関数として呼...
> a = cos
> a
.Primitive("cos")
> a(0)
[1] 1
後でaの元の関数名である、"cos"という文字列を得る方法はあ...
> f2s(a)
[1] "cos"
のf2sようなことはできますでしょうか？~
//
-出来ないことは無いですが、まずなぜそんなことが必要なのか...
-理由がないと教えられないと言うことでもないようには思いま...
もっとも，私にはわかりませんでした。~
変数に関数定義を付値したら，その変数の class は function ...
> 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"
-答えるにはそれなりの手間がかかりますから、必要性が私には...
-勉強の機会を与えてくれて，ありがとうございます。後者の場...
> 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"
-うむ、あなたは初級者などでは決してない。何につかうのか秘...
f2a <- function(x) {capture.output(x ,file=(filename <- ...
-私は，原質問者じゃないよ。~
だから，原質問者は用途が秘密だとは言っていない。~
質問者をもてあそぶヒマはあっても，答えるヒマはないという...
-暇のある無しの問題ではありません。質問するからには、やは...
-回答ありがとうございます。おかげさまでうまくいきそうです...
funcs = c(cos,sin,tan)
のようにして、ループでいろいろな関数のグラフを描こうと思...
-2006-06-22 (木) 23:33:51 は，もっと簡単になるでしょう？ ...
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() は関数中でつかわれ、実引数の数...
func <- function(...) print(nargs())
-なるほど。~
純粋に ... の個数を知りたい場合には，以下の方がよいのかな...
> func = function(a, ...){print(length(c(...)))}
> func(1,2,3)
[1] 2
-【「Rの関数定義の基本」を見るとlength(...)という書き方も...
-ショーンさんはlength(1:3)とlength(c(1,2,3))、length(1,2,...
-COLOR(blue){【「Rの関数定義の基本」を見るとlength(...)と...
> 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(...)) でし...
-ありがとうございました。よくわかりました。 -- [[ショーン...
**データ行列の行をまたぐレコード参照 [#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と比較演算子を組み合わせればできそうに思いますが -...
-例に挙げたデータについての例解が間違えているし。条件の記...
-まずは，val1 < 35.0 という条件なら，それを満たさないもの...
-しかしまあ，val1 < 35.0 は，第 i 行に関する条件なのか，...
-いろいろ面倒だが，骨格を作ったから直してみて？val1 につ...
たぶん遅いから，R じゃなきゃいけないとか R の機能をふんだ...
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() 調べてみます&ご提示頂いた...
その他の皆様も，アドバイスありがとうございます(たしかに例...
-上のコードと本質的に同じコードで50万行の人工例を（一部実...
-一昼夜コンプータを動かしておけば解が出るなら，安いモンで...
**Williamsの多重比較のパーセント点 [#p8461227]
>[[tm]] (2006-06-20 (火) 22:23:00)~
~
Rとは直接関係なく、大変申し訳ありませんが、Williamsの多重...
身近で入手できる本をあさったのですが、5%点くらいしか手に...
できれば、1%点の表が欲しいのですが、どなたかご存知ありま...
//
-その本に，生成アルゴリズムなどは書いてなかったですか。。...
**.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>
の実際のコードを覗いてみることはできないのでしょうか。~
//
-過去に同じ類の質問が沢山あります。ちゃんと調べましょう。...
-.Internal については，あったかな？ -- &new{2006-06-20 (...
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の値が同じになると思っていたのですが、違います。ど...
//
-主成分負荷量を求めるのでしょうか？だとしたら，貴方は，説...
> d <- structure(c(0.1991, 0.2003, 0.2169, 0.1429, 0.184...
+ 0.1516, 0.1478, 0.2071, 0.1613, 0.2288, -15.0277, -16....
+ -17.3694, -16.7424, -15.0885, -15.7661, -15.1187, -14....
+ -13.5147, -15.0654, 0.279, 0.3462, 0.299, 0.2586, 0.25...
+ 0.155, 0.294, 0.3411, 0.2789, 0.2592, 0.2165), .Dim = ...
+ ), .Dimnames = list(c("1", "2", "3", "4", "5", "6", "7...
+ "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/PC...
res2$rotation
が行列Lですよね。それから、
res2$x
が主成分zですよね。投稿した質問のように、元の変数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
-ありがとうございました。よくわかりました。 -- [[ショーン...
**Sweave環境でsink()は使えないのでしょうか？ [#fda5d2d5]
>[[akira]] (2006-06-16 (金) 14:11:29)~
~
Sweave環境でxtableを書き換えたい場合、~
(1)xtableの出力先をtexファイルからsink()で別のファイル（...
(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(...
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がないからxta...
'RweaveLatex', 'Rtangle'当たりを見てもそれらしい記述が見...
//
-[[print.xtable2:http://sapporo.cool.ne.jp/matsut/r02.txt...
**multcomp パッケージでの Williams 検定について [#pa7f07ff]
>[[tm]] (2006-06-15 (木) 21:10:37)~
~
multcomp パッケージ中の simtest の Williams 検定について...
Rのバージョンは 2.0.1、multcomp のバージョンは 0.4-8、mvt...
> library(mvtnorm)
> library(multcomp)
> data(angina)
> summary(simtest(response ~ dose, angina, type="William...
Simultaneous tests: Williams contrasts
Call:
simtest.formula(formula = response ~ dose, data = angina...
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...
Simultaneous 95% confidence intervals: Williams ...
Call:
simint.formula(formula = response ~ dose, data = angina,...
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 検定については情報...
使用方法や結果の見方などご教授いただけないでしょうか？~
よろしくお願いいたします。~
//
-あるプログラムにバグがあるというとき，確証がないといけま...
プログラムの使い方などについては，教科書などに載っている...
> library(mvtnorm)
> library(multcomp)
> set.seed(12345678)
> df <- data.frame(g = factor(rep(1:6, each=20)), y = rn...
> summary(simtest(y ~ g, df, type="Williams", alternativ...
Simultaneous tests: Williams contrasts
Call:
simtest.formula(formula = y ~ g, data = df, type = "Will...
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...
Simultaneous 95% confidence intervals: Williams
contrasts
Call:
simint.formula(formula = y ~ g, data = df, type = "Willi...
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』に見えるのは出力の有効桁数だけの問題では？全ての...
-ご返信、アドバイスありがとうございます。 ~
まず、いいたかったのは「バグがあるよ。」ということではな...
おっしゃられるように、正解例がないとよくわかりませんので...
> library(mvtnorm)
> library(multcomp)
> set.seed(12345)
> df <- data.frame(
+ Dose = ordered(rep(c("Cont", "Low", "MidLow", "MidHigh...
+ 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="...
Simultaneous tests: Williams contrasts
Call:
simtest.formula(formula = Response ~ Dose, data = df, ty...
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 の値が...
この違いは、私の R 環境の問題でしょうか？それとも、実装さ...
ご教授いただけますよう、よろしくお願いいたします。
-- [[tm]] &new{2006-06-20 (火) 10:32:31};
-Contrast マトリックスと件の本を比べると違いがわかるので...
-contrast マトリックスの意味がなんとなく理解できました。...
ただ、おっしゃられるように反応の並びや greater、 less あ...
実装が、本と異なるのは理解できましたが、contrast マトリッ...
> df <- data.frame(
+ Dose = ordered(rep(c("Cont", "Low", "MidLow", "MidHigh...
+ 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="...
Simultaneous tests: Williams contrasts
Call:
simtest.formula(formula = Response ~ Dose, data = df, ty...
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 の中を調べる必要があ...
**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" [#ic...
>[[kaz]] (2006-06-12 (月) 16:00:47)~
~
glmを用いておよそ100万件のデータをロジスティック回帰して...
~
インターネットで検索したところ、perfect separation、つま...
~
他にこのエラーの生じる原因をご存知であれば教えていただけ...
//
-また、5つの説明変数は全て連続量です。Rと回帰分析の初心者...
-Rのバージョンは2.2.1（日本語化版）です。 -- [[kaz]] &new...
-Quasi-complete separation とか Complete separation とい...
簡単に説明すれば，ある変数（または複数の変数の線形結合）...
データが何百万，何千万あろうとも，ある変数の分布が重なら...
共用パソコンに入っているなら仕方ないけど，最新版を使うの...
> 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 M...
-5.647e-05 -2.107e-08 0.000e+00 2.107e-08 5.560e-...
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...
2: 数値的に 0 か 1 である確率が生じました in: glm.fit(x ...
-大変わかり易い具体例を有難うございます。頭がすっきりしま...
-私のやり方が間違っているのかもしれませんが、上記のように...
-先ほど追加した部分と関連しますが，複数の変数の線形結合値...
-set.seed も含めて。上の通りやってみましたか？set.seed が...
> 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 M...
-1.256e-04 -2.107e-08 0.000e+00 2.107e-08 1.282e-...
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...
2: 数値的に 0 か 1 である確率が生じました in: glm.fit(x ...
-また、（これは単にRのバージョンの問題かもしれませんが）w...
-なるほど、確かに単回帰のエラーは重回帰のエラーの十分条件...
-set.seed で生成される乱数列は，バージョンによっても違う...
-個々の変数の分布は重なるが，線形結合の分布が重ならない例...
> 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...
(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...
(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...
(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 M...
-9.418e-05 -2.107e-08 0.000e+00 2.107e-08 6.798e-...
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...
2: 数値的に 0 か 1 である確率が生じました in: glm.fit(x ...
-詳しい例示ありがとうございます。「個々の変数の分布は重な...
-> 「個々の変数の分布は重なるが，線形結合の分布が重ならな...
上の例を見ればおわかり頂けると思ったのですが。分布が重な...
その方が理解しやすいということならかまいませんが，必要条...
上の x から合成される z を，y 別に描いてみる
重なっていないことがわかる
boxplot(z ~ y, horizontal=T)
#ref(boxplot.png)
- ずいぶん時間が経ってから見させていただきましたが、参考...
-Firth bias-correctionという方法があるようです -- nu &new...
--参考
---http://cran.r-project.org/web/packages/brglm/index.html
---http://www.ats.ucla.edu/stat/mult_pkg/faq/general/comp...
**Sweaveでxtableをminipage環境に適用したいです [#x624f617]
>[[Akira]] (2006-06-09 (金) 10:40:42)~
~
latexコマンドをRスクリプトに埋め込んで、Sweaveで結果をま...
xtableで作成したtableをminipage環境で2つ横に並べたいので...
> 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のオプションを確認しましたが、出力しない設定が分かり...
何か他に良い方法がありますでしょうか？
R-2.3.0、WinXPSP1を使用しています。
//
-print.xtableをみましょう. カラクリ屋敷(ぼそっ)... -- [[...
?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は...
-なかまさまの前に記入されていたsinkを利用する方法がなくな...
-バックアップをたどればありますよ. (良く見ろとか怒られる...
-勝手にしていいものかと・・・でも、Tipsに乗せてしまいまし...
**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.yyyymmdd <- function(sdate, sep="/") {
require(survival) # もちろんシステムにあらかじめインス...
temp <- date.mdy(sdate) # date.mdy は survival の関数
ifelse(is.na(sdate), "NA", paste(temp$year, temp$month,...
> date.yyyymmdd(100:105)
[1] "1960/4/10" "1960/4/11" "1960/4/12" "1960/4/13" "196...
ふむ、何か変かな。次のようになるのだけれど。
> date.yyyymmdd(13149)
[1] "1996/1/1"
-パッケージ data にも survival にも全く同じ関数があります...
-dates関数って，どのライブラリにあるのか。。。。~
date.yyyymmdd は，以下のように微修正すればよいでしょう。 ...
> date.yyyymmdd <- function(sdate) {
+ require(survival) # もちろんシステムにあらかじめイン...
+ temp <- date.mdy(sdate) # date.mdy は survival の関数
+ ifelse(is.na(sdate), "NA", sprintf("%4i/%02i/%02i", ...
+ }
> date.yyyymmdd(100:105)
[1] "1960/04/10" "1960/04/11" "1960/04/12" "1960/04/13" ...
> date.yyyymmdd(13149)
[1] "1996/01/01"
-13149 て Julian day (1960/01/01) からの経過日数ですよね...
-変も何も，そういう前提で経過日数後の日付を yyyy/mm/dd の...
> dates2 <- function(x)
+ {
+ s <- as.character(dates(as.numeric(13149),out="yyyy/m...
+ sprintf("%s%02i%s",substring(s, 1, 5),which(month.abb...
+ }
> dates2(13149)
[1] "2006/01/01"
-うん？
> dates(0,out="yyyy/mm/dd")
[1] 1970/Jan/01
なんだから，dates 関数では 1970/01/01 からの日付では？~
結局は，関数の定義によるようで，別の関数で代替処理すると...
-datesって S-Plusの関数では? -- &new{2006-06-09 (金) 07:...
-google で少し調べたら Julian date とは本来紀元前4713/1/...
-元の発言にも記載がないので，いったいどのdates関数だと，...
> 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関数についての記述をきちんとせず、皆さんを混...
今回は自分で関数を書いてこの問題を解決しようと思います。 ...
**分散の違い [#hfeda6c6]
>[[kanako]] (2006-06-07 (水) 14:21:24)~
~
異なる平均と標準偏差を持つ分布を3種類用意して、それらを順...
> 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...
> 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
~
//
-乱数は文字通りランダムですから、毎回同じ物が作られるとは...
> rnorm(3,0,1)
> rnorm(3,0,1)
> set.seed(1234)
> rnorm(3,0,1)
> set.seed(1234)
> rnorm(3,0,1)
-コメントありがとうございました。説明が不足していました。...
-ちなみに、二つの例で分散（と言うよりは共分散行列ですね）...
**データから円グラフを作成 [#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 を見てく...
-ありがとうございます。頑張ってみます。 -- [[やきそばぱん...
**プロビットモデルでの疑似決定係数の計算 [#o080db1a]
>[[イマイ]] (2006-06-05 (月) 23:18:44)~
~
プロビットモデルの回帰分析を行っていますが、疑似決定係数...
超初心者なりに考え付く限り検索をかけてみたのですが、みつ...
//
-前にも同じ質問をどこかでしてませんでしたっけ？というのは...
「プロビットモデル 疑似決定係数」でググってみると，20件く...
「12）プロビットモデルで疑似決定係数を計算する方法はいく...
-尤度比インデックスに関しては存じ上げませんが，疑似決定係...
-検索がまだまだ不十分だったようです。アドバイス、本当にあ...
**Mac OSXでのR 2.3.1でパッケージアップデートができない [#...
>[[高井]] (2006-06-05 (月) 17:21:12)~
~
PowerPC G5（デュアル 2.3GHz）, OSX 10.4.6でR 2.3.1のGUI版...
以下にエラーupdate.packages(lib = ""/Library/Frameworks/...
オブジェクト "Library" は存在しません
とあります。ところが，ターミナルから起動してupdate.packag...
//
-補足です：エラーメッセージが出た後に，Rパッケージインス...
-R コンソールから，update.packages() すると，バージョンの...
-はい，そうですね。Rコンソールの結果もターミナルからの結...
-Libraryの前のダブルクォートが２つになっているのが怪しい...
-翻訳はアップデートされてたけどアップデートは直ってないで...
-間違いなくバグでした。PackageInstaller.mの中でupdate.pac...
-私の環境のせいではなくて一安心です。 -- [[高井]] &new{20...
**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.si...
2: Reached total allocation of 734Mb: see help(memory.si...
> summary(slm2)
以下にエラーsummary(slm2) : オブジェクト "slm2" は存在し...
となって計算が実行できません。~
ちなみにメモリーは
> memory.size(T)
[1] 61349888
です。~
Rではハードディスクを用いた仮想メモリを利用できる機能等は...
//
-help(memory.size) はしてみましたか？ この wiki 内を「メ...
-help(memory.size)も実行済みです。memory.limit(size=1024)...
-NULLが返ってきた後に、memory.limit()をすると上限が増えて...
-(A+B+C+D+E)^2を独立変数にした回帰で，しかも変数選択をし...
-memory.limit()で確かに上限が増えています。が、計算が上手...
--モデル自体の吟味をすることを，まずお薦めします~
そうかもしれません。実は回帰と言うより、実際は数量化I類で...
ところで｢起動コマンドライン｣とは何でしょうか?単語検索でも...
-Rを起動するアイコンの上でマウスを右クリックしてプロパテ...
-それを｢起動コマンドライン｣と言うのですね。実は--sdiにし...
-起動アイコンのコマンドラインパラメータなので，そう呼んで...
-なるほど。Tipsありがとうございます。数量化I類は結局どの...
-カテゴリー変数をダミーデータに変換して重回帰分析を行う場...
-Rのfactor関数を使ってダミー化したものですから、元締めの...
#ref(QQplot.jpg)
あれ、画像添付こんなのでいいのかな･･････?~
ご覧のように、正規Q-Qプロットを一応試みたのですが、大部分...
そんなわけで、交互作用の存在を疑ってみたのです。-- [[R初...
-step関数は大変便利な方法だと思いますが、強引な面も持って...
まずは、主効果のみモデルのAICを求めて、~
その後、2次の交互作用項を一つ入れたモデルを作って~
AICが改善されるかを見ていくのもいいと思います。-- [[MK]] ...
-なるほど。力技ですね。試してみます。 -- [[R初心者]] &new...
-バグってRが止まってしまいました。。。。。。○|￣|＿~
ただ、交互作用A:Bを入れたモデルだけは計算できたんですが、...
ところで、AIC関数を用いるとstep関数で言うRSSが表示される...
-その交互作用項のt値はどうでしたか？有意であれば、交互作...
-交互作用項A:BのP値は全て 2×10^-16以下を示しています。
#ref(QQplot2.jpg)
交互作用A:Bを考慮した場合の正規Q-Qプロットです。あんまり...
-交互作用を疑うよりは，81510番のデータを良く精査して，問...
-81510番は確かにイレギュラーに見えるので、次のようにして...
df <- df[-c(81510),]
それで元の交互作用項を含まないモデル式に戻して、正規Q-Qプ...
#ref(QQplot4.jpg)
今度は乖離が派手に見えます。指数関数的に増加していってま...
-これは，交互作用と言うよりは，独立変数と従属変数が直線相...
-大きすぎると思います。wiki上で図を小さくする方法が分から...
正規Q-Qプロットより残差プロットかS-Lプロットの方が宜しい...
-最初から小さな図を描けばいいのです。ま，ちょっと描いてみ...
残差は以下のような分布になることが期待されています。同じ...
#ref(temp.png)
#ref(temp2.png)
-大量のデータを全部使って分析しなければ得られない結果とい...
-確かにRとあんまり関係無くなってますね。・・・・・・どっ...
#ref(回帰診断.jpg)
Rで画像を小さくする、と言うのはこんな感じで宜しいのでしょ...
間瀬先生の本に従って、見よう見まねで回帰診断プロットを行...
--あるモデルで表現できるシステムならば，データ数は少なく...
--近代統計学はサンプルを用いて母集団を推測するという点を...
仰る事は分かります。その点も気になって、逆に、｢データが多...
場合に拠ってはWinBugsを試してみようかな、とは思っています...
-0を中心にしてばらついているというのではなく，予測値が70...
-外れ値の影響が大きいんですね。という事は・・・・ロバスト...
-とか言ってましたがダメでしたね。
エラー：lqs failed: all the samples were singular
ロバスト/抵抗回帰と数量化I類は相性が悪い?連続量じゃないと...
**windows版Rのインストール [#sae95e3a]
>[[ありんこ]] (2006-06-01 (木) 16:51:01)~
~
すみません、中級者の所にしゃしゃり出てしまって…。インスト...
//
-最新版の[[R-2.3.0:http://cran.md.tsukuba.ac.jp/]]を使い...
-と書いた矢先にR-2.3.1が。上のリンク先（筑波大ミラーサイ...
-確かにインストールの説明のところに新しい2.3.1が追加され...
-「インストールの説明のところに新しい2.3.1が追加されてい...
-実はみんな違う人が書いたんだったりして:-) -- [１つは書...
**ノンパラメトリックなモデルでの変数選択 [#i7bc81e2]
>[[金井]] (2006-06-01 (木) 01:35:15)~
~
【やろうとしていること】~
20個程度の説明変数があり，その説明変数を利用して，目的変...
~
【現状】~
現在，R2.2.1を，Windows XPで利用中です．~
線形判別分析，ロジスティック回帰分析の変数選択に関してはA...
~
【お聞きしたい点】~
お聞きしたのは，２点です．~
~
１．ニューラルネットや分類木で，線形判別分析で提供されて...
~
２．そもそもニューラルネットや分類木で0/1の２値を予測する...
//
-質問２は R プロパーの話題ではありませんね。統計プロパー...
-回答しないのに、わざわざへこませるようなこと書かなくても...
-ほっとくよりいいんじゃない？回答者の意気をくじくようなこ...
聞く人は多いけど，答える人は少ないよ。答える人がいなくな...
2. については，変数選択というか予測に重要な変数を選択する...
既存の仕組みがないと言うのは障害かもしれないが，やり方は...
そもそも，重回帰の変数選択の仕組み自体を評価しない人もい...
総当たり法も，変数が20もあると無理だけど，有望そうな変数...
理論根拠に基づき，試行錯誤的に変数選択するのも良いのでは...
**rglがOSXでコンパイルできない [#g1afb318]
>[[高井]] (2006-05-30 (火) 12:30:55)~
~
R 2.3.0のフルセットをMac OSX 10.4.6, Xcode 2.3，PowerMac ...
~
api.cpp: In function ‘void rgl_user2window(int*, int*, d...
api.cpp:608: error: invalid conversion from ‘int*’ to ‘c...
api.cpp:608: error: initializing argument 6 of ‘GLint ...
const GLdouble*, const GLdouble*, ...
api.cpp: In function ‘void rgl_window2user(int*, int*, d...
api.cpp:633: error: invalid conversion from ‘int*’ to ‘c...
api.cpp:633: error: initializing argument 6 of ‘GLint ...
const GLdouble*, const GLdouble*, ...
GLdouble*, GLdouble*)’
make: *** [api.o] Error 1
chmod: /Users/takai/Library/R/library/rgl/libs/ppc/*: No...
ERROR: compilation failed for package 'rgl'
//
-608行，633行を修正すればよいだけでしょう（誰でもソースを...
-CRANからソースをダウンロードして，ダブルクリックで解凍し...
-[[R-develメーリングリストのアーカイブの中:http://tolstoy...
-ちゃんとやれば，誰にでもできるはずなんですけどね。~
レスポンスがないようなのですが，以下は，ターミナルで操作...
中澤さんのコメントを参考に，まず解凍された rgl ディレクト...
601a602
> GLint viewport[4]; 追加
605a607
> for (int i=0; i < 4; i++) viewport[i] = view[i];...
607c609
< gluProject(point[0],point[1],point[2],mo...
---
> gluProject(point[0],point[1],point[2],mo...
624a627
> GLint viewport[4]; 追加
628a632
> for (int i=0; i < 4; i++) viewport[i] = view[i];...
632c636
< gluUnProject(pixel[0],pixel[1],pixel[2],...
---
> gluUnProject(pixel[0],pixel[1],pixel[2],...
その後
$ R CMD INSTALL rgl
でしょうか？？ -- &new{2006-05-30 (火) 18:24:22};
-上記のお二方のコメントに従って再トライしました。api.cpp...
collect2: ld シグナル 6 [Abort trap] で終了させられました
/usr/bin/ld: warning internal error: output_flush(offset...
flushed block(offset = 464972, size = 3068)
calling abort()
make: *** [rgl.so] Error 1
chmod: /Library/Frameworks/R.framework/Versions/2.3/Reso...
No such file or directory
ERROR: compilation failed for package 'rgl'
となってコンパイル失敗です。一歩前進しましたがエラーメッ...
-R CMD INSTALL rgl を実行したときの pwd が，rgl ディレク...
-ああ，「上位のディレクトリでR CMD INSTALL rgl」とあった...
-/usr/bin/ld の internal errorが出るのは, RかDeveloper To...
-岡田さん，rglのバイナリーありがとうございます。無事にイ...
-3時間くらいかかって，本家からDevelopper Tools 2.3をダウ...
-環境依存の問題であることはほぼ明らかになりましたから、あ...
-岡田さんに励まされ（？），ネジ捲かれ（？）ハードウェアチ...
-どうやらgccがおかしいんじゃないかと当たりをつけました。/...
**文字を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")) -- [[t...
-x が実際にどうなっているのかわからないのだけど。他の人が...
x <- letters
y <- rnorm(length(x))
plot(as.numeric(factor(x)), y, type="b", xaxt="n")
axis(1, at=factor(x), labels=x)
老婆心ながら，このような場合には折れ線グラフは不適切と，...
-ああ、これ苦労した記憶があります。結局出来なかったんです...
-これだとlabels内文字列でX軸の順番がソートされてしまいま...
-y<-c(1,2,3,4);plot(1:length(y),y,xaxt="n");axis(1,1:leng...
**スクリプトでread.delimを使うとエラーになってしまいます...
>[[KT]] (2006-05-25 (木) 17:52:50)~
~
コンソールから~
data <- read.delim("data.txt", header=T)
と入力・実行すると正しく実行されるのに、スクリプトファイ...
source("test.R")
と実行すると、~
以下にエラーeval.with.vis(expr, envir, enclos) :
関数 "raed.delim" を見つけることができませんでした
とエラーになってしまいます。対策をご教授いただければ幸い...
//
-x raed.delim / o read.delim -- [[takahashi]] &new{2006-0...
-takahashi様、コメントありがとうございます。ですが、超初...
-あのですね。「つづりが間違えてますよ」と指摘されたわけで...
-scriptの方がtypoなんじゃないですか？ x raed / o read -- ...
-まさにスペルミスでした。何たるケアレスミス。どうもありが...
**英語版R1.3の入手法は？ [#a767b2b7]
>[[Yabo]] (2006-05-25 (木) 08:32:16)~
~
R1.3を本家からダウンロードしても、インストールすると日本...
//
-R1.3って(@.@)何？ -- [[okinawa]] &new{2006-05-25 (木) 09...
-えっと
set LANG=C
set LC_ALL=C
"*****?jgr.exe" #JGRのパス
というバッチファイルを作って，そこから起動すればO.K.です...
-R3.1.0 のことでしょう。未来のお話です(^_^;) というのは...
-32bitかつUTF-8のロカールを持つOS上なら日本語でも動きます...
-ｇさん、ご教示の方法でうまくいきました。ひとりぼっちでや...
-R1.3はR2.3の誤りでした。すんません。 -- [[yabo]] &new{20...
-gさんのバッチファイルの set LC_ALL=C のCをUTF-8にしても...
**条件が長さが2以上なので [#wb1cb369]
>[[レッド]] (2006-05-23 (火) 04:30:21)~
~
次ののメッセージがでて条件式が効きません．
TSi <- amp$Si
> if (TSi < 8) {
+ TAl <- 8 - TSi
+ } else {
+ TAl <- 0
+ }
Warning message:
条件が長さが2以上なので，最初の一つだけが使われます in: ...
超初心者で，Win．Version 2.2.1 ~
//
-if()を調べましょう。~
cond: A length-one logical vector that is not 'NA'. C...
length greater than one are accepted with a war...
the first element is used.
if(x)ではis.logical(x)がTRUE、かつis.na(x)がFALSEかつ、le...
ただ、この条件を満たしているのでしたら、amp$Siがどういっ...
エラーが出る場合は、いきなり応用問題はやめて基礎からはじ...
-TSi がベクトルなんでは？ c(1,3) < 2 は論理ベクトル c(TR...
-このページの上の方の注意事項を読んでから質問・投稿しまし...
amp$Si がベクトルなんでしょうね。だとしたら，TAl <- ifels...
ifelse 関数を調べましょう。 -- &new{2006-05-23 (火) 10:4...
-amp$Siはベクトルです． ex.amp$Si <- c(8.123, 7.897, 7....
**平均と標準偏差を使ったboxplotは描けませんか？ [#r2b5ff6e]
>[[くに]] (2006-05-22 (月) 18:35:13)~
~
箱ひげ図には平均値と標準偏差を用いたものもありますが、こ...
//
-平均値・標準偏差を求める関数があり，plot, lines, points,...
-ありがとうございます。できたら、具体例か参考になるウェブ...
-可能ですが，平均値と標準偏差を使うなら，ストリップチャー...
-可能かどうか聞かれたので「可能ですと」お答えしましたが。...
何だってできますから，どこまでやらなきゃならないかを決め...
元のデータを記号で書き加えたり，そのほかなんでもあなたの...
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...
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", xl...
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 (探索的データ解析) の精神」とはどのような...
-「平均と標準偏差を使ったboxplotを描くにはどうすればよい...
-「EDA (探索的データ解析) の精神に反する」、というのはデ...
**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",k...
ここで、エラーメッセージが次のようにでます。~
以下にエラーmatch.arg(kernel, c("rbfdot", "polydot", "ta...
'arg' は以下の一つでなければなりません： rbfdot, ...
besseldo...
stringdotが選択肢にないようです。
パッケージのマニュアルには、stringdotは選択肢にあり、ソー...
ものがあることは確認しました。それともデータセットの与え...
//
-よく知らないのだけど，ヘルプを見る限り，
## S4 method for signature 'list':
ksvm(x, y = NULL, type = NULL, kernel = "stringdot",
:
x a symbolic description of the model to be fit. When n...
can be a matrix or vector containg the training data
or a kernel matrix of class kernelMatrix of the train...
character vectors (for use with the string kernel).
なんだから，x はリストでないといけないのでは？ -- &new{2...
-下記のように試してみました。 -- [[北陸人?]] &new{2006-05...
> 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....
ここで、エラーメッセージが次のようにでます。
以下にエラーas.double.default(t(x)) : (list)オブジェク...
ラベルの与え方がまだ理解できていませんが、ラベルyとすると...
指定しないということでしょうか？
-なかなか，コメントがつきませんね。「わたしは，その方面は...
-そうですね。string kernelが正確に使えるSVMは少ないですね...
**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, kern...
#SVMでの訓練データの学習
ここで
以下にエラーmodel.frame.default(data = list(V1 = c(33,19...
とエラーメッセージがでます。~
確認してみますと
> 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 ...
というスペース区切りの数値データです。~
どういう操作が足らないのかご教示いただけませんでしょうか？~
//
-問題を再現できる必要最小限のデータと共に，再度質問される...
-そのデータは１行しかなくて、行列型に旨く変換されなかった...
-説明不足で申し訳ありません。 -- [[北陸人]] &new{2006-05-...
データセット(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="rbfd...
以下にエラーmodel.frame.default(data = list(V1 = c(1, 0)...
:オブジェクトが行列ではありません
> dim(test1m)
[1] 3 3
> class(test1m.train)
[1] "matrix"
> dim(test1m.train)
[1] 2 3
上記のようになりました。ksvmにかけたときのオブジェクト(te...
-第一引数のtypeがデータの列名に見つからないからじゃないで...
-library(kernlab)であることを書いて頂いた方が早く返事が来...
-どんぴしゃり解決しました。アホでした。どうもありがとうご...
**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....
共有ライブラリ '/Library/Frameworks/R.framework/Resourc...
を読み込めません
dlopen(/Library/Frameworks/R.framework/Resources/libra...
Library not loaded: /usr/X11R6/lib/libX11.6.dylib
Referenced from: /Library/Frameworks/R.framework/Resou...
Reason: image not found
追加情報: Warning message:
ライブラリ ‘/Users/doima/Library/R/library’ はパッケージ...
エラー：.onLoad は 'tcltk' のための 'loadNamespace' に失...
エラー：パッケージ 'tcltk' をロードできませんでした
//
-Macを使っていないので予想ですが、Tcl/Tkがinstallされてい...
-R Commanderには"X11"が必須ですが、これは標準ではインスト...
**Rで使えるデータ一覧 [#y1977666]
>[[究極超具あーる]] (2006-05-12 (金) 12:57:54)~
~
Rで例えば、~
library(stats)
data(morley)
とすると、データmorleyが使えるようになりますよね。このよ...
//
-データセットmorleyはstatsパッケージではなくbaseパッケー...
-morley などは，data(morley) などとしなくても，使えます。...
-インストール済みの全パッケージの全データ一覧を知りたいで...
data(package="パッケージ名")
とするとパッケージ中のデータ一覧が得られることがわかりま...
ただ、この中のItem列を抜き出そうとすると、以下のようなエ...
datas = data(package="accuracy")[["results"]];
datas[["Item"]]
以下にエラーdatas[["Item"]] : 添え字が許される範囲外です
datasのデータ構造が良くわかりません。Item列を抜き出すには...
-datas[,"Item"]とするか、data(package="パッケージ名")[["r...
local({pkg <- select.list(sort(.packages(all.available =...
+ if(nchar(pkg)) library(pkg, character.only=TRUE)})
の部分を書き換えると対応できるのではないかと。 -- &new{2...
-単純ミスでした。ありがとうございました。 -- [[究極超具あ...
-次のようにすると、究極超具あーるさんの望みは叶うと思いま...
x<-sort(.packages(all.available = TRUE))
lspkg <- function(pkgname){
data(package=pkgname)[["results"]][,"Item"]}
sapply(x,lspkg)
データを含んでないパッケージまで一々character(0)と返すの...
**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,...
plot(y)
フィルタhがうまく動けば、yは(t,y) = (0,1), (10,9)を通る直...
アドバイスをいただけると幸いです。よろしくお願いいたしま...
//
-それぞれの行ごとに，x, y がどのようになっているか，確認...
x は全部1なんですけどいいのかなぁ。 x = ts(1:1000, freq=1...
-典型的なまずい質問の例ですね。自分のしたいことを人にもわ...
-まずいですか？「フィルタをfilterでどう実装するか」という...
-あなたの説明はフィルターやら積分器とやらをすでにある程度...
-お勉強の結果（1）~
filter の引数 method が "convolution" のときは，「（重み...
二番目の引数は項数ごとの重み，つまり，n 項移動平均なら re...
ただし，本当は移動和なので，重みは等しくなくても良いし和...
使用例 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", sid...
> a2 <- filter(x, rep(1/3, 3), method="convolution", sid...
> a3 <- filter(x, c(2,0.5,1), method="convolution", side...
> 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 と同じ結果になる（当然だ...
重みの長さが1で値が2なら，a4[1] = x[1], a4[2] = x[2]+2*a4...
重みの長さが2以上で，値が違っても同様。a5[4] = x[4]+1*a5[...
> 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
ということですか。~
どうも，書かれたプログラムはあなたがやろうとしたこととは...
-お勉強の結果（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
以下略
プログラムに書かれたとおりの動きをしているようだ。 -- &n...
-お勉強の結果（4）~
あの〜。もとの発言にあった z^-1 というのは，「遅延素子」...
そうだとすると，関数の定義上それはあなたの意図とは違った...
-お勉強の結果（5）~
遅延素子かどうかはともかくとして，あなたの最初の質問で，...
もしそうだとすると，
> y <- filter(filter(x, 0.01, method="con"), 1, method="...
でいいのでは。つまり，外側の filter の第2引数が -1 ではな...
#ref(filter.png)
-いろいろありがとうございました。おっしゃるとおり、ｚ^-1...
**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....
> 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 = ",", co...
> write.table(y, file = "svmsample_y.csv", sep = ",", co...
> rm(list=ls(all=TRUE))
> x<-read.table("svmsample_x.csv", header = TRUE, sep = ...
> y<-read.table("svmsample_y.csv", header = TRUE, sep = ...
> iris1 <- ksvm(x,y,kernel="anovadot",kpar=list(sigma=0....
以下にエラーksvm(x, y, kernel = "anovadot", kpar = list(...
no direct or inherited method for function 'ksvm...
> table(y,predict(iris1,x))
エラー：オブジェクト "iris1" は存在しません
以下にエラーpredict(iris1, x) : unable to find the argum...
selecting a method for function 'predict'
~
読み書きでデータ型が変化したわけでもないようです。~
申し訳ありませんがご教示宜しくお願いします。~
//
-自己解決しました。失礼しました。 >x<-as.matrix(read.tab...
-read.table 関数はデータフレームを返す一方、ksvm は第一引...
**成長曲線の当てはめ [#ja28eabb]
>[[たつ]] (2006-04-29 (土) 15:09:09)~
~
y(売上高）とx(キャンペーン額)の関係を成長曲線に当てはめた...
以上の目的に適当なパッケージはありますでしょうか？~
その際に売上高の上限と下限を指定できたりすると好都合なの...
~
アドバイスよろしくお願いします。~
//
-glm関数で出来るかもと思い色々と試行錯誤もしましたが、自...
-まず RSiteSearch("growth curve") を実行し、ヒットしたた...
-ヒットしました･･･覚束ない英語力でチェックしていたところ...
お陰さまで解決できたのですが最後に質問をさせてください。
上限と下限を指定しない場合はglmでも可能なのでしょうか？--...
-glm は非正規誤差にたいする線形回帰手法ですから、成長モデ...
~
SSlogis
SSgompertz(stats) Gompertz Growth Model
SSweibull(stats) Weibull growth curve model
-素人の私にはSSlogisが使い勝手が良さそうです。これから詳...
**最小絶対偏差回帰 [#r1f44837]
>[[たろう]] (2006-04-26 (水) 20:19:04)~
~
R初心者ですが、Splusをしばし使用しておりました。~
Splusでは、最小絶対偏差回帰は~
l1fitですが、Ｒで対応する関数を見出すことできずにおります...
ご存知のかた、ご教示宜しくお願いします。 ~
//
- r-help ML の過去記事を検索したらつぎのような(二件の別個...
> is there any package or method which enables to comput...
> 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 ...
-より頑健な回帰手法（おそらく単なる l1fit よりベター？）...
-ご教示ありがとうございます。おかげさまで無事できました。...
-ちなみに二件の記事は Rsitesearch でそれぞれキーワード "l...
-なるほど。RSiteSearch、初めてしり早速、確認いたしました...
-「より頑健な回帰手法」の実行方法についても、お教え頂きあ...
**部分一致検索による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)の内容も良く分かりま...
すみませんが、ご教示宜しくお願いします。~
//
-以下に示すように，
in1[grep('AA', as.character(in1$Name)),]
でできますが、これだと、"BAA" があれば、それも摘出してし...
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 (水) 0...
-ちなみに「'subset' は論理値として評価されなければなりま...
-さらさん、青木先生、どうもありがとうございました。as.cha...
-行頭を表す正規表現 ^ 。 そんな便利な物があったんですね...
-Rjpwikiにも正規表現の詳しい説明がある -- &new{2006-04-2...
**常微分方程式系のバラメータ推定をするには？ [#le048e5a]
>[[使用1ヶ月目]] (2006-04-19 (水) 11:44:48)~
~
常微分方程式系のバラメータ推定を考えています。常微分方程...
//
-それぞれの関数の help に example がありますでしょ？その...
-いくら教えてる立場だとは言え、全く関係ない人が見ても不愉...
**データファイルを読み込んでcor.testをする時のエラー [#f6...
>[[tensho]] (2006-04-17 (月) 23:26:15)~
~
基本的な質問で申し訳ありません。~
データファイルを読みこんでcor.testを行うと、「長いオブジ...
//
-文字とおりの意味では？これ以上の回答を期待するなら、読み...
-データ入力は以下の通りです。
> 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")
-データフレームの意味を誤解していませんか。行がケース、列...
-ご親切に教えていただき、どうもありがとうございました。 -...
**GLMでの決定係数 [#n25e8784]
>[[ごう]] (2006-04-11 (火) 06:30:08)~
~
GLMの基本的な部分かもしれないのですが、情報をいただけませ...
正規分布を前提とした普通の重回帰モデルでは決定係数が算出...
近似した値は「誰々（人名）」の決定係数といった雰囲気で算...
正規分布重回帰モデルでもサンプルサイズが大きければ有意な...
重複投稿ご容赦願います。他のサイトでも質問したのですが...
//
-＞ 近似した値は「誰々（人名）」の決定係数といった雰囲気...
できるのかもしれないのですがとは，できるんですか，できな...
誰々の決定係数というのは，あるんですか，ないんですか？~
あるのなら，具体的に何なんですか？~
＞ これはロジスティック回帰を見る限りではどうも低く算出さ...
気がするだけなんですか，実際に低いんですか？~
その判断は，あなたがしたものですか，そのように書いてあっ...
＞ モデルが成り立つかどうかで判断するというような記載を読...
気がするだけなのですか？ -- &new{2006-04-11 (火) 08:55:2...
-不明瞭な質問でした。決定係数の例：Nagelkerke R^2, Cox&Sn...
**bartlett.test関数の中身を見るには？ [#j47456dd]
>[[使用1ヶ月目]] (2006-04-10 (月) 17:17:00)~
~
通常、関数の中身を見るには、関数名をそのまま入力すればよ...
//
-このスレッドの下の方に同じような質問がありますよ。 -- &...
**尤度のプロットについて [#l09ecac0]
>[[にゃあ]] (2006-04-08 (土) 14:26:34)~
~
例えば二項分布binomの確率分布はdbinomとして計算できますが...
また尤度のプロットというのも可能なんでしょうか?~
//
-尤度の定義通りに計算すればよろしいのでは。ただし、普通は...
N <- 20
X <- rbinom(100, N, 0.55) # データ
loglikelihood <- function(p) {
log.prob <- dbinom(0:N, N, prob=p, lo...
LL <- 0
for (x in X) LL <- LL + log.prob[x+1]...
# もしくは
# for (i in 0:20)
# # LL <- LL + sum(X[X=i])*log.prob...
# LL <- LL + sum(X==i)*log.prob[i+1...
return(LL)
}
-後で消去しますが，というか，もし指摘が正しいなら，記事修...
-ベルヌーイ尤度にベータ分布を掛ける作業を視覚化してみたか...
関数を作ってみたんですが、以下の様になりました。~
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(likel...
plot(theta, posterior, ylim=c(0, 1), type="l")
distribution <- posterior
}
}
かなり雑なんですが、見よう見まねでやってみました。~
しかし、実行速度にかなり問題があります。~
どなたかいいアイディアがありますでしょうか?-- [[にゃあ]] ...
-関数を書いて，それに対する意見を求めて，回答があったとし...
-なるほど。言われてみればその通りかもしれません。 -- [[に...
-毎回プロットしていれば時間はかかるし、途中の結果も見られ...
-そうですね。ただし、根本的な考え方の間違いに気付いたので...
-関連して気になったこと。非負整数値データの度数分布ベクト...
> 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 ...
> 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-...
-参照先を示してあげるか，再度例示してあげると，おんぶにだ...
> 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)) ですか。ありが...
-ですよね(^_^;)要望してみればいかがでしょう。 -- &new{20...
**ライブラリ [#i1aeaca4]
>[[超初心者]] (2006-04-03 (月) 15:18:03)~
~
Package source: evd_2.1-7.tar.gz ~
Windows binary: evd_2.1-7.zip ~
ライブラリには、上記のように、拡張子がgzのものと、zipのも...
//
-拡張子より，その前の情報を見ましょう。「Package source」...
-R(最新2.2.1)では、CRANサイトにあるパッケーッジは動作する...
何卒よろしくお願いします。 -- [[超初心者]] &new{2006-04-...
-パッケージの添付ファイルに，その情報が書かれているものが...
あなたも,
[[パッケージの説明:http://cran.md.tsukuba.ac.jp/src/contr...
-度々申し訳ありません。~
ライブラリの拡張子gz「Package source」はソースプログラム...
（zipの場合は、Rメニューから呼び出し後インストール機能が...
-ソースパッケージをビルドするにはRをビルドするのとほぼ同...
-Windows版のRをお使いなら、パッケージのインストールはメニ...
-早速ご指導深謝します。舟尾様のページを見ました。 -- [[超...
-なかま様早速ご指導深謝します。~
状況としては、ライブラリの拡張子gz「Package source」がＰ...
ＲはＰＣにインストールされています。~
この場合、ソースをＲからインストールする場合
>install.packages("c:/***.gz",CRAN=NULL)
では、パッケージ***がインストールできないのでしょうか？~
誠に申し訳ありません。 -- [[超初心者]] &new{2006-04-07 (...
-Rtools.zipと言うのと,Mingw(C,Fortranコンパイラ)を入れて...
-[[超初心者]]さんがインストールしたいのは[[evdパッケージ:...
-超初心者さんなら，バイナリをインストールすればいいでしょ...
-礼儀作法がわからずに誠に申し訳ありません。具体的に申し上...
-上のなかまさんのコメントの「[[FAQのリンク先:http://www.m...
-[[ビルドのログ:http://cran.us.r-project.org/bin/windows/...
#include <sys/types.h>
/* --> int64_t ; if people don't have the above, they ca...
/* #include "int64.h" */
とか言ってるので, <sys/types.h>では無く<inttypes.h>にすれ...
-#include "int64.h" では？で，sys/types.h はある場合の方...
-[[データサイズ:http://www.unix.org/version2/whatsnew/dat...
-カレントディレクトリにある *.h は " でくくるんじゃなかっ...
-それはそうです. ただ int64.hの中で何を定義しているかと言...
-ライブラリ(zip)でインストール時にエラーが発生します。
何卒よろしくご指導をお願いします。
一例をあげますと、copula_0[1].3-5~
http://cran.r-project.org/bin/windows/contrib/r-release/...
エラーメッセージは
> utils:::menuInstallLocal()
以下にエラーgzfile(file, "r") : コネクションを開くことが...
追加情報: Warning message:
圧縮されたファイル 'copula_0[1].3-5/DESCRIPTION' を開く...
その他、e1071_1[1].5-13、やkernlab_0[1].6-2もインストール...
これらzip形式ライブラリのインストール方法をご教示下さい。...
-まだ質問の仕方がおわかりになっていないようですね。関係な...
-返す言葉もありません。正直なところ該当箇所を見ましたが、...
-ご指導頂いたように、アーカイブの「プロキシーサーバー」を...
-options 関数と Sys.putenv 関数は同じ行に入力したのですか...
-2段に別々に入力しました。本件に関連しているかどうかわか...
-何の ID とパスワードなんでしょう。どのページを見ても要求...
**数値の置き換え [#g6c16a4c]
>[[青木繁伸]] (2006-03-31 (金) 23:03:18)~
~
整数ベクトル x を，別の数値に置き換える関数を作る必要があ...
条件はn行3列の行列で与えようと思います。xの要素が1列目以...
> 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]...
xが次のようなとき~
> x <- c(7, 4, 3, 1, 8, 6, 5, 5, 8, 1, 1, 10, 7, 6, 3, 5...
答えが~
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()の定義の長さを考えると（あちら...
-コメント，ありがとうございます。私の悪い癖で，探すより作...
-青木さんの例のように置き換える数値が単調増加正整数であれ...
> x <- c(7, 4, 3, 1, 8, 6, 5, 5, 8, 1, 1, 10, 7, 6, 3, 5...
# 1,2,3 は第一の区間 [1,3.5)、4,5,6,7 は第5の区間 [4,7.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 ...
> 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 # これが何を出力するかとかも
-どうもありがとうございます。助かりました。 -- [[内藤]] &...
**図を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...
par(new=T)
plot (x,z, ylim=c(0, 10000), ann=F, pch=4)
上記の図を描いた後，ファイル−別名で保存より，pdfを選択す...
エラー：invalid character sent to 'PostScriptCIDMetricIn...
追加情報: Warning messages:
1: invalid string in 'PostScriptStringWidth'
2: invalid string in 'PostScriptStringWidth'
3: invalid string in 'PostScriptStringWidth'
何が原因でしょうか？基本的な質問で大変申し訳ありませんが...
//
-?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...
par(new=T)
plot (x,z, ylim=c(0, 10000), ann=F, pch=4)
dev.off()
初心者ならどんな基本的なことを質問しても良いということで...
-Akira様．初心者STです．準備不足でいろいろとご迷惑をおか...
-今回初めて pdf ファイルとして書き出そうとしたのでしょう...
-R version 2.2.1, 2005-12-20, i386-pc-mingw32 512M WXPSP2...
-R version 2.2.1, 2005-12-20, i386-pc-mingw32 256MB WXPSP...
-多分インストール時に東アジアを選択しなかったのでは？(こ...
-青木様，okinawa様，Akira様，なかま様．初心者STです．再イ...
**大量のグラフを描くためには [#a6a1ddb0]
>[[antt]] (2006-03-14 (火) 00:08:58)~
~
はじめまして。大量のグラフを描きたいのですが，以下に記し...
【環境】R version 2.2.1, i386-pc-mingw32, メモリ 1GB~
【問題】大量のグラフを描き続けると，オブジェクトは別に増...
【知りたいこと】大量のグラフを描きたい場合に，どのような...
【プログラム例】
alternative <- 1:5 #選択肢の数
nitem <- 100 #項目の数
nschool <- 10 #施設の数
n <- 10000 #標本サイズ
##仮想データを作成
data <- data.frame(matrix(alternative, ncol=nitem, nrow=...
data$school <- rep(1:nschool, each=n/nschool) #1施設 100...
##画像作成用の関数を定義する
plotitem <- function(i, j) {
dir.create(paste("school", i, sep=""), showWarnings=...
filename <- paste("./school", i,"/item", j, ".ps", s...
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)
}
}
//
-がんばって整形してくれましたが、無意味に長すぎるのも読む...
-R version 2.2.1, 2005-12-20, i386-pc-mingw32 メモリ512M...
-闇版ならFMの処理のリークです.setHook(packageEvent("grDev...
-日本語をどうしても使いたい場合は, Postscriptを吐くならTe...
-ご回答ありがとうございました。仰るとおり闇版を使っており...
-フォントのメトリック情報の処理によってリークが発生すると...
-リーク => メモリーリーク => 使われたまま回収再利用不可...
-はい。 -- [[なかま]] &new{2006-03-17 (金) 00:08:31};
**imageの色を絶対指定するには [#i140c8f9]
>[[子牛]] (2006-03-13 (月) 19:39:19)~
~
はじめまして。imageで得た画像にcolをつかって色を指定した...
数値にたいして絶対的な方法で色を指定するにはどうすればよ...
~
> 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" "uti...
[7] "base"
//
-100段階のパレットで一段階の違いを表現するという例題は，...
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...
#ref(pict2.png)
-ありがとうございました、色の使う範囲をこれで指定できるの...
ただ、教えていただいた方法ですと、オフセットを手で調整せ...
image(x,y,a, zlim=c(-2,2), col=cm.colors(110), main=past...
この方法で絶対的な値を固定できるようです。（すみません、...
**SVMのマージン距離出力方法について [#o5f5dbe6]
>[[ちょめ夫]] (2006-03-13 (月) 12:10:44)~
~
現在SVMに関するパッケージの検討を行っております。~
検討対象としているパッケージは以下の２つです。~
１）e1071~
２）kernlab~
なお、判別したいクラスは2クラスの分類となっております。上...
http://www.msi.co.jp/vmstudio/materials/tech/classificati...
「図 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" "ut...
other attached packages:
colorspace kernlab e1071 class
"0.9" "0.6-2" "1.5-11" "7.2-16"~
//
-投稿法についてご指摘いただきまして、誠に申し訳ありません...
-当てずっぽうですが、もし話題の量が SVM アルゴリズムにと...
終了行:
COLOR(green){SIZE(20){初心者のための R および RjpWiki に...
新規投稿はできません
----
-[[初級Ｑ＆Ａ アーカイブ(4)]] (元記事が 2005-11-09 より 2...
-[[初級Ｑ＆Ａ アーカイブ(3)]] (元記事が 2005-05-02 より 2...
-[[初級Ｑ＆Ａ アーカイブ(2)]] (元記事が 2004-12-13 より 2...
-[[初級Ｑ＆Ａ アーカイブ(1)]] (元記事が 2004-08-03 より 2...
----
#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...
-ありがとうございました。sinkでファイルに書き込むことがで...
-> 行列を一列にするプログラムを書きました~
プログラムを書くって，as.vector(行列オブジェクト）だけで...
それと，ファイルに取り込むの？書き込むの？「関数の値」っ...
-適当なエディタを開き，画面（コンソール）に表示されたもの...
**R-2.3.1のインストール手順10のファイルの上書き [#m9ece0ae]
>[[初心者]] (2006-06-24 (土) 12:03:31)~
~
インストールの手順の最後にRconsole，Rdevga，Rprofile.site...
//
-ファイル名のところにリンクが張ってありません？？ -- &ne...
-インターネットでは（そのほかでも），半角カタカナは使わな...
-ファイル名の所を押しても文字が出てくるだけなのですが？ -...
-MIMEタイプの設定がうまくできていないサーバーなんでしょう...
このような場合の対処法は，覚えておくと，今後得しますよ。~
そのような場合は，うんとね，ういんどーずのようですから，...
**関数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{...
**関数名を文字列に [#l2aaf4df]
>[[ショーン]] (2006-06-22 (木) 19:50:20)~
~
以下のように関数オブジェクトを変数に入れて、関数として呼...
> a = cos
> a
.Primitive("cos")
> a(0)
[1] 1
後でaの元の関数名である、"cos"という文字列を得る方法はあ...
> f2s(a)
[1] "cos"
のf2sようなことはできますでしょうか？~
//
-出来ないことは無いですが、まずなぜそんなことが必要なのか...
-理由がないと教えられないと言うことでもないようには思いま...
もっとも，私にはわかりませんでした。~
変数に関数定義を付値したら，その変数の class は function ...
> 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"
-答えるにはそれなりの手間がかかりますから、必要性が私には...
-勉強の機会を与えてくれて，ありがとうございます。後者の場...
> 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"
-うむ、あなたは初級者などでは決してない。何につかうのか秘...
f2a <- function(x) {capture.output(x ,file=(filename <- ...
-私は，原質問者じゃないよ。~
だから，原質問者は用途が秘密だとは言っていない。~
質問者をもてあそぶヒマはあっても，答えるヒマはないという...
-暇のある無しの問題ではありません。質問するからには、やは...
-回答ありがとうございます。おかげさまでうまくいきそうです...
funcs = c(cos,sin,tan)
のようにして、ループでいろいろな関数のグラフを描こうと思...
-2006-06-22 (木) 23:33:51 は，もっと簡単になるでしょう？ ...
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() は関数中でつかわれ、実引数の数...
func <- function(...) print(nargs())
-なるほど。~
純粋に ... の個数を知りたい場合には，以下の方がよいのかな...
> func = function(a, ...){print(length(c(...)))}
> func(1,2,3)
[1] 2
-【「Rの関数定義の基本」を見るとlength(...)という書き方も...
-ショーンさんはlength(1:3)とlength(c(1,2,3))、length(1,2,...
-COLOR(blue){【「Rの関数定義の基本」を見るとlength(...)と...
> 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(...)) でし...
-ありがとうございました。よくわかりました。 -- [[ショーン...
**データ行列の行をまたぐレコード参照 [#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と比較演算子を組み合わせればできそうに思いますが -...
-例に挙げたデータについての例解が間違えているし。条件の記...
-まずは，val1 < 35.0 という条件なら，それを満たさないもの...
-しかしまあ，val1 < 35.0 は，第 i 行に関する条件なのか，...
-いろいろ面倒だが，骨格を作ったから直してみて？val1 につ...
たぶん遅いから，R じゃなきゃいけないとか R の機能をふんだ...
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() 調べてみます&ご提示頂いた...
その他の皆様も，アドバイスありがとうございます(たしかに例...
-上のコードと本質的に同じコードで50万行の人工例を（一部実...
-一昼夜コンプータを動かしておけば解が出るなら，安いモンで...
**Williamsの多重比較のパーセント点 [#p8461227]
>[[tm]] (2006-06-20 (火) 22:23:00)~
~
Rとは直接関係なく、大変申し訳ありませんが、Williamsの多重...
身近で入手できる本をあさったのですが、5%点くらいしか手に...
できれば、1%点の表が欲しいのですが、どなたかご存知ありま...
//
-その本に，生成アルゴリズムなどは書いてなかったですか。。...
**.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>
の実際のコードを覗いてみることはできないのでしょうか。~
//
-過去に同じ類の質問が沢山あります。ちゃんと調べましょう。...
-.Internal については，あったかな？ -- &new{2006-06-20 (...
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の値が同じになると思っていたのですが、違います。ど...
//
-主成分負荷量を求めるのでしょうか？だとしたら，貴方は，説...
> d <- structure(c(0.1991, 0.2003, 0.2169, 0.1429, 0.184...
+ 0.1516, 0.1478, 0.2071, 0.1613, 0.2288, -15.0277, -16....
+ -17.3694, -16.7424, -15.0885, -15.7661, -15.1187, -14....
+ -13.5147, -15.0654, 0.279, 0.3462, 0.299, 0.2586, 0.25...
+ 0.155, 0.294, 0.3411, 0.2789, 0.2592, 0.2165), .Dim = ...
+ ), .Dimnames = list(c("1", "2", "3", "4", "5", "6", "7...
+ "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/PC...
res2$rotation
が行列Lですよね。それから、
res2$x
が主成分zですよね。投稿した質問のように、元の変数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
-ありがとうございました。よくわかりました。 -- [[ショーン...
**Sweave環境でsink()は使えないのでしょうか？ [#fda5d2d5]
>[[akira]] (2006-06-16 (金) 14:11:29)~
~
Sweave環境でxtableを書き換えたい場合、~
(1)xtableの出力先をtexファイルからsink()で別のファイル（...
(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(...
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がないからxta...
'RweaveLatex', 'Rtangle'当たりを見てもそれらしい記述が見...
//
-[[print.xtable2:http://sapporo.cool.ne.jp/matsut/r02.txt...
**multcomp パッケージでの Williams 検定について [#pa7f07ff]
>[[tm]] (2006-06-15 (木) 21:10:37)~
~
multcomp パッケージ中の simtest の Williams 検定について...
Rのバージョンは 2.0.1、multcomp のバージョンは 0.4-8、mvt...
> library(mvtnorm)
> library(multcomp)
> data(angina)
> summary(simtest(response ~ dose, angina, type="William...
Simultaneous tests: Williams contrasts
Call:
simtest.formula(formula = response ~ dose, data = angina...
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...
Simultaneous 95% confidence intervals: Williams ...
Call:
simint.formula(formula = response ~ dose, data = angina,...
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 検定については情報...
使用方法や結果の見方などご教授いただけないでしょうか？~
よろしくお願いいたします。~
//
-あるプログラムにバグがあるというとき，確証がないといけま...
プログラムの使い方などについては，教科書などに載っている...
> library(mvtnorm)
> library(multcomp)
> set.seed(12345678)
> df <- data.frame(g = factor(rep(1:6, each=20)), y = rn...
> summary(simtest(y ~ g, df, type="Williams", alternativ...
Simultaneous tests: Williams contrasts
Call:
simtest.formula(formula = y ~ g, data = df, type = "Will...
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...
Simultaneous 95% confidence intervals: Williams
contrasts
Call:
simint.formula(formula = y ~ g, data = df, type = "Willi...
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』に見えるのは出力の有効桁数だけの問題では？全ての...
-ご返信、アドバイスありがとうございます。 ~
まず、いいたかったのは「バグがあるよ。」ということではな...
おっしゃられるように、正解例がないとよくわかりませんので...
> library(mvtnorm)
> library(multcomp)
> set.seed(12345)
> df <- data.frame(
+ Dose = ordered(rep(c("Cont", "Low", "MidLow", "MidHigh...
+ 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="...
Simultaneous tests: Williams contrasts
Call:
simtest.formula(formula = Response ~ Dose, data = df, ty...
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 の値が...
この違いは、私の R 環境の問題でしょうか？それとも、実装さ...
ご教授いただけますよう、よろしくお願いいたします。
-- [[tm]] &new{2006-06-20 (火) 10:32:31};
-Contrast マトリックスと件の本を比べると違いがわかるので...
-contrast マトリックスの意味がなんとなく理解できました。...
ただ、おっしゃられるように反応の並びや greater、 less あ...
実装が、本と異なるのは理解できましたが、contrast マトリッ...
> df <- data.frame(
+ Dose = ordered(rep(c("Cont", "Low", "MidLow", "MidHigh...
+ 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="...
Simultaneous tests: Williams contrasts
Call:
simtest.formula(formula = Response ~ Dose, data = df, ty...
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 の中を調べる必要があ...
**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" [#ic...
>[[kaz]] (2006-06-12 (月) 16:00:47)~
~
glmを用いておよそ100万件のデータをロジスティック回帰して...
~
インターネットで検索したところ、perfect separation、つま...
~
他にこのエラーの生じる原因をご存知であれば教えていただけ...
//
-また、5つの説明変数は全て連続量です。Rと回帰分析の初心者...
-Rのバージョンは2.2.1（日本語化版）です。 -- [[kaz]] &new...
-Quasi-complete separation とか Complete separation とい...
簡単に説明すれば，ある変数（または複数の変数の線形結合）...
データが何百万，何千万あろうとも，ある変数の分布が重なら...
共用パソコンに入っているなら仕方ないけど，最新版を使うの...
> 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 M...
-5.647e-05 -2.107e-08 0.000e+00 2.107e-08 5.560e-...
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...
2: 数値的に 0 か 1 である確率が生じました in: glm.fit(x ...
-大変わかり易い具体例を有難うございます。頭がすっきりしま...
-私のやり方が間違っているのかもしれませんが、上記のように...
-先ほど追加した部分と関連しますが，複数の変数の線形結合値...
-set.seed も含めて。上の通りやってみましたか？set.seed が...
> 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 M...
-1.256e-04 -2.107e-08 0.000e+00 2.107e-08 1.282e-...
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...
2: 数値的に 0 か 1 である確率が生じました in: glm.fit(x ...
-また、（これは単にRのバージョンの問題かもしれませんが）w...
-なるほど、確かに単回帰のエラーは重回帰のエラーの十分条件...
-set.seed で生成される乱数列は，バージョンによっても違う...
-個々の変数の分布は重なるが，線形結合の分布が重ならない例...
> 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...
(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...
(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...
(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 M...
-9.418e-05 -2.107e-08 0.000e+00 2.107e-08 6.798e-...
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...
2: 数値的に 0 か 1 である確率が生じました in: glm.fit(x ...
-詳しい例示ありがとうございます。「個々の変数の分布は重な...
-> 「個々の変数の分布は重なるが，線形結合の分布が重ならな...
上の例を見ればおわかり頂けると思ったのですが。分布が重な...
その方が理解しやすいということならかまいませんが，必要条...
上の x から合成される z を，y 別に描いてみる
重なっていないことがわかる
boxplot(z ~ y, horizontal=T)
#ref(boxplot.png)
- ずいぶん時間が経ってから見させていただきましたが、参考...
-Firth bias-correctionという方法があるようです -- nu &new...
--参考
---http://cran.r-project.org/web/packages/brglm/index.html
---http://www.ats.ucla.edu/stat/mult_pkg/faq/general/comp...
**Sweaveでxtableをminipage環境に適用したいです [#x624f617]
>[[Akira]] (2006-06-09 (金) 10:40:42)~
~
latexコマンドをRスクリプトに埋め込んで、Sweaveで結果をま...
xtableで作成したtableをminipage環境で2つ横に並べたいので...
> 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のオプションを確認しましたが、出力しない設定が分かり...
何か他に良い方法がありますでしょうか？
R-2.3.0、WinXPSP1を使用しています。
//
-print.xtableをみましょう. カラクリ屋敷(ぼそっ)... -- [[...
?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は...
-なかまさまの前に記入されていたsinkを利用する方法がなくな...
-バックアップをたどればありますよ. (良く見ろとか怒られる...
-勝手にしていいものかと・・・でも、Tipsに乗せてしまいまし...
**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.yyyymmdd <- function(sdate, sep="/") {
require(survival) # もちろんシステムにあらかじめインス...
temp <- date.mdy(sdate) # date.mdy は survival の関数
ifelse(is.na(sdate), "NA", paste(temp$year, temp$month,...
> date.yyyymmdd(100:105)
[1] "1960/4/10" "1960/4/11" "1960/4/12" "1960/4/13" "196...
ふむ、何か変かな。次のようになるのだけれど。
> date.yyyymmdd(13149)
[1] "1996/1/1"
-パッケージ data にも survival にも全く同じ関数があります...
-dates関数って，どのライブラリにあるのか。。。。~
date.yyyymmdd は，以下のように微修正すればよいでしょう。 ...
> date.yyyymmdd <- function(sdate) {
+ require(survival) # もちろんシステムにあらかじめイン...
+ temp <- date.mdy(sdate) # date.mdy は survival の関数
+ ifelse(is.na(sdate), "NA", sprintf("%4i/%02i/%02i", ...
+ }
> date.yyyymmdd(100:105)
[1] "1960/04/10" "1960/04/11" "1960/04/12" "1960/04/13" ...
> date.yyyymmdd(13149)
[1] "1996/01/01"
-13149 て Julian day (1960/01/01) からの経過日数ですよね...
-変も何も，そういう前提で経過日数後の日付を yyyy/mm/dd の...
> dates2 <- function(x)
+ {
+ s <- as.character(dates(as.numeric(13149),out="yyyy/m...
+ sprintf("%s%02i%s",substring(s, 1, 5),which(month.abb...
+ }
> dates2(13149)
[1] "2006/01/01"
-うん？
> dates(0,out="yyyy/mm/dd")
[1] 1970/Jan/01
なんだから，dates 関数では 1970/01/01 からの日付では？~
結局は，関数の定義によるようで，別の関数で代替処理すると...
-datesって S-Plusの関数では? -- &new{2006-06-09 (金) 07:...
-google で少し調べたら Julian date とは本来紀元前4713/1/...
-元の発言にも記載がないので，いったいどのdates関数だと，...
> 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関数についての記述をきちんとせず、皆さんを混...
今回は自分で関数を書いてこの問題を解決しようと思います。 ...
**分散の違い [#hfeda6c6]
>[[kanako]] (2006-06-07 (水) 14:21:24)~
~
異なる平均と標準偏差を持つ分布を3種類用意して、それらを順...
> 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...
> 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
~
//
-乱数は文字通りランダムですから、毎回同じ物が作られるとは...
> rnorm(3,0,1)
> rnorm(3,0,1)
> set.seed(1234)
> rnorm(3,0,1)
> set.seed(1234)
> rnorm(3,0,1)
-コメントありがとうございました。説明が不足していました。...
-ちなみに、二つの例で分散（と言うよりは共分散行列ですね）...
**データから円グラフを作成 [#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 を見てく...
-ありがとうございます。頑張ってみます。 -- [[やきそばぱん...
**プロビットモデルでの疑似決定係数の計算 [#o080db1a]
>[[イマイ]] (2006-06-05 (月) 23:18:44)~
~
プロビットモデルの回帰分析を行っていますが、疑似決定係数...
超初心者なりに考え付く限り検索をかけてみたのですが、みつ...
//
-前にも同じ質問をどこかでしてませんでしたっけ？というのは...
「プロビットモデル 疑似決定係数」でググってみると，20件く...
「12）プロビットモデルで疑似決定係数を計算する方法はいく...
-尤度比インデックスに関しては存じ上げませんが，疑似決定係...
-検索がまだまだ不十分だったようです。アドバイス、本当にあ...
**Mac OSXでのR 2.3.1でパッケージアップデートができない [#...
>[[高井]] (2006-06-05 (月) 17:21:12)~
~
PowerPC G5（デュアル 2.3GHz）, OSX 10.4.6でR 2.3.1のGUI版...
以下にエラーupdate.packages(lib = ""/Library/Frameworks/...
オブジェクト "Library" は存在しません
とあります。ところが，ターミナルから起動してupdate.packag...
//
-補足です：エラーメッセージが出た後に，Rパッケージインス...
-R コンソールから，update.packages() すると，バージョンの...
-はい，そうですね。Rコンソールの結果もターミナルからの結...
-Libraryの前のダブルクォートが２つになっているのが怪しい...
-翻訳はアップデートされてたけどアップデートは直ってないで...
-間違いなくバグでした。PackageInstaller.mの中でupdate.pac...
-私の環境のせいではなくて一安心です。 -- [[高井]] &new{20...
**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.si...
2: Reached total allocation of 734Mb: see help(memory.si...
> summary(slm2)
以下にエラーsummary(slm2) : オブジェクト "slm2" は存在し...
となって計算が実行できません。~
ちなみにメモリーは
> memory.size(T)
[1] 61349888
です。~
Rではハードディスクを用いた仮想メモリを利用できる機能等は...
//
-help(memory.size) はしてみましたか？ この wiki 内を「メ...
-help(memory.size)も実行済みです。memory.limit(size=1024)...
-NULLが返ってきた後に、memory.limit()をすると上限が増えて...
-(A+B+C+D+E)^2を独立変数にした回帰で，しかも変数選択をし...
-memory.limit()で確かに上限が増えています。が、計算が上手...
--モデル自体の吟味をすることを，まずお薦めします~
そうかもしれません。実は回帰と言うより、実際は数量化I類で...
ところで｢起動コマンドライン｣とは何でしょうか?単語検索でも...
-Rを起動するアイコンの上でマウスを右クリックしてプロパテ...
-それを｢起動コマンドライン｣と言うのですね。実は--sdiにし...
-起動アイコンのコマンドラインパラメータなので，そう呼んで...
-なるほど。Tipsありがとうございます。数量化I類は結局どの...
-カテゴリー変数をダミーデータに変換して重回帰分析を行う場...
-Rのfactor関数を使ってダミー化したものですから、元締めの...
#ref(QQplot.jpg)
あれ、画像添付こんなのでいいのかな･･････?~
ご覧のように、正規Q-Qプロットを一応試みたのですが、大部分...
そんなわけで、交互作用の存在を疑ってみたのです。-- [[R初...
-step関数は大変便利な方法だと思いますが、強引な面も持って...
まずは、主効果のみモデルのAICを求めて、~
その後、2次の交互作用項を一つ入れたモデルを作って~
AICが改善されるかを見ていくのもいいと思います。-- [[MK]] ...
-なるほど。力技ですね。試してみます。 -- [[R初心者]] &new...
-バグってRが止まってしまいました。。。。。。○|￣|＿~
ただ、交互作用A:Bを入れたモデルだけは計算できたんですが、...
ところで、AIC関数を用いるとstep関数で言うRSSが表示される...
-その交互作用項のt値はどうでしたか？有意であれば、交互作...
-交互作用項A:BのP値は全て 2×10^-16以下を示しています。
#ref(QQplot2.jpg)
交互作用A:Bを考慮した場合の正規Q-Qプロットです。あんまり...
-交互作用を疑うよりは，81510番のデータを良く精査して，問...
-81510番は確かにイレギュラーに見えるので、次のようにして...
df <- df[-c(81510),]
それで元の交互作用項を含まないモデル式に戻して、正規Q-Qプ...
#ref(QQplot4.jpg)
今度は乖離が派手に見えます。指数関数的に増加していってま...
-これは，交互作用と言うよりは，独立変数と従属変数が直線相...
-大きすぎると思います。wiki上で図を小さくする方法が分から...
正規Q-Qプロットより残差プロットかS-Lプロットの方が宜しい...
-最初から小さな図を描けばいいのです。ま，ちょっと描いてみ...
残差は以下のような分布になることが期待されています。同じ...
#ref(temp.png)
#ref(temp2.png)
-大量のデータを全部使って分析しなければ得られない結果とい...
-確かにRとあんまり関係無くなってますね。・・・・・・どっ...
#ref(回帰診断.jpg)
Rで画像を小さくする、と言うのはこんな感じで宜しいのでしょ...
間瀬先生の本に従って、見よう見まねで回帰診断プロットを行...
--あるモデルで表現できるシステムならば，データ数は少なく...
--近代統計学はサンプルを用いて母集団を推測するという点を...
仰る事は分かります。その点も気になって、逆に、｢データが多...
場合に拠ってはWinBugsを試してみようかな、とは思っています...
-0を中心にしてばらついているというのではなく，予測値が70...
-外れ値の影響が大きいんですね。という事は・・・・ロバスト...
-とか言ってましたがダメでしたね。
エラー：lqs failed: all the samples were singular
ロバスト/抵抗回帰と数量化I類は相性が悪い?連続量じゃないと...
**windows版Rのインストール [#sae95e3a]
>[[ありんこ]] (2006-06-01 (木) 16:51:01)~
~
すみません、中級者の所にしゃしゃり出てしまって…。インスト...
//
-最新版の[[R-2.3.0:http://cran.md.tsukuba.ac.jp/]]を使い...
-と書いた矢先にR-2.3.1が。上のリンク先（筑波大ミラーサイ...
-確かにインストールの説明のところに新しい2.3.1が追加され...
-「インストールの説明のところに新しい2.3.1が追加されてい...
-実はみんな違う人が書いたんだったりして:-) -- [１つは書...
**ノンパラメトリックなモデルでの変数選択 [#i7bc81e2]
>[[金井]] (2006-06-01 (木) 01:35:15)~
~
【やろうとしていること】~
20個程度の説明変数があり，その説明変数を利用して，目的変...
~
【現状】~
現在，R2.2.1を，Windows XPで利用中です．~
線形判別分析，ロジスティック回帰分析の変数選択に関してはA...
~
【お聞きしたい点】~
お聞きしたのは，２点です．~
~
１．ニューラルネットや分類木で，線形判別分析で提供されて...
~
２．そもそもニューラルネットや分類木で0/1の２値を予測する...
//
-質問２は R プロパーの話題ではありませんね。統計プロパー...
-回答しないのに、わざわざへこませるようなこと書かなくても...
-ほっとくよりいいんじゃない？回答者の意気をくじくようなこ...
聞く人は多いけど，答える人は少ないよ。答える人がいなくな...
2. については，変数選択というか予測に重要な変数を選択する...
既存の仕組みがないと言うのは障害かもしれないが，やり方は...
そもそも，重回帰の変数選択の仕組み自体を評価しない人もい...
総当たり法も，変数が20もあると無理だけど，有望そうな変数...
理論根拠に基づき，試行錯誤的に変数選択するのも良いのでは...
**rglがOSXでコンパイルできない [#g1afb318]
>[[高井]] (2006-05-30 (火) 12:30:55)~
~
R 2.3.0のフルセットをMac OSX 10.4.6, Xcode 2.3，PowerMac ...
~
api.cpp: In function ‘void rgl_user2window(int*, int*, d...
api.cpp:608: error: invalid conversion from ‘int*’ to ‘c...
api.cpp:608: error: initializing argument 6 of ‘GLint ...
const GLdouble*, const GLdouble*, ...
api.cpp: In function ‘void rgl_window2user(int*, int*, d...
api.cpp:633: error: invalid conversion from ‘int*’ to ‘c...
api.cpp:633: error: initializing argument 6 of ‘GLint ...
const GLdouble*, const GLdouble*, ...
GLdouble*, GLdouble*)’
make: *** [api.o] Error 1
chmod: /Users/takai/Library/R/library/rgl/libs/ppc/*: No...
ERROR: compilation failed for package 'rgl'
//
-608行，633行を修正すればよいだけでしょう（誰でもソースを...
-CRANからソースをダウンロードして，ダブルクリックで解凍し...
-[[R-develメーリングリストのアーカイブの中:http://tolstoy...
-ちゃんとやれば，誰にでもできるはずなんですけどね。~
レスポンスがないようなのですが，以下は，ターミナルで操作...
中澤さんのコメントを参考に，まず解凍された rgl ディレクト...
601a602
> GLint viewport[4]; 追加
605a607
> for (int i=0; i < 4; i++) viewport[i] = view[i];...
607c609
< gluProject(point[0],point[1],point[2],mo...
---
> gluProject(point[0],point[1],point[2],mo...
624a627
> GLint viewport[4]; 追加
628a632
> for (int i=0; i < 4; i++) viewport[i] = view[i];...
632c636
< gluUnProject(pixel[0],pixel[1],pixel[2],...
---
> gluUnProject(pixel[0],pixel[1],pixel[2],...
その後
$ R CMD INSTALL rgl
でしょうか？？ -- &new{2006-05-30 (火) 18:24:22};
-上記のお二方のコメントに従って再トライしました。api.cpp...
collect2: ld シグナル 6 [Abort trap] で終了させられました
/usr/bin/ld: warning internal error: output_flush(offset...
flushed block(offset = 464972, size = 3068)
calling abort()
make: *** [rgl.so] Error 1
chmod: /Library/Frameworks/R.framework/Versions/2.3/Reso...
No such file or directory
ERROR: compilation failed for package 'rgl'
となってコンパイル失敗です。一歩前進しましたがエラーメッ...
-R CMD INSTALL rgl を実行したときの pwd が，rgl ディレク...
-ああ，「上位のディレクトリでR CMD INSTALL rgl」とあった...
-/usr/bin/ld の internal errorが出るのは, RかDeveloper To...
-岡田さん，rglのバイナリーありがとうございます。無事にイ...
-3時間くらいかかって，本家からDevelopper Tools 2.3をダウ...
-環境依存の問題であることはほぼ明らかになりましたから、あ...
-岡田さんに励まされ（？），ネジ捲かれ（？）ハードウェアチ...
-どうやらgccがおかしいんじゃないかと当たりをつけました。/...
**文字を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")) -- [[t...
-x が実際にどうなっているのかわからないのだけど。他の人が...
x <- letters
y <- rnorm(length(x))
plot(as.numeric(factor(x)), y, type="b", xaxt="n")
axis(1, at=factor(x), labels=x)
老婆心ながら，このような場合には折れ線グラフは不適切と，...
-ああ、これ苦労した記憶があります。結局出来なかったんです...
-これだとlabels内文字列でX軸の順番がソートされてしまいま...
-y<-c(1,2,3,4);plot(1:length(y),y,xaxt="n");axis(1,1:leng...
**スクリプトでread.delimを使うとエラーになってしまいます...
>[[KT]] (2006-05-25 (木) 17:52:50)~
~
コンソールから~
data <- read.delim("data.txt", header=T)
と入力・実行すると正しく実行されるのに、スクリプトファイ...
source("test.R")
と実行すると、~
以下にエラーeval.with.vis(expr, envir, enclos) :
関数 "raed.delim" を見つけることができませんでした
とエラーになってしまいます。対策をご教授いただければ幸い...
//
-x raed.delim / o read.delim -- [[takahashi]] &new{2006-0...
-takahashi様、コメントありがとうございます。ですが、超初...
-あのですね。「つづりが間違えてますよ」と指摘されたわけで...
-scriptの方がtypoなんじゃないですか？ x raed / o read -- ...
-まさにスペルミスでした。何たるケアレスミス。どうもありが...
**英語版R1.3の入手法は？ [#a767b2b7]
>[[Yabo]] (2006-05-25 (木) 08:32:16)~
~
R1.3を本家からダウンロードしても、インストールすると日本...
//
-R1.3って(@.@)何？ -- [[okinawa]] &new{2006-05-25 (木) 09...
-えっと
set LANG=C
set LC_ALL=C
"*****?jgr.exe" #JGRのパス
というバッチファイルを作って，そこから起動すればO.K.です...
-R3.1.0 のことでしょう。未来のお話です(^_^;) というのは...
-32bitかつUTF-8のロカールを持つOS上なら日本語でも動きます...
-ｇさん、ご教示の方法でうまくいきました。ひとりぼっちでや...
-R1.3はR2.3の誤りでした。すんません。 -- [[yabo]] &new{20...
-gさんのバッチファイルの set LC_ALL=C のCをUTF-8にしても...
**条件が長さが2以上なので [#wb1cb369]
>[[レッド]] (2006-05-23 (火) 04:30:21)~
~
次ののメッセージがでて条件式が効きません．
TSi <- amp$Si
> if (TSi < 8) {
+ TAl <- 8 - TSi
+ } else {
+ TAl <- 0
+ }
Warning message:
条件が長さが2以上なので，最初の一つだけが使われます in: ...
超初心者で，Win．Version 2.2.1 ~
//
-if()を調べましょう。~
cond: A length-one logical vector that is not 'NA'. C...
length greater than one are accepted with a war...
the first element is used.
if(x)ではis.logical(x)がTRUE、かつis.na(x)がFALSEかつ、le...
ただ、この条件を満たしているのでしたら、amp$Siがどういっ...
エラーが出る場合は、いきなり応用問題はやめて基礎からはじ...
-TSi がベクトルなんでは？ c(1,3) < 2 は論理ベクトル c(TR...
-このページの上の方の注意事項を読んでから質問・投稿しまし...
amp$Si がベクトルなんでしょうね。だとしたら，TAl <- ifels...
ifelse 関数を調べましょう。 -- &new{2006-05-23 (火) 10:4...
-amp$Siはベクトルです． ex.amp$Si <- c(8.123, 7.897, 7....
**平均と標準偏差を使ったboxplotは描けませんか？ [#r2b5ff6e]
>[[くに]] (2006-05-22 (月) 18:35:13)~
~
箱ひげ図には平均値と標準偏差を用いたものもありますが、こ...
//
-平均値・標準偏差を求める関数があり，plot, lines, points,...
-ありがとうございます。できたら、具体例か参考になるウェブ...
-可能ですが，平均値と標準偏差を使うなら，ストリップチャー...
-可能かどうか聞かれたので「可能ですと」お答えしましたが。...
何だってできますから，どこまでやらなきゃならないかを決め...
元のデータを記号で書き加えたり，そのほかなんでもあなたの...
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...
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", xl...
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 (探索的データ解析) の精神」とはどのような...
-「平均と標準偏差を使ったboxplotを描くにはどうすればよい...
-「EDA (探索的データ解析) の精神に反する」、というのはデ...
**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",k...
ここで、エラーメッセージが次のようにでます。~
以下にエラーmatch.arg(kernel, c("rbfdot", "polydot", "ta...
'arg' は以下の一つでなければなりません： rbfdot, ...
besseldo...
stringdotが選択肢にないようです。
パッケージのマニュアルには、stringdotは選択肢にあり、ソー...
ものがあることは確認しました。それともデータセットの与え...
//
-よく知らないのだけど，ヘルプを見る限り，
## S4 method for signature 'list':
ksvm(x, y = NULL, type = NULL, kernel = "stringdot",
:
x a symbolic description of the model to be fit. When n...
can be a matrix or vector containg the training data
or a kernel matrix of class kernelMatrix of the train...
character vectors (for use with the string kernel).
なんだから，x はリストでないといけないのでは？ -- &new{2...
-下記のように試してみました。 -- [[北陸人?]] &new{2006-05...
> 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....
ここで、エラーメッセージが次のようにでます。
以下にエラーas.double.default(t(x)) : (list)オブジェク...
ラベルの与え方がまだ理解できていませんが、ラベルyとすると...
指定しないということでしょうか？
-なかなか，コメントがつきませんね。「わたしは，その方面は...
-そうですね。string kernelが正確に使えるSVMは少ないですね...
**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, kern...
#SVMでの訓練データの学習
ここで
以下にエラーmodel.frame.default(data = list(V1 = c(33,19...
とエラーメッセージがでます。~
確認してみますと
> 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 ...
というスペース区切りの数値データです。~
どういう操作が足らないのかご教示いただけませんでしょうか？~
//
-問題を再現できる必要最小限のデータと共に，再度質問される...
-そのデータは１行しかなくて、行列型に旨く変換されなかった...
-説明不足で申し訳ありません。 -- [[北陸人]] &new{2006-05-...
データセット(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="rbfd...
以下にエラーmodel.frame.default(data = list(V1 = c(1, 0)...
:オブジェクトが行列ではありません
> dim(test1m)
[1] 3 3
> class(test1m.train)
[1] "matrix"
> dim(test1m.train)
[1] 2 3
上記のようになりました。ksvmにかけたときのオブジェクト(te...
-第一引数のtypeがデータの列名に見つからないからじゃないで...
-library(kernlab)であることを書いて頂いた方が早く返事が来...
-どんぴしゃり解決しました。アホでした。どうもありがとうご...
**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....
共有ライブラリ '/Library/Frameworks/R.framework/Resourc...
を読み込めません
dlopen(/Library/Frameworks/R.framework/Resources/libra...
Library not loaded: /usr/X11R6/lib/libX11.6.dylib
Referenced from: /Library/Frameworks/R.framework/Resou...
Reason: image not found
追加情報: Warning message:
ライブラリ ‘/Users/doima/Library/R/library’ はパッケージ...
エラー：.onLoad は 'tcltk' のための 'loadNamespace' に失...
エラー：パッケージ 'tcltk' をロードできませんでした
//
-Macを使っていないので予想ですが、Tcl/Tkがinstallされてい...
-R Commanderには"X11"が必須ですが、これは標準ではインスト...
**Rで使えるデータ一覧 [#y1977666]
>[[究極超具あーる]] (2006-05-12 (金) 12:57:54)~
~
Rで例えば、~
library(stats)
data(morley)
とすると、データmorleyが使えるようになりますよね。このよ...
//
-データセットmorleyはstatsパッケージではなくbaseパッケー...
-morley などは，data(morley) などとしなくても，使えます。...
-インストール済みの全パッケージの全データ一覧を知りたいで...
data(package="パッケージ名")
とするとパッケージ中のデータ一覧が得られることがわかりま...
ただ、この中のItem列を抜き出そうとすると、以下のようなエ...
datas = data(package="accuracy")[["results"]];
datas[["Item"]]
以下にエラーdatas[["Item"]] : 添え字が許される範囲外です
datasのデータ構造が良くわかりません。Item列を抜き出すには...
-datas[,"Item"]とするか、data(package="パッケージ名")[["r...
local({pkg <- select.list(sort(.packages(all.available =...
+ if(nchar(pkg)) library(pkg, character.only=TRUE)})
の部分を書き換えると対応できるのではないかと。 -- &new{2...
-単純ミスでした。ありがとうございました。 -- [[究極超具あ...
-次のようにすると、究極超具あーるさんの望みは叶うと思いま...
x<-sort(.packages(all.available = TRUE))
lspkg <- function(pkgname){
data(package=pkgname)[["results"]][,"Item"]}
sapply(x,lspkg)
データを含んでないパッケージまで一々character(0)と返すの...
**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,...
plot(y)
フィルタhがうまく動けば、yは(t,y) = (0,1), (10,9)を通る直...
アドバイスをいただけると幸いです。よろしくお願いいたしま...
//
-それぞれの行ごとに，x, y がどのようになっているか，確認...
x は全部1なんですけどいいのかなぁ。 x = ts(1:1000, freq=1...
-典型的なまずい質問の例ですね。自分のしたいことを人にもわ...
-まずいですか？「フィルタをfilterでどう実装するか」という...
-あなたの説明はフィルターやら積分器とやらをすでにある程度...
-お勉強の結果（1）~
filter の引数 method が "convolution" のときは，「（重み...
二番目の引数は項数ごとの重み，つまり，n 項移動平均なら re...
ただし，本当は移動和なので，重みは等しくなくても良いし和...
使用例 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", sid...
> a2 <- filter(x, rep(1/3, 3), method="convolution", sid...
> a3 <- filter(x, c(2,0.5,1), method="convolution", side...
> 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 と同じ結果になる（当然だ...
重みの長さが1で値が2なら，a4[1] = x[1], a4[2] = x[2]+2*a4...
重みの長さが2以上で，値が違っても同様。a5[4] = x[4]+1*a5[...
> 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
ということですか。~
どうも，書かれたプログラムはあなたがやろうとしたこととは...
-お勉強の結果（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
以下略
プログラムに書かれたとおりの動きをしているようだ。 -- &n...
-お勉強の結果（4）~
あの〜。もとの発言にあった z^-1 というのは，「遅延素子」...
そうだとすると，関数の定義上それはあなたの意図とは違った...
-お勉強の結果（5）~
遅延素子かどうかはともかくとして，あなたの最初の質問で，...
もしそうだとすると，
> y <- filter(filter(x, 0.01, method="con"), 1, method="...
でいいのでは。つまり，外側の filter の第2引数が -1 ではな...
#ref(filter.png)
-いろいろありがとうございました。おっしゃるとおり、ｚ^-1...
**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....
> 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 = ",", co...
> write.table(y, file = "svmsample_y.csv", sep = ",", co...
> rm(list=ls(all=TRUE))
> x<-read.table("svmsample_x.csv", header = TRUE, sep = ...
> y<-read.table("svmsample_y.csv", header = TRUE, sep = ...
> iris1 <- ksvm(x,y,kernel="anovadot",kpar=list(sigma=0....
以下にエラーksvm(x, y, kernel = "anovadot", kpar = list(...
no direct or inherited method for function 'ksvm...
> table(y,predict(iris1,x))
エラー：オブジェクト "iris1" は存在しません
以下にエラーpredict(iris1, x) : unable to find the argum...
selecting a method for function 'predict'
~
読み書きでデータ型が変化したわけでもないようです。~
申し訳ありませんがご教示宜しくお願いします。~
//
-自己解決しました。失礼しました。 >x<-as.matrix(read.tab...
-read.table 関数はデータフレームを返す一方、ksvm は第一引...
**成長曲線の当てはめ [#ja28eabb]
>[[たつ]] (2006-04-29 (土) 15:09:09)~
~
y(売上高）とx(キャンペーン額)の関係を成長曲線に当てはめた...
以上の目的に適当なパッケージはありますでしょうか？~
その際に売上高の上限と下限を指定できたりすると好都合なの...
~
アドバイスよろしくお願いします。~
//
-glm関数で出来るかもと思い色々と試行錯誤もしましたが、自...
-まず RSiteSearch("growth curve") を実行し、ヒットしたた...
-ヒットしました･･･覚束ない英語力でチェックしていたところ...
お陰さまで解決できたのですが最後に質問をさせてください。
上限と下限を指定しない場合はglmでも可能なのでしょうか？--...
-glm は非正規誤差にたいする線形回帰手法ですから、成長モデ...
~
SSlogis
SSgompertz(stats) Gompertz Growth Model
SSweibull(stats) Weibull growth curve model
-素人の私にはSSlogisが使い勝手が良さそうです。これから詳...
**最小絶対偏差回帰 [#r1f44837]
>[[たろう]] (2006-04-26 (水) 20:19:04)~
~
R初心者ですが、Splusをしばし使用しておりました。~
Splusでは、最小絶対偏差回帰は~
l1fitですが、Ｒで対応する関数を見出すことできずにおります...
ご存知のかた、ご教示宜しくお願いします。 ~
//
- r-help ML の過去記事を検索したらつぎのような(二件の別個...
> is there any package or method which enables to comput...
> 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 ...
-より頑健な回帰手法（おそらく単なる l1fit よりベター？）...
-ご教示ありがとうございます。おかげさまで無事できました。...
-ちなみに二件の記事は Rsitesearch でそれぞれキーワード "l...
-なるほど。RSiteSearch、初めてしり早速、確認いたしました...
-「より頑健な回帰手法」の実行方法についても、お教え頂きあ...
**部分一致検索による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)の内容も良く分かりま...
すみませんが、ご教示宜しくお願いします。~
//
-以下に示すように，
in1[grep('AA', as.character(in1$Name)),]
でできますが、これだと、"BAA" があれば、それも摘出してし...
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 (水) 0...
-ちなみに「'subset' は論理値として評価されなければなりま...
-さらさん、青木先生、どうもありがとうございました。as.cha...
-行頭を表す正規表現 ^ 。 そんな便利な物があったんですね...
-Rjpwikiにも正規表現の詳しい説明がある -- &new{2006-04-2...
**常微分方程式系のバラメータ推定をするには？ [#le048e5a]
>[[使用1ヶ月目]] (2006-04-19 (水) 11:44:48)~
~
常微分方程式系のバラメータ推定を考えています。常微分方程...
//
-それぞれの関数の help に example がありますでしょ？その...
-いくら教えてる立場だとは言え、全く関係ない人が見ても不愉...
**データファイルを読み込んでcor.testをする時のエラー [#f6...
>[[tensho]] (2006-04-17 (月) 23:26:15)~
~
基本的な質問で申し訳ありません。~
データファイルを読みこんでcor.testを行うと、「長いオブジ...
//
-文字とおりの意味では？これ以上の回答を期待するなら、読み...
-データ入力は以下の通りです。
> 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")
-データフレームの意味を誤解していませんか。行がケース、列...
-ご親切に教えていただき、どうもありがとうございました。 -...
**GLMでの決定係数 [#n25e8784]
>[[ごう]] (2006-04-11 (火) 06:30:08)~
~
GLMの基本的な部分かもしれないのですが、情報をいただけませ...
正規分布を前提とした普通の重回帰モデルでは決定係数が算出...
近似した値は「誰々（人名）」の決定係数といった雰囲気で算...
正規分布重回帰モデルでもサンプルサイズが大きければ有意な...
重複投稿ご容赦願います。他のサイトでも質問したのですが...
//
-＞ 近似した値は「誰々（人名）」の決定係数といった雰囲気...
できるのかもしれないのですがとは，できるんですか，できな...
誰々の決定係数というのは，あるんですか，ないんですか？~
あるのなら，具体的に何なんですか？~
＞ これはロジスティック回帰を見る限りではどうも低く算出さ...
気がするだけなんですか，実際に低いんですか？~
その判断は，あなたがしたものですか，そのように書いてあっ...
＞ モデルが成り立つかどうかで判断するというような記載を読...
気がするだけなのですか？ -- &new{2006-04-11 (火) 08:55:2...
-不明瞭な質問でした。決定係数の例：Nagelkerke R^2, Cox&Sn...
**bartlett.test関数の中身を見るには？ [#j47456dd]
>[[使用1ヶ月目]] (2006-04-10 (月) 17:17:00)~
~
通常、関数の中身を見るには、関数名をそのまま入力すればよ...
//
-このスレッドの下の方に同じような質問がありますよ。 -- &...
**尤度のプロットについて [#l09ecac0]
>[[にゃあ]] (2006-04-08 (土) 14:26:34)~
~
例えば二項分布binomの確率分布はdbinomとして計算できますが...
また尤度のプロットというのも可能なんでしょうか?~
//
-尤度の定義通りに計算すればよろしいのでは。ただし、普通は...
N <- 20
X <- rbinom(100, N, 0.55) # データ
loglikelihood <- function(p) {
log.prob <- dbinom(0:N, N, prob=p, lo...
LL <- 0
for (x in X) LL <- LL + log.prob[x+1]...
# もしくは
# for (i in 0:20)
# # LL <- LL + sum(X[X=i])*log.prob...
# LL <- LL + sum(X==i)*log.prob[i+1...
return(LL)
}
-後で消去しますが，というか，もし指摘が正しいなら，記事修...
-ベルヌーイ尤度にベータ分布を掛ける作業を視覚化してみたか...
関数を作ってみたんですが、以下の様になりました。~
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(likel...
plot(theta, posterior, ylim=c(0, 1), type="l")
distribution <- posterior
}
}
かなり雑なんですが、見よう見まねでやってみました。~
しかし、実行速度にかなり問題があります。~
どなたかいいアイディアがありますでしょうか?-- [[にゃあ]] ...
-関数を書いて，それに対する意見を求めて，回答があったとし...
-なるほど。言われてみればその通りかもしれません。 -- [[に...
-毎回プロットしていれば時間はかかるし、途中の結果も見られ...
-そうですね。ただし、根本的な考え方の間違いに気付いたので...
-関連して気になったこと。非負整数値データの度数分布ベクト...
> 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 ...
> 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-...
-参照先を示してあげるか，再度例示してあげると，おんぶにだ...
> 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)) ですか。ありが...
-ですよね(^_^;)要望してみればいかがでしょう。 -- &new{20...
**ライブラリ [#i1aeaca4]
>[[超初心者]] (2006-04-03 (月) 15:18:03)~
~
Package source: evd_2.1-7.tar.gz ~
Windows binary: evd_2.1-7.zip ~
ライブラリには、上記のように、拡張子がgzのものと、zipのも...
//
-拡張子より，その前の情報を見ましょう。「Package source」...
-R(最新2.2.1)では、CRANサイトにあるパッケーッジは動作する...
何卒よろしくお願いします。 -- [[超初心者]] &new{2006-04-...
-パッケージの添付ファイルに，その情報が書かれているものが...
あなたも,
[[パッケージの説明:http://cran.md.tsukuba.ac.jp/src/contr...
-度々申し訳ありません。~
ライブラリの拡張子gz「Package source」はソースプログラム...
（zipの場合は、Rメニューから呼び出し後インストール機能が...
-ソースパッケージをビルドするにはRをビルドするのとほぼ同...
-Windows版のRをお使いなら、パッケージのインストールはメニ...
-早速ご指導深謝します。舟尾様のページを見ました。 -- [[超...
-なかま様早速ご指導深謝します。~
状況としては、ライブラリの拡張子gz「Package source」がＰ...
ＲはＰＣにインストールされています。~
この場合、ソースをＲからインストールする場合
>install.packages("c:/***.gz",CRAN=NULL)
では、パッケージ***がインストールできないのでしょうか？~
誠に申し訳ありません。 -- [[超初心者]] &new{2006-04-07 (...
-Rtools.zipと言うのと,Mingw(C,Fortranコンパイラ)を入れて...
-[[超初心者]]さんがインストールしたいのは[[evdパッケージ:...
-超初心者さんなら，バイナリをインストールすればいいでしょ...
-礼儀作法がわからずに誠に申し訳ありません。具体的に申し上...
-上のなかまさんのコメントの「[[FAQのリンク先:http://www.m...
-[[ビルドのログ:http://cran.us.r-project.org/bin/windows/...
#include <sys/types.h>
/* --> int64_t ; if people don't have the above, they ca...
/* #include "int64.h" */
とか言ってるので, <sys/types.h>では無く<inttypes.h>にすれ...
-#include "int64.h" では？で，sys/types.h はある場合の方...
-[[データサイズ:http://www.unix.org/version2/whatsnew/dat...
-カレントディレクトリにある *.h は " でくくるんじゃなかっ...
-それはそうです. ただ int64.hの中で何を定義しているかと言...
-ライブラリ(zip)でインストール時にエラーが発生します。
何卒よろしくご指導をお願いします。
一例をあげますと、copula_0[1].3-5~
http://cran.r-project.org/bin/windows/contrib/r-release/...
エラーメッセージは
> utils:::menuInstallLocal()
以下にエラーgzfile(file, "r") : コネクションを開くことが...
追加情報: Warning message:
圧縮されたファイル 'copula_0[1].3-5/DESCRIPTION' を開く...
その他、e1071_1[1].5-13、やkernlab_0[1].6-2もインストール...
これらzip形式ライブラリのインストール方法をご教示下さい。...
-まだ質問の仕方がおわかりになっていないようですね。関係な...
-返す言葉もありません。正直なところ該当箇所を見ましたが、...
-ご指導頂いたように、アーカイブの「プロキシーサーバー」を...
-options 関数と Sys.putenv 関数は同じ行に入力したのですか...
-2段に別々に入力しました。本件に関連しているかどうかわか...
-何の ID とパスワードなんでしょう。どのページを見ても要求...
**数値の置き換え [#g6c16a4c]
>[[青木繁伸]] (2006-03-31 (金) 23:03:18)~
~
整数ベクトル x を，別の数値に置き換える関数を作る必要があ...
条件はn行3列の行列で与えようと思います。xの要素が1列目以...
> 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]...
xが次のようなとき~
> x <- c(7, 4, 3, 1, 8, 6, 5, 5, 8, 1, 1, 10, 7, 6, 3, 5...
答えが~
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()の定義の長さを考えると（あちら...
-コメント，ありがとうございます。私の悪い癖で，探すより作...
-青木さんの例のように置き換える数値が単調増加正整数であれ...
> x <- c(7, 4, 3, 1, 8, 6, 5, 5, 8, 1, 1, 10, 7, 6, 3, 5...
# 1,2,3 は第一の区間 [1,3.5)、4,5,6,7 は第5の区間 [4,7.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 ...
> 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 # これが何を出力するかとかも
-どうもありがとうございます。助かりました。 -- [[内藤]] &...
**図を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...
par(new=T)
plot (x,z, ylim=c(0, 10000), ann=F, pch=4)
上記の図を描いた後，ファイル−別名で保存より，pdfを選択す...
エラー：invalid character sent to 'PostScriptCIDMetricIn...
追加情報: Warning messages:
1: invalid string in 'PostScriptStringWidth'
2: invalid string in 'PostScriptStringWidth'
3: invalid string in 'PostScriptStringWidth'
何が原因でしょうか？基本的な質問で大変申し訳ありませんが...
//
-?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...
par(new=T)
plot (x,z, ylim=c(0, 10000), ann=F, pch=4)
dev.off()
初心者ならどんな基本的なことを質問しても良いということで...
-Akira様．初心者STです．準備不足でいろいろとご迷惑をおか...
-今回初めて pdf ファイルとして書き出そうとしたのでしょう...
-R version 2.2.1, 2005-12-20, i386-pc-mingw32 512M WXPSP2...
-R version 2.2.1, 2005-12-20, i386-pc-mingw32 256MB WXPSP...
-多分インストール時に東アジアを選択しなかったのでは？(こ...
-青木様，okinawa様，Akira様，なかま様．初心者STです．再イ...
**大量のグラフを描くためには [#a6a1ddb0]
>[[antt]] (2006-03-14 (火) 00:08:58)~
~
はじめまして。大量のグラフを描きたいのですが，以下に記し...
【環境】R version 2.2.1, i386-pc-mingw32, メモリ 1GB~
【問題】大量のグラフを描き続けると，オブジェクトは別に増...
【知りたいこと】大量のグラフを描きたい場合に，どのような...
【プログラム例】
alternative <- 1:5 #選択肢の数
nitem <- 100 #項目の数
nschool <- 10 #施設の数
n <- 10000 #標本サイズ
##仮想データを作成
data <- data.frame(matrix(alternative, ncol=nitem, nrow=...
data$school <- rep(1:nschool, each=n/nschool) #1施設 100...
##画像作成用の関数を定義する
plotitem <- function(i, j) {
dir.create(paste("school", i, sep=""), showWarnings=...
filename <- paste("./school", i,"/item", j, ".ps", s...
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)
}
}
//
-がんばって整形してくれましたが、無意味に長すぎるのも読む...
-R version 2.2.1, 2005-12-20, i386-pc-mingw32 メモリ512M...
-闇版ならFMの処理のリークです.setHook(packageEvent("grDev...
-日本語をどうしても使いたい場合は, Postscriptを吐くならTe...
-ご回答ありがとうございました。仰るとおり闇版を使っており...
-フォントのメトリック情報の処理によってリークが発生すると...
-リーク => メモリーリーク => 使われたまま回収再利用不可...
-はい。 -- [[なかま]] &new{2006-03-17 (金) 00:08:31};
**imageの色を絶対指定するには [#i140c8f9]
>[[子牛]] (2006-03-13 (月) 19:39:19)~
~
はじめまして。imageで得た画像にcolをつかって色を指定した...
数値にたいして絶対的な方法で色を指定するにはどうすればよ...
~
> 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" "uti...
[7] "base"
//
-100段階のパレットで一段階の違いを表現するという例題は，...
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...
#ref(pict2.png)
-ありがとうございました、色の使う範囲をこれで指定できるの...
ただ、教えていただいた方法ですと、オフセットを手で調整せ...
image(x,y,a, zlim=c(-2,2), col=cm.colors(110), main=past...
この方法で絶対的な値を固定できるようです。（すみません、...
**SVMのマージン距離出力方法について [#o5f5dbe6]
>[[ちょめ夫]] (2006-03-13 (月) 12:10:44)~
~
現在SVMに関するパッケージの検討を行っております。~
検討対象としているパッケージは以下の２つです。~
１）e1071~
２）kernlab~
なお、判別したいクラスは2クラスの分類となっております。上...
http://www.msi.co.jp/vmstudio/materials/tech/classificati...
「図 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" "ut...
other attached packages:
colorspace kernlab e1071 class
"0.9" "0.6-2" "1.5-11" "7.2-16"~
//
-投稿法についてご指摘いただきまして、誠に申し訳ありません...
-当てずっぽうですが、もし話題の量が SVM アルゴリズムにと...
ページ名:
![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
