Ｑ＆Ａ (初級者コース)/13

if (iris$Species[i] == "setosa") {
  iris$Sp[i] <- 1
} else {
  iris$Sp[i] <- 0
}
library(geoR)
library(gstat)
# データの読み込み
data(ca20)
#   geoR
vcloud <- variog(ca20)
(vcloud)
plot(vcloud)
#   gstat
ca20gstat <- cbind(ca20$coords[, "east"], ca20$coords[, "north"], ca20$data)
colnames(ca20gstat) <- c("east", "north", "data")
ca20gstat <- as.data.frame(ca20gstat)
vfgstat  <- variogram(data ~ 1, loc = ~east+north, data = ca20gstat,
                      cutoff = vcloud$max.dist, width = vcloud$max.dist / 13)
(vfgstat)
plot(vfgstat)
　　$n     $u       $v
1   543   43.77376  55.50552
2  1648   131.3213  78.09254
3  2364   218.8688  97.05838
4  2246   306.4163  117.1554
5  2531   393.9638  129.0646
6  2152   481.5114  146.4889
7  1708   569.0589  151.9713
8  1309   656.6064  153.6482
9   694   744.1539  144.7442
10  343   831.7014  146.8878
11  144    919.249  131.8958
12   59   1006.796  139.2288
13   11   1094.344  139.2727
     np       dist     gamma
1   543   60.35719  55.50552
2  1648  132.58886  78.09254
3  2364  220.90992  97.05838
4  2246  305.65494 117.15539
5  2531  391.80321 129.06460
6  2152  482.02453 146.48885
7  1708  568.46999 151.97131
8  1309  652.62969 153.64820
9   694  739.59074 144.74424
10  343  826.14460 146.88776
11  144  913.61963 131.89583
12   59  995.10190 139.22881
13   12 1082.94416 165.16667
R version 2.12.1 (2010-12-16)
Platform: i386-pc-mingw32/i386 (32-bit)
gstat_1.0-8
geoR_1.6-35
library(tcltk)
draw.graph <- function(...) {
  upper <- as.double(tclvalue(Slider))
  hist(x, freq = FALSE, ylim = c(0, upper),
       col = "cyan", border = "white")
  lines(den, lwd = 3)
}
x <- rnorm(200)
den <- density(x)
tt <- tktoplevel()
Slider <- tclVar("0.5")
slider <- tkscale(tt, from = 0.1, to = 1,
            resolution = 0.01, showvalue = TRUE,
            variable = Slider, orient = "horizontal",
            command = draw.graph)
tkgrid(tklabel(tt, text = "upper"), sticky = "w")
tkgrid(slider)
tkfocus(tt)
エラー：  サイズ 176.4 Mb のベクトルを割り当てることができません 
> memory.limit()
[1] 2047
> memory.size()
[1] 817.03
> data(meuse)
> variogram(log(zinc) ~ x + y + elev, meuse)
month <- c("1月", "2月", "3月", "4月", "5月", "6月",
           "7月", "8月", "9月", "10月", "11月", "12月")
sex    <- c("F", "F", "M", "M", "M")
height <- c(158, 162, 177, 173, 166)
weight <- c(50, 50, 60, 60, 70)
( x    <- data.frame(SEX=sex, HEIGHT=height, WEIGHT=weight) )
levels(factor(x$SEX))
levels(factor(x$HEIGHT))
levels(factor(x$WEIGHT))
colname <- names(x)
for(i in 1:3) {
  levels(factor(x$colname[i]))
}
> abc <- function() (print("ABC"))
> abc()
[1] "ABC"
[1] "ABC"
> abc <- function() {print("ABC")}
> abc()
[1] "ABC"
> abc <- function() print("ABC")
> abc()
[1] "ABC"
> (a <- 999)
[1] 999
> print(a <- 999)
[1] 999
Ｑ＆Ａ (初級者コース)/13
http://www.okada.jp.org/RWiki/?%A3%D1%A1%F5%A3%C1%20%28%BD%E9%B5%E9%BC%D4%A5%B3%A1%BC%A5%B9%29%2F13
[ トップ | Tips紹介 | 中級Ｑ＆Ａ | 初級Ｑ＆Ａ | R掲示板 | 日本語化掲示板 | リンク集 ]
[ リロード ]   [ 新規 | 編集 | 凍結 | 差分 | ファイル添付 ]   [ 一覧 | 検索 | 単語検索 | 最終更新 | バックアップ | ヘルプ ] 
> gonyu <- function(x, digits = 0) {
>	     ceiling(x * 10 ^ digits) * 10 ^ (-digits) - 
+           round(ceiling(x * 10 ^ digits) *
+           10 ^ (-digits) - x, digits)
>	}
> gonyu(1.45, 1)
> gonyu(2.45,1)
> x <- 1.45
> digits <- 1
> ceiling(x * 10 ^ digits) * 10 ^ (-digits) - x
> round(ceiling(x * 10 ^ digits) * 10 ^ (-digits) - x, digits)
> round(ceiling(1.45 * 10 ^ 1) * 10 ^ (-1) - 1.45, 1)
> round(0.05, 1)
> ceiling(1.45 * 10 ^ 1) * 10 ^ (-1) - 1.45
> ceiling(1.45 * 10 ^ 1) * 10 ^ (-1) - 1.45 - 0.05
> ceiling(1.45 * 10 ^ 1) * 10 ^ (-1) - (1.45 + 0.05)
> func <- function(x, n)
+ {
+ 	a <- x * 10 ^ (n + 1)
+ 	b <- a - floor(a / 10 ^ n) * 10 ^ n
+ 	return(round(a + 5 * (b >= 5), -1) / 10 ^ (n + 1))
+ }
> func(1.45, 1)
[1] 1.5
> func(1.449999999999999, 1)
[1] 1.4
> func(2.5, 1)
[1] 2.5
> func(2.449999999999999, 1)
[1] 2.4
「編集」→「コンソール画面を消去」
または
「Ctrl」＋「L」
y <- function(x,a,b){ log(a*x + x^-(a/(2*b)) } 
matplot(outer(1:5, 1:3), type="b", pch=1:3, lty=1:3)
legend("topleft", LETTERS[1:3], pch=1:3, lty=1:3, col=1:3, seg.len=4)
matplot(outer(1:5, 1:3), type="b", pch=1:3, lty=1:3)
a <- legend("topleft", LETTERS[1:3], pch=1:3, lty=1:3,
            col=1:3, seg.len=4, plot=FALSE)
i <- par("cxy")
sapply(1:3, function(j) {
  with(a, {
    segments(rect$left+i[1]/2, text$y[j], rect$left+i[1]*1.9,
             text$y[j], col=j, lty=j)
    segments(rect$left+i[1]*3.1, text$y[j], rect$left+i[1]*4.5,
             text$y[j], col=j, lty=j)})
})
points(rep(a$rect$left+i[1]*2.5, 3), a$text$y, pch=1:3, col=1:3)
text(a$text$x, a$text$y, LETTERS[1:3])
with(a$rect, lines(c(left+w, left+w, left), c(top, top-h, top-h)))
cat("\a") #音が鳴る
cat("\b") #不明
cat("\f") #フィード文字（半角スペース？）
cat("\n") #改行
cat("\r") #不明
cat("\t") #タブ
cat("\v") #不明（半角スペース？）
> test <- read.csv("sample.csv", header=TRUE, row.names="x")
> testa <- read.csv("sample.csv", header=TRUE, row.names="x")
> testb.rda <- rda(test~.testa, nangteng, scale=TRUE)
> testb.rda
Call: 
> plot(testb.rda)
rda(x ~ z$y)
rda(x ~ y, data = z)
dune.Manure <- rda(dune ~ Manure, data = dune.env)
dune.Manure
mean(<data.frame>) is deprecated
> mean(iris[,1:4])
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
    5.843333     3.057333     3.758000     1.199333 
 警告メッセージ： 
mean(<data.frame>) is deprecated.
 Use colMeans() or sapply(*, mean) instead. # これからは，こうしてね！という示唆
> colMeans(iris[,1:4]) # つまり，こんな風にとか
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
    5.843333     3.057333     3.758000     1.199333 
> sapply(iris[,1:4], mean) # こんな風にとか
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
    5.843333     3.057333     3.758000     1.199333 
> mean.data.frame
function (x, ...) 
{
    msg <- "mean(<data.frame>) is deprecated.\n Use colMeans() or
            sapply(*, mean) instead."
    warning(paste(msg, collapse = ""), call. = FALSE, domain = NA)
    sapply(X = x, FUN = mean, ...)
}
> mean.data.frame
function (x, ...) 
{
    sapply(X = x, FUN = mean, ...)
}
bi <- glm(y ~ x, family = binomial(link = "cloglog"))
lines(x, 1-exp(-exp(b + a*x)))
bioff <- glm(x ~ y + offset(log(z)),
             family = binomial(link ="cloglog"))
po <- glm(y ~ x, family = poisson)
lines(x, exp(b+ax))
pooff <- glm(y ~ x + offset(log(z)), family = poisson)
> data <- read.csv("read_sample.csv", header=TRUE, row.names=1)
data <- read.csv("sample.csv")
data[, -1] <- function(data[, -1])
write.csv(data)
R version 2.14.0 (2011-10-31)
Platform: i386-pc-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=Japanese_Japan.932  LC_CTYPE=Japanese_Japan.932
[3] LC_MONETARY=Japanese_Japan.932 LC_NUMERIC=C
[5] LC_TIME=Japanese_Japan.932

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base 
> x <- c(1, 1, 1, 0, 0, 0, 0, 1, 1, 1,
         1, 1, 0, 1, 1, 1, 1, 0, 1, 0,
         1, 0, 0, 1, 0, 1, 1, 0, 1, 0,
         0, 0, 0, 1, 1, 1, 0, 1, 1, 0,
         0, 1, 1, 0, 0, 1, 1, 0, 0, 1,
         0, 1, 0, 1, 1, 0, 1, 0, 0, 1,
         1, 0, 0, 0, 0, 1, 0, 0, 1, 0,
         1, 1, 0, 0, 1, 0, 0, 0, 1, 1,
         1, 0, 1, 0, 0, 1, 1, 0, 0, 1,
         1, 0, 0, 0, 0, 1, 1, 1, 0, 0)
> y <- matrix(x, nrow=10, ncol=10, byrow=TRUE)
> apply(y, 2, cumsum)
> apply(y, 1, cumsum)
> (y <- matrix(1:12, 3, 4))
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12
> apply(y, 1, cumsum)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    5    7    9
[3,]   12   15   18
[4,]   22   26   30
> t(apply(y, 1, cumsum))
     [,1] [,2] [,3] [,4]
[1,]    1    5   12   22
[2,]    2    7   15   26
[3,]    3    9   18   30
> conn <- file("sample.bin", "rb")
> readBin(conn, integer(), n=100*100, size=4)
> readBin(conn, integer(), n=100*100, size=2)
options(CRAN="http://cran.md.tsukuba.ac.jp/bin/windows/contrib/2.13")
Sys.setenv("http_proxy"="http://(アドレス):(ポート番号)")
A <- c("aaa","aaaa","aaaaa")
B <- c("bbb","aaa","bbbbb")
C <- c("aaa","ccc","cccc")
system.time({
	### テストデータを作る
	### 要素数 n=500 個のリスト
	### 各要素は 503 個の要素を持つ文字列ベクトル
	###   テストのため最後の 3 つ, foo, bar, baz は全ベクトルに含まれる
	### 各要素は 1 〜 5 文字からなる
	set.seed(1234567)
	n <- 500
	k <- 500
	d <- vector("list", n)
	for (i in 1:n) {
		str <- character(k)
		for (j in seq_along(str)) {
			str[j] <- paste(sample(letters,
                         sample(5, 1), replace=TRUE), sep="", collapse="")
		}
		d[[i]] <- c(str, "foo", "bar", "baz")
	}
})
system.time({
	### 重複する文字列の抽出
	ans <- intersect(d[[1]], d[[2]])
	for (i in 3:n) {
		ans <- intersect(ans, d[[i]])
	}
	print(ans)
})
テストデータの作成にかかった時間
    ユーザ   システム       経過  
     9.297      0.257      9.508

共通要素
[1] "foo" "bar" "baz"

共通要素を抽出するのにかかった時間
    ユーザ   システム       経過  
     0.049      0.001      0.055
ほげ（ほげ）
1
2
3
> read.csv("input/testCSV.csv", header=TRUE)
  ほげ.ほげ.
1          1
2          2
3          3
x <- data.frame( matrix(rnorm(500), ncol=10) )
names(x) <- c( sample(LETTERS[1:26], 9, FALSE), 'y' )
f <- paste( 'y ~', paste(names(x)[1:9], collapse='+') )
lm( as.formula(f), data=x )
> x <- data.frame( matrix(rnorm(500),ncol=10) )
> lmtmp <- lm(X10 ~ ., data=x)
> step(lmtmp)
Start:  AIC=-3.58
X10 ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 

       Df Sum of Sq    RSS     AIC
- X7    1    0.0288 31.227 -5.5375
- X1    1    0.0427 31.241 -5.5154
- X2    1    0.0650 31.263 -5.4796
- X9    1    0.0659 31.264 -5.4783
- X4    1    0.0936 31.291 -5.4340
- X5    1    0.1595 31.357 -5.3287
- X3    1    0.4099 31.608 -4.9310
<none>              31.198 -3.5837
- X6    1    1.3648 32.563 -3.4429
- X8    1    8.1977 39.395  6.0814
(途中省略)

Step:  AIC=-16.17    ### AIC基準で最終的に選ばれたモデル
X10 ~ X8 

       Df Sum of Sq    RSS      AIC
<none>              33.405 -16.1657
- X8    1    7.3309 40.736  -8.2455

Call:
lm(formula = X10 ~ X8, data = x)

Coefficients:
(Intercept)           X8  
   -0.04151      0.39186
> x <- data.frame(matrix(rnorm(10000),ncol=100))
> Formula <- unique(sapply(seq(10^4),FUN=function(i) paste("X100 ~ ",
             paste(sort(sample(names(x)[-100],9)),collapse="+"))))
> res <- lapply(seq(10^4), FUN=function(i) lm(Formula[i],data=x))
> Formula[9999]; res[[9999]]
[1] "X100 ~  X20+X25+X30+X41+X55+X6+X68+X80+X99"

Call:
lm(formula = Formula[i], data = x) 

Coefficients:
(Intercept)          X20          X25          X30          X41          X55  
   0.054830     0.020898     0.036853     0.025012    -0.007347    -0.138181  
         X6          X68          X80          X99  
   0.013766     0.071437    -0.010897     0.083487  
> library(scatterplot3d)
以下にエラー library(scatterplot3d) : 
  '‘scatterplot3d’' という名前のパッケージはありません
　　　：  関数 "scatterplot3d" を見つけることができませんでした
> f <- function(p) 3*p[1]^2-4*p[1]*p[2]+3*p[1]+6*p[2]^2+8*p[2]
> nlm(f, c(1,1))
$minimum
[1] -6.107143

$estimate
[1] -1.214287 -1.071429

$gradient
[1] 6.422047e-07 4.269174e-07

$code
[1] 1

$iterations
[1] 3
x <- c("A","B","C","D","E")
y <- c("D","E","B")

x[x %in% y]
> x <- c("B","A","B","D","C","D","E")
> y <- c("D","E","B")
> unlist(sapply(y, function(ch) x[x==ch]))
 D1  D2   E  B1  B2 
"D" "D" "E" "B" "B"
> (x <- data.frame(str.x=LETTERS[1:5], num.x=1:5))
  str.x num.x
1     A     1
2     B     2
3     C     3
4     D     4
5     E     5
> (y <- data.frame(str.y=c("D", "E", "B"), num.y=c(1,3,6)))
  str.y num.y
1     D     1
2     E     3
3     B     6
> x[sapply(y$str.y, grep, x$str.x),]
  str.x num.x
4     D     4
5     E     5
2     B     2
png("hoge.png", width=640, height=640, unit="px")
plot(1,1)
text(1,1,"ほげ")
dev.off()
png("test.png")
par(family="HiraMaruProN-W4")
hist(rnorm(10000), main="日本語を使う")
dev.off()
> object.size(1L)
32 bytes
> str(1L)
 int 1
> object.size(as.single(1L))
120 bytes
> str(as.single(1L))
 atomic [1:1] 1
 - attr(*, "Csingle")= logi TRUE
> object.size(as.double(1L))
32 bytes
> str(as.double(1L))
 num 1
> library(parallel)
> cl <- makeCluster(rep("localhost", 4), type="SOCK")
> y <- 3
> tmp <- function(x) x+y
> clusterApply(cl, 1:4, tmp)
 以下にエラー checkForRemoteErrors(val) : 
  4 nodes produced errors; first error:  オブジェクト 'y' がありません
> stopCluster(cl)
> library(parallel)
> cl <- makeCluster(rep("localhost", 4), type="SOCK") # ４つのコアで実行
> y <- 3
> tmp <- function(x) x+y
# すべての実行時環境オブジェクトを子プロセスにエクスポート。
# この例では clusterExport(cl,"y") で十分
> clusterExport(cl, ls())
> clusterApply(cl, 1:4, tmp)
[[1]]
[1] 4

[[2]]
[1] 5 

[[3]]
[1] 6

[[4]]
[1] 7

> stopCluster(cl)

a = 1
cat(a,"\n")
> R --vanilla --quiet < test.r
=====コンソール表示======
> a = 1
> cat(a)
1
=========================
> R --vanilla --slave < test.r
> library(ggplot2)
要求されたパッケージ reshape をロード中です 
要求されたパッケージ plyr をロード中です 
Error in dyn.load(file, DLLpath = DLLpath, ...) : 
  共有ライブラリ '/Library/Frameworks/R.framework/Versions/2.14/
    Resources/library/plyr/>libs/x86_64/plyr.so' を読み込めません： 
 dlopen(/Library/Frameworks/R.framework/Versions/2.14/Resources/
    library/plyr/libs/x86_64/plyr.so, 6): Library not loaded:
    /Library/Frameworks/R.framework/Versions/2.13/Resources/lib/libR.dylib
 Referenced from: /Library/Frameworks/R.framework/Versions/
    2.14/Resources/library/plyr/libs/x86_64/plyr.so
 Reason: image not found 
エラー：  パッケージ '‘plyr’' をロードできませんでした
set.seed(100) ; x <- sample(c(1,2,5,7), 200, TRUE) ; y <- rnorm(200)
boxplot(y~x, at=sort(unique(x)), xlim=c(0,8))
> x <- matrix(c(1:4), nrow=2, ncol=2, byrow=T)
> (x/apply(x, 1, sum))

y <- matrix(c(1:4), nrow=2, ncol=2, byrow=F)
t(y/apply(y, 2, sum))
> (x <- matrix(c(1:4), nrow=2, ncol=2, byrow=T))
     [,1] [,2]
[1,]    1    2
[2,]    3    4
> apply(x, 1, sum)
[1] 3 7
> x/apply(x, 1, sum)
          [,1]      [,2]
[1,] 0.3333333 0.6666667
[2,] 0.4285714 0.5714286
> (y <- matrix(c(1:4), nrow=2, ncol=2, byrow=F))
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> t(y/apply(y, 2, sum))
          [,1]      [,2]
[1,] 0.3333333 0.2857143
[2,] 1.0000000 0.5714286
> apply(y, 2, sum)
[1] 3 7
> y/apply(y, 2, sum)
          [,1]      [,2]
[1,] 0.3333333 1.0000000
[2,] 0.2857143 0.5714286
> t(y/apply(y, 2, sum))
          [,1]      [,2]
[1,] 0.3333333 0.2857143
[2,] 1.0000000 0.5714286
> t(y)/apply(y, 2, sum)
          [,1]      [,2]
[1,] 0.3333333 0.6666667
[2,] 0.4285714 0.5714286
> x <- matrix(c(1:4), nrow=2, ncol=2, byrow=T)
> x/apply(x, 1, sum)
          [,1]      [,2]
[1,] 0.3333333 0.6666667
[2,] 0.4285714 0.5714286
> prop.table(x, 1)
          [,1]      [,2]
[1,] 0.3333333 0.6666667
[2,] 0.4285714 0.5714286
(x <- stats::arima(lh, c(2, 0, 1)))
f <- deriv(~x^2*y, c("x", "y"), function(x, y) { }, hessian=TRUE)
attr[f, "hessian"][x, y]
> f12 <- function(x, y) attr(f(x, y), "hessian")
> f12(1, 2)
[1] 2
> f12(4, 3)
[1] 8
A <- c("a","a","a","b","c","a","b","a","a")
B <- rep(0,length(A))
A.str <- paste(A, collapse="")
#開始位置と長さの取得
C <- gregexpr("a+", A.str)
#C.mat[1,] 開始位置 C.mat[2, ] 長さ
C.mat <- rbind(C[[1]], attr(C[[1]], "match.length"))
B[A=="a"] <- unlist(apply(C.mat, 2, function(x){c(1:x[2])}))
> my.func <- function(A)
+ {
+ 	B <- integer(length(A))
+ 	C <- rle(A == "a")
+ 	B[A == "a"] <- unlist(lapply(C$lengths[C$values], function(x) 1:x))
+ 	B
+ }
> your.func <- function(A)
+ {
+ 	B <- rep(0, length(A))
+ 	A.str <- paste(A, collapse="")
+ 	# 開始位置と長さの取得
+ 	C <- gregexpr("a+", A.str)
+ 	# C.mat[1,] 開始位置 C.mat[2, ] 長さ
+ 	C.mat <- rbind(C[[1]], attr(C[[1]], "match.length"))
+ 	B[A == "a"] <- unlist(apply(C.mat, 2, function(x) {c(1:x[2])}))
+ 	B
+ }
> a <- c("a","a","a","b","c","a","b","a","a")
> my.func(a)
[1] 1 2 3 0 0 1 0 1 2
> your.func(a)
[1] 1 2 3 0 0 1 0 1 2
> a <- sample(letters[1:5], 1000000, replace=TRUE)
> system.time(my.func(a))
   ユーザ   システム       経過  
     0.615      0.035      0.646 
> system.time(your.func(a))
   ユーザ   システム       経過  
     2.448      0.040      2.453
End.index.all <- cumsum(C$lengths)
Start.index.all <- End.index.all - (C$lengths - 1)
End.index <- End.index.all[C$values]
Start.index <- Start.index.all[C$values]
my2.func <- function(A)
{
	B <- A == "a"
	C <- rle(B)
	B[B] <- unlist(lapply(C$lengths[C$values], function(x) 1:x))
	B
}
> A <- read.table("matrix.dat")
> B <- matrix(A[1,], 3, 3)
> kappa(B, exact=TRUE)
Error in svd(z, nu = 0, nv = 0) : infinite or missing values in 'x'
if (any(!is.finite(x))) 
   stop("infinite or missing values in 'x'")
> any(!is.finite(B))
> kappa(A, exact=TRUE)
> any(!is.finite(A))
> A <- matrix(scan("matrix.dat"), ncol=9)
> any(!is.finite(A))
> any(!is.finite(B))
$ uname -sr
Linux 2.6.18-92.el5
$ R --version
R version 2.10.0 (2009-10-26)
library(XML)
data <- c("変数1の値です。","value of variable 2.")
xml <- xmlTree()
xml$addTag('データセット', close=FALSE)
xml$addTag("データ", close=FALSE)
xml$addTag("変数1", data[1])
xml$addTag("変数2", data[2])
xml$closeTag();xml$closeTag();xml$closeTag()
cat(saveXML(xml))
> cat(saveXML(xml,encoding = "UTF-8"))
<?xml version="1.0"?> 

<データセット>
  <データ>
    <変数1>変数1の値です。</変数1>
    <変数2>value of variable 2.</変数2>
  </データ>
</データセット>
function() {
  nt <- 10000
  pt <- 0
  for(kk in 1:nt){
    set.seed(12+kk*144)
    xx1 <- abs(rnorm(50, mean = 0, sd = 1))
    xx2 <- abs(rnorm(50, mean = 0, sd = 1))
    xx <- c(xx1, - xx2)
    ttest1 <- t.test(xx, mu = 0)
    pp <- ttest1$p.value
    if (pp < 0.05){
      pt <- pt +1
    }
  }
  pt <- pt/nt
  print(pt)
}
plot(res$baseMeanA, res$baseMeanB, pch=20, cex=.1, col = ifelse( res$padj < .1, "red", "black" ))
install.packages(NULL, .libPaths()[1L], dependencies = NA, type = type) 中で警告がありました： 
 'lib = "C:/PROGRA~1/R/R-2.12.1/library"' は書き込み可能ではありません  
以下にエラー sprintf(msg, userdir) : 
  <83>ｼ<e3>Ν縺吶ｋ縺溘ａ縺ｫ蛟倶ｺｺ逧<84>縺ｪ繝ｩ繧､繝悶Λ繝ｪ 
'%s' 
繧剃ｽ懊ｊ縺溘＞縺ｮ縺ｧ縺吶°<ef>ｼ<9f> に不正なマルチバイト文字があります 
arima(時系列,order=(1,1,0))
> set.seed(20100409)
> d <- data.frame(y=rnorm(24),
+                 A=gl(2, 12, labels=c("A1", "A2")),
+                 B=gl(3, 4, 12, labels=c("B1", "B2", "B3")))
> d$temp <- as.integer(d$A)*10+as.integer(d$B) # 策を弄する
> d$C <- factor(d$temp,                        # 新たに factor 変数を作る
+               labels=c("A1B1", "A1B2", "A1B3", "A2B1", "A2B2", "A2B3"))
> d
              y  A  B temp    C
1   0.807748961 A1 B1   11 A1B1
2   0.513417877 A1 B1   11 A1B1
3  -0.156009578 A1 B1   11 A1B1
4  -0.192902834 A1 B1   11 A1B1
5   0.043651547 A1 B2   12 A1B2
6  -1.530954060 A1 B2   12 A1B2
7  -1.984433103 A1 B2   12 A1B2
8   0.416585172 A1 B2   12 A1B2
9  -0.001342267 A1 B3   13 A1B3
10  0.364838567 A1 B3   13 A1B3
11  1.182524866 A1 B3   13 A1B3
12 -1.087637431 A1 B3   13 A1B3
13 -1.552197067 A2 B1   21 A2B1
14  0.323540547 A2 B1   21 A2B1
15 -0.635204370 A2 B1   21 A2B1
16 -2.889526450 A2 B1   21 A2B1
17  0.223709445 A2 B2   22 A2B2
18  0.118526662 A2 B2   22 A2B2
19  0.172598117 A2 B2   22 A2B2
20  1.846951151 A2 B2   22 A2B2
21  0.673207028 A2 B3   23 A2B3
22  1.080020322 A2 B3   23 A2B3
23  1.383334771 A2 B3   23 A2B3
24 -0.260967601 A2 B3   23 A2B3
> bartlett.test(y~C, d)

	Bartlett test of homogeneity of variances

data:  y by C 
Bartlett's K-squared = 3.125, df = 5, p-value = 0.6807
> set.seed(20100409)
> d <- data.frame(y=rnorm(24),
+                 A=gl(2, 12, labels=c("A1", "A2")),
+                 B=gl(3, 4, 12, labels=c("B1", "B2", "B3")))
> d$C <- factor(paste(d$A, d$B, sep=""), # 今の例ではこの方がよいかな？
+               labels=c("A1B1", "A1B2", "A1B3", "A2B1", "A2B2", "A2B3"))
> d
              y  A  B    C
1   0.807748961 A1 B1 A1B1
2   0.513417877 A1 B1 A1B1
 途中省略
23  1.383334771 A2 B3 A2B3
24 -0.260967601 A2 B3 A2B3
> bartlett.test(y~C, d)

	Bartlett test of homogeneity of variances

data:  y by C 
Bartlett's K-squared = 3.125, df = 5, p-value = 0.6807
cat("Hello World")
R --vanilla < test.R > test.Rout
以下にエラー match.call() : 
  'match.call' がその中から呼び出されたクロージャを見付けることができません。
formula <- y ~ 1 + x
xx <- c(1.386294, 2.079442 ,2.772589, 3.465736)
yy <- c(0.06254298 ,0.71444544, 1.27103337, 1.28339677)
ww <- c( 1.0000000 ,0.5000000, 0.3333333, 0.2500000)

data <- data.frame(x = xx, y = yy, w = ww)
subset <- 0
weights <- ww
na.action <- 0 
method <- "qr"
model <- TRUE 
x <- FALSE
y <- FALSE
qr <- TRUE
singular.ok <- TRUE 
contrasts <- NULL
offset <- 0  

ret.x <- x
ret.y <- y
cl <- match.call()
mf <- match.call(expand.dots = FALSE)
> cl <- match.call()
 エラー： match.call() was called from outside a function
> mf <- match.call(expand.dots = FALSE)
 エラー： match.call() was called from outside a function
a<-rnorm(10000, sd=2, mean=10)
b<-rnorm(10000, sd=1, mean=0)
var(a)+var(b)
var(a-b)
Cochran's C test  A test, introduced by Cochran in a 1941 paper, for equality
of variance in m  normal populations having unknown means and having sample 
variances denoted by . Let be the ... 
> Map (function (x) x(), Map (function (y) {function () y}, c(1, 2, 3))) ## list(1, 2, 3) が返ることを期待
[[1]]
[1] 3

[[2]]
[1] 3

[[3]]
[1] 3
>
> Map (function (x) x(), Map (function (y) {y;function () y}, c(1, 2, 3)))
[[1]]
[1] 1

[[2]]
[1] 2

[[3]]
[1] 3

>
          :
          :
df<-data.frame(lat=c(35,38,37,38),
               long=c(156,165,150,148),
               year=c(2000,2000,2001,2001)
               )
cat(deparse(df),fill=T)

library (maps);
par(ask=T)
lapply(unique(df$year),function(y)
       {
         cat(y,fill=T)
         longlat<-df[df$year==y,c("long","lat")]
         plot (longlat,
               xlim = c(135, 155),
               ylim = c(32, 50),
               axes = F, main = y, pch = 20)
         map ("world2", xlim = c(135, 155), ylim = c(32, 50), add = T)
         map.axes()
         dev.copy(pdf, file = paste(y,".pdf",sep=""))
         dev.off()
       }
)
x <- c(1:10)
y <- c(10,10,9,8,5,3,2,1,0,0)
plot(y~x)
> f <- function(x) class(x)
> 
> setMethod("f", "B0", function(x) c(x@b0^2, callNextMethod()))
 以下にエラー sprintf(gettext(fmt, domain = domain), ...) : 
   引数が少なすぎます 
Calls: setMethod -> message -> gettextf -> sprintf
 実行が停止されました 
--- config.site ----
CFLAFGS="-g -O2 -wd188 -ip -std=c99"
F77=ifort
FFLAGS="-ieee-fp -g -O3"
CXX=icpc
CXXFLAGS="-g -O3 -ieee-fp"
FC=ifort
FCFLAGS="-g -O3 -ieee-fp"
ICC_LIBS=/opt/intel/Compiler/11.0/074/lib/intel64
IFC_LIBS=/opt/intel/Compiler/11.0/074/lib/intel64
LDFLAGS="-L$ICC_LIBS -L$IFC_LIBS -L/usr/local/lib64"
SHLIB_CXXLD=icpc

----- configure 実行引数 -----
OMP_NUM_THREADS=8

MKL_LIB_PATH=/opt/intel/Compiler/11.0/074/mkl/lib/em64t

MKL="   -L${MKL_LIB_PATH}                               \
            -Wl,--start-group                               \
                    ${MKL_LIB_PATH}/libmkl_intel_lp64.a        \
                    ${MKL_LIB_PATH}/libmkl_intel_thread.a     \
                    ${MKL_LIB_PATH}/libmkl_core.a           \
            -Wl,--end-group                                 \
            -liomp5 -lguide -lpthread -lgomp"

./configure --with-lapack="$MKL" --with-blas="$MKL" 
patch <  methods.patchのパス名　
(cd po; make remove-potcdate.sed; make en@quot.insert-header) 
(cd po; make methods.pkg-update) 
$ grep -r "co2" /usr/lib/R/* | grep "/help"
      for (j in 1:nrow(x$est)) {
      		temp <- as.integer(unlist(strsplit(rownames(x$est)[j], " ")))
      		temp <- ftime[fstatus == 0 & group == temp[1]]
      		temp2 <- timepoints(fit, temp)
     		mapply(function(x, y) segments(x, y, x, y+0.02, col=col[j]),
                       sort(temp), temp2$est[j,])
      }
      legend("topleft", legend =  rownames(x$est), x.intersp = 2, 
             bty = "n", xjust = 1, col = col, lty = lty, lwd = lwd)
      for (j in 1:nrow(x$est)) {
      		temp <- as.integer(unlist(strsplit(rownames(x$est)[j], " ")))
      		temp <- ftime[fstatus == 0 & group == temp[1]]
      		temp2 <- timepoints(fit, temp)
     		mapply(function(x, y) segments(x, y, x, y+0.02, col=col[j]), 
                      sort(temp), temp2$est[j,])
      }
      out <- list(test = tests, est = tfit$est, se = sqrt(tfit$var))
>  ( x <- read.csv("bmt.csv", header=TRUE, sep=";") )
   dis ftime status
1    0    13      2
2    0     1      1
3    0    72      0
4    0     7      2
　　途中省略
34   1    32      0
35   0    12      1
> CumIncidence(x$ftime, x$status, x$dis)

+-------------------------------------------------------------------+
| Cumulative incidence function estimates from competing risks data |
+-------------------------------------------------------------------+
Test equality across groups:
  Statistic  p-value df
1     1.302 0.253915  1
2     7.082 0.007785  1 

Estimates at time points:
          0     10     20     30     40     50     60     70
0 1 0.05882 0.1176 0.1765 0.1765 0.1765 0.1765 0.1765 0.1765
途中省略
0 2 0.05882 0.12631 0.12667 0.12077 0.12077 0.12077 0.12077 0.12077
1 2 0.00000 0.11558 0.11558 0.11558 0.11558 0.11558 0.11558 0.11558
以上で終了し，tick mark のついている図が描かれており，warning なんて出ないのですけど。
> DF <- data.frame(A=letters[1:4], B=runif(4), C=rnorm(4))
> DF
  A         B          C
1 a 0.1429617  0.8767259
2 b 0.7656706  0.2077312
3 c 0.5311202 -0.6043531
4 d 0.9375921  1.0804659
> rownames(DF) <- DF[,1]
> DF
  A         B          C
a a 0.1429617  0.8767259
b b 0.7656706  0.2077312
c c 0.5311202 -0.6043531
d d 0.9375921  1.0804659
> DF <- DF[,-1]
> DF
          B          C
a 0.1429617  0.8767259
b 0.7656706  0.2077312
c 0.5311202 -0.6043531
d 0.9375921  1.0804659
> DF.new <- data.frame(DF, row.names="A")
> DF.new
           B          C
a 0.79124741  1.2432360
b 0.03667521 -1.0248516
c 0.78015687  0.1808611
d 0.30211288  0.2649450
(plotのスクリプト省略)
par(xpd=NA)
legend(locator(1),c("a-pro","a-height","b-pro",
  "b-height"),col=c("black","black","Royalblue","Royalblue"),
  lwd=c(2,3,1,1.5),pch=c(1,NA,1,NA))
legend(locator(1),c("a=pro","a-height","b-pro",
  "b-height"),col=c("black","black","Royalblue","Royalblue"),
  lty=c(1,0,1,0),lwd=c(2,0,1,0),pch=c(1,15,1,15))
   school y x
 1      1 4 1
 2      1 3 2
 3      1 2 3
 4      5 4 2
 5      5 5 3
 6      5 3 4
 7      8 4 3
 8      8 4 3
 9      8 1 2
> A <- data.frame(school=sample(1:3, 10, 1),x =sample(1:4,10,1), y=sample(10:14,10,1))
> A
   school x  y
1       2 1 12
2       1 4 10
3       1 1 12
4       1 1 13
5       1 1 14
6       3 1 12
7       3 4 12
8       2 4 13
9       3 3 13
10      2 1 11
> A1 <- A[A$school==1,,]
> A1
  school x  y
2      1 4 10
3      1 1 12
4      1 1 13
5      1 1 14
> A1[sample(seq(ncol(A1)),3),,] # school=1 から3ケースを非復元ランダム抽出
  school x  y
2      1 4 10
4      1 1 13
3      1 1 12
> A1[order(sample(seq(ncol(A1)),3)),,]  # 行番号を整列させたければ
  school x  y
2      1 4 10
3      1 1 12
4      1 1 13
> install.packages("hoge",lib="~/R/library",
  configure.args = "CPPFLAGS = -I~/include LDFLAGS = -L~/lib")
> library()
 > .libPaths()
[1] "~/R/library"              "/usr/local/lib/R/library"
* Installing * source packages 'gsl'...
checking for gcc... gcc
checking for C compiler dafault output... configure: error C compiler cannot create executables
$ tar xzf gsl_1.8-8.tar.gz && cd $(tar tzf gsl_1.8-8.tar.gz |head -n1)
$ pwd
/tmp/gsl
$ ./configure LDFLAGS="/opt/local/lib" CPPFLAGS="-I/opt/local/include" 
checking for gcc... gcc
checking for C compiler default output... configure: error: C compiler cannot  create executables
See `config.log' for more details.
$ ./configure LDFLAGS="-L/opt/local/lib" CPPFLAGS="-I/opt/local/include" 
checking for gcc... gcc
checking for C compiler default output... a.out
checking whether the C compiler works... yes
[以下略]
#include <math.h> 
void test1(double *seq, int *seqn, double *sf)
{
	int i;
	for(i=0;i<=*seqn-1;i++){
	sf[i] = sin(seq[i]);
	}
}
Rcmd SHLIB test1.c
dyn.load("test1.dll")
a<-seq(0,pi,length=10)
.C("test1",as.double(a),length(a),sf=double(length(a)))$sf
         PC1         PC2          PC3
a -0.2138188 -0.09690295 -0.056538208
b -0.1983966 -0.06068939  0.108130244
c -0.3581248 -0.04182101  0.157024270
d -0.2755906 -0.07243737  0.004878913
e -0.4614072 -0.02512825  0.006534084
f -0.3580950 -0.11927931  0.036777825
g  0.2977400  0.60951743  0.180188932
h  0.3043281  0.62594929  0.290022187
i  0.1160385 -0.55632764  0.073677929
j  0.7941420 -0.39650995  0.289428417
k  0.4994019 -0.38653550 -0.234907993
l  0.3296465  0.32195986 -0.691875477
m -0.2982350 -0.06277819 -0.159138712
n -0.1776291  0.26098297 -0.004202413
pcalabel<-row.names(pca3)     #各プロットにつけたいラベルの名前。
rgl.bg(color=c("white","black"))
rgl.lines(c(-1,1),0,0,color="black")
rgl.lines(0,c(-1,1),0,color="black")
rgl.lines(0,0,c(-1,1),color="black")
rgl.texts(c(1,0,0),c(0,1,0),c(0,0,1),
          text=c("PC1","PC2","PC3"),color="blue")
rgl.points(pca3,color=heat.colors(1000),size=2)
rgl.texts(pca3[,1],pca3[,2],pca3[,3],pcalabel,adj = 1)  #ここでラベル付けにトライ
rp.plot3d {rpanel}	R Documentation
Interactive display of a plot of three variables
Description

This function produces a scatterplot of three variables, using the rgl package
for three-dimensional display.
Usage

 rp.plot3d(x, y, z, xlab  = NA, ylab = NA, zlab = NA, 
                  axes = TRUE, type = "p", size = 3, col = "red",
                  xlim = NA, ylim = NA, zlim = NA, ...)

Arguments
x,y,z 	vectors of observed values.
xlab 	a character variable used for the first axis label.
ylab 	a character variable used for the second axis label.
zlab 	a character variable used for the third axis label.
axes 	a logical variable determining whether the axes are shown.
type  a character variable controlling the type of plotting. 
       If the value is  set to "n", the points are not plotted.
size   the size of the plotted points.
col 	the colour of the plotted points.
xlim 	the plotting range for the first variable.
ylim 	the plotting range for the second variable.
zlim 	the plotting range for the third variable.
... 	other rgl parameters which control the appearance of the plotted points.
Details

The plot is produced by appropriate calls to the rgl package. This allows 
interactive control of the viewing position. Other objects may subsequently be 
added to the plot by using rgl functions and data which are centred and scaled 
by the returned values indicated below.
library(rgl)
open3d()
pcalabel<-row.names(pca3)     #各プロットにつけたいラベルの名前。
rgl.bg(color=c("white","black"))
rgl.lines(c(-1,1),0,0,color="black")
rgl.lines(0,c(-1,1),0,color="black")
rgl.lines(0,0,c(-1,1),color="black")
rgl.texts(c(1,0,0),c(0,1,0),c(0,0,1),
         text=c("PC1","PC2","PC3"),color="blue")
rgl.points(pca3,color=heat.colors(1000),size=2)
text3d(pca3[,1],pca3[,2],pca3[,3],text=pcalabel,adj =1)  #text3d()を利用
A1 <- c(9,8,19,8,2,23,2,10)
A1label<-c("AAAAAAA","BBBBBBBBBBBBBBBBBBBBBBB","CCCCCCCCCCCCCC",
           "DDDDDDDDDDDD","EEEEEEEEEEEEEEEEEE","FFFFFFFFFFFFFFFFFF",
           "GGGGGGGGG","HHHH")
barplot(A1,names=A1label,horiz=T)
> x
   StandardLength Feed Area
1             100   10    N
2             100   11    N
3             100   12    N
4             100   13    N
5             100   14    N
6             100   15    N
7             100   16    N
8             100   17    N
9             100   18    N
10            100   19    N
11            100   20    N
12            100   10    N
13            150   21    W
14            150   22    W
15            150   23    W
16            150   24    W
17            150   25    W
18            150   26    W
19            150   27    W
20            150   28    W
21            150   29    W
22            200   31    S
23            200   32    S
24            200   33    S
25            200   34    S
26            200   35    S
27            200   36    S
28            200   37    S
29            200   38    S
30            200   39    S
> result=lm(StandardLength~Feed+Area)　 #StandardLengthを独立変数、
                                         FeedとAreaを従属変数として重回帰
> summary(result)
Call:
lm(formula = StandardLength ~ Feed + Area)

Residuals:
       Min         1Q     Median         3Q        Max 
-4.422e-14 -1.798e-15  1.179e-15  4.160e-15  9.431e-15 

Coefficients:
             Estimate Std. Error   t value Pr(>|t|)    
(Intercept) 1.000e+02  9.435e-15 1.060e+16  < 2e-16 ***
Feed        2.209e-15  6.171e-16 3.579e+00  0.00139 ** 
AreaS       1.000e+02  1.332e-14 7.506e+15  < 2e-16 *** #問題2の該当箇所。
AreaW       5.000e+01  7.750e-15 6.452e+15  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 9.815e-15 on 26 degrees of freedom
Multiple R-Squared:     1,      Adjusted R-squared:     1 
F-statistic: 1.791e+32 on 3 and 26 DF,  p-value: < 2.2e-16 

> sd(Feed)　　#問題3の該当箇所。
               今となっては恥ずかしいかぎりです。
[1] 9.123986
> sd(AreaS)
 以下にエラーsd(AreaS) :  オブジェクト "AreaS" は存在しません
             #問題1の該当箇所。
              sd()が計算されないため、
              前に投稿したやり方では標準化できませんでした。
StandardLength の予測値 = 100 + 0*Feed + 0*AreaN + 50*AreaW + 100*AreaS
StandardLength の予測値バージョンII = 150 + 0*Feed - 50*AreaN + 0*AreaW + 50*AreaS
# 2 群のヒストグラム
hist2 <- function(     x1,        # 第一群のデータ
                       x2,        # 第二群のデータ
                       brks=NULL, # 階級分割点
                       ...)       # barplot に引き渡す任意の引数
{
       if (is.null(brks)) {       # 階級分割点が与えられないときには，適切に設定
               brks <- hist(c(x1, x2), right=FALSE, plot=FALSE)$breaks
       }
       c1 <- hist(x1, breaks=brks, right=FALSE, plot=FALSE)$counts # 度数1
       c2 <- hist(x2, breaks=brks, right=FALSE, plot=FALSE)$counts # 度数2
       barplot(rbind(c1, c2), beside=TRUE, space=c(0, 0.2),        # 棒を並べて描く
               names.arg=brks[-length(c1)],                        # 横軸の目盛りラベル等
               axisnames=TRUE, axis.lty=1, ...)
}
data(USArrests)
USArrests
str(USArrests)
x1<-USArrests$Murder
x2<-USArrests$Rape
str(x1)
str(x2)
hist2(x1,x2)
if (is.data.frame(x1)) x1 <- as.matrix(x1)
if (is.data.frame(x2)) x2 <- as.matrix(x2)
x <- data.frame(matrix(1:12,ncol=3,dimnames=list(1:4,c("AAA","BBB","CCC"))))
x$"AAA"
x$AAA
x$AA
x$A
> x <- data.frame(matrix(1:12,ncol=3,dimnames=list(1:4,c("AAA","BBB","CCC"))))
> x["AAA"]
  AAA
1   1
2   2
3   3
4   4
> x[,"A"]
NULL
> x[[1]]      # 成分をベクトルとして取り出す「古風な」しかた
[1] 1 2 3 4
> x[["AAA"]]
[1] 1 2 3 4
> x[["A"]]    # 二重鈎括弧演算子も補完機能を持つことも今回初めて知った　
[1] 1 2 3 4
R : Copyright 1997, Robert Gentleman and Ross Ihaka
Version 0.49 Beta (April 23, 1997)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type "license()" for details.

[Previously saved workspace restored]

Error: workspace file corrupted -- no data loaded
Error in device("X11", as.character(c(display[1], paper)), as.double(c(width,  : unable to start device X11
> x <- data.frame(matrix(1:12,ncol=3,dimnames=list(1:4,c("AAA","BBB","CCC"))))
> x$"AAA"
[1] 1 2 3 4
> x$AAA
[1] 1 2 3 4
> x$AA
[1] 1 2 3 4
> x$A
[1] 1 2 3 4
> xy <- expand.grid(1:100, 1:100)  # これが予測位置グリッド座標
> names(xy) <- c("x", "y")
> g.dummy <- gstat(formula = z ~ 1, locations = ~x + y, dummy = TRUE, 
                   beta = 0, model = vgm(1, "Exp", 15), nmax = 20)
> yy <- predict(g.dummy, newdata = xy, nsim = 4)
片方のサンプルサイズNc，もう一方のサンプルサイズNt，その比r(=Nc/Nt)
Pc, Pt は属性を持つものの割合，Qc=1-Pc, Qt=1-Pt
P =(r Pc + Pt) / (r+1)
Zα=qnorm(α/2, lower.tail=FALSE)
Zβ=qnorm(β)
Nt ={Zα*sqrt((r+1) P Q)+Zβ*sqrt(r Pt Qt +Pc Qc)}^2 /(r(Pt-Pc)^2)
> library(Hmisc)
> foo <- function(N) {
           power <- bpower(p1=P1, p2=P2, 
                           n1=N, n2=ceiling(Nratio*N), alpha=Alpha)
           if (power >= Power) { # 要求パワーを越えた時点でストップ
                                cat(N, ceiling(Nratio*N), power, fill=TRUE); stop()}
         }
> P1 <- 0.3; P2 <- 0.5; Nratio <- 0.2; Alpha=0.05; Power = 0.8
> sapply(100:10000, foo)
266 54 0.8003697
 以下にエラーFUN(100:10000[[167L]], ...) : # これはループを中断したことによる注意で無視可能
> P1 <- 0.3; P2 <- 0.5; Nratio <- 0.3; Alpha=0.05; Power = 0.8 
> sapply(100:10000, foo)
195 59 0.8002799
 以下にエラーFUN(100:10000[[96L]], ...) : 
> P1 <- 0.3; P2 <- 0.35; Nratio <- 0.3; Alpha=0.05; Power = 0.8
> sapply(100:10000, foo)
2954 887 0.8002141
 以下にエラーFUN(100:10000[[2855L]], ...) : 
> P1 <- 0.3; P2 <- 0.35; Nratio <- 1; Alpha=0.05; Power = 0.8
> sapply(100:10000, foo)
1377 1377 0.8002008
 以下にエラーFUN(100:10000[[1278L]], ...) :
power.prop.test2 <- function( p1=NULL, p2=NULL, n1=NULL, r=NULL, sig.level=0.05, power=NULL)
{
        if (sum(sapply(list(p1, p2, n1, r, sig.level, power), is.null)) != 1) {
                stop("p1, p2, n1, r, sig.level, power のどれか一つだけを NULL とする")
        }
        power.function <- quote(bpower(p1=p1, p2=p2, n1=n1, n2=ceiling(n1*r), alpha=sig.level))
        if (is.null(power)) {
                power <- eval(power.function)
        }
        else if (is.null(n1)) {
                n1 <- uniroot(function(n1) eval(power.function)-power, c(1, 1e7))$root
        }
        else if (is.null(p1)) {
                p1 <- uniroot(function(p1) eval(power.function)-power, c(1e-5, 1))$root
        }
        else if (is.null(p2)) {
                p1 <- uniroot(function(p2) eval(power.function)-power, c(1e-5, 1))$root
        }
        else if (is.null(r)) {
                r <- uniroot(function(r) eval(power.function)-power, c(1e-5, 1000))$root
        }
        else if (is.null(sig.level)) {
                sig.level <- uniroot(function(sig.level) eval(power.function)-power, c(1e-5, 0.99999))$root
        }
        else {
                stop("internal error")
        }
        METHOD <- "Power calculation of the tow-sample proportion test."
        structure(list(p1=p1, p2=p2, n1=ceiling(n1), n2=ceiling(n1*r), r=r, 
                       sig.level=sig.level,power=power, alternative="two.sided",
                       method=METHOD), class="power.htest")
}
> power.prop.test2(p1=P1, p2=P2, n1=NULL, r=0.2, power=0.8)

     Power calculation of the tow-sample proportion test. 

             p1 = 0.3
             p2 = 0.5
             n1 = 266
             n2 = 54
              r = 0.2
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

> power.prop.test2(p1=P1, p2=P2, n1=NULL, r=1, power=0.8)

     Power calculation of the tow-sample proportion test. 

             p1 = 0.3
             p2 = 0.5
             n1 = 93
             n2 = 93
              r = 1
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

> power.prop.test2(p1=P1, p2=P2, n1=266, r=NULL, power=0.8)

     Power calculation of the tow-sample proportion test. 

             p1 = 0.3
             p2 = 0.5
             n1 = 266
             n2 = 54
              r = 0.1992972
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

> power.prop.test2(p1=P1, p2=P2, n1=266, r=0.2, power=NULL)

     Power calculation of the tow-sample proportion test. 

             p1 = 0.3
             p2 = 0.5
             n1 = 266
             n2 = 54
              r = 0.2
      sig.level = 0.05
          power = 0.8003697
    alternative = two.sided
## 二項分布 B(n,p),B(m,p) に従う独立標本 N 個を有意水準 1-alpha で
## 両側検定したとき棄却される割合(標本検出力)を計算する関数
> foo <- function(N,n,p,m,q, alpha) {
    X <- rbinom(N,n,p)                 # B(n.p) に従う標本N個
    Y <- rbinom(N,m,q)                 # B(m,q) に従う標本N個
    phat <- (X+Y)/(n+m)                # 帰無仮説p=qの下での共通比率の推定値 
    S <- sqrt((1/n+1/m)*phat*(1-phat)) # 帰無仮説の下での(X/n-Y/m)の標準偏差推定値   
    z <- qnorm(1-alpha/2)              # 両側検定棄却限界値
    mean(abs(X/n - Y/m)/S > z)         # N 回の検定で棄却される標本比率
  } 
> bar <- function(n) { # m=0.2*n のケース
    x <- foo(2e6, n, p, ceiling(0.2*n), q, alpha)
    cat(k, ceiling(0.2*n), x, fill=TRUE)
  }
# 実行例
> p=0.3; q=0.5; alpha=0.05; beta=0.8
> set.seed(1); sapply(260:300, bar)
260 52 0.7882775
261 53 0.793533
262 53 0.795987
263 53 0.791743
264 53 0.7943355
265 53 0.797369
266 54 0.802263  # 標本検出力が初めて80%を越える標本数(bpower関数の場合,266 54 0.8003697,とほぼ一致)
267 54 0.799412  
268 54 0.801248
269 54 0.8035085
.................
>library(mgcv); data(airquality)
>airq.gam<-gam(Ozone~s(Solar.R)+s(Wind)+s(Temp),data=airquality)
>summary(airq.gam)
17
5
8
26
15
26
15
26
5
0
26
...
17　　サブグループ１
5　　サブグループ１
8　　サブグループ１
26　　サブグループ１
15　　サブグループ１
26　　サブグループ２
15 　　サブグループ２
26 　　サブグループ２
5 　　サブグループ２
0 　　サブグループ２
26 　　サブグループ３
...
サブグループ１　　17, 5, 8, 26, 15
サブグループ２　　26, 15, 26, 5, 0
サブグループ３　　26, ...
> x
 [1]  1 14 19  5  8  4  8 10  1  4 16  7  3 18  8 20  4  1  8 12  7  8  3 12 10
[26] 19  5 15  8  9
> xx <- matrix(x, 6, 5, byrow=TRUE)
> xx
     [,1] [,2] [,3] [,4] [,5]
[1,]    1   14   19    5    8
[2,]    4    8   10    1    4
[3,]   16    7    3   18    8
[4,]   20    4    1    8   12
[5,]    7    8    3   12   10
[6,]   19    5   15    8    9
> as.data.frame(xx)
  V1 V2 V3 V4 V5
1  1 14 19  5  8
2  4  8 10  1  4
3 16  7  3 18  8
4 20  4  1  8 12
5  7  8  3 12 10
6 19  5 15  8  9
x <- rnorm(1:20)
y <- c(rep(1, 4), rep(2, 6), rep(3, 10))
z <- data.frame("dat"=x, "sub.g"=as.factor(y))
> filled.contour(Time,Frequency,t(z),nlevel=50,color.palette=topo.colors)
以下にエラーdiff(x) : 既定引数の再帰的な参照です
#include <R.h>
void countnan(double *x, int *n, int *count)
 {
   int i;
   for (i=0,*count=0;i<*n;i++) 
     if (ISNAN(x[i])) (*count)++;
 }

> countnan <- function(x){ 
 n=length(x);
 .C("countnan",arg1=as.double(x),
  arg2=numeric(n),arg3=numeric(1))$arg3
}
> x<-c(1, NA , 3 , 4)
> countnan(x)
 以下にエラーcountnan(x) :  外部関数の呼び出し(引数 1) 中に NA/NaN/Inf があります
> countnan <- function(x){ 
 .C("countnan",arg1=as.double(x),
  arg2=length(x),arg3=integer(1),NAOK=TRUE)$arg3 }
> countnan(c(1,2,NA,4))
[1] 1
#include <R.h>
void countnan(double *x, int *n, int *count)
{
 int i;
 for (i=0,*count=0;i<*n;i++) 
 if (ISNAN(x[i])) (*count)++;
 (*count)--;
 if (ISNAN(R_NaReal)) (*count)++;
}
("leaf","leaf","node","node")
DNapplyTest <- function(d) {
	if (!is.leaf(d)) {
		for (j in seq_along(d)) Recall(d[[j]])
		cat("node:",attr(d,"height"),"?n")
	}else{
		cat("leaf:",labels(d),"?n")
	}
}

> DNapplyTest(as.dendrogram(hclust(dist(USArrests))))
leaf:Florida
leaf:North Carolina
node:38.52791
...
node:87.32634
node:168.6114
node:293.6228
## 文字列が長すぎるので適宜改行しています
> dhc <- as.dendrogram(hclust(dist(USArrests)))
> x <- capture.output(dendrapply(dhc, function(n) {cat(labels(n),"?n");  cat(attr(n,"height"),"?n")}))
> xx <- x[1:198] # 手作業で末尾につく無意味な部分を削除
> xx[5:20]
 [5] "Florida North Carolina "
 [6] "38.52791 "
 [7] "Florida "
 [8] "0 "
 [9] "North Carolina "
[10] "0 "
[11] "Delaware Alabama Louisiana Alaska Mississippi South Carolina Maryland
Arizona New Mexico California  Illinois New York Michigan Nevada "
[12] "64.99362 "
[13] "Delaware Alabama Louisiana Alaska Mississippi South Carolina "
[14] "48.72515 "  
lines(c(apx[c1], hc$height[i], hc$height[i], apx[c2]),
      c(apy[c1], apy[c1], apy[c2], apy[c2]))
f<-function(X){
  ret<-NULL
  f<-function(x){
    if(is.leaf(x)){
      attr(x,"label")
    }else{
      r<-c(sapply(x,f),recursive=T)
      ret<<-rbind(ret,list(leafs=r,size=attr(x,"members"),height=attr(x,"height")))
      r
    }
  }
  f(X)
  ret
}
r<-f(as.dendrogram(hclust(dist(USArrests), "ave")))
関数定義
dend2nodelist <- function(hc) {
	GetElement <- function(i, j) {
		if (i < 0) {
			if (j < 0) return(c(hc$labels[-i], hc$labels[-j]))
			else return(c(hc$labels[-i], GetElement(hc$merge[j,1], hc$merge[j,2])))
		} else {
			if (j < 0) return(c(GetElement(hc$merge[i,1], hc$merge[i,2]), hc$labels[-j]))
			else return(c(GetElement(hc$merge[i,1], hc$merge[i,2]), 
                                     GetElement(hc$merge[j,1], hc$merge[j,2])))
		}
	}
	n <- length(hc$height)
	leafs <- NULL
	size <- integer(n)
	for (i in 1:n) {
		labels <- GetElement(hc$merge[i,1], hc$merge[i,2])
		size[i] <- length(labels)
		leafs[i] <- list(labels)
	}
	return(cbind(leafs=leafs, size=size, height=hc$height))
}
実行例
> hc <- hclust(dist(USArrests[1:10,]))
> result <- dend2nodelist(hc)
> result
      leafs        size height  
 [1,] Character,2  2    16.80625
 [2,] Character,2  2    23.19418
 [3,] Character,2  2    25.09303
 [4,] Character,3  3    36.73486
 [5,] Character,3  3    45.18296
 [6,] Character,3  3    60.98073
 [7,] Character,6  6    77.19741
 [8,] Character,9  9    148.7357
 [9,] Character,10 10   226.303 
> result[1,]
$leafs
[1] "Alabama"  "Delaware"

$size
[1] 2

$height
[1] 16.80625

> result[5,]
$leafs
[1] "Alaska"   "Alabama"  "Delaware"

$size
[1] 3

$height
[1] 45.18296

> result[9,]
$leafs
 [1] "Connecticut" "Florida"     "Arizona"     "California" 
 [5] "Arkansas"    "Colorado"    "Georgia"     "Alaska"     
 [9] "Alabama"     "Delaware"   

$size
[1] 10

$height
[1] 226.303
dend2nodelist <- function(hc) {
 	temp <- function(hc, j) GetElement(hc$merge[j,1], hc$merge[j,2])
 	temp2 <- function(hc, j) if (j < 0) hc$labels[-j] else temp(hc, j)
 	GetElement <- function(i, j) if (i < 0) c(hc$labels[-i], temp2(hc, j)) 
                                    else c(temp(hc,i), temp2(hc, j))
 	leafs <- size <- integer(n <- length(hc$height))
 	for (i in 1:n) {size[i] <- length(labels <- temp(hc, i)); leafs[i] <- list(labels) }
 	return(cbind(leafs=leafs, size=size, height=hc$height))
}
   info: From hessian inversion. Should be 0. If not, no variances
         could be estimated. You could try fixing sigma at the
         estimated value and rerun.
> library(glmmML)
> data(anorexia, package="MASS")
> id <- rep(1:nrow(anorexia))
> class <- c(rep(1,10), rep(0,62))
> glmmML(class ~ Prewt, data=anorexia, cluster=id, family=binomial)
info = 3

Call:  glmmML(formula = class ~ Prewt, family = binomial, data = anorexia, cluster = id) 

                  coef se(coef)  z Pr(>|z|)
(Intercept) -12.124711       NA NA       NA
Prewt        -0.001919       NA NA       NA

Standard deviation in mixing distribution:  48.71 
Std. Error:                                 NA 

Residual deviance: 31.01 on 69 degrees of freedom       AIC: 37.01 
Warning message:
Hessian non-positive definite. No variance! 
> 
> glmmML(class ~ Prewt, data=anorexia, cluster=id, family=binomial, 
       fix.sigma=TRUE)

Call:  glmmML(formula = class ~ Prewt, family = binomial, data = anorexia,      
      cluster = id, fix.sigma = TRUE) 

                 coef se(coef)        z Pr(>|z|)
(Intercept) -1.772713  5.59883 -0.31662    0.752
Prewt       -0.001720  0.06783 -0.02536    0.980 

Standard deviation in mixing distribution:  0 
Std. Error:                                 0 

Residual deviance: 57.99 on 69 degrees of freedom 	AIC: 61.99 
> glmmML(class ~ Prewt, data=anorexia, cluster=id, family=binomial, 
      start.sigma=1e-4)

Call:  glmmML(formula = class ~ Prewt, family = binomial, data = anorexia,      
      cluster = id, start.sigma = 1e-04) 

                 coef se(coef)        z Pr(>|z|)
(Intercept) -1.686092  5.46507 -0.30852    0.758
Prewt       -0.001680  0.06621 -0.02538    0.980

Standard deviation in mixing distribution:  0.0002467 
Std. Error:                                 406.6 

Residual deviance: 58.02 on 69 degrees of freedom 	AIC: 64.02 
The optimization may not converge with the default value of start.sigma. In 
that case, try different start values for sigma. If still no convergence, 
consider the possibility to fix the value of sigma at several values and study 
the profile likelihood.
F=lm(Z~X)
G=lm(Z~X-1)
Z=c(1,2,3,4,5,6,7,8,9,10)
X=c(1,3,4,2,5,6,7,8,10,10)
F=lm(Z~X)
G=lm(Z~X-1)
summary(F)
Call:
lm(formula = Z ~ X)
Residuals:
Min     1Q       Median  3Q      Max
1.11283 -0.47400 0.09181 0.27600 1.80531
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.35841  0.59297    0.604   0.562
          X 0.91814  0.09329    9.842   9.56e-06 ***
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.887 on 8 degrees of freedom
Multiple R-Squared: 0.9237, Adjusted R-squared: 0.9142
F-statistic: 96.86 on 1 and 8 DF, p-value: 9.561e-06
summary(G)
Call:
lm(formula = Z ~ X - 1)
Residuals:
Min    1Q      Median 3Q     Max
0.9035 -0.5006 0.1770 0.2494 2.0644
Coefficients:
  Estimate Std. Error t value Pr(>|t|)
X 0.96782  0.04255    22.75   2.91e-09 ***
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8552 on 9 degrees of freedom
Multiple R-Squared: 0.9829, Adjusted R-squared: 0.981
F-statistic: 517.5 on 1 and 9 DF, p-value: 2.909e-09
> attach(learning)
> learning.glm<-glm(def~lastprice+volume+size+rtn+next1rtn+lsize+epr1
                    +greps1+dy+bpr+wcr+rrl+rrs+cfpr+roe+spr+ebitda+igr
                    +sgrowth+esr1+lev+cgrowth+vol+roa1,family=binomial)
> summary(learning.glm)
> detach(learning)
> test<-read.csv("property200006.csv",header=TRUE)
> attach(test)
> h0 <- read.delim("c:/2000R2.txt")
> fs <- h0$fru/h0$flo
> result <- glm(fs~h0$ftl+h0$ptl+h0$ptw+h0$sta+h0$pis+h0$blt+h0$wit, weights=h0$flo)
> result2 <- step(result)
Start:  AIC= -2819.3 
 fs ~ h0$ftl + h0$ptl + h0$ptw + h0$sta + h0$pis + h0$blt + h0$wit 

         Df Deviance      AIC
- h0$wit  1    22.05 -2821.19
- h0$ptw  1    22.05 -2821.16
- h0$sta  1    22.06 -2821.13
- h0$blt  1    22.12 -2820.58
- h0$ptl  1    22.14 -2820.44
<none>         22.04 -2819.30
- h0$pis  1    22.33 -2818.70
- h0$ftl  1    23.20 -2811.30
Error in step(result) : number of rows in use has changed: remove missing values?
「2006-11-21 10:48:31.215 R[261] CFLog (21): dyld returns 2 when trying to load /Users/（ユーザ名）
/Library/ScriptingAdditions/YouHelper.osax/Contents/MacOS/YouHelper」
Error in X11(paste("jpeg::", quality, ":", filename, sep = ""), width,  : 
        unable to start device JPEG
In addition: Warning message:
unable to open connection to X11 display ''
Execution halted
> a<-adult/colonysize
> yy<-glm(cbind(x,y)~a+b+c+d,binomial)
> sub("foo?.bar.baz", "match", c("foo.bar.baz", "foo.bar,baz", "foo,bar.baz", "foo,bar,baz"))
[1] "match" "match" "match" "match"
Warning messages:
1: '?.' is an unrecognized escape in a character string 
2: '?.' is an unrecognized escape in a character string 
> sub("foo.bar.baz", "match", c("foo.bar.baz", "foo.bar,baz", "foo,bar.baz", "foo,bar,baz"))
[1] "match" "match" "match" "match"
> sub("foo??.bar.baz", "match", c("foo.bar.baz", "foo.bar,baz", "foo,bar.baz", "foo,bar,baz"))
[1] "match"    "match"    "foo,bar.baz"    "foo,bar,baz"
以下にエラーodbcQuery(channel, query, rows_at_time) : 
Calloc がメモリー (263168 of 1) を割りあてられませんでした
追加情報: Warning message:
Reached total allocation of 1015Mb: see help(memory.size)
> m <- structure(c(1, 0.1, 0.8, 0.1, 1, -0.7, 0.8, -0.7, 1),
+ .Dim = c(3, 3))
> m
     [,1] [,2] [,3]
[1,]  1.0  0.1  0.8
[2,]  0.1  1.0 -0.7
[3,]  0.8 -0.7  1.0
> chol(m)
以下にエラーchol(m) : 次数 3 の主対角行列が正定値ではありません
> eigen(m)
$values
[1]  2.0147034  1.0990993 -0.1138027 # 3つめの固有値がマイナス

$vectors
           [,1]        [,2]       [,3]
[1,] -0.5262037 -0.65755945  0.5391894
[2,]  0.4471271 -0.75329549 -0.4823104
[3,] -0.7233167  0.01270732 -0.6903995

> m <- matrix(-1, 1,1) # 次数1というと，1行1列の行列ですが。。。？
> chol(m)
以下にエラーchol(m) : 次数 1 の主対角行列が正定値ではありません
cols=rainbow(101)
x <- rnorm(1000)
y <- rnorm(1000)*200
z <- sample(101, 1000, replace=TRUE)
plot(x, y, col=cols[z], type="p", pch=20)
for (i in 1:101) rect(1.5, 502-2*i, 2.5, 500-2*i,col=cols[i], border=cols[i])
text(2.5, seq(502, 300, length=3), c("lo", "med", "hi"), pos=4)
gammaRange <- 10^(-5:5)
costRange <- 10^(-2:2)
t <- tune.svm(Species ~ ., data = iris, gamma=gammaRange,
cost=costRange,tunecontrol = tune.control(sampling="cross", cross=8))
以下にエラー sample(n) : 
  .Random.seed は整数ベクトルではなくタイプ 'list' で
Dataset <- 
read.table("C:/Users/Documents/ワークショップ/Rデータ処理用/（ここで強制改行）
修論素材/加工後/R用/111024/tsushin_growth.csv",
header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE)
[2] エラー: NA
OS Windows7(64bit)
R version 2.14.0
Rcmdr version 1.7-0
Total Accuracy:94
Single Accuracies:
90 96 96
model <- svm(judge~., data=dataset, corss=3)
summary(model)

Parameters:
   SVM-Type:  eps-regression 
 SVM-Kernel:  radial 
       cost:  1 
      gamma:  0.01190476 
    epsilon:  0.1 

Number of Support Vectors:  51

3-fold cross-validation on training data:

Total Mean Squared Error: 0.6049732 
Squared Correlation Coefficient: 0.06871624 
Mean Squared Errors:
 0.4767778 0.5592772 0.7788645 
> Df <- deriv(~x^2, "x", func=TRUE)
> Df(2)
[1] 4
attr(,"gradient")
     x
[1,] 4
> f <- expression(x^2)
> Df <- deriv(f, "x", func=TRUE)
> Df(2)
expression(x^2)
attr(,"gradient")
     x
[1,] 0
> f <- expression(x^2)
> Df <- deriv(f, "x", func=TRUE)
> Df(2)
[1] 4
attr(,"gradient")
     x
[1,] 4
test <- c(0.8699, 0.9029, 0.9133, 0.9014, 0.9160)
plot(test, type="l", lwd=2, xlab="test", ylab="AUC", xaxt="n")
mtext(1:4, 1, 1, at=1:4)
mtext("Q", 1, 1, at=5)
sessionInfo()
R version 2.10.1 (2009-12-14) 
i386-pc-mingw32 
plot(test, type="l", lwd=2, xlab="test", ylab="AUC", xaxt="n")
axis(side=1, at=1:5, label=c(1:4, 'Q'))
x <- c("マツホ−ム", "タケハウス", "ウメコウムテン",
		"カブ）マツホ−ム", "タケハウスカブシキガイシャ", "ウメコウムテンカブ")
> grep("_*ウメコウムテン_*", x, value = TRUE)
[1] "ウメコウムテン"     "ウメコウムテンカブ"
> grep("_*タケハウス_*", x, value = TRUE)
character(0)
> grep("_*マツホ−ム_*", x, value = TRUE)
[1] "マツホ−ム"         "ウメコウムテン"     "カブ）マツホ−ム"   "ウメコウムテンカブ"
> x <- c("マツホ−ム", "タケハウス", "ウメコウムテン",
+ 		"カブ）マツホ−ム", "タケハウスカブシキガイシャ", "ウメコウムテンカブ")
> grep("_*ウメコウムテン_*", x, value = TRUE)
[1] "ウメコウムテン"     "ウメコウムテンカブ"
> grep("_*タケハウス_*", x, value = TRUE)
[1] "タケハウス"                 "タケハウスカブシキガイシャ"
>  grep("_*マツホ−ム_*", x, value = TRUE)
[1] "マツホ−ム"       "カブ）マツホ−ム"
data <- read.csv("csvfile.csv")
data <- data[which(data[,1] == "hoge test"):nrow(data),]
x <- "csvfile.csv"
table <- read.table(x, sep="\n", header=F)
c <- grep("Well", table$V1) - 1
read.table(x, sep = ",", header = T, skip = c)
x <- "csvfile.csv"
table <- read.table(x, sep="\n", header=F, blank.lines.skip = F, fill = T)
c <- grep("Well", table$V1) - 1
read.table(x, sep = ",", header = T, skip = c, blank.lines.skip = F, fill = T)
x <- "test.csv"
table <- readLines(x)
c <- grep("Well", table) - 1
read.csv(x, skip = c, blank.lines.skip = FALSE)
A <- c(3, 7, 4, 6)
(mlist <- list(diag(1, A[1]), diag(1, A[2]), diag(1, A[3]), diag(1, A[4])))
lapply(A, diag)
A <- c(3, 7, 4, 6)
mlist <- lapply(A, function(N) matrix(1, N, N)) # このようにする
library(Matrix)
bdiag(mlist)
> group1
    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
[1,]   15   12   12   14   14   12   16   11   19    18    15    11    10     9     7    15
[2,]   14   16   12    9   14   16   14   11   11     6     7    11    16    12     9     9
[3,]   16   10   11   16   20   20   16   18   17    16    12     8    14    12    17    10
> group2
    [,1] [,2] [,3] [,4]
[1,]   18   19   12   16
> group3
    [,1] [,2] [,3] [,4]
[1,]   14   17   17   15
> group4
    [,1] [,2] [,3] [,4]
[1,]   13   19   11   17
> t.test(group1,group2,var.equal=F)
       Welch Two Sample t-test
data:  group1 and group2 
t = -1.9197, df = 3.667, p-value = 0.1337
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
-7.810966  1.560966 
sample estimates:
mean of x mean of y 
  13.125    16.250 
> t.test(group1,group3,var.equal=F)
       Welch Two Sample t-test
data:  group1 and group3 
t = -2.9049, df = 6.241, p-value = 0.02597
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:~
-4.8156168 -0.4343832 
sample estimates:
mean of x mean of y 
  13.125    15.750 
> t.test(group1,group4,var.equal=F)
       Welch Two Sample t-test
data:  group1 and group4 
t = -0.9899, df = 3.474, p-value = 0.3861
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:~
-7.463495  3.713495 
sample estimates:
mean of x mean of y 
  13.125    15.000 
> data

1     にんじん,たまねぎ
2     にんじん
3     にんじん,たまねぎ,じゃがいも
にんじん　3
たまねぎ   2
じゃがいも 1
> table(unlist(strsplit(as.character(data[,1]), ",")))
 じゃがいも   たまねぎ   にんじん
          1          2          3
ageage(x)<-function(x) {
m<-sprintf("title\t%s\n",x)
cat(m,file="result.txt",append=T)
}
title  50
title  54
title  76
> d <- data.frame(age=c(50, 54, 76), sex=c("male", "female", NA))
> deparse(substitute(d$age))
[1] "d$age"
> attach(d)
> deparse(substitute(age))
[1] "age"
> ageage <- function(x) {
+     m <- sprintf("title\t%s\n", x)
+     cat(m, file="result.txt", append=T)
+ }
> ageage(d$age)
> ageage(age)
func <- function(data.frame, vector)
{
	for (i in vector) {
		variable.name <- colnames(data.frame)[i]
		values <- data.frame[,i]
		cat("variable=", variable.name, "  n=", length(values),   
                   "  mean=", mean=mean(values, na.rm=TRUE), 
                   "  s d=", sd=sd(values, na.rm=TRUE), "\n")
	} 
}
> func(iris, 1:4)
variable= Sepal.Length   n= 150   mean= 5.843333   sd= 0.8280661 
variable= Sepal.Width   n= 150   mean= 3.057333   sd= 0.4358663 
variable= Petal.Length   n= 150   mean= 3.758   sd= 1.765298 
variable= Petal.Width   n= 150   mean= 1.199333   sd= 0.7622377 
Plate	Repeat	Well	Type	Time	Luciferase (CPS)	Time	 Luciferase (CPS)
1	1	A01	M	00:00:06.80	53530	00:02:38.05	47900
1	1	A02	M	00:00:07.07	107230	00:02:38.31	94150
1	1	A03	M	00:00:07.34	103510	00:02:38.58	95070
> library(RODBC)
> sheet <- odbcConnect("Excel Files")
> tab <- sqlQuery(sheet, "select * from [List ; Plates 1 - 1$]")
> t <- as.numeric(difftime(tab[,7], tab[,5], units="secs"))
> t[1]
[1] 152
> (x <- difftime("2011-10-11 00:02:38.05 JST", "2011-10-11 00:00:06.80 JST", units="secs"))
Time difference of 151.25 secs
> as.numeric(x)
[1] 151.25
x <- c(rep(0, 30000), 1, rep(0, 19999))
m <- matrix(x, nrow=50000, ncol=2)
par(mfrow=c(2,1)) # plot関数との表示比較のため
image(m, col=c("white", "black"))
plot(1:nrow(m), m[,1], type="l")
x <- c(rep(0, 300), 1, rep(0, 199))
m <- matrix(x, nrow=500, ncol=2)
par(mfrow=c(2, 1))
image(m, col=c("white", "black"))
plot(1:nrow(m), m[,1], type="l")
plot(1:nrow(m), m[,1], type="l") # あなたのプログラムだとスパイク状に描画される
plot(1:nrow(m), m[,1], type="n")
segments(1:nrow(m), 0, 1:nrow(m), m[,1]) # 垂直線を描く
st <- c("2008/12/19", "2009/3/24")
st <- strptime(st, "%Y/%m/%d", tz="")
y.days <- seq(st[1], st[2], by="1 day")
x.hours <- seq(0, 48, length=3600*48/10) 

m <- matrix(sample(c(0, 1), (3600*48/10)*length(y.days), rep=T, prob=c(10, 1)),
            17280, length(y.days))

par(oma=c(2, 1.5, 1, 1))

y.lim <- c("2008/12/17", "2009/3/26")
y.lim <- strptime(y.lim, "%Y/%m/%d", tz="")
image(x.hours, y.days, m, col=c("white", "black"), yaxt="n", xaxp=c(0, 48, 4),
      xlab="hours", ylab=NA,
      xlim=c(-0.1, 48), ylim=c(as.numeric(y.lim[2]), as.numeric(y.lim[1])))

y.gr <- c("2008/12/18 12", "2009/3/24 12")
y.gr <- strptime(y.gr, "%Y/%m/%d %H", tz="")
y.gr <- seq(y.gr[1], y.gr[2], by="1 days")
abline(h=y.gr, col="gray")

axis(2, at=as.numeric(y.days[(1:4)*20]), labels=as.character(y.days[(1:4)*20],
     format="%Y/%m/%d"), las=2)
png(, length(x1), length(y1))
par(mar = c(0, 0, 0, 0))
image(x.hours, y.days, m, col = c("white", "black"))
dev.off()
x <- x.hours[matrix(1:length(x.hours), length(x.hours), length(y.days))[m == 1]]
y <- y.days[matrix(1:length(y.days), length(x.hours), length(y.days), byrow = TRUE)[m == 1]]
plot(x, y, pch = 20, cex = 1) # 真っ黒になるときはcex = 1の値を小さくしてみる。
for(i in 1:ncol(m)){
segments(x.hours[m[,i]==1], y.days[i], x.hours[m[,i]==1], y.days[i+1]) 
}
st <- c("2008/12/19", "2009/3/24")
st <- strptime(st, "%Y/%m/%d", tz="")
y.days <- seq(st[1], st[2], by="1 day")
x.hours <- seq(0, 48, length=3600*48/10) 

m <- matrix(sample(c(0, 1), (3600*48/10)*length(y.days), rep=T, prob=c(10, 1)),
            17280, length(y.days))

par(oma=c(2, 1.5, 1, 1))

y.lim <- c("2008/12/17", "2009/3/26")
y.lim <- strptime(y.lim, "%Y/%m/%d", tz="")
plot(c(-0.1, 48), c(as.numeric(y.lim[2]), as.numeric(y.lim[1])), type="n",
     yaxt="n", ylab="")
y.days2 <- rep(y.days, each=17280)[m==1]
x.hours2 <- rep(x.hours, 96)[m==1]
mk <- diff(y.days)[1]
segments(x.hours2, y.days2, x.hours2, y.days2+mk)
y.gr <- c("2008/12/18 12", "2009/3/24 12")
y.gr <- strptime(y.gr, "%Y/%m/%d %H", tz="")
y.gr <- seq(y.gr[1], y.gr[2], by="1 days")
abline(h=y.gr+mk/2, col="gray")

axis(2, at=as.numeric(y.days[(1:4)*20]), labels=as.character(y.days[(1:4)*20],
     format="%Y/%m/%d"), las=2)
a <- c(0, 0, 1, 1, 0, 0, 0, 1, 0, 0) 
a2 <- c(0, 1, 1, 0, 1, 0)
a2 <- 0
for (i in 1:9) {
   if (a[i]==0 && a[i+1]==0) next
   a2 <- c(a2, a[i+1])
}
> a3 <- c(FALSE,TRUE,TRUE,TRUE,FALSE,FALSE,TRUE,TRUE,FALSE,TRUE)
> a2 <- a[a3]
> a2
[1] 0 1 1 0 1 0
b <- rle(a)
b$lengths[which(b$values==0)] = 1
a2 <- inverse.rle(b)
set.seed(10); x <- rnorm(100)
hi <- hist(x, plot=F)
plot(0, 0, xlim=range(hi$counts), ylim=range(hi$breaks), type='n',
     xlab='', ylab='', las=1, bty='L')
for(i in seq(length(hi$counts))) {
  rect(0, hi$breaks[i], hi$counts[i], hi$breaks[i+1], col=grey(0.8),  
       border=grey(0.7))
}
box(bty='L')
set.seed(10); x <- rnorm(100)
h <- hist(x, plot = FALSE)
barplot(h$counts, space = 0, horiz = TRUE)
plot.window(xlim = range(h$counts), ylim = range(h$breaks))
axis(2)
> data
       data1 data2 data3 data4  data5
case1     NA    10   100  1000  10000
case2      2    NA   200    NA  20000
case3      3    30   300  3000  30000
case4     NA    40    NA  4000     NA
case5      5    NA   500  5000  50000
case6      6    60   600    NA     NA
case7      7    70   700  7000  70000
case8     NA    80   800  8000     NA
case9      9    90    NA  9000  90000
case10    10    NA  1000 10000 100000
> similarity
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,] 0.00 0.16 0.46 0.50 0.77 0.84 1.38 1.41 1.61  1.73
 [2,] 0.16 0.00 0.20 0.00 0.60 0.67 1.00 0.67 1.17  1.61
 [3,] 0.46 0.20 0.00 0.17 0.46 0.63 1.05 1.00 1.42  1.61
 [4,] 0.50 0.00 0.17 0.00 0.11 0.25 0.50 0.67 0.84  0.67
 [5,] 0.77 0.60 0.46 0.11 0.00 0.17 0.46 0.47 0.80  1.15
 [6,] 0.84 0.67 0.63 0.25 0.17 0.00 0.21 0.33 0.53  0.67
 [7,] 1.38 1.00 1.05 0.50 0.46 0.21 0.00 0.20 0.47  0.69
 [8,] 1.41 0.67 1.00 0.67 0.47 0.33 0.20 0.00 0.17  0.31
 [9,] 1.61 1.17 1.42 0.84 0.80 0.53 0.47 0.17 0.00  0.20
[10,] 1.73 1.61 1.61 0.67 1.15 0.67 0.69 0.31 0.20  0.00
n <- nrow(data)   # 列数
m <- ncol(data)   # 行数
suf <- integer(3) # 類似度が高いもの3つの添え字を入れる
for (i in 1:n) {  # 各 case について
  for (j in 1:m) { # 各 data について
    if (is.na(data[i, j])) { # もし，欠損値なら，類似度によりランク付け
      # ties.method は average 以外ならどれでもよい?
      hi <- rank(similarity[1,], ties.method="first")
      # 類似度が高く，かつ，データが欠損値でない case の添え字を
      # suf に格納（3つまで）
      l <- 0           
      for (k in n:1) {
        if (!is.na(data[hi[k], j])) { # 欠損値でないなら
          l <- l+1              # カウンタを増やす
          suf[l] <- k           # 添え字を格納
          if (l == 3) break     # 3つ見つかったら脱出
        }
      }
      if (l == 3) {    # 3つないときには NA のままにする?
                       # （平均値でよいと思うけど）
         data[i, j] <- mean(data[hi[suf], j])
      }
    }
  }
}
           data1 data2     data3 data4     data5
case1   8.666667    10  100.0000  1000  10000.00
case2   2.000000    80  200.0000  9000  20000.00
case3   3.000000    30  300.0000  3000  30000.00
case4   8.666667    40  833.3333  4000  86666.67
case5   5.000000    80  500.0000  5000  50000.00
case6   6.000000    60  600.0000  9000  86666.67
case7   7.000000    70  700.0000  7000  70000.00
case8   8.666667    80  800.0000  8000  86666.67
case9   9.000000    90  833.3333  9000  90000.00
case10 10.000000    80 1000.0000 10000 100000.00
hi <- rank(similarity[1,], ties.method="first") # 変更前
hi <- rank(similarity[i,], ties.method="first") # 変更後
hi  <- rank(similarity[i,], ties.method="first")
hi2 <- hi
for (a in 1:max(hi)) {
	for (b in 1:max(hi)) {
		if (hi[b] == max(hi) + 1 - a) {
			hi2[a] <- b
			break
		}
	}
}
ifelse(is.na(data), t(apply(similarity, 2,
       function(i) apply(data[order(i, decreasing = TRUE), ], 2,
       function(x) mean(x[head(which(!is.na(x)), 3)])))), data)
          data1    data2     data3     data4     data5
case1   8.666667 10.00000  100.0000  1000.000  10000.00
case2   2.000000 73.33333  200.0000  8666.667  20000.00
case3   3.000000 30.00000  300.0000  3000.000  30000.00
case4   8.666667 40.00000  633.3333  4000.000  66666.67
case5   5.000000 60.00000  500.0000  5000.000  50000.00
case6   6.000000 60.00000  600.0000  4666.667  43333.33
case7   7.000000 70.00000  700.0000  7000.000  70000.00
case8   3.333333 80.00000  800.0000  8000.000  20000.00
case9   9.000000 90.00000  200.0000  9000.000  90000.00
case10 10.000000 36.66667 1000.0000 10000.000 100000.00
while (1) {
    ...
    if (1) break
    ...
}
...　#ここへジャンプ
> strsplit("abcdedcba","")[[1]]
[1] "a" "b" "c" "d" "e" "d" "c" "b" "a"
> which(strsplit("abcdedcba","")[[1]]=="d")
[1] 4 6
> AB <- c(names(A), names(B))   # 名前を取り出して一つのベクトルにします。
> CD <- AB[duplicated(AB)]   # 重複している名前を一つにする。
ここで以下のような行列xがあります。
V1には名前が入っていて、ベクトルAとベクトルBの名前は
この中に全て含まれています。
> x
      V1   V2
1  名前1   値
2  名前2   値
3  名前3   値
:
> x[x$V1 == CD, ] # CDの要素と行列xの名前とが
                    一致する行を取り出します。
> CD <- c("foo", "bar", "baz")
> (x <- data.frame(V1=c("abc", "foo", "zoo", "boo",
                        "baz", "yah", "bar"), V2=1:7))
   V1 V2
1 abc  1
2 foo  2
3 zoo  3
4 boo  4
5 baz  5
6 yah  6
7 bar  7
> x[x$V1 %in% CD,]
   V1 V2
2 foo  2
5 baz  5
7 bar  7
a <- c("水", "金", "地", "火", "木", "土", "天", "海")
b <- c("水", "地", "土")
a[!a %in% b]
spin <- function(x0, y0, factor=1, rotation=2)
{
	r <- seq(0.1, 1.0, length=100)*factor
	theta <- seq(0, -rotation*2*pi, length=100)
	x <- r*cos(theta)+x0
	y <- r*sin(theta)+y0
	lines(x, y)
}
plot(0:10, 0:10, type="n", axes=FALSE, xlab="", ylab="")
spin(3, 5, 0.5)
spin(5, 5, 1)
spin(7, 5, 0.2)
spin(3, 7, 0.5, rotation=5)
spin(5, 7, 1, rotation=4)
spin(7, 7, 0.2, rotation=1)
X	Y
4	1.1
4.1	0.6
4.2	0
4.3	0.9
4.4	0.9
4.5	0
4.6	0
4.7	0
4.8	0.7
4.9	0.9
5	0
5.1	1
5.2	1.1
5.3	0
5.4	0.4
5.5	1.5
5.6	0.9
5.7	1.5
5.8	2.1
5.9	2.6
6	3.2
6.1	4.3
6.2	6.1
6.3	8.3
6.4	10
6.5	11.9
6.6	15.3

data <- read.table("dataXY.txt", header=T)
data1 <- data[1:16,]
c <- mean(data1$Y)
c
data2 <- data[20:27,]
L <- lm(Y~X, data= data2)
L
plot(data)
abline(h=c, col=2)
abline(L, col=4)
(x <- (c-coefficients(L)[1]) / coefficients(L)[2])
points(x, c, pch=19)
text(x, c, sprintf("(%.3f, %.3f)", x, c), pos=4, xpd=TRUE)
library(maptools)
map <- readShapeLines("tokyok.shp")

lon1 <- 138   ### 経度138E-141E
lon2 <- 141
lat1 <- 35    ### 緯度35N-35.5N   
lat2 <- 35.5

par(usr=c(lon1, lon2, lat1, lat2))
plot(map, xlim=c(lon1, lon2), ylim=c(lat1, lat2))
> par("usr")
[1] 138.00000 141.00000  33.90625   36.59375
> sessionInfo()
R version 2.13.1 Patched (2011-08-31 r56868)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] ja_JP.UTF-8/ja_JP.UTF-8/C/C/ja_JP.UTF-8/ja_JP.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets
    methods   base     

other attached packages:
[1] maptools_0.8-9  lattice_0.19-33 sp_0.9-83
    foreign_0.8-45 

loaded via a namespace (and not attached):
[1] grid_2.13.1  tools_2.13.1
> par("usr")
[1] 137.88000 141.12000  32.23003  38.26997
library("maptools")
map <- readShapeLines("tokyok.shp")
lon1 <- 138   ### 経度138E-141E
lon2 <- 141
lat1 <- 35    ### 緯度35N-35.5N   
lat2 <- 36 # 35.5だと何もプロットされないので変更してます。
plot(map, xlim = c(lon1, lon2), ylim = c(lat1, lat2), setParUsrBB = TRUE)
pdf("test.pdf")
plot(1:10)
hist(1:10)
dev.off()
x <- xmean + σ # 親集合の重心周辺に点を取っているのだと思う(σは変数)
x <- 10 * ( x / sum(x) )      #xの合計が5なら5で割ってその後10倍
...##色々計算　よくなっていたら更新
...##くり返し
> n <- 30
> x <- seq(-n, n)
> y <- 0.04*x^3+0.03*x^2+0.02*x+rnorm(2*n+1, 0, 1) # 真の関係
> plot(x, y)
> 
> res <- lm(y ~ x) # 直線回帰
> func <- function(x){
+ coef(res)[1] + coef(res)[2]*x # これが直線回帰の回帰式
+ }
> points(x, func(x=x))
> 
> res2 <- gam(y~s(x, k=3)) # 加法モデル（
                           # デフォルトは薄板スプライン）による回帰
> points(x, predict(res2), type="l", col=2, lwd=3)
                           # predictによる予測で何故か直線回帰と一致
> coef(res2) #回帰係数が何か出力されてはいるがヘルプを見ても分からない
  (Intercept)        s(x).1        s(x).2 
 9.309025e+00 -1.199402e-09  3.932092e+02 
> 
> res2 <- gam(y~s(x, k=4)) #節点の数を増やすと曲線となる
> points(x, predict(res2), type="l", col=2, lwd=3)
> coef(res2) 
(Intercept)      s(x).1      s(x).2      s(x).3 
   9.309025 -843.683785 -129.386077 1162.199292 
> set.seed(1)
> x <- seq(0, 100)
> y <- 0.02*x+rnorm(101, 0, 1)
> plot(x, y)
> 
> res <- lm(y ~ x)
> func <- function(x){
+ coef(res)[1] + coef(res)[2]*x 
+ }
> points(x, func(x=x), type="l", lwd=3, col=2)
> x2 <- seq(0, 80)
> points(x2, func(x=x2), type="l", lwd=3, col=3)
lines(x.range, coefficients(a)[1]+x.range*coefficients(a)[2], col=2)
> a <- matrix(sample(c(1:9,NA),1000,replace=TRUE),100)
> sum(is.na(a[,3]))/nrow(a)
[1] 0.12
> b <- list(a1=a,a2=a)
> is.na(b)
   a1    a2 
FALSE FALSE
curve(dnorm(x, sd=sqrt(1/10)), -3, 3)
curve(dnorm(x, sd=sqrt(1/5)), -3, 3, add=TRUE)
curve(dnorm(x, sd=sqrt(1/3)), -3, 3, add=TRUE)
curve(dnorm(x, sd=sqrt(1/1)), -3, 3, add=TRUE)
> d <- table(x)
        b
a         1  2  3
  female 13  2  6
  male    6 14  9
> data.frame(d)
       a b Freq
1 female 1   13
2   male 1    6
3 female 2    2
4   male 2   14
5 female 3    6
6   male 3    9
st <- c("2008/12/19", "2009/3/24")
st2 <- strptime(st, "%Y/%m/%d", tz="")
y.days <- seq(st2[1],st2[2], by="1 day")
x.hours <- seq(0, 24, by=1)

m <- matrix(sample(c(0, 1), 24*length(y.days), rep=T),
            nrow=24, ncol=length(y.days))

image(x.hours, y.days, m, col=c("white", "black"), yaxt="n",
      xaxp=c(0,24,4), xlab="hours", ylab=NA)

r <- as.POSIXct(round(range(y.days), "days"))
axis.POSIXct(2, at=seq(r[1],r[2], by="20 day"), 
             format="%Y/%m/%d", las=2)
image(...,  ylim=c(as.numeric(y.days[96]), as.numeric(y.days[1])))
r <- as.POSIXct(round(range(y.days), "days"))
axis.POSIXct(2, at=seq(r[1],r[2], by="20 day"), 
             format="%Y/%m/%d", las=2)
image(x.hours, y.days, m, col=c("white", "black"),
      yaxt="n", xaxp=c(0,24,4),xlab="hours", ylab=NA,
      ylim=c(as.numeric(y.days[96]), as.numeric(y.days[1])))
axis(2,at=as.numeric(y.days[(1:4)*20]),
     labels=as.character(y.days[(1:4)*20]), las=2)
>x
	1	2	3	4	.	.
A	20	50	60	29    .       .
B	49	74	90	38    .       .
C	57	49	20	39    .       .
>Y
[1] "C"	"A"	.....
result <- NULL
for (i in 1:length(Y)) {
    result <- rbind(result, x[rownames(x) == Y[i], ])
}
t(sapply(Y, function(a) x[a,]))
result <- matrix(0, length(Y), 2000)
for (i in 1:length(Y)) {
    result[i,] <- x[rownames(x) == Y[i], ]
}
#簡単な例
> set.seed(1234)
> n <- 6; m <- 4
> x <- matrix(rnorm(n^2),nrow=n,ncol=n)
> colnames(x) <- seq(n)
> rownames(x) <- sample(LETTERS,n,rep=TRUE)

> B <- rownames(x); uB <- unique(B)
> A <- sample(uB,m,rep=TRUE)

> z <- factor(B,levels=uB)
> zx <- cbind(z,x) # 行ラベルを数値（因子）化し行列の先頭列に加える

> Ax <- zx[factor(A,levels=uB),][,seq(2,n+1)]

> x
           1          2           3          4          5          6
A -1.2070657 -0.5747400 -0.77625389 -0.8371717 -0.6937202  1.1022975
U  0.2774292 -0.5466319  0.06445882  2.4158352 -1.4482049 -0.4755931
C  1.0844412 -0.5644520  0.95949406  0.1340882  0.5747557 -0.7094400
N -2.3456977 -0.8900378 -0.11028549 -0.4906859 -1.0236557 -0.5012581
J  0.4291247 -0.4771927 -0.51100951 -0.4405479 -0.0151383 -1.6290935
B  0.5060559 -0.9983864 -0.91119542  0.4595894 -0.9359486 -1.1676193
> B
[1] "A" "U" "C" "N" "J" "B"
> A
[1] "U" "J" "B" "C"
> z
[1] A U C N J B
Levels: A U C N J B
> zx
  z          1          2           3          4          5          6
A 1 -1.2070657 -0.5747400 -0.77625389 -0.8371717 -0.6937202  1.1022975
U 2  0.2774292 -0.5466319  0.06445882  2.4158352 -1.4482049 -0.4755931
C 3  1.0844412 -0.5644520  0.95949406  0.1340882  0.5747557 -0.7094400
N 4 -2.3456977 -0.8900378 -0.11028549 -0.4906859 -1.0236557 -0.5012581
J 5  0.4291247 -0.4771927 -0.51100951 -0.4405479 -0.0151383 -1.6290935
B 6  0.5060559 -0.9983864 -0.91119542  0.4595894 -0.9359486 -1.1676193
> Ax
          1          2           3          4          5          6
U 0.2774292 -0.5466319  0.06445882  2.4158352 -1.4482049 -0.4755931
J 0.4291247 -0.4771927 -0.51100951 -0.4405479 -0.0151383 -1.6290935
B 0.5060559 -0.9983864 -0.91119542  0.4595894 -0.9359486 -1.1676193
C 1.0844412 -0.5644520  0.95949406  0.1340882  0.5747557 -0.7094400 

# n=10000, m=3500の場合の実行速度
> system.time(Ax <- zx[factor(A,levels=uB),][,seq(2,n+1)])
   ユーザ   システム       経過  
     1.804      0.156      1.965
> x <- matrix(round(runif(12)*100),ncol=4)
> colnames(x) <- as.character(1:4)
> rownames(x) <- LETTERS[1:3]
> x
   1  2  3  4
A 93 68 20  5
B 26 20 17 74
C 86 35 59 64
> Y <- sample(LETTERS[1:3],10000,replace=TRUE)
> head(Y)
[1] "C" "C" "A" "C" "C" "B"
> system.time(x[Y,])
   ユーザ   システム       経過  
     0.000      0.000      0.001 
> system.time({result <- matrix(0, length(Y), 4);
+             for (i in 1:length(Y)) {
+               result[i,] <- x[rownames(x) == Y[i], ]
+             }})
   ユーザ   システム       経過  
     0.080      0.000      0.087
>  system.time(Ax <- x[A,])
   ユーザ   システム       経過  
     1.048      0.080      1.132
> chisq.test(c(1,2,1,7))
Chi-squared test for given probabilities
data:  c(1, 2, 1, 7)
X-squared = 9, df = 3, p-value = 0.02929
警告メッセージ： 
In chisq.test(c(1, 2, 1, 7)) :  カイ自乗近似は不正確かもしれません
     ID class Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q11 Q15
1     1     a  5  4  4  4  4  2  4   4   5
2     2     a  4  4  3  3  4  4  3   4   2
3     3     a  4  4  3  5  4  4  3   4   5
4     4     a  4  4  1  3  3  4  5   3   5
5     5     a  5  4  1  5  5  3  2   4   4
6     6     a  4  4  4  4  4  2  1   2   1
         　　　-- 省略--
94   94     d  4  4  3  4  4  2  3   4   4
95   95     d  4  3  4  3  3  4  4   4   4
96   96     d  4  4  1  4  3  3  2   4   3
97   97     d  4  5  2  5  5  5  4   4   4
98   98     d  5  4  4  4  3  4  5   4   4
99   99     d  4  4  2  4  4  3  3   4   4
100 100     d  4  4  2  3  4  4  5   4   4
101 101     d  5  4  2  4  4  5  4   4   1
102 102     d  5  4  4  4  3  5  2   4   2
$Q1
   x
     3  4  5
  a  2 15  8
  b  0 21  6
  c  2 19  4
  d  0 14 11

$Q2
   x
     2  3  4  5
  a  1  5 19  0
  b  1  6 20  0
  c  0  8 14  3
  d  0  4 17  4
 -- 以下省略--
sapply(d01[,3:11],function(x) { tbl <- table(d01$class,x); print(round(tbl/rowSums(tbl)*100,1)) }))
  x
    3  4  5
 a  2 15  8
 b  0 19  3

  x
    2  3  4
 a  1  5 19
 b  1  5 16
sapply(split(d, d[,2]),
       function(x) colSums(x[3:9])/nrow(x)*100)
sapply(split(d, d[,2]), function(x) colMeans(x[3:9])*100)
sapply(d[3:9], function(x) tapply(x, d$class, mean)*100)
by(d[3:9], d[,2], colMeans)
(Intercept)      s(x).1      s(x).2      s(x).3 
   9.309025 -843.683785 -129.386077 1162.199292 
  ID class R Wr L Spk Vocab Gr Spd
   1     a 0  0 0   0     0  0   0
   2     a 0  0 0   0     0  0   0
   3     a 0  0 0   0     0  0   1
   4     a 0  0 0   0     0  0   0
   5     a 1  0 0   0     0  1   1
      :  省略
 100     d 1  0 0   0     1  0   0
 101     d 1  0 0   0     1  0   0
 102     d 0  0 0   0     0  0   0
 103     d 1  0 0   0     1  0   0
sapply(split(d, d[,2]), function(x) colSums(x[3:9]))
by(d[3:9], d[,2], colSums)
sapply(d[3:9], function(x) tapply(x, d$class, sum))
set.seed(123456789)
n <- 500
y <- rnorm(n)
x <- sample(10000, n, replace=TRUE)
for (i in 1:n) {
	xi <- x[i]
	if (xi < 4500) x[i] <- xi/4500
	else if (xi < 5500) x[i] <- (xi-4500)/500+1
	else x[i] <- (xi-5500)/4500+3
}
x <- x*2500
plot(x, y, xaxt="n")
axis(1, at=0:4*2500,
     labels=c(0, 4500, 5000, 5500, 10000))
ex)
     0002jh 0043sss
  1         2           4
    X0002jh X0043sss
   1       2        4
> f <- function(x) exp(1i*x)
> integrate(function(x) Re(f(x)),0,1)$value+
            (1i)*integrate(function(x) Im(f(x)),0,1)$value 
[1] 0.841471+0.4596977i
要素1,要素2,要素3
りんご,バナナ,みかん
バナナ,なし,
パイン,もも,みかん
data1 <- read.table("C:/lab/kankou.txt",sep=",",
                    header=TRUE)
　要素1  要素2  要素3
1 りんご バナナ みかん
2 バナナ なし
3 パイン もも   みかん
data.tran <- as(data1, "transactions")
data.frame <- as(data.tran, "data.frame")
1{要素1=りんご,要素2=バナナ,要素3=みかん}
2{要素1=バナナ,要素2=なし,要素3=}
3{要素1=パイン,要素2=もも,要素3=みかん}
1{りんご,バナナ,みかん}
2{バナナ,なし,}
3{パイン,もも,みかん}
data.frame[,2] <- gsub("要素[1-3]=", "", data.frame[,2])
1{りんご,バナナ,みかん}
2{バナナ,なし,}
3{パイン,もも,みかん}
> library(plm)
 以下にエラー library(plm) :
    'plm' という名前のパッケージはありません 
> install.packages(plm)
 以下にエラー install.packages(plm) :
    オブジェクト 'plm' がありません 
> utils:::menuInstallLocal()
 パッケージ 'plm' は無事に開封され、
 MD5 サムもチェックされました 
> library(plm)
 要求されたパッケージ bdsmatrix をロード中です 
 エラー：  パッケージ 'bdsmatrix' を
           ロードできませんでした 
 追加情報:   警告メッセージ： 
1:  パッケージ 'plm' はバージョン 2.13.1 の R の下で
    造られました  
2: In library(pkg, character.only = TRUE, 
              logical.return = TRUE, lib.loc = lib.loc) :
   'bdsmatrix' という名前のパッケージはありません 
utils:::menuInstallLocal()
 以下にエラー gzfile(file, "r") :
    コネクションを開くことができません 
 追加情報:   警告メッセージ： 
In gzfile(file, "r") :
   圧縮されたファイル 'bdsmatrix_1.0[1]/DE木SCRIPTION'
     を開くことができません,
   理由は 'No such file or directory' です 
> utils:::menuInstallLocal()
 以下にエラー gzfile(file, "r") :
    コネクションを開くことができません 
 追加情報:   警告メッセージ： 
In gzfile(file, "r") :
   圧縮されたファイル 'bdsmatrix_1.0[1]/DESCRIPTION' を
   開くことができません,
   理由は 'No such file or directory' です 
以下にエラー as.vector(x, mode) : 
 cannot coerce type 'closure' to vector of type 'any'
     ID class Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q11 Q15
1     1     a  5  4  4  4  4  2  4   4   5
2     2     a  4  4  3  3  4  4  3   4   2
3     3     a  4  4  3  5  4  4  3   4   5
4     4     a  4  4  1  3  3  4  5   3   5
5     5     a  5  4  1  5  5  3  2   4   4
6     6     a  4  4  4  4  4  2  1   2   1
7     7     a  5  4  2  2  4  3  4   4   2
8     8     a  4  4  2  5  4  3  3   4   3
9     9     a  4  4  2  3  3  4  3   5   5
10   10     a  5  4  4  4  5  3  3   4   4
11   11     a  4  3  2  3  2  3  1   3   3
12   12     a  3  2  3  4  3  5  4   3   2
13   13     a  5  4  2  4  4  4  5   4   4
14   14     a  3  3  4  3  4  3  1   3   2
15   15     a  4  4  3  4  4  4  4   4   3
16   16     a  5  3  2  4  4  3  3   5   5
17   17     a  4  4  2  4  3  3  4   3   3
18   18     a  4  4  2  2  4  3  2   3   4
19   19     a  5  4  4  4  4  2  5   4   5
20   20     a  4  4  3  4  4  5  4   4   5
21   21     a  4  4  1  4  3  4  4   4   2
22   22     a  4  3  2  3  3  5  3   3   3
23   23     a  4  3  3  3  4  3  3   4   3
24   24     a  5  4  1  4  4  4  4   4   2
25   25     a  4  4  3  4  4  5  4   3   3
26   26     b  4  4  1  4  4  5  3   5   5
27   27     b  5  4  3  4  3  3  3   4   3
28   28     b  4  3  2  4  3  5  4   3   2
29   29     b  4  4  2  3  3  4  4   4   2
30   30     b  4  4  2  3  4  4  3   4   3
31   31     b  4  4  4  3  4  4  2   4   4
32   32     b  4  4  4  4  3  4  3   4   3
33   33     b  4  4  4  3  3  4  4   3   3
34   34     b  4  4  5  3  3  4  3   3   3
35   35     b  5  4  3  4  4  3  3   4   4
36   36     b  4  3  1  3  3  4  5   3   1
37   37     b  4  4  5  4  3  4  4   4   2
38   38     b  4  4  5  2  2  4  2   3   3
39   39     b  4  4  4  4  3  4  2   4   3
40   40     b  4  4  3  4  4  3  3   4   4
41   41     b  4  3  2  3  3  4  4   4   2
42   42     b  5  3  3  4  3  4  5   3   2
43   43     b  4  3  2  3  4  5  5   4   5
44   44     b  4  2  3  3  3  4  3   3   3
45   45     b  4  4  3  5  3  3  5   4   1
46   46     b  4  4  3  3  5  4  5   4   4
47   47     b  4  4  3  4  3  4  4   4   2
> dat01 <- data.frame(ID=1:47,class=c(rep('a',25),rep('b',22)),
                             Q1=sample.int(5,size=47,replace=TRUE),
                             Q2=sample.int(5,size=47,replace=TRUE))
> apply(dat01[,3:4],2,function(x){tapply(x,dat01$class,table)})                  
$Q1
$Q1$a

 1  2  3  4  5 
 3 10  1  5  6 

$Q1$b

1 2 3 4 5 
4 7 3 4 4 

$Q2
$Q2$a

1 2 3 4 5 
6 4 2 7 6 

$Q2$b 

1 2 3 4 5 
2 5 4 4 7 
> d <- read.table("data01.dat", header=TRUE)
> sapply(d[3:11], function (x) table(d$class, x))
$Q1
   x
     3  4  5
  a  2 15  8
  b  0 19  3

$Q2
   x
     2  3  4
  a  1  5 19
  b  1  5 16
   ： 以下略
以下にエラー grep("リー", "あ", "い", "う", "え") : 
  "リー" は不正な正則表現です, 理由は 'Missing ']''
> grep("リー", c("スリー","い","リード","え"))
[1] 1 3
#サンプル
x <- 2
x

Ｒ　Version:  2.13.1
ＲコマンダーVersion:  1.6-4
In getDependencies(pkgs, dependencies, available, lib) :
  package ‘adapt’ is not available (for R version 2.13.1
> sessionInfo()
R version 2.13.1 (2011-07-08)
Platform: i386-pc-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=Japanese_Japan.932
    LC_CTYPE=Japanese_Japan.932   
[3] LC_MONETARY=Japanese_Japan.932 LC_NUMERIC=C                  
[5] LC_TIME=Japanese_Japan.932    

attached base packages:
[1] stats     graphics  grDevices utils     
    datasets  methods   base     

other attached packages:
[1] rgl_0.92.798

loaded via a namespace (and not attached):
[1] tools_2.13.1
x <- 0:4*0.25 # 等間隔であるとする
y <- c(0, 4, 10, 7, 2)
plot(x, y, type="l", yaxs="i")
polygon(c(x, rev(x)), c(y, rep(0, 5)),
        col="#00FF0040")
abline(h=y, lty=3)
segments(x, y, x, rep(0, 5))
(y[1]+y[2])*0.25/2 + (y[2]+y[3])*0.25/2 + (y[3]+y[4])*0.25/2 + (y[4]+y[5])*0.25/2
(2*sum(y)-y[1]-y[5])*0.25/2
> delta.x <- x[2]-x[1] # x 座標の間隔
> sum(y)*delta.x
[1] 3.803799
> system.time(for (i in 1:10000) sum(y)*delta.x)
   ユーザ   システム       経過  
     0.068      0.001      0.069 
> polygon <- cbind(x, y) # 開曲線
> areapl(polygon)
[1] 3.803799
> system.time(for (i in 1:10000) areapl(polygon))
   ユーザ   システム       経過  
     4.943      0.211      4.643 
> polygon2 <- cbind(c(x, 0), c(y, 0)) # 閉曲線
> areapl(polygon2)
[1] 3.803799
> system.time(for (i in 1:10000) areapl(polygon2))
   ユーザ   システム       経過  
     2.586      0.123      2.337
> chisq.test(c(4, 5, 6, 17))

	Chi-squared test for given probabilities

data:  c(4, 5, 6, 17) 
X-squared = 13.75, df = 3, p-value = 0.003266
> library(Matching)
> data(lalonde)
> Y <- lalonde$re78
> Tr <- lalonde$treat
> glm1 <- glm(Tr~age+educ+black+hisp+married+nodegr+
              re74+re75, family=binomial, data=lalonde)
> rr1 <- Match(Y = Y, Tr = Tr, X = glm1$fitted)
> summary(rr1)

Estimate...  2624.3 
AI SE......  802.19 
T-stat.....  3.2714 
p.val......  0.0010702 

Original number of observations..............  445
Original number of treated obs...............  185 
Matched number of observations...............  185 
Matched number of observations  (unweighted).  344 
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    0    0    0    0    0    0    0    0     0
 [2,]    0    0    0    0    0    0    0    0    0     0
 [3,]    0    0    0    1    0    0    0    0    0     0
 [4,]    0    0    0  999    0    0    0    0    0     0
 [5,]    0    0    0  999    0    0    0    0    0     0
 [6,]    0    1    0  999    0    0    0    0    0     0
 [7,]    0    1    0  999    0    0    0    0    0     0
 [8,]    0    1    0  999    0    0    0    0    0     0
 [9,]    0    1    0  999    0    0    0    0    0     0
[10,]    0  999    0  999    0    0    0    0    0     0
[11,]    0  999    0  999    0    0    0    0    0     0
[12,]    1  999    0  999    0    0    0    0    0     0
[13,]    1  999    0  999    0    0    0    0    0     0
x[4:13, 4] <- rep(2:4, c(3, 5, 2)) # x[4:13, 4] <- c(rep(2, 3), rep(3, 5), rep(4, 2)) と同じ
# 乱数発生数
k <- 100
# 解析年数
n <- 50
###### データ入力(塩化物イオン拡散予測) ########
C0 <- rweibull (k,   2.600, 9.631)
D  <- rlnorm   (k, -17.454, 0.995)
d  <- 6.85
###### データ入力(腐食環境) ####################
a   <- 1500
c   <- 0.005
alf <- 220
fai <- 1.3
Fe  <- 7850
################################################
# 誤差関数の設定
erf <- function(x) 2 * pnorm(x * sqrt(2)) - 1
# 塩化物配列の設定
C <- matrix(0, n, k)
###### 塩化物拡散予測 ##########################
for (i in 1:n) {      
  for (j in 1:k) {  
     C[i,j] <- C0[j]*(1-erf(d/(2*sqrt(D[j]*
                      i*31536000))))
  } 
}
###### 鉄筋腐食予測（年数予測) #################
t1 <- round((a/alf)*log(0.01/((4*a*c)/(fai*Fe))))
t2 <- round((a/alf)*log(0.05/((4*a*c)/(fai*Fe))))
t3 <- round((a/alf)*log(0.20/((4*a*c)/(fai*Fe))))

# 劣化度の配列設定
CORR<- matrix(0, n, k)
###### 劣化度設定 ##############################
for (j in 1:k) {
  for (i in 1:n) {
    if (C[i, j] < 1.2) CORR[i, j] <- 0
      else if(C[i, j] < 2.0) CORR[i, j] <- 1
       else CORR[i, j] <- 999
  }
}
################################################
set.seed(123456789) # ★（テスト用）
# 乱数発生数
k <- 100
# 解析年数
n <- 50
###### データ入力(塩化物イオン拡散予測) ########
C0 <- rweibull (k,   2.600, 9.631)
D  <- rlnorm   (k, -17.454, 0.995)
d  <- 6.85
###### データ入力(腐食環境) ####################
a   <- 1500
c   <- 0.005
alf <- 220
fai <- 1.3
Fe  <- 7850
################################################
# 誤差関数の設定
erf <- function(x) 2 * pnorm(x * sqrt(2)) - 1
# 塩化物配列の設定
C <- matrix(0, n, k)
###### 塩化物拡散予測 ##########################
for (i in 1:n) {      
  for (j in 1:k) {  
     C[i,j] <- C0[j]*(1-erf(d/(2*sqrt(D[j]*i*
                            31536000))))
  } 
}
###### 鉄筋腐食予測（年数予測) #################
t1 <- round((a/alf)*log(0.01/((4*a*c)/(fai*Fe))))
t2 <- round((a/alf)*log(0.05/((4*a*c)/(fai*Fe))))
t3 <- round((a/alf)*log(0.20/((4*a*c)/(fai*Fe))))
t12 <- t1+t2   # ★
t123 <- t12+t3 # ★

# 劣化度の配列設定
CORR<- matrix(0, n, k)
###### 劣化度設定 ##############################
for (j in 1:k) {
  nn <- 0 # ★
  for (i in 1:n) {
    if (C[i, j] < 1.2) CORR[i, j] <- 0
    else if(C[i, j] < 2.0) CORR[i, j] <- 1
    else { # ★ ここから
      nn <- nn+1
      CORR[i, j] <- ifelse(nn <= t1, 2,
                    ifelse(nn <= t12, 3,
                    ifelse(nn <= t123, 4, 5)))
    } # ★ ここまで
  }
}
################################################
C <- d/(2*sqrt(D*t(row(C)*31536000)))
C <- t(C0*(1-erf(C)))
result <- t(apply(CORR, 2, function(x) {
	a <- table(x)
	b <- numeric(6)
	b[as.integer(names(a))+1] <- a
	return(b)
}))
result2 <- t(apply(result, 1, cumsum)) # 面グラフを描くために数値を積み上げ
result2 <- cbind(rep(0, nrow(result2)), result2) # 基準線（0）を追加
plot(c(1, nrow(result2)), range(result2), type="n") # plot 枠
for (i in 2:ncol(result2)) { 隣り合う2列が構成する多角形を塗りつぶす
	polygon(c(1:nrow(result2), nrow(result2):1), c(result2[, i-1], rev(result2[, i])), col=i)
}
result2 <- t(apply(result, 1, cumsum))
result2 <- cbind(rep(0, nrow(result2)), result2) 
color   <- gray(seq(1.0, 0.5, length=7))
plot(c(1, nrow(result2)), range(result2),
     xlab="経過年数(年)",
     ylab="各劣化状態の面積割合(%)", type="n") 
for (i in 2:ncol(result2)) { 
  polygon(c(1:nrow(result2), nrow(result2):1),
          c(result2[, i-1],
          rev(result2[, i])), col=color[i])
}   
legend(2, 40, c("潜伏期", "進展期", "加速期",
                "劣化期", "耐力比1.0以下"),
       col=1, fill=color, cex=0.7)
#1
N <- 4 ; TIME <- seq(6)
TIME2 <- rep(TIME, N)
ID <- rep(seq(N), each=length(TIME))
VAR1 <- rep(seq(3), length(TIME)/3*N)
VAR2 <- round(runif(N*length(TIME), min=1, max=9))
(xx <-  data.frame(ID, TIME2, VAR1, VAR2))
layout(matrix(1:4, byrow=TRUE, ncol=2))
old <- par(mar=c(2, 2.7, 1.3, 1), mgp=c(1.6, 0.7, 0))
for (i in 1:N) {
	plot(xx[xx[,"ID"]==i, "TIME2"], xx[xx[, "ID"]==i,
            "VAR1"], xlim=c(0, 7), ylim=c(0, 10),
	     xlab="", ylab="", main=sprintf("ID%02d", i),
            col="red", type="b")
	points(xx[xx[, "ID"]==i, "TIME2"],
              xx[xx[, "ID"]==i, "VAR2"],
	xlim=c(0, 7), ylim=c(0, 10), col="blue", type="b")
	}
par(old)
layout(1)
#2
#library(lattice)
#xyplot(VAR1 ~ TIME2 | ID, data=xx, ylim=c(0, 10),
        type="b", col="red")
#xyplot(VAR2 ~ TIME2 | ID, data=xx, ylim=c(0, 10),
        type="b", col="blue")
Learn more about GTK+ at http://www.gtk.org
If the package still does not load, please ensure
that GTK+ is installed
and that it is on your PATH environment variable
IN ANY CASE, RESTART R BEFORE TRYING TO LOAD THE
PACKAGE AGAIN
 以下にエラー as.GType(type) : Cannot convert
              RGtkBuilder to GType
 追加情報:   警告メッセージ： 
In unzip(path, exdir = .windows_gtk_path) :
   zip ファイルから抽出中に書き込みエラーが生じました
library(RthroughExcelWorkbooksInstaller)
installRthroughExcel()  
library(scatterplot3d)
x <- runif(100, 0, 10)
y <- runif(100, 0, 20)
z <- runif(100, 0, 5)
scatterplot3d(x, y, z, pch=20)
old <- par(mar=c(3, 7, 1, 2), las=3, xpd=TRUE, crt=90)
a <- hist(rnorm(1000), xlab="", main="")
mtext("Rotated Histogram", side=2, line=5, cex=1.5)
text(max(a$breaks)+15*strheight("H", units="figure"),
     -strwidth("H", units="figure")*10, "Value", pos=1)
par(old)
hist2 <- function(x, breaks = "Sturges",
                  xlab="Counts", ylab="Value",
                  main="Horizontal Histogram")
{
	a <- hist(x, breaks=breaks, plot=FALSE)
	class <- length(a$counts)
	plot(range(a$count), range(a$breaks),
            type="n", xlab=xlab,
            ylab=ylab, main=main)
	for (i in 1:class) {
		rect(0, a$breaks[i], a$counts[i],
                    a$breaks[i+1])
	}
}
str <- ""
for (i in c("a", "b", "c")) {
    str <- paste(str, i, sep="")
}
hist(c(0.1*1:10), probability=TRUE)
> x # データ
 [1] 0.1 0.2 0.2 0.3 0.3 0.3 0.4 0.4 0.4 0.4
     0.5 0.5 0.5 0.5 0.5 0.6 0.6
[18] 0.6 0.6 0.6 0.7 0.7 0.7 0.7 0.8 0.8 0.8
     0.9 0.9 1.0
> a <- hist(x,probability=TRUE) # 元のデータのヒストグラム
> a$breaks # 階級幅は 0.2
[1] 0.0 0.2 0.4 0.6 0.8 1.0
> a$density # 各階級の density
[1] 0.500000 1.166667 1.666667 1.166667 0.500000
> sum(diff(a$breaks)[1]*a$density) # 確率の合計
[1] 1
> b <- hist(x*10,probability=TRUE) # 元のデータを 10 倍すると
> b$breaks # 階級幅は 2（元のデータを使った場合の 10倍）
[1]  0  2  4  6  8 10
> b$density # 各階級の density（元のデータを使った場合の 1/10）
[1] 0.0500000 0.1166667 0.1666667 0.1166667 0.0500000
> sum(diff(b$breaks)[1]*b$density)
[1] 1 # 確率の合計はやはり 1
  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.0     3.5   100.0   466.4  1002.0  1005.0 
  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.00    3.50   10.00   49.09  102.50  105.00 
> library(RColorBrewer)
エラー：  パッケージ 'RColorBrewer' は R 2.10.0 以前に造らました。新しくインストールして下さい 
library(nlme)
res.lme <- lme(distance ~ age + Sex,
               random= ~ 1|Subject, data=Orthodont)
anova(res.lme)

library(lme4)
res.lmer <- lmer(distance ~ age + Sex + (1|Subject),
                 data=Orthodont)
anova(res.lmer)
Orthodont$X1<-as.factor(rep(letters[1:3],36))
res.lme2 <- lme(distance ~ age + X1,
                random= ~ 1|Subject, data=Orthodont)  
summary(res.lme2)  
anova(res.lme2) 
   Acid Strange Hard Acid1 Strange2 Hard2
C1  2.0       2  2.0     2        1   2.0
C2  4.0       2  1.0     4        3   4.0
C3  7.0       1  2.5    11        1   7.0
C4  4.5       2  1.5     1        2   2.8
C5  5.0       3  4.0     5        6   4.0
library(FactoMineR)
res.GPA <- GPA(test, group=c(3, 3))
> traceback()
2: coeffRV(Xi, Xj)
1: GPA(test, group = c(3, 3))
   Acid Strange Hard Acid1 Strange2 Hard2
C1  2.0       2  2.0     2        1   2.0
C2  4.0       2  1.0     4        3   4.0
C3  7.0       1  2.5    11        1   7.0
C4  4.5       2  1.5     1        2   2.8
C5  5.0       4  6.0     4        5   4.0
> a <- NaN
> if (a >= 0) 1
 以下にエラー if (a >= 0) 1 : 
      TRUE/FALSE が必要なところが欠損値です 
> a <- -Inf
> if (a >= 0) 1
> if (a < 0) 2
[1] 2
> test
   Acid Strange Hard Acid1 Strange2 Hard2
C1  2.0       2  2.0     2        1   2.0
C2  4.0       2  1.0     4        3   4.0
C3  7.0       1  2.5    11        1   7.0
C4  4.5       2  1.5     1        2   2.8
C5  5.0       4  6.0     4        5   4.0
> res.GPA <- GPA(test, group=c(3, 3))
 以下にエラー
  if (a >= 0) prob <- pgamma(rvstd - (-2/a),
              shape = (4/a^2), scale = (a/2),  : 
   TRUE/FALSE が必要なところが欠損値です 
> test
    Acid Strange Hard Acid1 Strange2 Hard2
C1  2.0       2  2.0     2        1   2.0
C2  4.0       2  1.0     4        3   4.0
C3  7.0       1  2.5    11        1   7.0
C4  4.5       2  1.5     1        2   2.8
C5  5.0       4  6.0     4        5   4.0
> (res.GPA <- GPA(test, group=c(3, 3)))
**Results of the Generalized Procrustes Analysis (GPA)**
There are individuals, characterized by variables
*Results are available in the following objects :
  name           description                                                   
1 "$RV"          "RV Coefficients between partial
                 configurations"              
2 "$RVs"         "standardized RV Coefficients
                 between partial configurations" 
3 "$simi"        "procrustes similarity indexes
                 between partial configurations"
4 "$scaling"     "isotropic scaling factors"                                   
5 "$dep"         "PCA of initial configuration "                               
6 "$consensus"   "coordinates of the consensus
                 configuration"                  
7 "$Xfin"        "coordinates of partial
                 configurations"                       
8 "$PANOVA"      "list of Procrustes Analysis of
                 Variance tables"              
9 "$correlation" "Correlations by sets"                                        

> utils:::menuInstallLocal()
 以下にエラー gzfile(file, "r") :  コネクションを開くことができません 
 追加情報:   警告メッセージ： 
In gzfile(file, "r") :
  圧縮されたファイル 'neural/DESCRIPTION'
  を開くことができません,
  理由は 'No such file or directory' です 
> install.packages("C/neural.zip",repose=NULL)
ファイル名から 'repos = NULL' を推測
 以下にエラー zip.unpack(pkg, tmpDir) : 
   zip ファイル 'C/neural.zip' が見付かりません 
> install.packages("C/neural.zip",contriburl=NULL)
 以下にエラー zip.unpack(pkg, tmpDir) : 
   zip ファイル 'C/neural.zip' が見付かりません 
plot.glm <- function(a)
{
	x <- a$data$x
	y <- a$data$y
	plot(x, y)
	points(x, a$fitted.values, col="red")
	d <- data.frame(x=seq(x[1], x[length(x)],
                       length=500))
	lines(d$x, predict(a, d, type="response"),
             col="red")
}
# 使用例
d <- data.frame(x=1:10, y=c(0,0,0,1,0,1,1,1,1,1))
a <- glm(y~x, data=d, family=binomial)
plot(a)
# 当然ながら，logit 以外にも対応
plot(glm(y~x, data=d,
     family=binomial(link="probit")))
plot(glm(y~x, data=d, family=gaussian))
a<-iris[,1:4]
b<-cor(a)
write.table(b,"C:/b.txt",append=F,
            quote=F,col.names=T,row.names=T)
subset(data, data["variable"] == "2*")
set.seed(123456)
x <- data.frame(before=rnorm(10),
                after=rnorm(10, mean=0.5))
plot(c(1, 2), range(x), type="n", xlab="",
     xaxt="n", ylab="value",
     xlim=c(0.8, 2.2))
junc <- sapply(1:nrow(x),
  function(i) segments(1, x[i, 1], 2, x[i, 2], col=i))
axis(1, at=1:2, labels=colnames(x))
text(2, x$after, 1:nrow(x), col=1:nrow(x),
     pos=4, xpd=TRUE)
for (i in 1:25) {
	eval(parse(text = paste("fv <- fv", i)))
                      # 変数fvにi番目のfvを代入
	#処理
}
以下にエラー parse(text = paste("fv <- fv", i)) : 
   <text>:1:10:  予想外の 数値定数 です  
 1: fv <- fv 1 
             ^
環境：
R version 2.12.1 (2010-12-16)
Platform: x86_64-pc-mingw32/x64 (64-bit)
# データが個々のファイルに用意されており，ファイル名が連番を含む規則的なものであるとき
# 入力すべきファイル名を生成する
file.name <- sprintf("file%03d.dat", 1:200)
   # file001.dat 〜 file200.dat まで
for (fn in file.name) {
  cat(sprintf("file: %s\n"))
     # どのファイルを処理した結果かを明記しておく
  data <- read.csv(fn)
     # いつも同じデータフレームに読み込む
  # data というデータフレームについて処理
}
z <- volcano       
x <- 10*(1:nrow(z))
y <- 10*(1:ncol(z))
brks <- quantile(z, seq(0,1, len = 21))
zi <- z[-nrow(z),-ncol(z)]
cols <- terrain.colors(20)[cut(zi,brks,
        include.lowest = TRUE)]
persp(x, y, 2*z, theta = 120, phi = 15,
      col = cols, scale = FALSE, 
      ltheta = -120, border = NA, box = FALSE) 
警告：  パッケージ 'urca' の前のインストールを取り除くことが出来ませんでした  
パッケージ 'zoo' は無事に開封され、MD5 サムもチェックされました 
ダウンロードされたパッケージは、以下にあります 
C:\Users\***\AppData\Local\Temp\RtmpVFYyw4\downloaded_packages 
install.packages(update[instlib == l, "Package"], l,
contriburl = contriburl,  中で警告がありました： 
 'lib = "C:/PROGRA~1/R/R-212~1.1/library"' は
 書き込み可能ではありません  
以下にエラー install.packages(update[instlib == l, "Package"], 
l, contriburl = contriburl,  : 
パッケージをインストール出来ませんでした
str(ip <- installed.packages(priority = "high"))
ip[, c(1,3:5)]
plic <- installed.packages(priority = "high", fields="License")
## what licenses are there:
table( plic[,"License"] )
> library(help=zoo)　により表示されるウインドウの一部
Package:            zoo
Version:            1.6-5

> library(help=MASS)　により表示されるウインドウの一部
Package:            MASS
Priority:           recommended
Version:            7.3-13
;;essのロード
(load "c:/Meadow/site-lisp/ess/Lisp/ess-site.el")
(setq-default inferior-R-program-name "c:/r/bin/i386/Rterm")
Editor:Meadow 3.01-dev
OS:Windows XP
R:Version 2-13-0
(set-language-environment "Japanese")
(set-default-coding-systems 'shift_jis)
(set-terminal-coding-system 'shift_jis)
(set-keyboard-coding-system 'shift_jis)
(set-buffer-file-coding-system 'shift_jis)
(require 'ess-site)
(setq ess-pre-run-hook
 '((lambda () (setq S-directory default-directory)
              (setq default-process-coding-system 
              '(shift_jis .  shift_jis))
  )))
;;essのロード
(load "c:/Meadow/site-lisp/ess/Lisp/ess-site.el")
&#65533;1つ目に使用したスクリプト
> df <- read.table("clipboard", header=TRUE)
> df
> library(MASS)
> corresp(df, nf=2)
&#65533;2つ目のエラー
> df <- read.table("clipboard", header=TRUE)
> df
> library(MASS)
> corresp(df, nf=2)
以下にエラー corresp.matrix(as.matrix(x), ...) :
empty row or column in table←これが表示がされます
> library(MASS)
> (x <- matrix(c(1,0,0, 0,0,0, 0,0,1, 1,0,1), 4))
    # 第2列が全部 0 です
     [,1] [,2] [,3]
[1,]    1    0    1
[2,]    0    0    1
[3,]    0    0    0
[4,]    0    0    1
> corresp(x)
 以下にエラー corresp.matrix(x) :
   empty row or column in table
> (x <- matrix(c(1,0,0, 0,1,0, 0,0,0, 1,0,1), 4))
    # 第3行が全部 0 です
     [,1] [,2] [,3]
[1,]    1    1    0
[2,]    0    0    1
[3,]    0    0    0
[4,]    0    0    1
> corresp(x)
 以下にエラー corresp.matrix(x) : 
 empty row or column in table
*	Added menu item 'Show Help for current
       Function' CTRL+H
	to context menu of RConsole and each Script
       Editor. It
	parses backwards from the current cursor's
       position to
	to find the current function and opens the
       help page.
	If the user did select a text chunk CTRL+H
       will try
	to come up with help page of the selected
       text.
	If nothing was found it sets the focus to
       the Help 
	Search Field.
	Apoly <- array(c(1,-0.5,0.3,  0,0.2,0.1,
                      0,-0.2,0.7,  1,0.5,-0.3),
                      c(3,2,2))
	Apoly[1,,]
	Apoly[2,,]
	Apoly[3,,]
> Apoly[1,,]
     [,1] [,2]
[1,]    1    0
[2,]    0    1
> Apoly[2,,]
     [,1] [,2]
[1,] -0.5 -0.2
[2,]  0.2  0.5
> Apoly[3,,]
     [,1] [,2]
[1,]  0.3  0.7
[2,]  0.1 -0.3
> Apoly[1:3,,]
, , 1

     [,1] [,2]
[1,]  1.0  0.0
[2,] -0.5  0.2
[3,]  0.3  0.1

, , 2

     [,1] [,2]
[1,]  0.0  1.0
[2,] -0.2  0.5
[3,]  0.7 -0.3
> lapply(1:3, function(x) Apoly[x,,])
[[1]]
     [,1] [,2]
[1,]    1    0
[2,]    0    1 

[[2]]
     [,1] [,2]
[1,] -0.5 -0.2
[2,]  0.2  0.5

[[3]]
     [,1] [,2]
[1,]  0.3  0.7
[2,]  0.1 -0.3
> junk <- apply(Apoly, 1, print)
     [,1] [,2]
[1,]    1    0
[2,]    0    1
     [,1] [,2]
[1,] -0.5 -0.2
[2,]  0.2  0.5
     [,1] [,2]
[1,]  0.3  0.7
[2,]  0.1 -0.3
> aperm(Apoly,c(2,3,1))
, , 1

    [,1] [,2]
[1,]    1    0
[2,]    0    1

, , 2

     [,1] [,2]
[1,] -0.5 -0.2
[2,]  0.2  0.5

, , 3 

     [,1] [,2]
[1,]  0.3  0.7
[2,]  0.1 -0.3
> x <- c("foo", "bar", "baz", "bar")
> (a1 <- factor(x))
[1] foo bar baz bar
Levels: bar baz foo
> as.integer(a1)
[1] 3 1 2 1
> (a2 <- factor(x,
  levels=c("bar", "baz", "foo")))
[1] foo bar baz bar
Levels: bar baz foo
> as.integer(a2)
[1] 3 1 2 1
> (a3 <- factor(x,
  levels=c("foo", "baz", "bar")))
[1] foo bar baz bar
Levels: foo baz bar
> as.integer(a3)
[1] 1 3 2 3
Error in structure(.External("dotTclObjv",
objv, PACKAGE = "tcltk"), class = "tclObj") : 
[tcl] invalid command name "font".
Error :  .onAttach は 'attachNamespace' で
失敗しました 
エラー：  'Rcmdr' に対するパッケージ
もしくは名前空間のロードが失敗しました
要求されたパッケージ tcltk をロード中です 
Tcl/Tkインターフェースのロード中  終了済 
要求されたパッケージ car をロード中です 
要求されたパッケージ MASS をロード中です 
要求されたパッケージ nnet をロード中です 
要求されたパッケージ survival をロード中です 
要求されたパッケージ splines をロード中です
library(misc3d)
## ptsは30点のデータフレームです
## range(pts[,1])
## range(pts[,2])
a <- -0.001709210
b <- 0.001085967
c <- -0.009890907
d <- 1
v <- surfaceTriangles(seq(55,90,len=30),
     seq(45,75,len=30),
     function(x, y){-(a/c)*x-(b/c)*y-(d/c)},
     color2 = "green")
p <- pointsTetrahedra(pts[,1],pts[,2],
     pts[,3],size = 0.1)
## drawScene.rgl(list(v, p))
drawScene(list(v, p))
library(geometry)
library(rgl)
tetramesh(delaunayn(as.matrix(B)),B)
    # Bは曲面Bのデータフレーム
> data
   x  y
1   1  2
2   2  3
3   3  5
4   4  9
5   5  4
6   6  4
7   7 10
8   8 12
9   9 15
10 10 18
11 11 19
12 12 19
13 13 20
14 14 22
15 15 25
16 16 29
17 17 30
18 18 35
19 19 36
20 20 35

> m=glm(y~x, data,
  family=gaussian(link=identity))
> predict(m, newdata=data.frame(x=x),
  se.fit=T)

$fit
        1          2          3          4          5          6 
-0.3428571  1.5458647  3.4345865  5.3233083  7.2120301  9.1007519 
        7          8          9         10         11         12 
10.9894737 12.8781955 14.7669173 16.6556391 18.5443609 20.4330827
       13         14         15         16         17         18 
22.3218045 24.2105263 26.0992481 27.9879699 29.8766917 31.7654135 
       19         20 
33.6541353 35.5428571 
$se.fit
       1         2         3         4         5         6 
0.9877626 0.9129465 0.8408734 0.7723116 0.7082816 0.6501237 
       7         8         9        10        11        12 
0.5995491 0.5586215 0.5295821 0.5144480 0.5144480 0.5295821 
      13        14        15        16        17        18 
0.5586215 0.5995491 0.6501237 0.7082816 0.7723116 0.8408734 
      19        20 
0.9129465 0.9877626 
$residual.scale
[1] 2.292081
> sessionInfo() 
R version 2.12.1 (2010-12-16)
Platform: i386-pc-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=Japanese_Japan.932 
    LC_CTYPE=Japanese_Japan.932   
[3] LC_MONETARY=Japanese_Japan.932
    LC_NUMERIC=C                  
[5] LC_TIME=Japanese_Japan.932    

attached base packages:
[1] stats     graphics  grDevices utils
    datasets  methods  
[7] base     

loaded via a namespace (and not attached):
[1] MASS_7.3-9
jpeg(file.name, quality = 50)
split.screen(c(8,4) )
screen(1)
plot(x,y,type="n")
screen(2)
......
screen(1, new = FALSE)
liens(x,y) 
........
dev.off()
pdf("test.pdf",  width=10, height=5)
split.screen(c(1, 2))
x <- 1:10
y <- x^2
screen(1)
plot(x,y,type="n")
screen(2)
plot(x,y,type="n")
screen(1, new=FALSE)
lines(x,y) 
close.screen(1) # これを入れると期待どおり
screen(2, new=FALSE)
points(x, y)
close.screen(all=TRUE)
dev.off()
require(km.ci)
data(rectum.dat)
require(survival)
fit <- survfit(Surv(time, status) ~ 1,
               data=rectum.dat)
plot(fit$time, 1-fit$surv, type="s")
text(fit$time, 1-fit$surv, fit$n.risk,
     pos=3, xpd=TRUE)
> install.packages("Rcmdr")
 パッケージを ‘C:\Documents and Settings\
                Administrator\Local
                Settings\Application Data\
                R-core/R/win-library/2.12
   中にインストールします 
 (‘lib’ が指定されていないので) 
 以下にエラー gzfile(file, mode) :
 コネクションを開くことができません 
 追加情報:   警告メッセージ： 
In gzfile(file, mode) :
  圧縮されたファイル 'C:\DOCUME~1\ADMINI~1\LOCALS~1\
                       Temp\RtmpNubfVC/libloc_C
%3a%5cDocuments%20and%20Settings%5cAdministrator
%5cLocal%20Settings%5cApplication%20Data%5cR-core%2fR%2fwin-
library%2f2.12Version,Priority,Depends,Imports,
LinkingTo,Suggests,Enhances,OS_type,License,Archs,Built.rds'
を開くことができません, 
理由は 'No such file or directory' です
x <- read.delim("clipboard")
以下にエラー file(file, "rt") :  コネクションを開くことができません :
追加情報:   警告メッセージ： 
In file(file, "rt") :
クリップボードを開くことができないか，中身がありません 
> AllL[2,3]
[1] 87.5242
47 Levels:  84.7952 86.1777 86.4596 86.4788
  86.5069 86.6655 ... Z
> AllL[3,3]
[1] 88.8609
47 Levels:  84.7952 86.1777 86.4596 86.4788
  86.5069 86.6655 ... Z
> AllL[2,3]- AllL[3,3]
[1] NA
 警告メッセージ： 
In Ops.factor(AllL[2, 3], AllL[3, 3]) :
 -  因子に対しては無意味です
library(lattice)
library(grid)
levels <- levels(as.factor(ethanol$C))
xyplot(NOx ~ E|C, data=ethanol,
	panel = function(x, y) {
	panel.xyplot(x, y)
	panel.loess(x, y)
       # 各パネルごとにCの水準を示したい。
	grid.text(label=paste("C=", levels, sep=""), 
	   x=unit(0.05, "npc"), y=unit(0.95, "npc"),
          just="left")})
library(lattice)
library(grid)
levels <- levels(as.factor(ethanol$C))
p <- xyplot(NOx ~ E|C, data=ethanol,
        panel = function(x, y) {
        panel.xyplot(x, y)
        panel.loess(x, y)})
print(p)
update(p, 
 index.cond = list(rep(1, length(levels))))
panel.locs <- trellis.currentLayout()

i <- 1

for (row in 1:nrow(panel.locs))
        for (column in 1:ncol(panel.locs))
              if (panel.locs[row, column] > 0)
  {
        trellis.focus("panel", row = row,
                       column = column,
                       highlight = FALSE)
        grid.text(label=paste("C=",
                  levels[i], sep=""),
                  x=unit(0.05, "npc"),
                  y=unit(0.95, "npc"),
                  just="left")
        trellis.unfocus()
        i <- i + 1
}
> library(mgcv)
> data(airquality)
> airq.gam <- gam(Ozone~s(Solar.R)+s(Wind)+s(Temp),
                  data=airquality)
> airq.gam2 <- gam(Ozone~s(Solar.R, df=10)+
      s(Wind, df=10)+s(Temp, df=10), data=airquality)
> library(gam)
> data(airquality)
> airq.gam2 <- gam(Ozone~s(Solar.R, df=10)+
       s(Wind, df=10)+s(Temp, df=10), data=airquality)
> library(mgcv)
> data(airquality)
> airq.gam <- gam(Ozone~s(Solar.R)+s(Wind)+s(Temp),  
  sp=c(1,3,7),data=airquality)
> print(airq.gam)
 Estimated degrees of freedom:
1.5182 1.4512 1.1420  total = 5.111303 
> plot(pam(data, 3), ask=TRUE)
> a <- pam(iris[,1:4],3)
> cbind(iris[,5], a$clustering)
    [,1] [,2]
[1,]    1    1
[2,]    1    1
[3,]    1    1
       :
[149,]    3    3
[150,]    3    2

> table(iris[,5], a$clustering)

              1  2  3
  setosa     50  0  0
  versicolor  0 48  2
  virginica   0 14 36
> FullModel <- glm(Outcome ~ a+b+c+d+e+f+g+h+i+J+K+L+M+N+O, family=binomial, data=mydata)
summary(stepModel)
Call:
glm(formula = Outcome ~ a+b+c+d+ J + K,
  family = binomial, data = mydata)
Deviance Residuals: 
    Min        1Q    Median        3Q       Max  
-1.53608  -0.65721  -0.36439  -0.09402   2.83204  
Coefficients:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -6.7502     1.8614  -3.626 0.000287 ***
a                      0.8752     0.6204   1.411 0.158350    
b                     -1.2636     0.5959  -2.120 0.033974 *  
c                      1.4079     0.8331   1.690 0.091028 .  
d                      0.8834     0.6045   1.461 0.143914    
T1                    3.6828     0.9318   3.952 7.74e-05 ***
T2                    -0.8963     0.4643  -1.930 0.053555 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’
  0.05 ‘.’ 0.1 ‘ ’ 1 
(Dispersion parameter for binomial family taken to be 1)
   Null deviance: 114.717  on 103  degrees of freedom
Residual deviance:  81.622  on  97  degrees of freedom
AIC: 95.622
GeneticModel <- glmulti(Outcome ~ a+b+c+d+e+f+g+h+i+J+K+L+M+N+O, family=binomial, data=mydata,
method="g", report = FALSE, marginality = TRUE, deltaB = 0, deltaM = 0.01, conseq = 6,
sexrate = 0.15, imm = 0.2, level = 1)
> summary(GeneticModel.glm)
Call:
glm(formula = Outcome ~ a+b+J+K,
  family = binomial, data = mydata)
Deviance Residuals: 
   Min       1Q   Median       3Q      Max  
-1.6986  -0.5952  -0.4402  -0.1794   2.3542  
Coefficients:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)         -5.3365     1.5652  -3.410 0.000651 ***
a                   -1.1344     0.5692  -1.993 0.046268 *  
b                    1.1870     0.7798   1.522 0.127945    
J                    3.4808     0.8931   3.897 9.72e-05 ***
K                   -0.8159     0.4454  -1.832 0.066993 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’
  0.05 ‘.’ 0.1 ‘ ’ 1 
(Dispersion parameter for binomial family taken to be 1)
   Null deviance: 114.717  on 103  degrees of freedom
Residual deviance:  85.879  on  99  degrees of freedom
AIC: 95.879
Number of Fisher Scoring iterations: 5
変数bのp値が0.05をかなり超えてます（0.128）．~
exp(cbind(OR=coef(GeneticModel.glm),
   confint(GeneticModel.glm)))
Waiting for profiling to be done...
                    OR                2.5 %         97.5 %
(Intercept)    0.00481246   0.0001591826    0.08237885
a                   0.32160075   0.1006644837    0.96164214
b                   3.27732742   0.8119375358   18.53407469
J                  32.48671705   6.5308248197   224.16425104
K                   0.44222572   0.1700117573   0.98476474
x <- c(1,2,3,4,5)
y <- c(10,20,30,40,50)
plot(x,y)
>load("./test.R")
エラー：  ファイルのマジック・ナンバーが不正です
（ファイルが壊れているかもしれません）データはロードされませんでした 
追加情報:  Warning message:
ファイル 'test3.R' はマジックナンバー 'x <- ' を持ちます 
二度を越えるセーブバージョンの使用は廃止予定です  
set.seed(12345)
# 元の行列
(a <- matrix(sample(12), 3))                       
# 二重中心化行列   
(b <- t(t(a-rowMeans(a))-colMeans(a))+mean(a))  
# これも可   
# b <- a-outer(rowMeans(a), colMeans(a), "+")+mean(a)
実行結果
> set.seed(12345)
> (a <- matrix(sample(12), 3))
     [,1] [,2] [,3] [,4]
[1,]    9   11    7   12
[2,]   10    4    3    1
[3,]    8    2    5    6
> rowMeans(a)
[1] 9.75 4.50 5.25
> colMeans(a)
[1] 9.000000 5.666667 5.000000 6.333333
> mean(a)
[1] 6.5
# 二重中心化やりま〜〜〜す
> (b <- t(t(a-rowMeans(a))-colMeans(a))+mean(a))
      [,1]       [,2]  [,3]       [,4]
[1,] -3.25  2.0833333 -1.25  2.4166667
[2,]  3.00  0.3333333  0.00 -3.3333333
[3,]  0.25 -2.4166667  1.25  0.9166667
# 二重中心化の効果
> rowMeans(b) # 行平均
[1]  0.000000e+00  2.220446e-16 -2.220446e-16 
> colMeans(b) # 列平均
[1]  0.00000e+00 -5.92119e-16  0.00000e+00  5.92119e-16
> mean(b) # 全体の平均
[1] 0
> 以下にエラー anova(ret.lm, X = anvX[[j]], M = anvM[[j]], idata = lvI, test = "Spherical") : 
  使われていない引数 (X = ~1, M = ~b, idata = list(b = 1:4), test = "Spherical")
OS:Windows7
R version 2.8.1 (2008-12-22) 
i386-pc-mingw32 
locale:
LC_COLLATE=Japanese_Japan.932;
LC_CTYPE=Japanese_Japan.932;
LC_MONETARY=Japanese_Japan.932;
LC_NUMERIC=C;
LC_TIME=Japanese_Japan.932
attached base packages:
[1] datasets  graphics  stats     utils
    grDevices tcltk     
    methods   base     
other attached packages:
[1] Rcmdr_1.4-8 car_1.2-12
R version 2.12.0 (2010-10-15)
Platform: i386-pc-mingw32/i386 (32-bit)
データフレームDF
   group data
1      1 77.4
2      1 78.2
3      1 78.1
4      1 77.8
5      1 77.9
6      2 78.3
7      2 78.2
8      2 78.4
9      2 77.3
10     2 79.1
11     3 79.2
12     3 79.3
13     3 79.1
14     3 78.2
15     3 79.3
16     4 78.9
17     4 78.8
18     4 78.1
19     4 78.1
20     4 78.9 

DF$group <- factor(DF$group)
plotmeans(DF$data ~ DF$group) 

result <- aov(DF$data ~ DF$group)
result.pre <- predict(result,
               interval="confidence")

points(1:4, unique(result.pre[,"lwr"]),
       pch=2, col=2)
points(1:4, unique(result.pre[,"upr"]),
       pch=2, col=2)

m <- sapply(1:4,
   function(x) mean(DF$data[DF$group == x]) )
s <- sapply(1:4,
   function(x) sd(DF$data[DF$group == x]) )
points(1:4, m+s, pch=3, col=3)
points(1:4, m-s, pch=3, col=3)
> tapply(iris$Sepal.Length,iris$Species,var)
    setosa versicolor  virginica 
 0.1242490  0.2664327  0.4043429 
　　setosa versicolor  virginica
0.01635510  0.1828980  0.3032898
> sapply(levels(iris$Species),
         function(x) cov(as.matrix(iris[iris[5]==x,
            c("Sepal.Length","Petal.Length")]))[1,2])
    setosa versicolor  virginica 
0.01635510 0.18289796 0.30328980
> sapply(split(iris, iris$Species),
         function(x)
           cov(x$Sepal.Length, x$Petal.Length))
    setosa versicolor  virginica 
0.01635510 0.18289796 0.30328980 
# 変数名で指定した場合
> head(iris[Petal.Length])
 以下にエラー `[.data.frame`(iris, Petal.Length) : 
   オブジェクト 'Petal.Length' がありません 
# 変数名の文字列で指定した場合
> head(iris["Petal.Length"])
  Petal.Length
1          1.4
2          1.4
3          1.3
4          1.5
5          1.4
6          1.7
> X <- data.frame(x=1:3, y=10:12)
> X["x"]
  x
1 1
2 2
3 3
> x <- 1
> X[x]
  x
1 1
2 2
3 3
> x <- 2
> X[x]
   y
1 10
2 11
3 12
> x <- 3
> X[x]
 以下にエラー `[.data.frame`(X, x) :
  undefined columns selected
day,1,2,3,4,5,6,7,8,9,10,11,12,13,14,
    15,16,17,18,19,20
temp,25,20,18,22,23,24,28,26,21,23,24,
    27,26,22,25,26,24,20,23,25
y <- read.table("temp.txt", header=T)
day<-y[, 1]
temp<-y[, 2]
library(fields)
library(spam)
y.sr<-sreg(day, temp, df=5, method = "GCV")
***** FALSE CONVERGENCE *****
Coefficient(s):
    Estimate  Std. Error  t value Pr(>|t|)    
a0 3.569e-06   5.868e-07    6.082 1.19e-09 ***
a1 8.174e-02   6.977e-03   11.716  < 2e-16 ***
b1 8.949e-01   9.361e-03   95.603  < 2e-16 ***
   name month    x
1     A     1 43.0
2     A     2 41.9
3     A     3 43.6
4     B     1 67.4
5     B     2 27.4
6     B     3 33.9
7     C     1 56.5
8     C     2 59.0
9     C     3 57.5
10    D     1 47.7
11    D     2 43.0
12    D     3 55.1
w <- reshape(data, idvar="name", timevar="month",
             v.names="x", direction="wide")
> w
   name  x.1  x.2  x.3
1     A 43.0 41.9 43.6
4     B 67.4 27.4 33.9
7     C 56.5 59.0 57.5
10    D 47.7 43.0 55.1
l <- reshape(w, direction="long")
> l
    name month    x
A.1    A     1 43.0
B.1    B     1 67.4
C.1    C     1 56.5
D.1    D     1 47.7
A.2    A     2 41.9
B.2    B     2 27.4
C.2    C     2 59.0
D.2    D     2 43.0
A.3    A     3 43.6
B.3    B     3 33.9
C.3    C     3 57.5
D.3    D     3 55.1
エラー：  関数 "garchFit" を見つけることができませんでした 
>  help(package="fSeries")
R version 2.12.1 (2010-12-16)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
> library(kernlab) - spamデータセットのみ利用しました
> library(e1071)
> data(spam)
> m1 <- svm(type~., data=spam, probability = TRUE)
> p1 <- predict(m1, spam)
> p2 <- predict(m1, spam, probability = TRUE)
> table(p1, p2)
         p2
p1        nonspam spam
  nonspam    2834   14
  spam          0 1753
by(measurement[-1], measurement[1], summary)
measurement <- matrix(1:110, nrow=10, ncol=11, byrow=T)
data
factor  mean lower upper
 A       50   40     55
 B       30   19     45
 C       55   25     60
> local({pkg <- 
+ select.list(sort(.packages(all.available = TRUE)),
+ graphics=TRUE)
+ if(nchar(pkg)) library(pkg, character.only=TRUE)})
 要求されたパッケージ MASS をロード中です 
 要求されたパッケージ strucchange をロード中です 
 要求されたパッケージ zoo をロード中です 
 要求されたパッケージ sandwich をロード中です 
Error in gzfile(file, "rb") :
 コネクションを開くことができません 
 追加情報:   警告メッセージ： 
In gzfile(file, "rb") :
  圧縮されたファイル 'C:/Users/fujimot/Documents/R/win-library/
  2.12/sandwich/data/Rdata.rdx' を
  開くことができません, 理由は 'No such file or directory' です 
 エラー：  パッケージ 'sandwich' をロードできませんでした 
> sessionInfo()
R version 2.12.1 (2010-12-16)
Platform: i386-pc-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=Japanese_Japan.932  
    LC_CTYPE=Japanese_Japan.932
[3] LC_MONETARY=Japanese_Japan.932 
    LC_NUMERIC=C
[5] LC_TIME=Japanese_Japan.932

attached base packages:
[1] stats     graphics  grDevices utils
    datasets  methods   base
library(rgdal)
x<-readOGR("c:/GISdata/gunma.kml","gunma")
writeOGR(x2,"c:/GISdata/gunmak","gunmak",
 "ESRI Shapefile")
R Version 2.12.1(2010-12-16)
$ ogr2ogr  -f "ESRI Shapefile" "hoge" KML_Samples.kml
$ ls hoge
Absolute and Relative.dbf  Google Campus.dbf    
Paths.dbf
Absolute and Relative.prj  Google Campus.prj    
Paths.prj
Absolute and Relative.shp  Google Campus.shp    
Paths.shp
Absolute and Relative.shx  Google Campus.shx    
Paths.shx
Extruded Polygon.dbf       Highlighted Icon.dbf 
Placemarks.dbf
Extruded Polygon.prj       Highlighted Icon.prj 
Placemarks.prj
Extruded Polygon.shp       Highlighted Icon.shp 
Placemarks.shp
Extruded Polygon.shx       Highlighted Icon.shx 
Placemarks.shx
> system("ogr2ogr  -f "ESRI Shapefile" "hoge
" KML_Samples.kml")
> M <- readOGR("KML_Samples.kml",
               "Absolute and Relative")
> class(M)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
二度を越えるセーブバージョンの使用は廃止予定です  
R version 2.11.1 (2010-05-31) 
i386-pc-mingw32
> library(biOps)
> data(logo)      # 画像例
> str(logo)
 int [1:275, 1:362, 1:3] 255 255 255 255 255
      255 255 255 255 255 ...
 - attr(*, "type")= chr "rgb"
 - attr(*, "class")= chr [1:2]
     "imagedata" "array"
> # RGB値を輝度値に変換（RGB値がこの順に並んでいるとして）
> lum <- 0.298912 * logo[,,1] + 0.586611 *
         logo[,,2] + 0.114478 * logo[,,3]
> str(lum)
 num [1:275, 1:362] 255 255 255 255 255 ...
> fivenum(lum) # 5数要約
[1]   0.0000 124.1494 165.8013 255.0003 255.0003
> hist(lum) # ヒストグラム
> plot.ecdf(lum) # 経験分布関数
data <- matrix(1:9, 3, 3)
prefix <- c("1:", "2:", "3:")
・・・
処理結果
1:1 2:4 3:7
1:2 2:5 3:8
1:3 2:6 3:9
> (Data <- matrix(as.character(data), 3, 3))
     [,1] [,2] [,3]
[1,] "1"  "4"  "7" 
[2,] "2"  "5"  "8" 
[3,] "3"  "6"  "9" 
> (Prefix <- matrix(prefix, 3, 3, byrow=TRUE))
     [,1] [,2] [,3]
[1,] "1:" "2:" "3:"
[2,] "1:" "2:" "3:"
[3,] "1:" "2:" "3:"
> matrix(paste(Prefix, Data, sep=""), 3, 3)
     [,1]  [,2]  [,3] 
[1,] "1:1" "2:4" "3:7"
[2,] "1:2" "2:5" "3:8"
[3,] "1:3" "2:6" "3:9"
> t(matrix(paste(prefix, t(data), sep=""),
    ncol(data)))
     [,1]  [,2]  [,3] 
[1,] "1:1" "2:4" "3:7"
[2,] "1:2" "2:5" "3:8"
[3,] "1:3" "2:6" "3:9"
> t(apply(data, 1, function(x) paste(prefix, x, sep=""))) # なども
> matrix(paste(prefix, t(data), sep=""),
         ncol=ncol(data), byrow=TRUE)
     [,1]  [,2]  [,3] 
[1,] "1:1" "2:4" "3:7"
[2,] "1:2" "2:5" "3:8"
[3,] "1:3" "2:6" "3:9"
> data <- subset(iris,subset=Species!="versicolor")
> data <- data.frame(class=factor(data$Species),
                     data[,1:4])
> model <- svm(class~.,data=data[-1,],lernel="linear")
> predict(model,data[1,])
     1 
setosa 
Levels: setosa virginica
> predict(model,data[1,],decision.values=TRUE)
     1 
setosa 
attr(,"decision.values")
  setosa/virginica
1         1.080803
Levels: setosa virginica
      X       Y       Z
56.8458 56.3972 95.2184
58.1262 55.2657 96.9028
61.0080 53.8846 97.6925
63.7117 52.5191 98.1147
65.9630 51.1603 98.0884
68.0598 49.8240 97.4626
70.3595 48.6390 96.2383
73.0480 47.8221 94.6007
76.0860 47.6060 92.8562
79.2539 48.1658 91.3380
82.2386 49.5663 90.3180
84.7249 51.7420 89.9494
86.4678 54.5090 90.2453
87.3370 57.6021 91.0912
87.3282 60.7262 92.2840
86.5478 63.6101 93.5838
85.1780 66.0516 94.7681
83.4311 67.9458 95.6757
81.5042 69.2908 96.2331
79.5429 70.1705 96.4590
77.6191 70.7182 96.4472
75.7289 71.0677 96.3317
73.8088 71.3037 96.2427
71.7678 71.4224 96.2643
69.5302 71.3165 96.4051
67.0781 70.7936 96.5927
64.4841 69.6327 96.6992
61.9234 67.6774 96.5962
59.6531 64.9476 96.2303
57.9525 61.7413 95.6913
      X       Y        Z
66.5214 52.3693 100.1890
69.0172 55.2805  98.6580
71.1465 57.6134  98.0792
73.4542 60.2416  96.9707
74.2816 61.4516  95.2830
75.8534 63.4528  93.5104
78.1356 65.5883  94.6485
79.3441 66.4813  96.3968
81.8057 68.7355  97.8613
69.4359 50.9051  99.5420
71.1025 54.5844  99.2688
71.9337 56.6633  98.3382
72.8589 58.7782  97.9511
73.5783 60.6862  96.7916
74.3058 62.7413  95.2100
74.8671 64.3807  93.8145
76.4086 67.2760  95.2155
78.2037 70.3165  97.9255
72.8598 48.6686  97.5556
73.0553 53.8495  97.0382
73.1433 56.6261  96.5826
73.3167 59.6391  96.9136
73.3231 61.7117  95.9461
73.3418 63.9996  95.0022
73.3420 65.9646  94.0229
73.5266 67.9228  94.9986
73.9282 71.3346  97.5466
76.6417 48.3653  95.4308
75.5219 52.4555  96.4391
74.9038 54.7740  97.3086
74.0257 57.8910  97.6356
73.1674 60.8503  97.5045
72.5119 63.0188  96.9356
71.6926 65.7293  96.2245
71.1617 67.6709  96.7142
70.0804 71.6452  97.8153
81.0857 49.4645  93.7371
77.7238 54.0901  94.3492
75.9176 56.5834  95.7828
74.4364 58.6303  97.2458
73.4941 59.9329  98.2465
72.1696 61.7559  98.5740
71.0152 63.3431  98.6402
69.0276 66.0708  98.0351
66.2486 69.8901  97.9546
85.2473 52.2150  93.0255
81.2928 54.9562  92.2683
79.3713 56.2273  93.0238
76.9566 57.8602  93.3170
75.3852 58.7892  95.9743
72.8658 60.3924  98.1361
70.5020 62.0146  97.9852
67.0597 64.3694  97.9051
62.4689 67.4921  98.1250
88.0380 56.8771  94.4088
84.4938 57.8743  92.7572
81.0642 58.7234  92.3931
78.9147 59.1500  93.2910
75.6953 59.7277  95.2876
72.8867 60.1105  98.3210
68.8648 61.0533  98.4587
65.0531 61.9484  98.5723
59.3881 63.4470  96.9480
87.7385 61.3387  95.9346
84.2589 61.3849  93.4009
80.7455 61.1783  92.9229
78.1388 60.8384  94.1021
74.5410 60.2960  96.3323
70.9528 59.7412  98.6698
67.4871 59.4730  98.7284
63.3083 59.2362  98.0871
58.8063 58.9720  97.4714
85.0653 65.4926  97.8556
83.9668 65.5278  94.3416
80.8889 64.4879  91.7567
77.4482 62.5331  93.9201
75.0270 61.1638  95.4027
73.0548 59.9064  97.5156
70.2444 58.5210  97.9352
66.1804 56.5065  98.6132
62.4675 54.5575  99.9249
p <- persp(range(A$x),range(A$y),
           matrix(rep(0,4),ncol=2),
           zlim=range(A$z),theta=10,
           phi=50,r=10)
xyz <- trans3d(A$x,A$y,A$z,p)
polygon(xyz$x,xyz$y,col="red")
> S <- function(p1,p2,p3) {
         w1 <- p1-p3; w2 <- p2-p3
         area <- as.vector(0.5*sqrt((w1%*%w1)*
                        (w2%*%w2)-(w1%*%w2)^2))
         return(area)
       }
> S(c(1,2,3),c(1,1,1),c(5,0,2))
[1] 4.716991
> D <- c(56.8458,56.3972,95.2184,58.1262,...)
> D <- matrix(D, ncol=3, nrow=30, byrow=TRUE)
> F <- function(x) sum((x[1]*D[,1]+x[2]*D[,2]+
                      x[3]*D[,3]+1)^2)/sum(x^2)
> (x <- optim(c(0.1,0.1,0.1),F)) # 適当な初期値で最適化
$par
[1]  0.09792275  0.37648968 -0.32677544
$value
[1] 1366.468
$counts
function gradient 
     116       NA   
$convergence
[1] 0 
$message
NULL
> (x <- optim(x$par,F)) # 上で求めた解を初期値として繰り返す
$par
[1] -0.001709210  0.001085967 -0.009890907  # これが解 a,b,c ?
$value
[1] 64.76301 
$counts
function gradient 
     158       NA  
$convergence
[1] 0 
$message
NULL
> library(geometry)
> eps <- 10^{-14} #ほとんどゼロの幅
# xy平面の単位正方形にepsの厚みを加えた三次元立方体の8個の頂点の座標を並べた行列
> E <- rbind(c(0,0,0),c(1,0,0),c(0,1,0),c(1,1,0),
             c(0,0,eps),c(1,0,eps),c(0,1,eps),c(1,1,eps))
> x <- convhulln(E, options = "FA") # オプションFAで実行
> x
$hull # Eの８個の点の凸包の表面を形成する12組の三角形分割
      [,1] [,2] [,3]
 [1,]    4    3    1 # 3点 E[4,],E[3,],E[1,] からなる
 [2,]    4    2    1 # 三角形という意味
 [3,]    6    2    1
 [4,]    6    5    1
 [5,]    6    4    2
 [6,]    6    4    8
 [7,]    7    3    1
 [8,]    7    5    1
 [9,]    7    4    3
[10,]    7    4    8
[11,]    7    6    5
[12,]    7    6    8

$area　# 凸包の表面積
       #（eps を極小に取ったため、上下の正方形の面積）
[1] 2

$vol　# 凸包の体積 1x1xeps　
[1] 1e-14

> D # 質問者の例の点群
         [,1]    [,2]    [,3]
 [1,] 56.8458 56.3972 95.2184
 [2,] 58.1262 55.2657 96.9028
 [3,] 61.0080 53.8846 97.6925
（途中省略） 
[29,] 59.6531 64.9476 96.2303
[30,] 57.9525 61.7413 95.6913
> x <- convhulln(D, options = "FA")
> x # その凸包は56組の三角形からなる
$hull
      [,1] [,2] [,3]
 [1,]   11   10    1
 [2,]   13   24   14
 [3,]   26   13   30
(途中省略)
[55,]    8    9    1
[56,]    8    7    1 

$area # 凸包の表面積
[1] 1194.405

$vol # 体積
[1] 1264.764

> length(unique(as.vector(x$hull)))
[1] 30 # この場合凸包は30個の点すべてを通る

# 81個の点群の場合
> x <- convhulln(DD, options = "FA")
> x$area
[1] 1116.008
> x$vol
[1] 1644.317
> length(unique(as.vector(x$hull)))
[1] 24 # 凸包は24個の点だけで構成されている！

> # 超平面 ax+by+cz+d=0 に
  # 点 p=c(X,Y,Z) から下ろした垂線の足の座標
> # v=c(a,b,c) と置く。
  # sum(v*ft(v,d,p))+d=0 となる
> ft <- function(v,d,p)
        p - (sum(v*p)+d)/(sum(v*v))*v
> # p=(X,Y,Z) に対する aX+bY+cZ+d の値
> # v=c(a,b,c)。c > 0 の時
  # fft(v,d,p) > 0 (< 0)
  #  ならば点 p は超平面の上(下)側
> fft <- function(v,d,p) sum(v*ft(v,d,p))+d
> # 先に計算した超平面Aに対する a,b,c,dの値。
  # cの値が正になるように調整
> v <- c(0.001709210,-0.001085967,0.009890907)
> d <- -1 
# v,d を固定した関数 fft
> FFT <- function(p) fft(v,d,p)
# 81個の点に対する ax+by+cz+d の値
> x <- apply(DD,1,FFT)
# 超平面Aの上及び上側にある点のインデックス
> xupper <- which(x >= 0)     
# 超平面Aの下側にある点のインデックス
> xlower <- which(x < 0)     
# v,d を固定した関数 ft
> FT <- function(p) ft(v,d,p)
# 超平面Aの上側にある点の平面Aへの射影  の座標行列
> xyzupper <- t(apply(DD[xupper,],1,FT))
# 超平面Aの下側にある点の平面Aへの射影 の座標行列
> xyzlower <- t(apply(DD[xlower,],1,FT))
> 
> # 例：平面Aの上側にある４つの点から定められる面
  #(二つの三角形からなる)
> #     と平面Aで囲まれる図形の表面積と体積を求める
> #     二つの三角形ごとに処理、共通する面の面積
  #     (の二倍)を最後に引く
> D0 <- rbind(c(60,60,100),c(70,60,120),
              c(60,70,110),c(70,70,100))
> tri1 <- c(1,2,4) # その三角形分割
> D1 <- D0[tri1,]  # 対応座標だけを取り出す
> # D1 と D1 の点から平面Aに下ろした点の
  # 凸包(三角柱)の表面積と体積を求める
> x1 <- convhulln(rbind(D1,t(apply(D1,1,FT))),"FA")
> c(x1$area,x1$vol)
[1] 419.3777 211.3755
> tri2 <- c(1,3,4) # もうひとつの三角形分割
> D2 <- D0[tri2,]  # 対応座標だけを取り出す
> x2 <- convhulln(rbind(D2,t(apply(D2,1,FT))),"FA")
> c(x2$area,x2$vol)
[1] 353.1260 360.3836
> eps <- 10^{-10}  # 極小の幅
# 二つの三角形分割に共通する二点 c(1,4)
> tri3 <- intersect(tri1,tri2)   
> # 共通辺をxy方向に少し平行移動してできる
  # 平行四辺形の4頂点の座標行列
> D3 <- rbind(D0[tri3,],D0[tri3,]+c(eps,eps,0)) 
> x3 <- convhulln(rbind(D3,t(apply(D3,1,FT))),"FA")
> c(x3$area,x3$vol)        # 体積はほとんどゼロ
[1] 8.284987e+01 9.942196e-10
# D0と平面Aで囲まれる図形の体積
> x1$vol+x2$vol            
[1] 571.7591
# D0と平面Aで囲まれる図形の面積(共通側面面積の2倍を引く)
> x1$area+x2$area-x3$area 
[1] 689.6538
> # Dは空間点群のx,y,z座標を並べた3列行列
> # TはDの点をある方法で三角形分割したものの
  # 点番号を並べた3列の整数行列
> # (任意の二つの三角形は辺を除いて重ならないとする)
> # ax+by+cz+d=0が参照超平面Aを定義(c=0でないとする)
> # Dの点はすべてAの上側か下側にあるとする
> Poly.av <- function(D,T,a,b,c,d) { 
+   require(geometry) # 必要なパッケージを読み込み
+   nD <- nrow(D)     # 点の総数
+   nT <- nrow(T)             # 三角形分割の総数
+   eps <- 10^{-10}*c(1,1,0)  # 僅かな平行移動用ベクトル
+   Area <- 0; Volume <- 0    # 面積、体積初期値
+   # 超平面 ax+by+cz+d=0 に
+   # 点 p=c(X,Y,Z) から下ろした垂線の足の座標
+   FT <- function(p)
+      p - (sum(c(a,b,c)*p)+d)/(a^2+b^2+c^2)*c(a,b,c)
+   for (i in seq(nT)) {       # 各三角形について繰返し
+     cat("番号 ",i," 対応三角形 (",T[i,],")",fill=TRUE) # デバッグ用
+     D0 <- D[T[i,],]          # 対応座標行列
+     D1 <- t(apply(D0,1,FT))  # D0の各点から超平面に下ろした足の座標行列
+     x <- convhulln(rbind(D0,D1),"FA")
+     Area <- Area + x$area    # 三角柱の共通側面の面積の2倍分余計に加算
+     Volume <- Volume + x$vol
+   }
+   for (i in 1:(nT-1))
      # すべての三角形の組 i,j について繰返し
+     for (j in (i+1):nT) {       
+       itsc <- intersect(T[i,],T[j,]) # 共有頂点番号
+       if (length(itsc)==2) {         # 二つの三角形が一辺を共有すれば
+         cat("共通辺 (", itsc, ")",fill=TRUE)　 # デバッグ用
# その辺をxy方向で少し平行移動した平行四辺形
+         D0 <- rbind(D[itsc,],D[itsc,]+eps)
+         x <- convhulln(rbind(D0,t(apply(D0,1,FT))),"FA")
+         Area <- Area - x$area       # 共有側面の面積の(約)2倍を引く
+       }}
+   return(list(area=Area,volume=Volume))
+ }
> 
> # 簡単な例：6点からなる点群
> D <- rbind(c(60,60,120), #P1
+            c(70,60,110), #P2
+            c(80,60,160), #P3
+            c(70,70,100), #P4
+            c(80,70,140), #P5
+            c(80,80,100)) #P6
> T <- rbind(c(1,2,4),c(2,3,4),
       c(3,4,5),c(4,5,6)) #その三角形分割(の例)       
> v <- c(0.001709210,-0.001085967,
       0.009890907,-1) # 参照超平面の係数 
> (x <- Poly.av(D,T,v[1],v[2],v[3],v[4]))
番号  1  対応三角形 ( 1 2 4 )
番号  2  対応三角形 ( 2 3 4 )
番号  3  対応三角形 ( 3 4 5 )
番号  4  対応三角形 ( 4 5 6 )
共通辺 ( 2 4 )
共通辺 ( 3 4 )
共通辺 ( 4 5 )
$area
[1] 2285.054 

$volume
[1] 994.3007
> D <- c(56.8458,56.3972,95.2184,58.1262,...)
> D <- matrix(D, nrow=30,ncol=3,byrow=TRUE)
> (ctr <- apply(D,2,mean))# 30個の点群の重心座標 
[1] 73.21028 60.43731 95.08736
> DD <- c(66.5214, 52.3693, 100.1890,69.0172, ...)  
> DD <- matrix(DD,nrow=length(DD)/3,ncol=3, byrow=TRUE)
> v <- c(0.001709210,-0.001085967,
         0.009890907,-1) # 参照超平面の係数 
> vv <- v[-4]
  vv <- vv/sqrt(sum(vv^2)) # (a,b,c)を長さ1に

> # 座標行列 D を点 cg を原点、
  # 方向ベクトル vv が (0,0,1) となるように変換
> Rot <- function(D,vv,cg) {
+   D0 <- sweep(D,2,cg)   # Dの各座標からcgを引く
# オイラー角 x[1],x[2],x[3]に対する回転行列
+   tmp1 <- function(x) { 
+     a <- x[1]; b <- x[2]; c <- x[3]
+     R1 <- rbind(c(1,0     ,0      ),
+                 c(0,cos(a),-sin(a)),
+                 c(0,sin(a),cos(a) ))
+     R2 <- rbind(c(cos(b) ,0,sin(b)),
+                 c(0      ,1,0     ),
+                 c(-sin(b),0,cos(b)))
+     R3 <- rbind(c(cos(c),-sin(c),0),
+                 c(sin(c),cos(c) ,0),
+                 c(0     ,0      ,1))
+     R1%*%R2%*%R3 }
 # vvをtmp1(x)で回転したベクトルとc(0,0,1)とのズレ
+   tmp2 <- function(x) {
+     W <- tmp1(x)%*%vv-c(0,0,1)
+     sum(W*W) }
+   # 方向ベクトル vv を (0,0,1) に移す
    # 回転行列のオイラー角を求める
+   # 初期値 (0,0,0) から始め念の為に
    # 最適化を三回繰り返す
+   angle <- optim(optim(optim(c(0,0,0),
             tmp2)$par,tmp2)$par,tmp2)$par
+   t(tmp1(angle)%*%t(D0)) # 変換した座標行列を返す
+ }
> DX <- Rot(DD,vv,apply(D,2,mean))
> DX[1:2,]
          [,1]      [,2]     [,3]
[1,] -5.098094 -9.352273 4.733394
[2,] -3.240725 -5.936903 3.342893
# tripackパッケージを読みこみ
> library(tripack)                              
# z座標が非負の点だけを取り出す
> DXupper <- DX[which(DX[,3]>=0),1:2]           
# そのドロネー三角分割を計算
> tri.upper <- tri.mesh(DXupper[,1],DXupper[,2])
# 三角分割のプロット
> plot(tri.upper)             
# 三角分割要約                  
> summary(tri.upper)                            
triangulation:
Call: tri.mesh(x = DXupper[, 1],
   y = DXupper[, 2]) 
number of nodes: 25 
number of arcs: 65 
number of boundary nodes: 7 
number of triangles: 41 
number of constraints: 0 
# 三角分割(各列が三角形のノード 番号)
> triangles(tri.upper)[,1:3]               
      node1 node2 node3
 [1,]     1    21     9
 [2,]     1     9     8
 [3,]     1     8    17
  (途中省略)
[40,]    22    24    23
[41,]    23    24    25
> 
# z座標が非負の点だけを取り出す
> DXlower <- DX[which(DX[,3]<=0),1:2]          
# そのドロネー三角分割を計算
> tri.lower <- tri.mesh(DXlower[,1],DXlower[,2])
# 三角分割のプロット
> plot(tri.lower)                               
# 三角分割要約
> summary(tri.lower)                            
triangulation:
Call: tri.mesh(x = DXlower[, 1], y = DXlower[, 2]) 
number of nodes: 56 
number of arcs: 156 
number of boundary nodes: 9 
number of triangles: 101 
number of constraints: 0 
# 三角分割(各列が三角形のノード 番号)
> triangles(tri.lower)[,1:3]                    
       node1 node2 node3
  [1,]     1    12    51
  [2,]     1    51    50
  [3,]     1    50    48
  (途中省略)
[100,]    49    51    52
[101,]    52    53    54 

setwd("D:\\")
d1<-read.table("data3.csv",
   header=TRUE, sep=",")
head(d1)
n.d1 <- nrow(d1)
mean.d1 <- apply(d1, 2, mean)
sd.d1 <- sd(d1)
cor.d1 <- cor(d1)
data.frame(n.d1, mean.d1, sd.d1, cor.d1)
library(sem)
model.1 <- specify.model()
x1 <- f1, b11, NA
x2 <- f1, b21, NA
x3 <- f1, b31, NA
x4 <- f1, b41, NA
x5 <- f1, b51, NA
x6 <- f1, b61, NA
x7 <- f1, b71, NA
x8 <- f1, b81, NA
x9 <- f1, b91, NA
x10 <- f1, b101, NA
x11 <- f1, b111, NA
x12 <- f1, b121, NA
x13 <- f1, b131, NA
x14 <- f2, b142, NA
x15 <- f2, b152, NA
x16 <- f2, b162, NA
x17 <- f2, b172, NA
x18 <- f2, b182, NA
f1 <-> f1, NA, 1
f2 <-> f2, NA, 1
f1 <-> f2, c12, NA
sem.model.1 <- sem(model.1,
  cor.d1, N=nrow(d1))
summary(sem.model.1)
以下にエラー scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings,  :
   scan 関数は 'a real' を期待したのに、得られたのは 'N=nrow(d1))'
v <- rep( c(1, 2, 3, 5, 6),
  c(4, 3, 4, 2, 1) )
require(gtools)
combinations( length(v), 3, v,
  set=FALSE, repeats=FALSE )
unique( combinations( length(v), 3, v,
  set=FALSE, repeats=FALSE ) )
getComb <- function(v, n)
{
	search <- function(x, y)
	{
		if (y == 1) {
			a <<- rbind(a, rep(v1, um[1,]))
		}
		else if (x == 1) {
			if ((t <- um[1, 1]-um[y, 1]) >= 0) {
				um[1, 1] <<- t
				search(nc, y-1)
				um[1, 1] <<- um[1, 1]+um[y, 1]
			}
		}
		else {
			search(x-1, y)
			while (um[y, 1] && um[1, x]) {
				um[y, x] <<- um[y, x]+1
				um[y, 1] <<- um[y, 1]-1
				um[1, x] <<- um[1, x]-1
				search(x-1, y)
			}
			um[y, 1] <<- um[y, 1]+um[y, x]
			um[1, x] <<- um[1, x]+um[y, x]
			um[y, x] <<- 0
		}
	}
	f1 <- table(v)
	v1 <- as.numeric(names(f1))
	um <- matrix(0, 2, length(v1))
	um[,1] <- c(n, sum(f1)-n)
	um[1,] <- f1
	a <- NULL
	search(nc, 2)
	return(a)
}
v <- rep( c(1, 2, 3, 5, 6), c(4, 3, 4, 2, 1) )
n <- 7
system.time(b <- unique(combinations(
  length(v), n, v, set=FALSE, repeats=FALSE )))

system.time(a <- getComb(v, n))

実行結果は以下の通り
unique(combinations(... による b とは，出力順が異なる
> head(a)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    3    3    3    3    5    5    6
[2,]    2    3    3    3    5    5    6
[3,]    1    3    3    3    5    5    6
[4,]    2    2    3    3    5    5    6
[5,]    1    2    3    3    5    5    6
[6,]    1    1    3    3    5    5    6
> tail(a)
      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[87,]    1    1    2    2    2    3    3
[88,]    1    1    1    2    2    3    3
[89,]    1    1    1    1    2    3    3
[90,]    1    1    1    2    2    2    3
[91,]    1    1    1    1    2    2    3
[92,]    1    1    1    1    2    2    2
--- Please select a CRAN mirror for
  use in this session ---~
警告：unable to access index for repository
      http://www.stats.ox.ac.uk/pub/RWin/bin/
      windows/contrib/2.4
URL 'http://essrc.hyogo-u.ac.jp/cran/bin/windows/
    contrib/2.4/sem_0.9-8.zip'
    を試しています
Content type 'application/zip' length 245050 bytes
開かれた URL
downloaded 239Kb

以下にエラーzip.unpack(pkg, tmpDir) : ファイル
  'C:/Program Files/R/R-2.4.1/library/fileef56f57/
   sem/chtml/sem.chm'
  を開くことができません
> library(sem)
以下にエラーlibrary(sem) :
  'sem' という名前のパッケージはありません
R version 2.12.1 (2010-12-16)
Platform: i386-pc-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=Japanese_Japan.932  
    LC_CTYPE=Japanese_Japan.932
[3] LC_MONETARY=Japanese_Japan.932 
    LC_NUMERIC=C
[5] LC_
[2,] TIME=Japanese_Japan.932    

attached base packages:
[1] stats     graphics  grDevices utils     
    datasets  methods   base     

other attached packages:
[1] urca_1.2-4

loaded via a namespace (and not attached):
[1] grid_2.12.1     lattice_0.19-13 nlme_3.1-97     tools_2.12.1   
> A <- 1:5
> A
[1] 1 2 3 4 5
> rev(A)
[1] 5 4 3 2 1
> A <- "ABCDE"
>  system(paste("echo",A,"|rev"),intern=TRUE)
[1] "EDCBA"
> a <- "abcdefghij"
> paste(rev(unlist(strsplit(a, ""))),
   collapse="")
[1] "jihgfedcba"
> Napoleon <- "Able was I ere I saw Elba"
> paste(rev(unlist(strsplit(Napoleon, ""))),
   collapse="")
[1] "ablE was I ere I saw elbA"
pdf <- function (n, u=1, a=1) {
  x <- runif(n, min=a, max=1000)
  p <- x^-u
  y <- s  vv ample(x, n, replace=T, prob=p)
  return(y)
}
rpower <- function(n, m=1, M=100, b=0.1)
{
	((M^-b-m^-b)*runif(n)+m^-b)^(1/-b)
}
flb <- function(x)
   { p <- length(x); sum(c(1, rep(4, p-1)) *
          (x - c(1, x[-p])^2)^2)
}
## 25-dimensional box constrained
optim(rep(3, 25), flb, NULL,
method = "L-BFGS-B",
     lower=rep(2, 25), upper=rep(4, 25))
# par[24] is *not* at  boundary
sum(c(1, rep(4, p-1)) * (x - c(1, x[-p])^2)^2)
c(1, rep(4, p-1))
(x - c(1, x[-p])^2)^2
sum(c(1,2,3)*(1,2,3))
lower=rep(2, 25), upper=rep(4, 25)
> (x <- c(3, 5, 2)) # 例えば，こんな例を
[1] 3 5 2
> (p <- length(x)) # x の長さ（要素数）
[1] 3
# 1 と 4 を p-1 個繰り返した要素を持つベクトル
> c(1, rep(4, p-1))
[1] 1 4 4
# 1 と x の p 番目の要素を除いたものを要素とするベクトル
> c(1, x[-p])
[1] 1 3 5
> c(1, x[-p])^2 # その二乗
[1]  1  9 25 
> x - c(1, x[-p])^2 # x から引き算
[1]   2  -4 -23
> (x - c(1, x[-p])^2)^2 # さらにその二乗
[1]   4  16 529
> c(1, rep(4, p-1)) *
   (x - c(1, x[-p])^2)^2 # ベクトルの積
[1]    4   64 2116
> sum(c(1, rep(4, p-1)) *
   (x - c(1, x[-p])^2)^2) # その要素の和
[1] 2184
a <- function(b, c) {
    d <- e+f
    if (g > 0) {
        h <- i+8
    }
    else {
        h <- i+90
    }
    return(h*900)
}
optim(par=c(0.927, 0.982, 0.091, 0.058, 0.604,
            1.129, 0.368, 0.326),
      fn=Likelihood(ε, R),
      control=list(fnscale=-1, trace=TRUE))
nlm(Likelihood, c(0.927, 0.982, 0.091, 0.058, 0.604,
                  1.129, 0.368, 0.326), hessian=TRUE)
> f <- function(x) sum((x-1:length(x))^2)
> nlm(f, c(10,10))
$minimum
[1] 4.303458e-26

$estimate
[1] 1 2

$gradient
[1]  2.757794e-13 -3.099743e-13

$code
[1] 1

$iterations
[1] 2

> optim(c(10, 10), f,
    control=list(reltol=1e-14))
$par
[1] 0.9999998 1.9999994

$value
[1] 4.218356e-13

$counts
function gradient 
     103       NA 

$convergence
[1] 0

$message
NULL
ε<-x[,10]
R<-x[,5]
y<-c(1:305)
Likelihood<-function(ε,R){
return(function(par){
p<-par[1]
q<-par[2]
ω_0<-par[3]
α_0<-par[4]
β_0<-par[5]
μ_0<-0.1
ω_1<-par[6]
α_1<-par[7]
β_1<-par[8]
μ_1<--0.2

LogL<-0
E<-(1-p)/(2-p-q) # t=0の定常確率(スタート)
F<-(1-q)/(2-p-q)
ε_0<-ε[1] # t=1の確率密度関数計算
R_1<-R[2]
A<- q*E+(1-p)*F
B<- (1-q)*E+p*F
sigma_0_1<-sqrt(ω_0+α_0*ε_0^2+β_0*ε_0^2)
sigma_1_1<-sqrt(ω_1+α_1*ε_0^2+β_1*ε_0^2)
f1<-1/sqrt(2*pi*sigma_1_1^2)*exp(-(R_1-μ_1)^2/
   (2*sigma_1_1^2))
f0<-1/sqrt(2*pi*sigma_0_1^2)*exp(-(R_1-μ_0)^2/
   (2*sigma_0_1^2))
f<- f0*A+f1*B
LogL<-LogL+log(f)

E<-f0*A/(f0*A+f1*B)
F<-f1*B/(f0*A+f1*B)
G<-E
H<-F
C<-A
D<-B

ε_1<-ε[2]      #t=2
R_2<-R[3]
A<-q*E+(1-p)*F
B<-(1-q)*E+p*F
sigma_0_2<-sqrt(ω_0+α_0*ε_1^2+β_0*
  (sigma_0_1^2*C+sigma_1_1^2*D))
sigma_1_2<-sqrt(ω_1+α_1*ε_1^2+β_1*
  (sigma_0_1^2*C+sigma_1_1^2*D))
f1<-1/sqrt(2*pi*sigma_1_2^2)*
  exp(-(R_2-μ_1)^2/(2*(sigma_1_2)^2))
f0<-1/sqrt(2*pi*sigma_0_2^2)*
  exp(-(R_2-μ_0)^2/(2*(sigma_0_2)^2))
f<-f0*A+f1*B
LogL<-LogL+log(f)

E<-f0*A/(f0*A+f1*B)
F<-f1*B/(f0*A+f1*B)
A=q*E+(1-p)*F
B=(1-q)*E+p*F
C=q*G+(1-p)*H
D=(1-q)*G+p*H
sigma0t<-sigma_0_2
sigma1t<-sigma_1_2

z<-2　　　　　　　　　　　#t=3以降
ε_t<-0
R_t<-0
for(i in y){          
z<-z+1
ε_t<-ε[c(z)]
R_t<-R[c(z+1)]
A<-q*E+(1-p)*F
B<-(1-q)*E+p*F
sigma0t<-sqrt(ω_0+α_0*ε_t^2+
  β_0*((sigma0t)^2*C+(sigma1t)^2*D))
sigma1t<-sqrt(ω_1+α_1*ε_t^2+
  β_1*((sigma0t)^2*C+(sigma1t)^2*D))
f1<-1/sqrt(2*pi*(sigma1t)^2)*
  exp(-(R_t-μ_1)^2/(2*(sigma1t)^2))
f0<-1/sqrt(2*pi*(sigma0t)^2)*
  exp(-(R_t-μ_0)^2/(2*(sigma0t)^2))
f<-f0*A+f1*B
LogL<-LogL+log(f)

E<-f0*A/(f0*A+f1*B)
F<-f1*B/(f0*A+f1*B)
C<-A
D<-B
}
return(LogL)
})
}
optim(par=c(0.927,0.982,0.091,0.058,
      0.604,0.9,0.368,0.326),
      fn=Likelihood(ε,R),
      control=list(fnscale=-1, trace=TRUE))
Likelihood<-function(ε,R){
return(function(par){
return(LogL)
})
}
> ε<-x[,10]
 エラー：  オブジェクト 'x' がありません 
> R<-x[,5]
 エラー：  オブジェクト 'x' がありません 
  :
 エラー：  オブジェクト 'z' がありません 
> return(LogL)
 エラー：  オブジェクト 'LogL' がありません 
> })
 エラー：   予想外の '}' です  ( "}" の) 
> }
 エラー：   予想外の '}' です  ( "}" の) 
> optim(par=c(0.927,0.982,0.091,0.058,
  0.604,0.9,0.368,0.326),fn=Likelihood(ε,R),
        control=list(fnscale=-1, trace=TRUE))
 Nelder-Mead direct search function minimizer
 以下にエラー function (par)  : 
   関数 "Likelihood" を見つけることができませんでした 
> x <- matrix(runif(305*2), 305)
> ε<-x[,10]
 エラー：  添え字が許される範囲外です 
> R<-x[,5]
 エラー：  添え字が許される範囲外です 
set.seed(12345)
ε<-runif(308)
R<-runif(309)
44: In log(f) :  計算結果が NaN になりました 
45: In sqrt(ω_0 + α_0 * ε_t^2 + β_0 *
  ((sigma0t)^2 * C +  ... :
   計算結果が NaN になりました 
NewFunctionName <- function(par) { 関数名は何でも良い
p<-par[1]
以下同じ（この中で εや R を参照すると，大域変数が参照されるので，引数で渡す必要がない
return(LogL)
}
nlm(NewFunctionName, c(0.927,0.982,0.091,
  0.058,0.604,0.9,0.368,0.326), hessian=TRUE)
y <- c(1:7)
ε <- c(0.104794415, 0.042511256, 0.125477647,
        0.314732511, -0.411718232,
        0.670690239, -0.006636651, -1.099916039,
       -0.380745473, 1.415399685)
R <- c(0.043765052, 0.019836787, 0.103638308,
       0.289968381, -0.442985542,
       0.665004049, - 0.050993233, -1.11928331,
      -0.362827922, 1.406915882)

Likelihood <- function(ε, R) {
    return(function(par) {
        p <- par[1]
        q <- par[2]
        ω_0 <- par[3]
        α_0 <- par[4]
        β_0 <- par[5]
        ω_1 <- par[6]
        α_1 <- par[7]
        β_1 <- par[8]
        μ_0 <- 0.1
        μ_1 <- -0.2

        LogL <- 0
        E <- (1-p)/(2-p-q) # t=0 の定常確率(スタート)
        F <- (1-q)/(2-p-q)

        ε_0 <- ε[1] # t=1 の確率密度関数計算
        R_1 <- R[2]
        A <-  q*E+(1-p)*F
        B <-  (1-q)*E+p*F
        sigma_0_1 <- sqrt(ω_0+α_0*
           ε_0^2+β_0*ε_0^2)
        sigma_1_1 <- sqrt(ω_1+α_1*
           ε_0^2+β_1*ε_0^2)
        f1 <- 1/sqrt(2*pi*sigma_1_1^2)*
            exp(-(R_1-μ_1)^2/(2*sigma_1_1^2))
        f0 <- 1/sqrt(2*pi*sigma_0_1^2)*
            exp(-(R_1-μ_0)^2/(2*sigma_0_1^2))
        f <-  f0*A+f1*B
        LogL <- LogL+log(f)

        E <- f0*A/(f0*A+f1*B)
        F <- f1*B/(f0*A+f1*B)
        G <- E
        H <- F
        C <- A
        D <- B

        ε_1 <- ε[2]      # t=2
        R_2 <- R[3]
        A <- q*E+(1-p)*F
        B <- (1-q)*E+p*F
        sigma_0_2 <- sqrt(ω_0+α_0*ε_1^2+
          β_0*(sigma_0_1^2*C+sigma_1_1^2*D))
        sigma_1_2 <- sqrt(ω_1+α_1*ε_1^2+
          β_1*(sigma_0_1^2*C+sigma_1_1^2*D))

        f1 <- 1/sqrt(2*pi*sigma_1_2^2)*
           exp(-(R_2-μ_1)^2/(2*(sigma_1_2)^2))
        f0 <- 1/sqrt(2*pi*sigma_0_2^2)*
           exp(-(R_2-μ_0)^2/(2*(sigma_0_2)^2))
        f <- f0*A+f1*B
        LogL <- LogL+log(f)

        E <- f0*A/(f0*A+f1*B)
        F <- f1*B/(f0*A+f1*B)
        A <- q*E+(1-p)*F
        B <- (1-q)*E+p*F
        C <- q*G+(1-p)*H
        D <- (1-q)*G+p*H

        sigma0t <- sigma_0_2
        sigma1t <- sigma_1_2

        z <- 2
        ε_t <- 0
        R_t <- 0

        for (i in y) {          # t=3 以降
            z <- z+1
            ε_t <- ε[c(z)]
            R_t <- R[c(z+1)]

            A <- q*E+(1-p)*F
            B <- (1-q)*E+p*F
            sigma0t <- sqrt(ω_0+α_0*ε_t^2+
              β_0*((sigma0t)^2*C+(sigma1t)^2*D))
            sigma1t <- sqrt(ω_1+α_1*ε_t^2+
              β_1*((sigma0t)^2*C+(sigma1t)^2*D))
            f1 <- 1/sqrt(2*pi*(sigma1t)^2)*
              exp(-(R_t-μ_1)^2/(2*(sigma1t)^2))
            f0 <- 1/sqrt(2*pi*(sigma0t)^2)*
              exp(-(R_t-μ_0)^2/(2*(sigma0t)^2))
            f <- f0*A+f1*B
            LogL <- LogL+log(f)

            E <- f0*A/(f0*A+f1*B)
            F <- f1*B/(f0*A+f1*B)
            C <- A
            D <- B
        }

        return(LogL)
    })
}
optim(par=c(0.9, 0.9, 0.3, 0.3, 0.4,
      0.9, 0.2, 0.5), fn=Likelihood(ε, R),
      control=list(fnscale=-1, trace=TRUE))
optim(rep(3, 25), flb, NULL, method = "L-BFGS-B",
      lower=rep(2, 25), upper=rep(4, 25))
      # par[24] is *not* at boundary
splom (~ USArrests[c(3, 1, 2, 4)] | state.region, pscales = 0, types = c("g", "p", "smooth"))
barchart (VADeaths, groups = F, layout = c(2,2))
dat.long <- stack (data.wide)
histogram( ~ values | ind, data = data.long,
          layout = c(2,2))
library( ggplot2)
p <- ggplot(diamonds, aes(carat)) + 
     geom_histogram(binwidth = 1)
p + facet_grid(cut ~ .)
library(lattice)
x <- as.data.frame(matrix(round(rlnorm(100*10,
    log(25), 0.4)), 100))
# xの中身を確認しておく
names(x) <- sprintf("Q%02d", 1:10)
histogram( ~ values | ind, data = stack(x),
    layout = c(5, 2))
set.seed(123)
x <- as.data.frame(matrix(round(rlnorm(
   100*10, log(25), 0.4)), 100))
# xの中身を確認しておく
names(x) <- sprintf("Q%02d", 1:10)
layout(matrix(1:10, byrow=TRUE, ncol=5))
old <- par(mar=c(2, 2.7, 1.3, 1),
    mgp=c(1.6, 0.7, 0))
for (i in 1:ncol(x)) {
	hist(x[, i], breaks=seq(-5, 100, by=5),
            xlim=c(0, 80), ylim=c(0, 25),
	     xlab="", main=sprintf("Q%02d", i),
            col="darkorchid1")
	}
par(old)
layout(1)
a <- c("a.1", "a.2", "b.1", "b.2",
       "c.1", "c.2")
b <- array(a, dim=6)
y <- seq(1, 6, 2)
split(b, y)
$`1`
[1] "a.1" "b.2"

$`3`
[1] "a.2" "c.1"

$`5`
[1] "b.1" "c.2"
> a <- c("a.1", "a.2", "b.1", "b.2",
         "c.1", "c.2")
> split(a, rep(1:3, each=2))
$`1`
[1] "a.1" "a.2"

$`2`
[1] "b.1" "b.2"

$`3`
[1] "c.1" "c.2"
[Workspace restored from /Users/foo/.RData]
[History restored from /Users/foo/.Rhistory]
library(pls)
data(yarn)
yarn.plsr <- plsr(density ~ NIR, 6,
   data = yarn, validation = "CV") 
plot(yarn.plsr, "prediction",
    which = "train", 6, line = T)
ss.t <- sum((yarn$density-
   mean(yarn$density))^2)
ss.e <- sum(yarn.plsr$residuals[,,
   yarn.plsr$ncomp]^2)
1-ss.e/ss.t
a <- c(0, 1, 1, 1, 0, 0, 1, 1)
b <- c(1, 1, 1, 1, 1, 0, 1, 0)
> sum(unlist(strsplit(gsub("01", "@",
    paste(a, collapse="")), "")) == "@")
[1] 2
> sum(unlist(strsplit(gsub("01", "@",
    paste(b, collapse="")), "")) == "@")
[1] 1
> a <- c(1, 1, 1, 0, 0, 1)
> aa <- a[-length(a)]
> sum( aa[which(aa==0)+1] == 1 )
[1] NA
> a <- c(1, 1, 1, 0, 0, 0)
> aa <- a[-length(a)]
> sum( aa[which(aa==0)+1] == 1 )
[1] NA
> sum(a[-length(a)]*10^8+a[-1]==0*10^8+1)
[1] 2
　path worst cost pathpattern throughput
stddev Numpulse Nummeasure AveF
1 □ 0 1 0 217.7 0.682583 90441
140519 0.985497
2 ○ 1 1.83 1 118.447 0.255297
163628 254691 0.984438
3 □□ 0 2 0 94.7204 0.29984
392181 609967 0.994922
4 △ 2 2.93 2 74.3604 0.188577
258852 404117 0.985847
5 □□□ 0 3 0 98.3713 0.170828
476040 739740 0.980546
6 □×○□ 0 4 0 97.9405 0.171691
753243 1172179 0.98989
df <- read.table("ex_data.txt", header=T)
sortlist <- order(df$throughput)
df <- df[sortlist,]

cs <- as.character(df$path)

asdf <- barplot(df$throughput, 0.1,
ylim=c(0, 250), xlab="", ylab="", axes=F, ann=F)
axis(1, at=asdf, labels=FALSE)
axis(2, at=seq(0, 250, 10))

library(gridBase)
vps <- baseViewports()
pushViewport(vps$inner, vps$figure, vps$plot)
grid.text(cs, x=unit(asdf, "native"),
y=unit(-1, "lines"), just="right", rot=90)
popViewport(3)
警告メッセージ： 
1: In grid.Call.graphics("L_texta.1",
as.graphicsAnnot(x$label),  ... :
   'mbcsToSbcs' 中の 'ﾃ療療療冷味笆｡笆｡笆｡
' で変換に失敗: <e2>をドットで置き換えました 
dev.copy2eps(file="throughput+cost.eps")
dev.copy(png, "throughput+cost.png")
dev.off()
# dataオリジナル設定
> data <- data.frame(x=1:10, y=1:10, z=1:10)
> names(data) # 変数名の確認
[1] "x" "y" "z"
> data$newx=data$x # 新しい変数名を定義
# データフレームから変数xを削除
> data <- data[setdiff(colnames(data), "x")]　
> names(data) # 変数名の確認
[1] "y"    "z"    "newx"
> library(gdata)
> data <- rename.vars(data, c("x", "y", "z"),
  c("first", "second", "third"))
> data <- data.frame(x=1:10,
                     y=1:10, z=1:10)
> names(data)
[1] "x" "y" "z"
> names(data)[2] <- "Y"
> names(data)
[1] "x" "Y" "z"

> names(data)[names(data)=="x"] <- "X"
> names(data)
[1] "X" "Y" "z"

> names(data)[2:3] <- c("A","B")
> names(data)
[1] "X" "A" "B"

> names(data)[names(data) %in% c("X", "B")] <-
    c("foo", "bar")
> names(data)
[1] "foo" "A"   "bar"
> data <- data.frame(x=1:10, y=1:10, z=1:10)
# x, z のみを foo, bar に変更（付値を忘れないように）
> data <- rename.vars(data, c("x", "z"),
  c("foo", "bar")) 

Changing in data             
From: x   z  
To:   foo bar

> names(data)
[1] "foo" "y"   "bar"

> data <- rename.vars(data, "y", "ほげほげ")

Changing in data              
From: y       
To:   ほげほげ

> names(data)
[1] "foo"      "ほげほげ" "bar"
data <- 
> ?glm
 警告メッセージ： 
In file.show(temp, title = gettextf("R Help on '%s'",
   topic), delete.file = TRUE) :
   file.show()：ファイル 'C:\DOCUME~1\蜿､蟾晏忽蠢予
                LOCALS~1\Temp\RtmpIDu7yb\Rtxt678418be'
                は存在しません 
> help("glm")
 警告メッセージ： 
In file.show(temp, title = gettextf("R Help on '%s'",
   topic), delete.file = TRUE) :
   file.show()：ファイル 'C:\DOCUME~1\蜿､蟾晏忽蠢予
                LOCALS~1\Temp\RtmpIDu7yb\Rtxt3d6c4ae1'
                は存在しません
Microsoft Windows XP Professional
Version 2002
Service Pack 3
Intel(R) Core(TM)2 Duo CPU
E8400 @ 3.00GHz
2.99 GHz、976 MB RAM
R version 2.12.1 (2010-12-16)
> utils:::menuInstallPkgs()
 --- このセッションで使うために、CRANのミラーサイトを選んでください --- 
 install.packages(NULL, .libPaths()[1L],
  dependencies = NA, type = type) 中で警告がありました： 
  'lib = "C:/PROGRA~1/R/R-212~1.1/library"' 
  は書き込み可能ではありません  
 以下にエラー sprintf(msg, userdir) : 
   <83>ｼ<e3>Ν縺吶ｋ縺溘ａ縺ｫ蛟倶ｺｺ逧<84>縺ｪ繝ｩ繧､繝悶Λ繝ｪ 
 '%s' 
繧剃ｽ懊ｊ縺溘＞縺ｮ縺ｧ縺吶°<ef>ｼ<9f> に
不正なマルチバイト文字があります 
> F <- function(x) sum((x[1]*DD[,1]+
       x[2]*DD[,2]+x[3]*DD[,3]+1)^2)/sum(x^2)
# 適当なパラメータ初期値で始める
> (x <- optim(c(0.1,0.1,0.1),F))  
$par
[1]  0.09792275  0.37648968 -0.32677544

$value
[1] 1366.468

$counts
function gradient 
     116       NA 

$convergence
[1] 0

$message
NULL

#　上で求めた解を初期値として繰り返す
> (x <- optim(x$par,F))　
$par
# これが求めたい a,b,c の値？
[1] -0.001709210  0.001085967 -0.009890907

$value
[1] 64.76301

$counts
function gradient 
     158       NA 

$convergence
[1] 0

$message
NULL

Ｑ＆Ａ (初級者コース)/13

場合分けに応じた新しい列を追加†

geoRパッケージの標本バリオグラムについて†

R-2.14以降、Tcl/Tkのスライダーを使うとグラフが一部消える†

メモリの割り当てについて†

gstatライブラリーのvariogramについて†

繰り返し処理について†

変数を使って，levels()などの関数を使う方法†

function定義と代入の違い†

なぜか発生した誤差で、四捨五入が妨げられる†

scatterplot3dについて†

スクリプトからのコンソール画面消去は実行することが出来ますか？†

例えば dgamma の lse を求める†

type="b"のplotと整合するlegend()の使い方†

cat関数の特殊文字について†

結果の解釈†

RDA解析の解釈†

deprecatedとは？†

一般化線形モデルの推定結果の図示（offset項を入れた時）†

read.csvで、header=TRUE, row.names=1として読み込んだとき†

行列の行に対してcumsumを適用することは出来るでしょうか？†

readBin関数のサイズ指定†

ccfでp-value†

プロキシサーバーの設定について†

文字列ベクトルからの共通要素の抽出†

ヘッダに括弧が含まれるCSVから正しくヘッダを取得したい†

column の名前を入力せずに lm()†

ライブラリーscatterplot3について†

Rserveを使ってload()をする回数を抑える方法について†

Rで出した結果の見方について†

非線形最適化について†

集合演算について†

Mac OS X上のpng()で文字化け†

データの型を指定して出力†

並列処理したい関数における外部変数の扱い†

コマンドのコンソール出力を抑えるには†

persp()で軸のタイトルの文字サイズの変更†

plyrのロードができない†

Boxプロットで、項目軸を数値にする方法†

行列の割合表示†

ARIMAモデルでの解析結果を詳細に表示したい場合…†

2階の導関数のもとめかた†

連続の検出†

read.tableで読み込んだ行列がis.finiteでFALSE??†

LSMEANを求めながら，そのデータのsdを求める†

csvからxmlに変換したときの日本語文字化けについて†

座標値の等間隔プロット†

Excelから１セルの読込†

Rで他言語の展開†

t.test()の機能†

順序ロジットモデルの適合度について†

散布図のプロットの中で条件によって色を変える方法†

マルチコアCPUでの、Rの効率的な利用について†

Rで作図をする際に、レイヤー表示のように手前に持ってくるキャラクターの指定はできないのでしょうか？†

Rコマンダーがインストールできない†

GAM(一般化加法モデル)における変数選択について†

コネクションエラー†

パネルデータ分析 固定効果†

interp1?†

ARIMAモデルの定数項について†

splancsのkhatについて†

クラスター分析での項目名†

リッジ回帰について†

パッケージFEARのdeaでNDRSを実行したいのですが。†

線型等式制約Ax=bにおける最適化問題を解く方法†

２元配置分散分析の場合の等分散性の検定†

plot()出力の部分拡大†

Help†

R CMD BATCH < test.R > test.Rout†

パッケージSOMの仕様†

lm について†

関数の追加†

多項ロジット†

グンベル分布推定†

hclustでのward計算の原理†

2次元正規分布関数のパラメータ推定について†

コクランのC検定†

平滑化された値の出力について†

constrOptimのhessianは表示可能でしょうか？†

Map (mapply)から関数のリストを返すにあたって†

パネルデータ分析　固定効果†

順序ロジット　関数 polr() について†

JACKNIFE　という統計用語はどんなことですか？†

Rで重回帰　再投稿†