R_fromOldHtml3 のバックアップ(No.4) - アールメカブ

アールメカブ


R_fromOldHtml3 のバックアップ(No.4)


個人メモを整理中です.不適切な記述が多々あるかと思います.お気づきの際は,ishida-m.at.ias.tokushima-u.ac.jp. までご連絡ください.

_ scatterplot.matrix 散布図行列の作成 2007 01 18

plot だと対角要素は変数名

対角要素をヒストグラムにしたければ

car パッケージの

 scatterplot.matrix

例えば エヴェリット, RとS-PLUSによる多変量解析

 usair.dat<-source(
 "/everitt/Data/chap3usair.dat")$value
 # install.packages("car")
 library(car)
 scatterplot.matrix(usair.dat[,-1],diagonal= "histo" )
 plot(usair.dat)

_ エヴェリット, RとS-PLUSによる多変量解析 の累積分布密度グラフ 2007 04 05

 pdf("p9.pdf")
 par(mar=c(2,2,2,2))
 plot(seq(-10, 10, 0.01), seq(0, 1, 0.0005), type = "n",
   axes =F, xlab = "", ylab = "", 
   main = "Cumulative distribution function")
 lines(x[order(x)], xp, xlim=c(-10, 10), lwd = 3 )
 lines(y[order(y)], yp,  xlim=c(-10, 10), lwd = 3 )
 lines(c(-9,9), c(1,1))
 lines(c(-9,-9), c(1,0))
 lines(c(-9,9), c(0,0))
 lines(c(9,9), c(1,0))
 arrows(-9, 0.8, x.p0.8, 0.8, length = 0.2, 
   lty = 2, lwd = 2)
 arrows(-9, 0.8, y.p0.8, 0.8, length = 0.2, 
   lty = 1, lwd = 2)
 arrows( x.p0.8, 0.8,  x.p0.8, 0, length = 0.2, 
   lty = 2, lwd = 2)
 arrows(y.p0.8, 0.8, y.p0.8, 0, length = 0.2, 
   lty = 1, lwd = 2)
 text(3.8, 0.85, expression(paste("c.d.f ", F[x], "( )")), 
 cex = 1.1)
 text(-2.0, 0.85,  expression( paste("c.d.f ", F[y], 
   "( )")), cex = 1.1)
 text(-10, 0.8, "p", cex = 1.1)
 text(x.p0.8+0.5, -.02, expression(q [x] (p)), cex = 1.1)
 text(y.p0.8, -.02, expression(q [y] (p)), cex = 1.1)
 arrows(x.p0.2, 0.0, x.p0.2, 0.2, length = 0.2, 
   lty = 1, lwd = 2)
 arrows(x.p0.2, 0.0, x.p0.2,  pnorm( x.p0.2, 
   mean = -1.5, sd = 2,    length = 0.2, lty = 2, lwd = 2)
 arrows(x.p0.2, 0.2, -9, 0.2, length = 0.2, 
    lty = 1, lwd = 2)
 arrows(x.p0.2,  pnorm( x.p0.2, mean = -1.5, sd = 2)  , -9,
   pnorm( x.p0.2, mean = -1.5, sd = 2), length = 0.2, 
    lty = 2, lwd = 2)
 text(x.p0.2, -.02, "q", cex = 1.1)
 text(-10,  pnorm( x.p0.2, mean = -1.5, sd = 2),
   expression(p [x] (q)), cex = 1.1)
 text(-10, 0.2, expression(p [y] (q)), cex = 1.1)
 text(-10,0, "0", cex = 1.2)
 text(-10,1, "1", cex = 1.2)
 ## text(-10, 1, pch = 64, cex = 1.2)# for test
 dev.off()

_ R の日本語化 2007 04 15

Rjp Wiki を参照

_ 正規分布の理論曲線の引きかた 2007 05 23

 x<- rnorm(1000, mean = 50, sd = 10)
 riron<- dnorm(sort(x), mean = 50, sd = 10) * sum(x) / sd(x)
 hist(x)
 curve(sort(x), riron)

他にリゲス,Rの基礎とプログラミング技法Crawley Statistics Using R p. 56 - 57.

_ 画像内の日本語フォント Illustrator 2006 06 06

Rjp Wiki より

_ pdf の作成 2007 06 06

 options(X11fonts = 
  c("-alias-mincho-%s-%s-*-*-%d-*-*-*-*-*-*-*",
     "-adobe-symbol-*-*-*-*-%d-*-*-*-*-*-*-*"))
 ps.options(family= "Japan1Ryumin") # これを指定しない場合,
 # 例えば Illustrator で文字化けする(開くことはできる)
 pdf("/home/ishida/test.pdf")
 plot(sin,-pi,pi,main="正弦")
 dev.off()

_ フォントの埋め込み embedFonts 2006 06 06

 options(X11fonts = 
 c("-alias-mincho-%s-%s-*-*-%d-*-*-*-*-*-*-*",
    "-adobe-symbol-*-*-*-*-%d-*-*-*-*-*-*-*"))
 ps.options(family= "Japan1Ryumin")
 plot(sin,-pi,pi,main="正弦")
 dev.copy2eps(file = "/home/ishida/test.eps")

このファイルは,例えばIllustrator では開けない.

 embedFonts(file = "/home/ishida/test.eps", 
   outfile =  "/home/ishida/test2.eps")

このファイルは,Illustrator では開けるが,文字はビットマップ化されて美しくない.

ちなみに

 options(X11fonts = 
 c("-alias-mincho-%s-%s-*-*-%d-*-*-*-*-*-*-*",
    "-adobe-symbol-*-*-*-*-%d-*-*-*-*-*-*-*"))
 # ps.options(family= "Japan1Ryumin")
 plot(sin,-pi,pi,main="正弦")
 dev.copy2eps(file="/home/ishida/test3.eps")

では日本語がそのまま入るが,例えばIllustratorでは開けない.

  ps.options(family= "Japan1Ryumin")

を実行すると,日本語部分はコード化されて記号となるが,Illustrator では開けない.

_ Ghostscript 用の環境変数 R_GSCMD 2007 06 12

 system("which gs")
 /usr/local/TeX/bin/gs
 R_GSCMD=${R_GSCMD='/usr/local/TeX/bin/gs'}
 Sys.putenv("R_GSCMD"="/usr/local/TeX/bin/gs")

_ CarbonEL for Mac Carbon Emacs 2007 06 12

How to install:

 install.packages("CarbonEL",,"http://rforge.net/")

to load:

 library(CarbonEL)

 Make sure CarbonEL package gets loaded by ESS such that
 Quartz is always responsive. This also works for R 
 run from a shell such that Terminal or xterm.

 Cheers,
 Simon

 PS: The above doesn't fix other problems with the old
 Quartz device (R.app uses a different, newer Quartz device!
 E.g. quartz.save(..) works only in R.app). It is just a
 hack for people that don't mind using the old Quartz
 device ;).

http://login.stat.ucla.edu/pipermail/statcompute/2003-November/000629.html

_ 文字ベクトルの一部を取り除く 2007 06 14

  writers[writers!="Kajii"]
  writers[c(writers!= "Kajii" & writers != "Koda")]
  writers[writers!= "Kajii" | writers != "Koda"] 

_ プロットを重ねる par new = T 2007 06 15

 par(new = T)

_ テーブルの区間ラベルの操作.行列名を取り出すなど 2007 06 16

 as.numeric(dimnames(table(buntyo))$buntyo)
 x<- c(1,1,1,1,2, 5,5,5, 6,6)
 table(cut(x, breaks = 0:max(x)))
 table(cut(x, breaks = 0:max(x), lab = F))
# 頻度表を作る
 as.numeric(labels(table(cut(x, breaks = 0:max(x),
    lab = F)))[[1]])
 table(cut(x, breaks = 0:max(x), lab = 0:max(x)))
# 頻度表を作る

_ file open close ファイルを一行ずつ読み込む 2007 07 10

 con<- file("test.txt")
 open(con)
 # 最初の 2 行は con<- file("file", "r") でも良い
 while(1){
   x<- readLines(con, n=1)
   if(length(x)> 0){
     print(x)
   }else{
     break;
   }
 }
 close(con)

_ UTF-8 で R で Rstem を実行する方法 2007 07 11

始めに

 [ishida@amd64 tmp]$ export LANG=ja_JP.UTF-8

ちなみにこの時,コンソールのメニュー設定も変更すること

コンソールから

 sessionInfo()  
 R version 2.5.1 (2007-06-27)
 x86_64-unknown-linux-gnu

 locale:
   LC_CTYPE=ja_JP.UTF-8;LC_NUMERIC=C;
   LC_TIME=ja_JP.UTF-8;LC_COLLATE=ja_JP.UTF-8;
   LC_MONETARY=ja_JP.UTF-8;LC_MESSAGES=ja_JP.UTF-8;
   LC_PAPER=ja_JP.UTF-8;LC_NAME=C;
   LC_ADDRESS=C;LC_TELEPHONE=C;
   LC_MEASUREMENT=ja_JP.UTF-8;LC_IDENTIFICATION=C

 kafka <- readLines("tmp/kafkaUTF8.txt")  
 library(Rstem)  
 kafka.vec<- unlist(strsplit(kafka[5], "[[:space:]]"))
  kafka.vec
  [1] ""           "Zu"         "diesem"     
 [16] "Gesetz."
 Warning message:
 Currently, only 'english' is tested. 
 You will need support for UTF characters. in:
 wordStem(kafka.vec, language = "german")

_ 正規表現 使用例 2007 07 11

日本語文字コードチェック perl での日本語処理は以下を参照いたしました. http://www.din.or.jp/~ohzaki/perl.htm#

 # or ~/progSource/Perl/perl.html regex.html
     # Perl で間違ってマッチしてしまう例
     $str = 'これはテストです';
     $pattern = '好';
     if ($str =~ /$pattern/) {
       print "マッチした\n";
     }

日本語文字コードチェック

 str.jp<- "これはテストです"
 length(str.jp)
 gsub("です","--", str.jp)  # "これはテスト--"
 gsub("^(\\w)","J", str.jp)#"Jれはテストです"
 regexpr("です",str.jp)
 # [1] 7
 # attr(,"match.length")
 # [1] 2
 gregexpr("です",str.jp)
 # [[1]]
 # [1] 7
 # attr(,"match.length")
 # [1] 2

 gsub("好","--", str.jp)#"これはテストです"
 regexpr("好",str.jp)
 # [1] -1
 # attr(,"match.length")
 #[1] -1
 gregexpr("好",str.jp)
 # [[1]]
 # [1] -1
 # attr(,"match.length")
 # [1] -1

茶筅の解析結果を取り込んで読み込むたびにスペースで区切る

 con open(con)
 while(1){
   if(length(x)> 0){
     for(j in 1:length(z)){
       print(z[j])
     }
   }else{
     break;
   }
 }
 close(con)

読み込むたびにスペースで区切る

 con<- file("chasen.txt")
 x<- readLines(con, n=1)
 close(con)
 x
 [1] "ここ\tココ\tここ\t名詞-代名詞-一般\t\t"
 y <- 
 gsub("^(.+)[[:space:]]+(.+)
   [[:space:]]+(.+)[[:space:]]+(.+)[[:space:]]+", 
   "\\1_\\2_\\3_\\4", x)
 y
 [1] "ここ_ココ_ここ_名詞-代名詞-一般\t"

 #  y<- 
gsub("^(.*)[[:space:]]+(.*)
  [[:space:]]+(.*)[[:space:]]+(.*)[[:space:]]+", 
  "\\1_\\2_\\3_\\4", x)

 #  y<- 
gsub("([[:alpha:]]*)
  [[:space:]]+([[:alpha:]]*)
  [[:space:]]+([[:alpha:]]*)
  [[:space:]]+([[:alpha:]]*)[[:space:]]+", 
  "\\1_\\2_\\3_\\4", x)
#  y<- gsub("(.*)\\t(.*)\\t(.*)\\t(.*)", 
"\\1_\\2_\\3_\\4", x, perl = T)

文字列の長さ

 mixed.text<- "花a子"
 nchar(mixed.text )                      # 5
 nchar(mixed.text, type = "chars")       # 3
 sub('[[:alnum:]]', '', mixed.text)      #  "a子"
 gsub('[[:alnum:]]', '', mixed.text)     #   ""
 gsub('[[:lower:]]', '', mixed.text)     #  "花子"
 gsub('[[:upper:]]', '', mixed.text)     #   "花a子"

_ 箱鬚図のノッチについて 2007 07 30

以下からの引用です!(石田) http://d.hatena.ne.jp/syou6162/20070702/1183349416 より

 ボックスプロットは非常に分かりやすいプロットであるため、
 使い易いツールである。
 しかし、medianが統計的に有意に違うかどうかなどを
 顧客に説明するには足りない部分がある。
 そこでボックスプロットを用いてmedianが
統計的に有意に異なるかを調べるためのツールが、
 このノッチの付いたボックスプロットである。
 ノッチのついたボックスプロットは以下のようにして
使うことができる。

 boxplot(lograifall.control,lograinfall.seeded,notch=TRUE)

 図にてへこんだ部分がノッチである。
 これがかぶっていなければ統計的に有意にmedianが
 違いということが言える。

 また、ノッチのついたボックスプロットはふたつの違いが
有意かどうかしか検定できないことに注意した。
 みっつ以上のノッチについては使えないということである
 (おそらくt検定だけではなく、F検定が必要な理由と同じと思う)。

_ モデル診断 model.diagnostics (cook, residuals) プロットの見方 2007 07 31

行動計量学会チュートリアルより 最後に紹介したプロットは,モデルの適合度を確認する補助となるプロットを作成した例である. 詳細は Faraway\citep[pp.14 -- 17]{faraway05}, Crawley\citep[pp.357 -- 359]{crawley07} Wood, p.33 参照されたい. ここでは左上のプロットは,残差と適合値を対応させたもので, データの変動になんらかの規則性がなければ, この散布図上で各ポイントはランダムに散っているはずである. この図の場合,体重の変化には何らかのシステマティックな要因が働いている, あるいは分散に変動があると思われる.左下の図は残渣を標準化し (平均 0,標準偏差 1 とする), その絶対値の平方根の取った図である(絶対値のままだと skew が生じるため).

後者の図では点が三角形の形になる場合,heteroscedascity が疑われる. この種の図が得られた場合は,例えば変数の二次の項を導入するなどの変換を検討することになる.

右上の図は QQ プロットと呼ばれ,残差が正規分布に従っているかをチェックするのに使われる. 最後の図はクックの距離で,モデルの当てはめに影響のあるはずれ値を検出するのに利用される. ここで点線がクックの距離を表し,0.5 を越えるデータは影響が少なくない, 1 を越える場合はかなり影響があると見なす.

それぞれのグラフで,問題となるデータ点にはデータ番号が振られる.