初心者のための R および RjpWiki に関する質問コーナー
新規投稿はできません
maechan (2006-06-26 (月) 14:04:44)
ある行列データがありそれを一つのヒストグラムで見ようと思い行列を一列にするプログラムを書きました。それを画面上に表示することはできるのですがファイルに取り込む方法がわかりません。
関数の値をファイルに書き込む方法が知りたいです。
初心者 (2006-06-24 (土) 12:03:31)
インストールの手順の最後にRconsole,Rdevga,Rprofile.siteをダウンロードし、同名ファイルに上書きとありますが、これらのファイルはどこからダウンロードしたらよいのでしょうか?
チョコボール (2006-06-23 (金) 07:21:00)
現在、perlでのスクリプトをバッチモードで使用しています。
変数を用いているのですが、c()の中では使用できません。
(こんなかんじ↓)Y <- X[c(args[5])]実行は
# R --vanilla --quiet --args test.pdf 3 5 3 < test.Rとやっています。
他の変数は得られています。
どなたかおしえてください。
ショーン (2006-06-22 (木) 19:50:20)
以下のように関数オブジェクトを変数に入れて、関数として呼び出すことができると思いますが、> a = cos > a .Primitive("cos") > a(0) [1] 1後でaの元の関数名である、"cos"という文字列を得る方法はありますでしょうか?つまり、
> f2s(a) [1] "cos"のf2sようなことはできますでしょうか?
> 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 <- tempfile())); unlist(strsplit(readLines(filename),"?""))[2]}
funcs = c(cos,sin,tan)のようにして、ループでいろいろな関数のグラフを描こうと思っています。複数のグラフファイルができますが、そのときのファイル名を関数名にするためにf2sを使います。 -- ショーン 2006-06-23 (金) 10:47:09
f2a <- function(x) { unlist(strsplit(capture.output(x),'"'))[2] }
ショーン (2006-06-22 (木) 12:37:36)
可変長引数...の長さを知りたいのですが、以下のようにうまくいきません。> func = function(...){print(length(...))} > func(1,2,3) 以下にエラーprint(length(...)) : 'length' に対する引数の個数が正しくありません「Rの関数定義の基本」を見るとlength(...)という書き方もあるように思えるのですが、何が間違っているのでしょうか?
func <- function(...) print(nargs())
> func = function(a, ...){print(length(c(...)))} > func(1,2,3) [1] 2
> 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(...)) でしょうか。
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向きの考え方でないのであれば改めて別言語で検討します.
よろしくお願いします
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)) # 度数分布を求めてその長さを見れば,カテゴリーが何種類あったかわかる })
tm (2006-06-20 (火) 22:23:00)
Rとは直接関係なく、大変申し訳ありませんが、Williamsの多重比較のパーセント点の表を、web経由で入手できるところはありませんでしょうか?
身近で入手できる本をあさったのですが、5%点くらいしか手に入りませんでした。
できれば、1%点の表が欲しいのですが、どなたかご存知ありませんか?
のぞき (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>の実際のコードを覗いてみることはできないのでしょうか。
R-2.3.x/src/nmath ディレクトリに dnorm.c pnorm.c qnorm.c rnorm.c という C ソースがある
ショーン (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.779837prとpr2の値が同じになると思っていたのですが、違います。どこかで計算ミスをしているのでしょうか?
> d <- structure(c(0.1991, 0.2003, 0.2169, 0.1429, 0.184, 0.2231, 0.1834, + 0.1516, 0.1478, 0.2071, 0.1613, 0.2288, -15.0277, -16.809, -16.6008, + -17.3694, -16.7424, -15.0885, -15.7661, -15.1187, -14.4152, -17.1567, + -13.5147, -15.0654, 0.279, 0.3462, 0.299, 0.2586, 0.2597, 0.3075, + 0.155, 0.294, 0.3411, 0.2789, 0.2592, 0.2165), .Dim = c(12, 3 + ), .Dimnames = list(c("1", "2", "3", "4", "5", "6", "7", "8", + "9", "10", "11", "12"), c("x1", "x2", "x3"))) > # eigen 関数から求める > res <- eigen(cor(d)) > print(t(sqrt(res$values)*t(res$vectors))) [,1] [,2] [,3] [1,] 0.7725894 -0.04950733 -0.6329728 [2,] -0.7142710 -0.37685279 -0.5897448 [3,] -0.2484622 0.92942168 -0.2728404 > # prcomp 関数から求める > res2 <- prcomp(d, scale.=TRUE) > print(t(res2$sdev*t(res2$rotation))) PC1 PC2 PC3 x1 0.7725894 0.04950733 -0.6329728 x2 -0.7142710 0.37685279 -0.5897448 x3 -0.2484622 -0.92942168 -0.2728404
t(res2$sdev*t(res2$rotation))は、ここのサイト(http://aoki2.si.gunma-u.ac.jp/lecture/PCA/pca1.html )でいう相関係数(因子負荷量)aですよね。それから、
res2$rotationが行列Lですよね。それから、
res2$xが主成分zですよね。投稿した質問のように、元の変数xから主成分zを求めるには、xにaではなくLをかけるのではないのでしょうか?ですので、res2$xとL%*%xが一致するはずだと思っているのですが、どこで勘違いをしているのでしょうか? -- ショーン 2006-06-19 (月) 17:39:26
> 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
akira (2006-06-16 (金) 14:11:29)
Sweave環境でxtableを書き換えたい場合、
(1)xtableの出力先をtexファイルからsink()で別のファイル(例えば"temp")へ切り替えて、
(2)readLinesでtempを読み込んで修正して、
(3)texファイルへ書き出す
ということを試してみました。しかし、sinkは受け付けてくれませんでした。
何かオプションがあるのでしょうか??documentclass{jarticle} % ?begin{document} ?section{sinkでtestが書き出せない} <<dat1, echo=F>>= sink("test") cat("ここはsinkでtestに出るはず?n") cat("testに出てますよ") cat("ここまで?n") sink() cat("test2?n") x <- readLines("test") cat("読み出したtestを書き出すよ?n") cat(x,"?n") cat("ここまで?n") @ ?section{table出力もできない} <<dat2, echo=F, results=tex>>= library(xtable) x <- data.frame(a=1:3, b=5:7) sink("test") cat("ここはsinkでtestに出るはず?n") print(xtable(x, caption="変更前のtable")) cat("????", "?n") cat("ここまで?n") sink() try(x <- readLines("test")) # ここがエラーでとまる cat("????", "?n") x[grep("caption", x)] <- sub("変更前", "変更後", x[grep("caption", x)]) x <- paste(x, "?n") cat("読み出したtestを書き出すよ?n") cat(x,"?n") cat("????", "?n") cat("ここまで?n") @ ?end{document}そこで、xtableの出力をxに入れなおして、readLinesの部分を書き換えると上手く動きます。
# xtableをxに入れる x <- xtable(x, caption="変更前のtable") # readLinesの部分はこうする x <- unlist(strsplit(x, split="?n")) # testがないからxtableの出力を変更する'RweaveLatex', 'Rtangle'当たりを見てもそれらしい記述が見当たりません。
- print.xtable2を発見。しかし、これを修正するのは大変そうだ。 -- akira 2006-06-19 (月) 16:36:45
tm (2006-06-15 (木) 21:10:37)
multcomp パッケージ中の simtest の Williams 検定についてですが、うまく動作していない気がします。
Rのバージョンは 2.0.1、multcomp のバージョンは 0.4-8、mvtnorm のバージョンは 0.7-2 です。> library(mvtnorm) > library(multcomp) > data(angina) > summary(simtest(response ~ dose, angina, type="Williams", alternative="greater")) Simultaneous tests: Williams contrasts Call: simtest.formula(formula = response ~ dose, data = angina,type = "Williams", alternative = "greater") Williams contrasts for factor dose Contrast matrix: dose0 dose1 dose2 dose3 dose4 C1 0 -1 0.00 0.0000000 0.0000000 1.0000000 C2 0 -1 0.00 0.0000000 0.5000000 0.5000000 C3 0 -1 0.00 0.3333333 0.3333333 0.3333333 C4 0 -1 0.25 0.2500000 0.2500000 0.2500000 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj C1 10.499 6.778 1.549 0 0 0 C2 7.747 5.775 1.341 0 0 0 C3 6.297 4.979 1.265 0 0 0 C4 5.247 4.284 1.225 0 0 0といった結果です。また、Interval の計算結果も
> summary(simint(response ~ dose, angina, type="Williams", alternative="greater")) Simultaneous 95% confidence intervals: Williams contrasts Call: simint.formula(formula = response ~ dose, data = angina, type = "Williams", alternative = "greater") Williams contrasts for factor dose Contrast matrix: dose0 dose1 dose2 dose3 dose4 C1 0 -1 0.00 0.0000000 0.0000000 1.0000000 C2 0 -1 0.00 0.0000000 0.5000000 0.5000000 C3 0 -1 0.00 0.3333333 0.3333333 0.3333333 C4 0 -1 0.25 0.2500000 0.2500000 0.2500000 Absolute Error Tolerance: 0.001 95 % quantile: 1.978 Coefficients: Estimate 5 % -- t value Std.Err. p raw p Bonf p adj C1 10.499 7.435 Inf 6.778 1.549 0 0 0 C2 7.747 5.093 Inf 5.775 1.341 0 0 0 C3 6.297 3.795 Inf 4.979 1.265 0 0 0 C4 5.247 2.824 Inf 4.284 1.225 0 0 0となってしまいます。p値はすべて 0 ですし、Interval に Inf が含まれていますし、結果がおかしい気がするのですが…
WEBで検索しても、multcomp の Williams 検定については情報が見つけられませんでした。
使用方法や結果の見方などご教授いただけないでしょうか?
よろしくお願いいたします。
> library(mvtnorm) > library(multcomp) > set.seed(12345678) > df <- data.frame(g = factor(rep(1:6, each=20)), y = rnorm(120)+rep(1:6/6, each=20)) > summary(simtest(y ~ g, df, type="Williams", alternative="greater")) Simultaneous tests: Williams contrasts Call: simtest.formula(formula = y ~ g, data = df, type = "Williams", alternative = "greater") Williams contrasts for factor g Contrast matrix: g1 g2 g3 g4 g5 g6 C1 0 -1 0.0 0.00 0.0000000 0.0000000 1.0000000 C2 0 -1 0.0 0.00 0.0000000 0.5000000 0.5000000 C3 0 -1 0.0 0.00 0.3333333 0.3333333 0.3333333 C4 0 -1 0.0 0.25 0.2500000 0.2500000 0.2500000 C5 0 -1 0.2 0.20 0.2000000 0.2000000 0.2000000 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj C1 0.771 2.751 0.280 0.003 0.017 0.008 C2 0.464 1.911 0.243 0.029 0.117 0.047 C3 0.299 1.309 0.229 0.097 0.290 0.124 C4 0.242 1.094 0.221 0.138 0.290 0.156 C5 0.092 0.426 0.217 0.335 0.335 0.335 Warning message: full precision was not achieved in 'pnt' > summary(simint(y ~ g, df, type="Williams", alternative="greater")) Simultaneous 95% confidence intervals: Williams contrasts Call: simint.formula(formula = y ~ g, data = df, type = "Williams", alternative = "greater") Williams contrasts for factor g Contrast matrix: g1 g2 g3 g4 g5 g6 C1 0 -1 0.0 0.00 0.0000000 0.0000000 1.0000000 C2 0 -1 0.0 0.00 0.0000000 0.5000000 0.5000000 C3 0 -1 0.0 0.00 0.3333333 0.3333333 0.3333333 C4 0 -1 0.0 0.25 0.2500000 0.2500000 0.2500000 C5 0 -1 0.2 0.20 0.2000000 0.2000000 0.2000000 Absolute Error Tolerance: 0.001 95 % quantile: 1.976 Coefficients: Estimate 5 % -- t value Std.Err. p raw p Bonf p adj C1 0.771 0.217 Inf 2.751 0.280 0.003 0.017 0.008 C2 0.464 -0.016 Inf 1.911 0.243 0.029 0.146 0.057 C3 0.299 -0.153 Inf 1.309 0.229 0.097 0.483 0.166 C4 0.242 -0.195 Inf 1.094 0.221 0.138 0.690 0.227 C5 0.092 -0.336 Inf 0.426 0.217 0.335 1.000 0.476
> library(mvtnorm) > library(multcomp) > set.seed(12345) > df <- data.frame( + Dose = ordered(rep(c("Cont", "Low", "MidLow", "MidHigh", "High"), each=7), c("Cont", "Low", "MidLow", "MidHigh", "High")), + Response = c( + 415, 380, 391, 413, 372, 359, 401, # Cont + 387, 378, 359, 391, 362, 351, 348, # Low + 357, 379, 401, 412, 392, 356, 366, # MidLow + 361, 351, 378, 332, 318, 344, 315, # MidHigh + 299, 308, 323, 351, 311, 285, 297 # High + )) > summary(simtest(formula=Response~Dose, data=df, type="Williams", alternative="less")) Simultaneous tests: Williams contrasts Call: simtest.formula(formula = Response ~ Dose, data = df, type = "Williams", alternative = "less") Williams contrasts for factor Dose Contrast matrix: DoseCont DoseLow DoseMidLow DoseMidHigh DoseHigh C1 0 -1 0.00 0.0000000 0.0000000 1.0000000 C2 0 -1 0.00 0.0000000 0.5000000 0.5000000 C3 0 -1 0.00 0.3333333 0.3333333 0.3333333 C4 0 -1 0.25 0.2500000 0.2500000 0.2500000 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj C1 -79.571 -7.076 11.245 0 0 0 C2 -63.500 -6.520 9.739 0 0 0 C3 -45.571 -4.963 9.182 0 0 0 C4 -39.714 -4.467 8.890 0 0 0となります。前述の本によると、
> df <- data.frame( + Dose = ordered(rep(c("Cont", "Low", "MidLow", "MidHigh", "High"), each=7), c("Cont", "High", "MidHigh", "MidLow", "Low")), + Response = c( + 415, 380, 391, 413, 372, 359, 401, # Cont + 387, 378, 359, 391, 362, 351, 348, # Low + 357, 379, 401, 412, 392, 356, 366, # MidLow + 361, 351, 378, 332, 318, 344, 315, # MidHigh + 299, 308, 323, 351, 311, 285, 297 # High + )) > summary(simtest(formula=Response~Dose, data=df, type="Williams", alternative="less")) Simultaneous tests: Williams contrasts Call: simtest.formula(formula = Response ~ Dose, data = df, type = "Williams", alternative = "less") Williams contrasts for factor Dose Contrast matrix: DoseCont DoseHigh DoseMidHigh DoseMidLow DoseLow C1 0 -1 0.00 0.0000000 0.0000000 1.0000000 C2 0 -1 0.00 0.0000000 0.5000000 0.5000000 C3 0 -1 0.00 0.3333333 0.3333333 0.3333333 C4 0 -1 0.25 0.2500000 0.2500000 0.2500000 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj C4 -39.714 -4.467 11.245 0.000 0.000 0.000 C3 -26.429 -2.878 9.739 0.004 0.011 0.007 C1 -22.143 -1.969 9.182 0.029 0.058 0.042 C2 -15.929 -1.636 8.890 0.056 0.058 0.056とするのが近いのかも。
Goeppingen (2006-06-14 (水) 10:42:34)
R内より、WebにアクセスしてHTMLやXMLを取得できるのでしょうか?
kaz (2006-06-12 (月) 16:00:47)
glmを用いておよそ100万件のデータをロジスティック回帰しています。説明変数は5個です。たまに首記のエラーが発生してしまいます。日本語版のエラーは「数値的に 0 か 1 である確率が生じました」です。
インターネットで検索したところ、perfect separation、つまりある説明変数XによってYが完全に説明されてしまう(Y=1 if and only if X<X0 for an X0)場合に当該エラーが出るようですが、100万件もありますし、またいくつかデータをチェックしましたが、そのようなことは起こっていないようです。(perfect separation参考:https://stat.ethz.ch/pipermail/r-help/2005-July/074693.html)
他にこのエラーの生じる原因をご存知であれば教えていただけますでしょうか。なお、100万件もデータを用いることの是非については別問題としてください(データが多すぎることがこのエラーを引き起こしているのならば別ですが、同じ数のデータで回帰してもエラーが起こらないこともあります)
> require(MASS) [1] TRUE > set.seed(12345) > y <- rep(0:1, each=25) > x <- matrix(rnorm(150), ncol=3) > cat(min(x[1:25,1]), max(x[1:25,1]),"?n") -1.817956 1.817312 ここと > cat(min(x[26:50,1]), max(x[26:50,1]),"?n") -2.380358 2.196834 ここの数値を見ると分布は重なっていない > summary(res <- glm(y ~ x, family=binomial)) Call: glm(formula = y ~ x, family = binomial) Deviance Residuals: Min 1Q Median 3Q Max -1.65228 -1.06530 -0.06326 1.03977 1.68695 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.1365 0.3115 -0.438 0.661 x1 0.3894 0.2823 1.379 0.168 x2 0.2135 0.2689 0.794 0.427 x3 -0.4449 0.2768 -1.607 0.108 (Dispersion parameter for binomial family taken to be 1) Null deviance: 69.315 on 49 degrees of freedom Residual deviance: 64.623 on 46 degrees of freedom AIC: 72.623 Number of Fisher Scoring iterations: 4 > x[1:25,1] <- x[1:25,1]+4.2 定数を加えて分布が重ならないようにしてやる 言い換えると,x[,1]の分布を描くと,二群が完全に分離しているのだ > cat(min(x[1:25,1]), max(x[1:25,1]),"?n") 2.382044 6.017312 ここと > cat(min(x[26:50,1]), max(x[26:50,1]),"?n") -2.380358 2.196834 ここの数値を見れば,分布が全く重なっていないことがわかる > summary(res <- glm(y ~ x, family=binomial)) Call: glm(formula = y ~ x, family = binomial) Deviance Residuals: Min 1Q Median 3Q Max -5.647e-05 -2.107e-08 0.000e+00 2.107e-08 5.560e-05 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 229.574 100367.677 0.002 0.998 x1 -99.764 42439.569 -0.002 0.998 x2 -6.802 19068.857 -0.000357 1.000 x3 19.210 15815.467 0.001 0.999 (Dispersion parameter for binomial family taken to be 1) Null deviance: 6.9315e+01 on 49 degrees of freedom Residual deviance: 9.7282e-09 on 46 degrees of freedom AIC: 8 Number of Fisher Scoring iterations: 25 Warning messages: 1: アルゴリズムは収束しませんでした in: glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart, 2: 数値的に 0 か 1 である確率が生じました in: glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,
> set.seed(12345) > y <- rep(0:1, each=25) > x <- matrix(rnorm(150), ncol=3) > x[1:25,1] <- x[1:25,1]+4.2 > summary(res <- glm(y ~ x[,1], family=binomial)) Call: glm(formula = y ~ x[, 1], family = binomial) Deviance Residuals: Min 1Q Median 3Q Max -1.256e-04 -2.107e-08 0.000e+00 2.107e-08 1.282e-04 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 460.8 118212.0 0.004 0.997 x[, 1] -201.3 51633.8 -0.004 0.997 (Dispersion parameter for binomial family taken to be 1) Null deviance: 6.9315e+01 on 49 degrees of freedom Residual deviance: 3.2194e-08 on 48 degrees of freedom AIC: 4 Number of Fisher Scoring iterations: 25 Warning messages: 1: アルゴリズムは収束しませんでした in: glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart, 2: 数値的に 0 か 1 である確率が生じました in: glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,
> set.seed(12345) > y <- rep(0:1, each=25) > x <- round(matrix(rnorm(150)*10+50, ncol=3), 1) > z <- apply(c(1,2,3)*t(x), 2, sum) # 線形結合を作って > oz <- order(z) > x <- x[oz,] # その大きい順に並べ替えた > z <- z[oz] # Z は確認のために表示するだけ > cbind(x, z) # 確認 z [1,] 26.2 41.4 37.1 220.3 [2,] 35.9 44.9 34.5 229.2 [3,] 68.2 36.6 36.8 251.8 [4,] 56.1 54.8 28.8 252.1 [5,] 40.8 26.5 53.2 253.4 中略 [47,] 55.2 65.9 63.2 376.6 [48,] 34.5 73.3 65.3 377.0 [49,] 61.2 55.2 71.8 387.0 [50,] 52.5 59.7 76.6 401.7 > # それぞれの変数を使って分析するとうまくいくが > summary(res <- glm(y ~ x[,1], family=binomial)) Call: glm(formula = y ~ x[, 1], family = binomial) Deviance Residuals: Min 1Q Median 3Q Max -1.66515 -1.06756 0.05372 1.06560 1.69950 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.50275 1.60011 -2.189 0.0286 * x[, 1] 0.06747 0.03016 2.237 0.0253 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 69.315 on 49 degrees of freedom Residual deviance: 63.530 on 48 degrees of freedom AIC: 67.53 Number of Fisher Scoring iterations: 4 > summary(res <- glm(y ~ x[,2], family=binomial)) Call: glm(formula = y ~ x[, 2], family = binomial) Deviance Residuals: Min 1Q Median 3Q Max -1.79526 -0.77141 0.04637 0.90898 1.81508 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -6.78401 2.06861 -3.280 0.001040 ** x[, 2] 0.12770 0.03841 3.324 0.000886 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 69.315 on 49 degrees of freedom Residual deviance: 52.668 on 48 degrees of freedom AIC: 56.668 Number of Fisher Scoring iterations: 4 > summary(res <- glm(y ~ x[,3], family=binomial)) Call: glm(formula = y ~ x[, 3], family = binomial) Deviance Residuals: Min 1Q Median 3Q Max -1.65849 -0.81491 -0.05309 0.67276 1.95259 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -7.58755 2.15363 -3.523 0.000426 *** x[, 3] 0.15294 0.04298 3.558 0.000373 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 69.315 on 49 degrees of freedom Residual deviance: 48.393 on 48 degrees of freedom AIC: 52.393 Number of Fisher Scoring iterations: 4 > # 全部を使ったらこける > summary(res <- glm(y ~ x, family=binomial)) Call: glm(formula = y ~ x, family = binomial) Deviance Residuals: Min 1Q Median 3Q Max -9.418e-05 -2.107e-08 0.000e+00 2.107e-08 6.798e-05 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2164.803 719087.033 -0.003 0.998 x1 5.672 2238.150 0.003 0.998 x2 14.675 5396.431 0.003 0.998 x3 21.969 7225.664 0.003 0.998 (Dispersion parameter for binomial family taken to be 1) Null deviance: 6.9315e+01 on 49 degrees of freedom Residual deviance: 1.7709e-08 on 46 degrees of freedom AIC: 8 Number of Fisher Scoring iterations: 25 Warning messages: 1: アルゴリズムは収束しませんでした in: glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart, 2: 数値的に 0 か 1 である確率が生じました in: glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,
上の x から合成される z を,y 別に描いてみる 重なっていないことがわかる boxplot(z ~ y, horizontal=T)
Akira (2006-06-09 (金) 10:40:42)
latexコマンドをRスクリプトに埋め込んで、Sweaveで結果をまとめています。
xtableで作成したtableをminipage環境で2つ横に並べたいのですが、Sweave()が出力するtexファイルの該当箇所を探して修正するのは面倒なので、minipage環境でxtableを出力したいと思っています。> library(xtable) > x <- data.frame(matrix(1:12, 2)) > xtable(x) % latex table generated in R 2.3.0 by xtable 1.3-2 package % Fri Jun 09 10:27:57 2006 ?begin{table}[ht] ?begin{center} ?begin{tabular}{rrrrrrr} ?hline & X1 & X2 & X3 & X4 & X5 & X6 ?? ?hline 1 & 1.00 & 3.00 & 5.00 & 7.00 & 9.00 & 11.00 ?? 2 & 2.00 & 4.00 & 6.00 & 8.00 & 10.00 & 12.00 ?? ?hline ?end{tabular} ?end{center} ?end{table}この出力にminipage環境のコマンドを埋め込みたいです。
% latex table generated in R 2.3.0 by xtable 1.3-2 package % Fri Jun 09 10:27:57 2006 ?begin{table}[ht] ?begin{minipage}{0.5?hsize} ?begin{center} ?begin{tabular}{rrrrrrr} ?hline & X1 & X2 & X3 & X4 & X5 & X6 ?? ?hline 1 & 1.00 & 3.00 & 5.00 & 7.00 & 9.00 & 11.00 ?? 2 & 2.00 & 4.00 & 6.00 & 8.00 & 10.00 & 12.00 ?? ?hline ?end{tabular} ?end{center} ?end{minipage} #2列の場合は?begin{minipage}〜?end{minipage}が繰り返される ?end{table}こんな感じです。
汎用性を持たせる技術がないので、今は> cat("??begin{table}[ht]?n?n") > y <- print(xtable(x)) #(1) > y <- sub("????begin??{table??}??[ht??]?n", "", y) > y <- sub("????end??{table??}?n", "", y) > print(cat(y)) > cat("??end{table}?n?n")としていますが、この場合(1)のところの出力がtexファイルに書き込まれてしまうので、結局修正しなればなりません。
printのオプションを確認しましたが、出力しない設定が分かりませんでした。RSiteSearchでもminipage環境の情報が少なく、困っています。
何か他に良い方法がありますでしょうか? R-2.3.0、WinXPSP1を使用しています。
?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}
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
> date.yyyymmdd <- function(sdate, sep="/") { require(survival) # もちろんシステムにあらかじめインストールしておきましょう temp <- date.mdy(sdate) # date.mdy は survival の関数 ifelse(is.na(sdate), "NA", paste(temp$year, temp$month, temp$day, sep=sep) )} > date.yyyymmdd(100:105) [1] "1960/4/10" "1960/4/11" "1960/4/12" "1960/4/13" "1960/4/14" "1960/4/15"ふむ、何か変かな。次のようになるのだけれど。
> date.yyyymmdd(13149) [1] "1996/1/1"
> date.yyyymmdd <- function(sdate) { + require(survival) # もちろんシステムにあらかじめインストールしておきましょう + temp <- date.mdy(sdate) # date.mdy は survival の関数 + ifelse(is.na(sdate), "NA", sprintf("%4i/%02i/%02i", temp$year, temp$month, temp$day )) + } > date.yyyymmdd(100:105) [1] "1960/04/10" "1960/04/11" "1960/04/12" "1960/04/13" "1960/04/14" "1960/04/15" > date.yyyymmdd(13149) [1] "1996/01/01"
> dates2 <- function(x) + { + s <- as.character(dates(as.numeric(13149),out="yyyy/mm/dd")) + sprintf("%s%02i%s",substring(s, 1, 5),which(month.abb==substr(s,6,8)),substr(s,9,11)) + } > dates2(13149) [1] "2006/01/01"
> dates(0,out="yyyy/mm/dd") [1] 1970/Jan/01なんだから,dates 関数では 1970/01/01 からの日付では?
> 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
kanako (2006-06-07 (水) 14:21:24)
異なる平均と標準偏差を持つ分布を3種類用意して、それらを順番に足したものを3列作成する際、別々に分けて作成したものと、一度に作成したものでは分散が違うようなのですが、何が違うのでしょうか?Rは2.3を使っていますが、2.1でも同じでした。いろいろ調べてみたのですが、何が悪いのか良く分かりませんでした。またどちらが正しいのでしょうか?> a <- rnorm(100,0,1) > b <- rnorm(100,3,3) > c <- rnorm(100,5,5) > abc <- cbind(a,a+b,a+b+c) > var(abc) a a 1.0695965 1.092034 0.7527298 1.0920344 9.323825 7.0911109 0.7527298 7.091111 30.6758695 > test <- cbind(rnorm(100,0,1), rnorm(100,0,1)+rnorm(100,3,3), rnorm(100,0,1)+rnorm(100,3,3)+rnorm(100,5,5)) > var(test) [,1] [,2] [,3] [1,] 0.9833900 0.10720269 0.56137488 [2,] 0.1072027 10.47080201 -0.03146996 [3,] 0.5613749 -0.03146996 34.98877571
> rnorm(3,0,1) > rnorm(3,0,1) > set.seed(1234) > rnorm(3,0,1) > set.seed(1234) > rnorm(3,0,1)
ヤキソバパンマン (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の割合が大きいですが…
この計算などを行って、円グラフで出力したいです。
どなたか教えてください。
イマイ (2006-06-05 (月) 23:18:44)
プロビットモデルの回帰分析を行っていますが、疑似決定係数、あるいは尤度比インデックスを計算させるにはどうすればよいのか、教えていただけませんか。
超初心者なりに考え付く限り検索をかけてみたのですが、みつけられませんでした。よろしくお願い致します。
高井 (2006-06-05 (月) 17:21:12)
PowerPC G5(デュアル 2.3GHz), OSX 10.4.6でR 2.3.1のGUI版の「パッケージとデータ」-->「パッケージインストーラ」を選んでバイナリーを「すべてアップデート」しようとすると,「パッケージのアップデートに失敗しました。詳細についてはコンソールを参照してください」となります。コンソールには以下にエラーupdate.packages(lib = ""/Library/Frameworks/R.framework/Resources/library/"", : オブジェクト "Library" は存在しませんとあります。ところが,ターミナルから起動してupdate.packages()するとエラーはでません。バグでしょうか?
R初心者 (2006-06-01 (木) 19:56:59)
Windows XPでR2.2.1を用いて回帰分析を行っています。
step関数でAICによる交互作用も含めた変数選択を行おうとしたのですが、> slm2 <- step(lm(y~(A+B+C+D+E)^2)) エラー:サイズ 1027089 Kb のベクトルを割り当てることができません 追加情報: Warning messages: 1: Reached total allocation of 734Mb: see help(memory.size) 2: Reached total allocation of 734Mb: see help(memory.size) > summary(slm2) 以下にエラーsummary(slm2) : オブジェクト "slm2" は存在しませんとなって計算が実行できません。
ちなみにメモリーは> memory.size(T) [1] 61349888です。
Rではハードディスクを用いた仮想メモリを利用できる機能等はないのでしょうか?
交互作用A:Bを考慮した場合の正規Q-Qプロットです。あんまり見た目が変わりません。と言うことは、依然とまだどっかに交互作用がある、と言う事でしょうね。 -- R初心者 2006-06-05 (月) 16:25:19
df <- df[-c(81510),]それで元の交互作用項を含まないモデル式に戻して、正規Q-Qプロットを試みてみました。 今度は乖離が派手に見えます。指数関数的に増加していってますね。-- R初心者 2006-06-05 (月) 18:30:46
エラー:lqs failed: all the samples were singularロバスト/抵抗回帰と数量化I類は相性が悪い?連続量じゃないとダメなんですね。一つ賢くなりました。 -- R初心者 2006-06-05 (月) 21:01:57
ありんこ (2006-06-01 (木) 16:51:01)
すみません、中級者の所にしゃしゃり出てしまって…。インストールの説明のところに、2.2.1、2,2,0、2,1,0、2,1,1の四つがあったのですが…。初心者にはどれが良いのでしょう。
金井 (2006-06-01 (木) 01:35:15)
【やろうとしていること】
20個程度の説明変数があり,その説明変数を利用して,目的変数(0/1)を予測しようとしています.
【現状】
現在,R2.2.1を,Windows XPで利用中です.
線形判別分析,ロジスティック回帰分析の変数選択に関してはAICに基づく変数選択を過去ログから発見できたのですが,分類木(treeパッケージ),ニューラルネット(nnetパッケージ)に関してはここの過去ログを含め,Webを調べても情報がありませんでした.
【お聞きしたい点】
お聞きしたのは,2点です.
1.ニューラルネットや分類木で,線形判別分析で提供されている関数 step() のような簡単な仕組みで変数選択を行う方法はありますか?
2.そもそもニューラルネットや分類木で0/1の2値を予測する際に,変数選択を行うことが間違いでしょうか.
高井 (2006-05-30 (火) 12:30:55)
R 2.3.0のフルセットをMac OSX 10.4.6, Xcode 2.3,PowerMac G5で使っていますが,rgl 0.66のコンパイルに失敗します。エラーメッセージは以下のようなものです。R 2.2.1のときはwarningはでるものの,rgl 0.66はコンパイルできて動きました。
api.cpp: In function ‘void rgl_user2window(int*, int*, double*, double*, double*, double*, int*)’: api.cpp:608: error: invalid conversion from ‘int*’ to ‘const GLint*’ api.cpp:608: error: initializing argument 6 of ‘GLint gluProject(GLdouble, GLdouble, GLdouble, const GLdouble*, const GLdouble*, const GLint*, GLdouble*, GLdouble*, GLdouble*)’ api.cpp: In function ‘void rgl_window2user(int*, int*, double*, double*, double*, double*, int*)’: api.cpp:633: error: invalid conversion from ‘int*’ to ‘const GLint*’ api.cpp:633: error: initializing argument 6 of ‘GLint gluUnProject(GLdouble, GLdouble, GLdouble, const GLdouble*, const GLdouble*, const GLint*, GLdouble*, GLdouble*, GLdouble*)’ make: *** [api.o] Error 1 chmod: /Users/takai/Library/R/library/rgl/libs/ppc/*: No such file or directory ERROR: compilation failed for package 'rgl'
601a602 > GLint viewport[4]; 追加 605a607 > for (int i=0; i < 4; i++) viewport[i] = view[i]; 追加 607c609 < gluProject(point[0],point[1],point[2],model,proj,view, 次行のように変更 --- > gluProject(point[0],point[1],point[2],model,proj,viewport, 624a627 > GLint viewport[4]; 追加 628a632 > for (int i=0; i < 4; i++) viewport[i] = view[i]; 追加 632c636 < gluUnProject(pixel[0],pixel[1],pixel[2],model,proj,view, 次行のように変更 --- > gluUnProject(pixel[0],pixel[1],pixel[2],model,proj,viewport,その後
$ R CMD INSTALL rglでしょうか?? -- 2006-05-30 (火) 18:24:22
collect2: ld シグナル 6 [Abort trap] で終了させられました /usr/bin/ld: warning internal error: output_flush(offset = 462632, size = 2352) overlaps with flushed block(offset = 464972, size = 3068) calling abort() make: *** [rgl.so] Error 1 chmod: /Library/Frameworks/R.framework/Versions/2.3/Resources/library/rgl/libs/ppc/*: No such file or directory ERROR: compilation failed for package 'rgl'となってコンパイル失敗です。一歩前進しましたがエラーメッセージの意味が分からず(リンクローダーでのエラーでしょうか?)立ち往生です。-- 高井 2006-05-31 (水) 11:32:56
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")軸の値は上のケース、プロットは下のケースのようにするにはどうしたらいいのでしょうか?
x <- letters y <- rnorm(length(x)) plot(as.numeric(factor(x)), y, type="b", xaxt="n") axis(1, at=factor(x), labels=x)老婆心ながら,このような場合には折れ線グラフは不適切と,統計学の教科書に書いてないか? -- 2006-05-26 (金) 18:18:26
KT (2006-05-25 (木) 17:52:50)
コンソールからdata <- read.delim("data.txt", header=T)と入力・実行すると正しく実行されるのに、スクリプトファイル("test.R")に同じ内容を記述して、
source("test.R")と実行すると、
以下にエラーeval.with.vis(expr, envir, enclos) : 関数 "raed.delim" を見つけることができませんでしたとエラーになってしまいます。対策をご教授いただければ幸いです。よろしくお願いします。
Yabo (2006-05-25 (木) 08:32:16)
R1.3を本家からダウンロードしても、インストールすると日本語表示になってしまいます。すでに"RとJava"のところでR1.2について報告があるのと同じく、JGRを動かすと文字化けします。残念ですが日本語が使えなくてもかまわないことにして、JGRを動かしたいです。英語版Rの導入法をどなたかご教示下さい。
set LANG=C set LC_ALL=C "*****?jgr.exe" #JGRのパスというバッチファイルを作って,そこから起動すればO.K.です. -- g 2006-05-25 (木) 10:29:28
レッド (2006-05-23 (火) 04:30:21)
次ののメッセージがでて条件式が効きません.TSi <- amp$Si > if (TSi < 8) { + TAl <- 8 - TSi + } else { + TAl <- 0 + } Warning message: 条件が長さが2以上なので,最初の一つだけが使われます in: if (TSi < 8) {超初心者で,Win.Version 2.2.1
cond: A length-one logical vector that is not 'NA'. Conditions of length greater than one are accepted with a warning, but only the first element is used.if(x)ではis.logical(x)がTRUE、かつis.na(x)がFALSEかつ、length(x)が1(ただし2以上は一つ目が使われるようです)である必要があります。レッドさんのTSiはこのどれかを満たしていないのです。
くに (2006-05-22 (月) 18:35:13)
箱ひげ図には平均値と標準偏差を用いたものもありますが、これをRで作図することは可能でしょうか?
bp <- function(x, px, wx=0.25) { m <- mean(x) s <- sd(x) mx <- max(x) mn <- min(x) x1 <- px-wx y1 <- m-s x2 <- px+wx y2 <- m+s rect(x1, y1, x2, y2) lines(c(px-wx, px+wx), c(m, m)) arrows(px, m+s, px, mx, angle=90, length=2*wx) arrows(px, m-s, px, mn, angle=90, length=2*wx) } x <- rnorm(100) plot(c(0,1), c(min(x), max(x)), type="n", xaxt="n", xlab="", bty="n", ylim=c(min(x), max(x))) bp(x, 0.5, 0.4) z <- matrix(rnorm(300), nr=100, nc=3) plot(c(0,3), c(min(z), max(z)), , type="n", xaxt="n", xlab="", bty="n", xlim=c(0,4), ylim=c(min(z), max(z)), ylab="foo bar baz") for (i in 1:3) { bp(z[,i], i, wx=0.12) }fig 1. example1 fig 2. example2
北陸人? (2006-05-18 (木) 22:42:37)
ksvmでのstring kernelを使いたいと考えています。
データセットを次のようにセットしました。(sktest.txt)AA 1 AB -1 BB 1コマンドは次のようにしました。
> library(kernlab) > sktest<-read.table("C:ファイルの場所") > sktest V1 V2 1 AA 1 2 AB -1 3 BB 1 > set.seed(50) > sktest.num<-sample(3,2) > sktest.train<-sktest[sktest.num,] > sktest.test<-sktest[-sktest.num,] > sktest.svm<-ksvm(V2~.,data=sktest,kernel="stringdot",kpar=list (lambda=0.5))ここで、エラーメッセージが次のようにでます。
以下にエラーmatch.arg(kernel, c("rbfdot", "polydot", "tanhdot", "vanilladot", : 'arg' は以下の一つでなければなりません: rbfdot, polydot, tanhdot, vanilladot, laplacedot, besseldot, anovadot, splinedot, matrixstringdotが選択肢にないようです。 パッケージのマニュアルには、stringdotは選択肢にあり、ソース(stringk.c)という ものがあることは確認しました。それともデータセットの与え方がまずいのでしょうか?ご教示いただければ幸いです。
## S4 method for signature 'list': ksvm(x, y = NULL, type = NULL, kernel = "stringdot", : x a symbolic description of the model to be fit. When not using a formula x can be a matrix or vector containg the training data or a kernel matrix of class kernelMatrix of the training data or a list of character vectors (for use with the string kernel).なんだから,x はリストでないといけないのでは? -- 2006-05-18 (木) 23:16:13
> x<-c("aaa","abb","bbb") > xl<-list(x) > class(xl) [1] "list" > xl [[1]] [1] "aaa" "abb" "bbb" > s.svm<-ksvm(xl,kernel="stringdot",kpar=list(lambda=0.5))ここで、エラーメッセージが次のようにでます。
以下にエラーas.double.default(t(x)) : (list)オブジェクトは 'double' に変換できませんラベルの与え方がまだ理解できていませんが、ラベルyとすると、y=NULLというのは 指定しないということでしょうか?
北陸人 (2006-05-17 (水) 13:40:12)
ksvmを使いたいと思っていますが、
入力するデータの変換についてお伺いします。testdata <- read.table("C:ファイルの場所") #データ読み込み testdatam <- data.matrix(testdata) #list(data.frame)から行列に変換 set.seed(50) #乱数の設定 testdatam.num <- sample(14965, 9000) #訓練データの数の設定 testdatam.train <- testdatam[testdatam.num,] #訓練データ testdatam.test <- testdatam[-testdatam.num,] #テストデータの設定 testdatam.svm <- ksvm(type~., data=testdatam.train, kernel="rbfdot", kpar=list(sigma=0.001)) #SVMでの訓練データの学習ここで
以下にエラーmodel.frame.default(data = list(V1 = c(33,19,13,24, 21,23,: オブジェクトが行列ではありませんとエラーメッセージがでます。
確認してみますと> class(testdatam.train) [1] "matrix" > is.list(testdatam.train) [1] FALSE > is.matrix(testdatam.train) [1] TRUEとでます。
データは7 11 4 1 2 3 0 2 (略) 9 4 3 8 5 3 1 2 0 3 1 -1というスペース区切りの数値データです。
どういう操作が足らないのかご教示いただけませんでしょうか?
データセット(test1.txt)データ数を少なくしました。 0 1 -1 2 3 -1 1 5 +1コマンドは次のようにしました。
> test1<-read.table("C:ファイルの場所",header=F) > test1 V1 V2 V3 1 0 1 -1 2 2 3 -1 3 1 5 1 > test1m<-data.matrix(test1) > class(test1m) [1] "matrix" > set.seed(50) > test1m.num<-sample(3,2) > test1m.train<-test1m[test1m.num,] > test1m.test<-test1m[-test1m.num,] > test1m.svm<-ksvm(type~.,data=test1m.train,kernel="rbfdot",kpar=list (sigma=0.001)) 以下にエラーmodel.frame.default(data = list(V1 = c(1, 0), V2 = c(5, 1), V3 = c(1, :オブジェクトが行列ではありません > dim(test1m) [1] 3 3 > class(test1m.train) [1] "matrix" > dim(test1m.train) [1] 2 3上記のようになりました。ksvmにかけたときのオブジェクト(test1m.train?)が行列ではないというメッセージでした。dim(test1m.train)で2*3の行列だというメッセージがでていますので、行列型になっていると考えました。初歩的なミスだと思いますがどの点が間違っているでしょうか?長文になりご迷惑をお掛けして申し訳ありません。
どいま (2006-05-16 (火) 14:42:13)
MacOS XでのRコマンダーのインストールが出来ず困っています。
Rコマンダーをロードしようとすると、以下のようなメッセージが表示されます。
Rのバージョンは2.3.0です。要求されたパッケージ tcltk をロード中です Loading Tcl/Tk interface ... 以下にエラーdyn.load(x, as.logical(local), as.logical(now)) : 共有ライブラリ '/Library/Frameworks/R.framework/Resources/library/tcltk/libs/i386/tcltk.so' を読み込めません dlopen(/Library/Frameworks/R.framework/Resources/library/tcltk/libs/i386/tcltk.so, 10): Library not loaded: /usr/X11R6/lib/libX11.6.dylib Referenced from: /Library/Frameworks/R.framework/Resources/library/tcltk/libs/i386/tcltk.so Reason: image not found 追加情報: Warning message: ライブラリ ‘/Users/doima/Library/R/library’ はパッケージを含んでいません in: library() エラー:.onLoad は 'tcltk' のための 'loadNamespace' に失敗しました エラー:パッケージ 'tcltk' をロードできませんでした
究極超具あーる (2006-05-12 (金) 12:57:54)
Rで例えば、library(stats) data(morley)とすると、データmorleyが使えるようになりますよね。このような既にRで利用可能なデータの一覧と、その内容を知りたいのですが、どこかにまとめられたページはありますか?もしくは調べ方がわかればいいのですが、教えてください。
data(package="パッケージ名")とするとパッケージ中のデータ一覧が得られることがわかりました。 ただ、この中のItem列を抜き出そうとすると、以下のようなエラーが出ます。
datas = data(package="accuracy")[["results"]]; datas[["Item"]] 以下にエラーdatas[["Item"]] : 添え字が許される範囲外ですdatasのデータ構造が良くわかりません。Item列を抜き出すにはどうしたらいいでしょうか? -- 究極超具あーる 2006-05-12 (金) 14:03:22
local({pkg <- select.list(sort(.packages(all.available = TRUE))) + if(nchar(pkg)) library(pkg, character.only=TRUE)})の部分を書き換えると対応できるのではないかと。 -- 2006-05-12 (金) 14:20:37
x<-sort(.packages(all.available = TRUE)) lspkg <- function(pkgname){ data(package=pkgname)[["results"]][,"Item"]} sapply(x,lspkg)データを含んでないパッケージまで一々character(0)と返すのはご愛嬌。 -- 2006-05-12 (金) 14:51:45
ワース (2006-05-09 (火) 19:59:58)
下のフィルタをfilterで実装したいと思っています。h = 0.01/(1-z^-1)以下のようにスクリプトを書きました。
x = ts(rep(1, 1000), freq=100, start=0) x = lag(x,k=-100) y = filter(filter(x, 0.01, method="con", sides = 1), -1, method="rec") plot(y)フィルタhがうまく動けば、yは(t,y) = (0,1), (10,9)を通る直線となるはずなのですが、上記の例ではそうなりません。何が間違っているのでしょうか?
アドバイスをいただけると幸いです。よろしくお願いいたします。
> x <- ts(c(1,3,4,2,1,4,6,2), freq=2, start=0) > a1 <- filter(x, rep(1/3, 3), method="convolution", sides=1) > a2 <- filter(x, rep(1/3, 3), method="convolution", sides=2) > a3 <- filter(x, c(2,0.5,1), method="convolution", sides=2) > a4 <- rep(NA, 8) > for (i in 2:7) a4[i] <- 1*x[i-1]+0.5*x[i]+2*x[i+1] > cbind(x, a1, a2, a3, a4) Time Series: Start = c(0, 1) End = c(3, 2) Frequency = 2 x a1 a2 a3 a4 0.0 1 NA NA NA NA 0.5 3 NA 2.666667 10.5 10.5 1.0 4 2.666667 3.000000 9.0 9.0 1.5 2 3.000000 2.333333 7.0 7.0 2.0 1 2.333333 2.333333 10.5 10.5 2.5 4 2.333333 3.666667 15.0 15.0 3.0 6 3.666667 4.000000 11.0 11.0 3.5 2 4.000000 NA NA NAということですか。 -- 2006-05-11 (木) 11:15:05
> 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ということですか。
> 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 以下略プログラムに書かれたとおりの動きをしているようだ。 -- 2006-05-11 (木) 12:01:24
> y <- filter(filter(x, 0.01, method="con"), 1, method="rec")でいいのでは。つまり,外側の filter の第2引数が -1 ではなくて 1 って,あっけない。 -- 2006-05-11 (木) 12:22:22
こねこ (2006-04-29 (土) 17:36:31)
ksvmに関連してcsvファイルの読み書を行おうとしているのですが、詰ってしまいました。
環境はWindowsXP + R(2.3.0)です。
アドバイスをいただけますと幸いです。
まずlibrary(kernlab)の後、下記コマンドを実行すると正常に処理が行われます。
> x<-as.matrix(iris[51:150,3:4]) > y<-as.matrix(iris[51:150,5]) > iris1 <- ksvm(x,y,kernel="anovadot",kpar=list(sigma=0.5)) > table(y,predict(iris1,x))
ところが下記のように一度データをcsvファイルに保存してから、読み込みなおして実行するとエラーが発生してしまいます。
> x<-as.matrix(iris[51:150,3:4]) > y<-as.matrix(iris[51:150,5]) > write.table(x, file = "svmsample_x.csv", sep = ",", col.names = NA) > write.table(y, file = "svmsample_y.csv", sep = ",", col.names = NA) > rm(list=ls(all=TRUE)) > x<-read.table("svmsample_x.csv", header = TRUE, sep = ",", row.names=1) > y<-read.table("svmsample_y.csv", header = TRUE, sep = ",", row.names=1) > iris1 <- ksvm(x,y,kernel="anovadot",kpar=list(sigma=0.5)) 以下にエラーksvm(x, y, kernel = "anovadot", kpar = list(sigma = 0.5)) : no direct or inherited method for function 'ksvm' for this call > table(y,predict(iris1,x)) エラー:オブジェクト "iris1" は存在しません 以下にエラーpredict(iris1, x) : unable to find the argument 'object' in selecting a method for function 'predict'
読み書きでデータ型が変化したわけでもないようです。
申し訳ありませんがご教示宜しくお願いします。
たつ (2006-04-29 (土) 15:09:09)
y(売上高)とx(キャンペーン額)の関係を成長曲線に当てはめたいと思っています。
以上の目的に適当なパッケージはありますでしょうか?
その際に売上高の上限と下限を指定できたりすると好都合なのですが・・・・
アドバイスよろしくお願いします。
SSlogis SSgompertz(stats) Gompertz Growth Model SSweibull(stats) Weibull growth curve model
たろう (2006-04-26 (水) 20:19:04)
R初心者ですが、Splusをしばし使用しておりました。
Splusでは、最小絶対偏差回帰は
l1fitですが、Rで対応する関数を見出すことできずにおります。
ご存知のかた、ご教示宜しくお願いします。
> is there any package or method which enables to compute a linear > regression with the leas absolute value fit-criterion? Try rq() from the quantreg package.
You can also look at package nlrq, if you intend to do non linear LAD regression.
hogehoge (2006-04-25 (火) 23:21:19)
R初心者です。
以下のようなtableからNameに関する部分一致検索で行を抽出したいのですが、どのように書けば良いのでしょうか?>in1 Name Score AAA 25 AAB 50 ABB 75 BBB 100完全一致の場合は以下のように書けるのを見てなるほどと思い
>in2 <- subset(in1, Name=="AAB") >in2 Name Score 2 AAB 50次のように書いたのですが
>in2 <- subset(in1, match("AA", in4$Name)) >in2 <- sunbset(in1, match("AA", Name)) etc... 'subset' は論理値として評価されなければなりません。と怒られてしまいます。またhelp(match)の内容も良く分かりませんでした。RはWin用2.3.0を使用しています。
すみませんが、ご教示宜しくお願いします。
in1[grep('AA', as.character(in1$Name)),]でできますが、これだと、"BAA" があれば、それも摘出してしまいます。もし"AA" で始まる物だけ、ということなら、
in1[grep('AA', substr(as.character(in1$Name),1,2)), ]でできます。 -- さら 2006-04-25 (火) 23:50:27
in1[grep('^AA', as.character(in1$Name)),]で良いと思います。 -- 青木繁伸 2006-04-26 (水) 00:19:45
使用1ヶ月目 (2006-04-19 (水) 11:44:48)
常微分方程式系のバラメータ推定を考えています。常微分方程式系のバラメータ推定をするには、おそらくodesolveの中の関数rk4とnlmを組み合わせる必要があると思います。まだ、Rの操作に慣れていないので、うまく記述することができません。どこかにコーディング例がありましたらご紹介ください。よろしくお願いします。
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")
ごう (2006-04-11 (火) 06:30:08)
GLMの基本的な部分かもしれないのですが、情報をいただけませんでしょうか。
正規分布を前提とした普通の重回帰モデルでは決定係数が算出できますが、GLMでのポアソン回帰などでは普通の決定係数は算出できませんよね。
近似した値は「誰々(人名)」の決定係数といった雰囲気で算出できるのかもしれないのですが、これはロジスティック回帰を見る限りではどうも低く算出されているような気がします。
正規分布重回帰モデルでもサンプルサイズが大きければ有意なモデルとして評価されますが、R^2=0.1ではあまり意味がありません。GLMではこのあたりをどのように解釈すれば良いのでしょうか。Rで step(model)とすればAICで評価されたモデルが手に入りますが。特にGLMMではもっとめんどうで、モデルが成り立つかどうかで判断するというような記載を読んだ気がします。GLM、GLMMではサンプルサイズ至上主義なのでしょうか?
重複投稿ご容赦願います。他のサイトでも質問したのですが返答頂けなかったもので、こちらで再度質問させて頂きました。
使用1ヶ月目 (2006-04-10 (月) 17:17:00)
通常、関数の中身を見るには、関数名をそのまま入力すればよいのですが、bartlett.test関数は、この方法で中身を見ることができません。S-plusには、bartlett検定の関数がないので、Rのbartlett.test関数を移植しようと思っています。bartlett.test関数の中身を見る方法をご存知の方がいらっしゃいましたら、よろしくお願いします。
にゃあ (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, log=TRUE) # 対数二項確率 LL <- 0 for (x in X) LL <- LL + log.prob[x+1] # x=n に対する確率は log.prob[n+1] であることに注意 # もしくは # for (i in 0:20) # # LL <- LL + sum(X[X=i])*log.prob[i+1] #ではなくて(冷汗) # LL <- LL + sum(X==i)*log.prob[i+1] # ですね(失礼) return(LL) }
mybayes1 <- function(n, prob){ theta <- seq(0, 1, length=100) distribution <- dbeta(theta, 1, 1) plot(theta, distribution, ylim=c(0, 1), type="l") count <- 0 for(i in 1:n){ .x <- runif(1) if (.x <= prob) cointoss <- 1 else cointoss <- 0 if (cointoss == 1) count <- count + 1 likelihood <- dbinom(count, n, theta) posterior <- likelihood*distribution/sum(likelihood*distribution) plot(theta, posterior, ylim=c(0, 1), type="l") distribution <- posterior } }かなり雑なんですが、見よう見まねでやってみました。
> x <- sample(0:10, 20, replace=TRUE) > x [1] 0 2 0 10 0 6 0 8 10 3 10 3 5 6 2 4 3 1 4 6 > table(x) # 欠損した値 7,9 は無視される x 0 1 2 3 4 5 6 8 10 4 1 2 3 2 1 3 1 3 > sapply(0:max(x), function(i) sum(x==i)) # 例えばこうすればよいがもっと直接的な方法があってもよさそうな? [1] 4 1 2 3 2 1 3 0 1 0 3
> 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かな? -- 2006-04-10 (月) 22:45:32
超初心者 (2006-04-03 (月) 15:18:03)
Package source: evd_2.1-7.tar.gz
Windows binary: evd_2.1-7.zip
ライブラリには、上記のように、拡張子がgzのものと、zipのものがあります。当方はWINDOWS環境でRを使用することを考えていますが、2種類の拡張子は使い分けが必要なのでしょうか?あるいは、zipだけで十分なのでしょうか???
>install.packages("c:/***.gz",CRAN=NULL)では、パッケージ***がインストールできないのでしょうか?
#include <sys/types.h> /* --> int64_t ; if people don't have the above, they can forget about it.. */ /* #include "int64.h" */とか言ってるので, <sys/types.h>では無く<inttypes.h>にすれば, あなたも私も幸せになれるとお伝え下さい. -- なかま 2006-04-07 (金) 16:40:02
http://cran.r-project.org/bin/windows/contrib/r-release/copula_0.3-6.zipエラーメッセージは
> utils:::menuInstallLocal() 以下にエラーgzfile(file, "r") : コネクションを開くことができません 追加情報: Warning message: 圧縮されたファイル 'copula_0[1].3-5/DESCRIPTION' を開くことができませんその他、e1071_1[1].5-13、やkernlab_0[1].6-2もインストールしたいのですが、同様の現象です。
青木繁伸 (2006-03-31 (金) 23:03:18)
整数ベクトル x を,別の数値に置き換える関数を作る必要がありました。
条件はn行3列の行列で与えようと思います。xの要素が1列目以上で2列目以下であるときに3列目の数値に変換しようと言うことになります。> y <- matrix(c(1,3,1, 4,7,5, 8,10,9), byrow=TRUE, nc=3) > y [,1] [,2] [,3] [1,] 1 3 1 [2,] 4 7 5 [3,] 8 10 9xの要素が,2のときはy[1,1] <= x <= y[1,2] ですので y[1,3] つまり1に置き換える,8のときは9になるなどです。
xが次のようなとき> x <- c(7, 4, 3, 1, 8, 6, 5, 5, 8, 1, 1, 10, 7, 6, 3, 5, 10, 8, 3, 3)答えが
5, 5, 1, 1, 9, 5, 5, 5, 9, 1, 1, 9, 5, 5, 1, 5, 9, 9, 1, 1になることを期待されます。
そこで作ったのが,recode <- function(x, y) { sapply(x, function(z) y[,3][y[,1] <= z & z <= y[,2]]) }ですが,もっとスマートな定義があるでしょうか...
> x <- c(7, 4, 3, 1, 8, 6, 5, 5, 8, 1, 1, 10, 7, 6, 3, 5, 10, 8, 3, 3) # 1,2,3 は第一の区間 [1,3.5)、4,5,6,7 は第5の区間 [4,7.5)、8,9,10 は第9の区間 [8,10.5) に入る > y <- c(1,3.5, 3.6, 3.7, 4,7.5, 7.6,7.7, 8,10.5) > findInterval(x,y) [1] 5 5 1 1 9 5 5 5 9 1 1 9 5 5 1 5 9 9 1 1
内藤 (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を抽出するには,
どのようにすればよいでしょうか?
> 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 # これが何を出力するかとかも
初心者ST (2006-03-20 (月) 16:50:24)
皆様.初めて質問させていただきます.初心者STです.宜しくお願いいたします.x <- c(0:100) y <- 10000 / (1+0.05*x) z <- 10000 / (1+0.10*x) plot (x,y, ylim=c(0, 10000), main="hypothetical function", sub="Fig 1",pch=2) par(new=T) plot (x,z, ylim=c(0, 10000), ann=F, pch=4)上記の図を描いた後,ファイル−別名で保存より,pdfを選択すると,次のようなエラーが生じて,pdfファイルとして正常に保存できません.
エラー:invalid character sent to 'PostScriptCIDMetricInfo' in a single-byte locale 追加情報: Warning messages: 1: invalid string in 'PostScriptStringWidth' 2: invalid string in 'PostScriptStringWidth' 3: invalid string in 'PostScriptStringWidth'何が原因でしょうか?基本的な質問で大変申し訳ありませんが,もしお分かりの方がいらっしゃれば,ご教示頂ければ幸いです.なお,使用環境は,R version 2.2.1, 2005-12-20, i386-pc-mingw32です.以上,宜しくお願い申し上げます.
pdf() x <- c(0:100) y <- 10000 / (1+0.05*x) z <- 10000 / (1+0.10*x) plot (x,y, ylim=c(0, 10000), main="hypothetical function", sub="Fig 1", pch=2) par(new=T) plot (x,z, ylim=c(0, 10000), ann=F, pch=4) dev.off()初心者ならどんな基本的なことを質問しても良いということではありません。まずは検索とヘルプで調べましょう。-- Akira 2006-03-20 (月) 17:28:36
antt (2006-03-14 (火) 00:08:58)
はじめまして。大量のグラフを描きたいのですが,以下に記したようにエラーが出て解決できません。どなたか解決方法をご存知の方,お教え頂けると幸いです。
【環境】R version 2.2.1, i386-pc-mingw32, メモリ 1GB
【問題】大量のグラフを描き続けると,オブジェクトは別に増えていないのに,メモリを消費したままで,その内,Rがフリーズする
【知りたいこと】大量のグラフを描きたい場合に,どのようなプログラムを書けば,メモリの消費を抑えられるか知りたい
【プログラム例】alternative <- 1:5 #選択肢の数 nitem <- 100 #項目の数 nschool <- 10 #施設の数 n <- 10000 #標本サイズ ##仮想データを作成 data <- data.frame(matrix(alternative, ncol=nitem, nrow=n)) data$school <- rep(1:nschool, each=n/nschool) #1施設 1000人 ##画像作成用の関数を定義する plotitem <- function(i, j) { dir.create(paste("school", i, sep=""), showWarnings=F) #施設ごとにディレクトリを作成 filename <- paste("./school", i,"/item", j, ".ps", sep="") #施設ごとに,各項目のファイルを出す postscript(filename) #保存する subdata <- which(data$school==i) #施設ごとのデータを抽出する hist(data[subdata,i], breaks=c(.5 + 0:5)) dev.off() } ##まずは,1施設,10項目を描画してみる for (i in 1) { for (j in 1:10){ plotitem(i,j) } } ##次に,1施設,100項目を描画してみる ##メモリサイズにより,この時点でRがフリーズする ##postscript(); dev.off() が終了した時点でも,メモリが消費しているようである for (i in 1) { for (j in 1:100) { plotitem(i,j) } }
子牛 (2006-03-13 (月) 19:39:19)
はじめまして。imageで得た画像にcolをつかって色を指定したいのですが、たとえば次の例にあるような二つのグラフの違いを明らかにしたいのです。この例ではグラフaはグラフa-1よりもいつも1だけ大きくなります。
数値にたいして絶対的な方法で色を指定するにはどうすればよろしいのでしょうか?
> a <- rnorm(100) > dim(a) <- c(10,10) > b<- a-1 > > x <- 1:10; y <- 1:10 > > op <- par(mfrow=c(1,2)) > image(x,y,a, col=cm.colors(100), main=paste("a")) > image(x,y,b, col=cm.colors(100), main=paste("a-1")) > > par(op)
以下は使用環境のメモです。> sessionInfo() R version 2.1.1, 2005-06-20, i386-pc-mingw32 attached base packages: [1] "methods" "stats" "graphics" "grDevices" "utils" "datasets" [7] "base"
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)これでいかが? -- 青木繁伸 2006-03-13 (月) 20:45:28
image(x,y,a, zlim=c(-2,2), col=cm.colors(110), main=paste("a"))この方法で絶対的な値を固定できるようです。(すみません、昨夜の段階ではHelpの説明が理解できてませんでした)。ちなみに上下の指定をはみだした値は、透明のイメージになるようです。はみ出さないようにしておかないと、落としてしまいそうですね。-- 子牛 2006-03-14 (火) 14:47:55
ちょめ夫 (2006-03-13 (月) 12:10:44)
現在SVMに関するパッケージの検討を行っております。
検討対象としているパッケージは以下の2つです。
1)e1071
2)kernlab
なお、判別したいクラスは2クラスの分類となっております。上記二つのパッケージのどちらでも良いのですが、SVMの予測結果において、与えた各サンプルと各クラスのマージンまでの距離を知りたいと考えておりますが、e1071,kernlabのpredict関数ともに、予測されたクラスやProbability、vote結果は出力可能ですがマージンまでの距離はどのように求めればよいか分かりません。イメージとしては、以下ウェブページの
http://www.msi.co.jp/vmstudio/materials/tech/classification.html#1-5
「図 Support Vector Machine 出力画面」のSVM.出力.1、2のような 数値情報です。もしご存じの方いらっしゃったらアドバイスいただければ幸いです。 どうぞよろしくお願いいたします。> sessionInfo() R version 2.1.1, 2005-06-20, i386-pc-mingw32 attached base packages: [1] "methods" "stats" "graphics" "grDevices" "utils" "datasets" [7] "base" other attached packages: colorspace kernlab e1071 class "0.9" "0.6-2" "1.5-11" "7.2-16"~