Site admin: motohiro PukiWiki 1.5.1 © 2001-2016 PukiWiki Development Team. Powered by PHP 7.4.33
HTML convert time: 0.003 sec.
個人メモを整理中です．不適切な記述が多々あるかと思います．お気づきの際は，ishida-mこの部分を"@"に変更下さいias.tokushima-u.ac.jp までご連絡ください． 個人メモを整理中です．不適切な記述が多々あるかと思います．お気づきの際は，ishida-m(この部分を"@"に変更下さい)ias.tokushima-u.ac.jp までご連絡ください． [[R_fromOldHtml1]] #contents * LaTeX でテーブル 1 [#h6be20ef] library(xtable) x<- matrix(1:24,4,6) mytable<- xtable(x) print(mytable) * LaTeX でテーブル 2 [#s319ea2e] [[RjpWikiを参照:http://www.okada.jp.org/RWiki/index.php?cmd=read&page=LaTeX&word=TeX%20#content_1_0]] * 行列のランクを求める。 [#s806caeb] qr(X)$rank # 舟尾 p.61 * ベクトルのノルム [#b92c3db0] 渋谷 柴田 p.37 norm if(p=="inf") return (max(abs(x))) sum(abs(x)^p)^(1/p) } これは sqrt(t(x) %*% x) に等しい。 * read.csv 読み込み時に変数名を付加する [#o54b5cb7] 例えば，以下のようにする場合 #x<- read.csv("anborgia.nr.csv",header=F,row.names=1);x # colnames(x)<- c("line");x 以下のように実行しても構わない。ただし元データが２列で，その一列目を行名と指定している。 その関係で，col.namesの引数は二ついる。 x<- read.csv("anborgia.nr.csv", header=F,row.names=1,col.names=c("","line1"));x * 生データから頻度表を作成し，適合度の検定を行なう [#zaf0a874] まず jena freq 1 15 2 10 3 21 4 20 5 30 ... kanzel freq 1 19 2 38 3 4 4 23 5 29 ... table( factor(cut(buntyo$V1, breaks=c(0, 5, 10, 15, 25, 35, 40))) ) などとする． ここでは 0 より大きく 5 以下，5 より大きく10以下，などを区間としている． としてデータを読み込み，これに頻度カテゴリを表す列を追加する。 jena$fac jena freq fac 1 15 (11,16] 2 10 (6,11] 3 21 (16,22] 4 20 (16,22] 5 30 (22,31] ... kanzel freq fac 1 19 (16,22] 2 38 (31,92] 3 4 (0,6] 4 23 (22,31] 5 29 (22,31] ここで table(kanzel$jena)とすれば，各カテゴリの頻度が求める。 さて，ここでこの二つの頻度について適合度の検定を行なうには， test chisq.test(test) とする。 * 図に凡例を付記する [#wa609321] legend(x軸, y軸, ラベルの指定, fill= 色の指定, pch = シンボル | lty = 線種) labs cols mypch legend(7, 10, labs, pch = mypch, col=cols) x,y 軸は locator(1)として調べれば良い。あるいは引数として与える。 [[John Fox:http//socserv.mcmaster.ca/jfox/]], p.138 legend(locator(1), c('high', 'low'), lty = c(1,3), lwd =c(3,3), pch = c(19,1)) [[Wood, Generalized Additive Models:http://www.amazon.co.jp/Generalized-Additive-Models-Introduction-Statistical/dp/1584884746/]] p.36 legend("topright", legend = levels(water$location), pch = c(1,2), bty = "n") S-Plus であれば，marks 。 * 重複を取る [#l940a327] 例:ishida/research/corpus/file.old/bak/file/writers19.R を参照 ここで unique は重複を取る。 legend(-4.3, 4.5, unique(writers.name), fill=(as.integer(unique(writers.name))), text.width=1) * 行列・配列の特定の次元を軸にあるベクトルを順に操作する、関数 sweep() [#x7d39075] x <- matrix(1:16, 4,4) x [,1] [,2] [,3] [,4] [1,] 1 5 9 13 [2,] 2 6 10 14 [3,] 3 7 11 15 [4,] 4 8 12 16 y= 1:4 # sweep されるベクトル sweep(x,1,y) # x の各列から y を引き去る(既定動作) [,1] [,2] [,3] [,4] [1,] 0 4 8 12 [2,] 0 4 8 12 [3,] 0 4 8 12 [4,] 0 4 8 12 sweep(x,2,y) # x の各行から y を引き去る(既定動作) [,1] [,2] [,3] [,4] [1,] 0 3 6 9 [2,] 1 4 7 10 [3,] 2 5 8 11 [4,] 3 6 9 12 sweep(x,1,y, "*") # x の各列に y を掛ける [,1] [,2] [,3] [,4] [1,] 1 5 9 13 [2,] 4 12 20 28 [3,] 9 21 33 45 [4,] 16 32 48 64 * updateの使い方 2005 /06 /04 [#l829ce4d] 例えば writers.lda table(writers$writer,predict(writers.lda)$class) と解析した結果と，und を取り除いた結果 writers.lda.und table(writers$writer, predict(writers.lda.und)$class) を比較する。なお以下はエラーとなる #writers.lda.step #anova(writers.lda , writers.lda.und) * 判別分析における交差妥当化 2005/06/04 [#ub22b09e] 例えば writers.lda #table(writers$writer,writers.lda$class) エラー table(writers$writer,predict(writers.lda)$class) #cross validation writers.lda.cv CV=T) table(writers$writer,writers.lda.cv$class) #table(writers$writer,predict(writers.lda.cv)$class) エラー * グラフに直線を引く 2005 / 06 / 14 [#j93a5c1d] plot(x ~ y) kaiseki.lm(x ~ y) abline(kaiseki.lm[1,1], kaiseki.lm[3,1]) として切片，傾きを指定しても，また abline(kaiseli.lm) とするのも可能。 * 分散分析の結果に summary.lm を実行 2005/06/15 [#ff0ccd41] Crawely p.163 分散分析で，平均値に差があるかを確認し，さらに，要因レベルの効果差を見るには， summary.lm(aov(ozone ~ garden)) などを実行する。 * quartile と quantile percentile の使い分け 2005/06/15 [#odf1c97d] Dalgaard p.58 //0.25, 0.5, 0.75 点の三つを three quartile といい， //これに min, max を含めてそれぞれを quantile という。 Crawely p.197 //50% - 75% quartile と範囲を指定 //75% quantile などと分位点を指定 // //quantile は，四分位点を一般化した p分位点，すなわち，100pパーセント点(0 * Boxplot での notch＝Ｔ オプション 2005/06/15 [#aee66ae1] [[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.167 //ノッチは，メディアンの中心線から，上下に同じ長さの裾を伸ばすが，それぞれの限界値で折り曲がることによって，幅を示す。 * treatment.control の見方 2005/06/16 [#i78e2ed7] [[Crawley Statistics Using R:http://www.amazon.co.jp/Statistics-Introduction-Michael-J-Crawley/dp/0470022981/]] p.174, [[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.287 treatment.control では，各要因の各レベルの効果を， その要因の切片(基軸)であるレベルと比較して判断する。 したがって，基軸以外のレベルの相互の有意差は表示されていない。 基軸以外のレベル相互の比較 要因Aと他要因 BCD全体との差，BCD 内の差を一度に検定 2006 03 20 [[John Fox:http://socserv.mcmaster.ca/jfox/]], p.142 - 143. 他に[[ここ>#ContrastsMatrix]] * loglin2 2*2 分割表の Likelihood ratio を求める [#u3b5389d] &aname(loglin2); 始めに ここを参照 こことここも またここも meyer p.132 の例より lessWords moreWords USA 57 2 Phil 71 12 NZ 66 16 GB 23 0 meyerP.132 dimnames(meyerP.132) model0 summary(model0) # anova(model0) # これは用意されていない model1 summary(model1) model1$lrt #P値は 1- pchisq(model1$lrt, model1$df) 1- pchisq(model1$pearson, model1$df) * コレスポンデンス分析 [#q882e6d5] R multiv パッケージのコレスポンデンス分析で出力される 行座標 $rproj, 列座標 $cproj はデフォルトでは，求められた x,y に単相関係数を乗じている。(言い方を変えると，固有値を平方根を乗じている) 参照 高橋信「Excelによるコレスポンデンス分析」, p.127 research/statistics/correspondence.R を参照 #高橋信「Excelによるコレスポンデンス分析」 Z 13,8,15,16, 18,11,14,8), ncol=4, byrow=T) #行ラベルとして使用する第1列をrowvarに代入 Z.rowvar<- c("M","H","C") #rowvarの値を行ラベルとして読み込む Z.colvar rownames(Z) colnames(Z) Z A B C D M 10 19 13 5 H 13 8 15 16 C 18 11 14 8 Z.diagC #> Z.diagC # [,1] [,2] [,3] [,4] #[1,] 41 0 0 0 #[2,] 0 38 0 0 #[3,] 0 0 42 0 #[4,] 0 0 0 29 #Z.diagR #> Z.diagR # [,1] [,2] [,3] #[1,] 47 0 0 #[2,] 0 52 0 #[3,] 0 0 51 #Z.diagR.sqrt # [,1] [,2] [,3] #[1,] 6.855655 0.000000 0.000000 #[2,] 0.000000 7.211103 0.000000 #[3,] 0.000000 0.000000 7.141428 # #Z.solve Z.eigen #$values #[1] 1.00000000 0.07629421 0.01828369 #$vectors # [,1] [,2] [,3] #[1,] -0.5597619 0.74374060 0.3653992 #[2,] -0.5887841 -0.66725782 0.4561802 #[3,] -0.5830952 -0.04021102 -0.8114081 # # [,1] [,2] [,3] # [1,] 0.3579767 0.2947642 0.3186919 # [2,] 0.2947642 0.3844402 0.3385965 # [3,] 0.3186919 0.3385965 0.3521610 # 固有値 １ と，対応する固有ベクトルは仮定 # 平均０，分散１に反する # p.121 # x, y が求まる，以下はxのみ Z.kai.X [,1] [,2] [,3] [1,] -1 1.32867325 -0.6527762 [2,] -1 -1.13328105 -0.7747835 [3,] -1 -0.06896133 1.3915534 # 以下は Y のみ Z.kai.Y # [,1] # [1,] -0.23728726 # [2,] 1.46910930 # [3,] -0.05964337 # [4,] -1.50318462 Z.kai.Y # [,1] # [1,] 1.5238380 # [2,] -0.6410599 # [3,] -0.1102451 # [4,] -1.1547168 # 固有値の平方根 は単相関係数 sqrt(Z.eigen$values) # [1] 1.0000000 0.2762141 0.1352172 # p.127 Z.kai[,2] * sqrt(Z.eigen$values)[2] # [1] 0.36699824 -0.31302816 -0.01904809 Z.kai[,3] * sqrt(Z.eigen$values)[3] # [1] -0.08826657 -0.10476405 0.18816195 #(solve(Z.diagR.sqrt, sqrt(150) * Z.eigen$vectors)) %o% sqrt(Z.eigen$values) library(multiv) Z.results> Z.results $evals [1] 0.07629421 0.01828369 $rproj Factor1 Factor2 [1,] 0.36699824 -0.08826657 [2,] -0.31302816 -0.10476405 [3,] -0.01904809 0.18816195 $cproj Factor1 Factor2 [1,] -0.06554208 0.20604911 [2,] 0.40578865 -0.08668233 [3,] -0.01647434 -0.01490704 [4,] -0.41520073 -0.15613757 $rcntr Factor1 Factor2 [1,] 0.553150080 0.1335166 [2,] 0.445232993 0.2081003 [3,] 0.001616926 0.6583831 $ccntr * 図の保存 eps 形式 2005 06 24 [#o1d0f24b] dev.copy2eps(file="log.eps") * 対数尤度比をコーパスで求める場合 2005 07 01 [#j65238fa] Dunning の式は，Oakes p.172 にまとめられている。必要となるのは a,b,c,d，と 当該コーパスにおける総語数 N (結局 a+b+c+d に等しい)。 a = あるキーワードにたいして設定したスパンの内部に，その共起語が現れた回数 b = キーワードにたいして，その共起語以外の語がスパンに現れた回数 c = キーワード以外の語にたいして，その共起語がスパンに現れた回数 d = キーワード以外の語にたいして，その共起語以外の語が現れた回数 * 東大「人文社会科学の統計学」，飽和モデル，独立モデル 2005 07 13 [#rc7b5e81] 始めに ここを参照 こことここも またここも &aname(social.mobility2); データは東大「人文社会科学の統計学」 p.277 social.mobility2<- array(c( 517, 244, 250, 152, 351, 328, 7, 13, 140, 441, 159, 272, 130, 227, 351, 24, 21, 268 ), dim = c(3,3, 2), dimnames = list( father = c("white","blue","agri" ), # = c(1) child = c( "white","blue","agri"), # = c(2) year= c("1985", "1975") # = c(3) ) ) class(social.mobility2 )<- "table" ftable(social.mobility2) library(MASS) # 以下は，3次の効果を外した同時独立モデルだが，ここに有意差は無い。すなわち year の効果は無い。 # 3次の効果を含めた飽和モデルは，自由度ゼロなので必ず適合する。p. 271 social.mobility2.loglm3<- loglm(~father + child + year + father:child + father:year + child:year, data=social.mobility2) anova(social.mobility2.loglm3, test="Chisq") 同じことだが，以下でもよい social.mobility2.loglm3<- loglm(social.mobility2, list(c(1), c(2), c(1,2), c(2,3), c(1,3))) 1 -pchisq(kekka$lrt,kekka$df)# 1 -pchisq(kekka$pearson,kekka$df)# # 東大「人文社会の統計学」 p.276 # log m = u... + ui.. + u.j. + u..k + uij. + ui.k + u.jk (9.33) # ここで uij. ui.k u.jk がそれぞれ，father:child, father:year, child:year # を表しているが # uijk (p.276下から12行目), すなわち father:child:year # は含まれていないモデル # なお u... は全体平均 grandmean は暗黙のうちに含まれる # ui.. 行平均 u.j. 列平均 u..k 軸平均は それぞれ要因 father, child, year # を指定して含まれる * 分散分析での母数モデル，変量モデル，混合モデル(豊田秀樹) 2005 07 27 [#h688e93f] S-Plus での混合モデルの構築 [[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.673 - ちなみに # p.678 lmedata3.model2 で。 random=~1|Genotype は Genotype の intercepts をランダム要因と見なす指定, p.678 # 豊田秀樹 違いを見抜く統計学 p.120, p.127, toyoda.Bunsan.R を参照 p120<- data.frame(type = rep(c("sea","mount"), each = 8), develop = rep(rep(c("old","new"), c(4,4)),2) , region = rep(c("izu","kuju","karui","yatsu"), each=4), price = c(18.0, 15.0, 12.0, 20.7, 7.6, 12.3, 15.0, 15.2, 11.7, 11.6, 14.5, 15.4, 7.8, 7.5, 4.2, 6.4)) #p120<- data.frame(type = rep(c("sea","mount"), each = 8), develop = rep(rep(c("old","new"), c(4,4)),2) , price = c(18.0, 15.0, 12.0, 20.7, 7.6, 12.3, 15.0, 15.2, 11.7, 11.6, 14.5, 15.4, 7.8, 7.5, 4.2, 6.4)) #> p120 # type develop region price #1 sea old izu 18.0 #2 sea old izu 15.0 #3 sea old izu 12.0 #4 sea old izu 20.7 #5 sea new kuju 7.6 #6 sea new kuju 12.3 #7 sea new kuju 15.0 #8 sea new kuju 15.2 #9 mount old karui 11.7 #10 mount old karui 11.6 #11 mount old karui 14.5 #12 mount old karui 15.4 #13 mount new yatsu 7.8 #14 mount new yatsu 7.5 #15 mount new yatsu 4.2 #16 mount new yatsu 6.4 attach(p120) 交互作用を含んだモデル，交互作用棄却後のプーリングモデル p120.aov<- aov(price ~ type * develop) summary(p120.aov) p.123と同じ出力# 豊田秀樹 違いを見抜く統計学 p.120, p.127 type develop の結果だけ見れば，プーリングした結果と同じ(p.123下) Df Sum Sq Mean Sq F value Pr(>F) type 1 84.181 84.181 10.1788 0.007769 ** develop 1 115.026 115.026 13.9084 0.002877 ** type:develop 1 8.556 8.556 1.0345 0.329170 Residuals 12 99.242 8.270 summary.lm(p.120.aov) > 84.181 / 8.270 # AB の平均平方ではなく，誤差の平均平方で割っている, p.108 [1] 10.17908 > 115.026 / 8.270 [1] 13.90883 > 1 - pf(10.1788, 1, 12) [1] 0.007769413 交互作用を変量としたモデル，反復測定 p121.aov = aov(price ~ type + develop + Error(type/develop)) summary(p121.aov) p.121 と同じ出力 母数因子だが，一つの標本から繰り返し測定した すなわち，交互作用は変量モデル Error: type:develop Df Sum Sq Mean Sq F value Pr(>F) type 1 84.181 84.181 9.8392 0.1965 develop 1 115.026 115.026 13.4444 0.1695 Residuals 1 8.556 8.556 Error: Within Df Sum Sq Mean Sq F value Pr(>F) Residuals 12 99.242 8.270 ネストを仮定したモデル p127.aov summary(p127.aov) p.127 と同じ出力 すなわち，地域はタイプにネストしている Error: region:type Df Sum Sq Mean Sq F value Pr(>F) type 1 84.181 84.181 1.3624 0.3635 Residuals 2 123.581 61.791 Error: Within Df Sum Sq Mean Sq F value Pr(>F) Residuals 12 99.242 8.270 ここで > 61.791 / 8.270 > 1 - pf(7.471705 , 2, 12) [1] 0.007804986 三段のネストを仮定したモデル これは上の２段ネストと同じ結果になっている p.157 とは異なり，残差が二つに分離されていないx p.157.aov <- summary( p.157.aov ) Error: region:type Df Sum Sq Mean Sq F value Pr(>F) type 1 84.181 84.181 1.3624 0.3635 Residuals 2 123.581 61.791 Error: Within Df Sum Sq Mean Sq F value Pr(>F) Residuals 12 99.242 8.270 また crawley, S-Plus, p.213下を参照。 p120.aov # = aov(price ~ type/develop) # = aov(price ~ type/develop + Error( type/develop)) summary(p120.aov) # # p.127 と同じ出力 # # Df Sum Sq Mean Sq F value Pr(>F) #type 1 84.181 84.181 10.1788 0.007769 ** #type:region 2 123.581 61.791 7.4715 0.007806 ** #Residuals 12 99.242 8.270 #anova(p120.aov) #Error # # Df Sum Sq Mean Sq F value Pr(>F) #type 1 84.181 84.181 10.1788 0.007769 ** #develop 1 115.026 115.026 13.9084 0.002877 ** #type:develop 1 8.556 8.556 1.0345 0.329170 #Residuals 12 99.242 8.270 #summary.lm(p120.aov) #> 84.181 / 8.270 # AB の平均平方ではなく，誤差の平均平方で割っている, p.108 #[1] 10.17908 #> 115.026 / 8.270 #[1] 13.90883 #> 1 - pf(10.1788, 1, 12) #[1] 0.007769413 # エラ−モデル 1 p120.aov summary(p120.aov) # # p.121 と同じ出力 # #Error: type:develop # Df Sum Sq Mean Sq F value Pr(>F) #type 1 84.181 84.181 9.8392 0.1965 #develop 1 115.026 115.026 13.4444 0.1695 #Residuals 1 8.556 8.556 #Error: Within # Df Sum Sq Mean Sq F value Pr(>F) #Residuals 12 99.242 8.270 # エラ−モデル 2 p120.aov summary(p120.aov) # # p.127 と同じ出力 # #Error: region:type # Df Sum Sq Mean Sq F value Pr(>F) #type 1 84.181 84.181 1.3624 0.3635 #Residuals 2 123.581 61.791 #Error: Within # Df Sum Sq Mean Sq F value Pr(>F) #Residuals 12 99.242 8.270 # ネストモデル p120.aov summary(p120.aov) # # p.127と同じ出力 # # Df Sum Sq Mean Sq F value Pr(>F) #type 1 84.181 84.181 10.1788 0.007769 ** #type:region 2 123.581 61.791 7.4715 0.007806 ** #Residuals 12 99.242 8.270 #anova(p120.aov) #Error * 対数尤度を求める関数 2005 08 01 [#o1feec7a] logLik, (crawley, S-Plus, p.121) * リストへのアクセス 2005 08 03 [#r8e750c9] test Error: factory Df Sum Sq Mean Sq F value Pr(>F) Residuals 6 260.70 43.45 Error: Within Df Sum Sq Mean Sq F value Pr(>F) treats 4 8.324 2.081 0.4637 0.7617 Residuals 24 107.704 4.488 以上のリストで，Error の F value にアクセスするには summary(aov(rate ~ treats + Error(factory)))[[2]][[1]]$"F value"[1] と実行する。 Crawely.S.plus.R を参照 * グラフィックス表示を一枚ずつ見る 2005 08 05 [#d0392771] par(ask=TRUE) * 共分散分析, Ancova モデル式におけるアスタリスク 2005 08 08, 2006 01 11 [#hc1a7d93] difference in the slop は 交互作用の項を見る。 [[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.295 covariate については p.293 [[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.216, p.341 w が sex カテゴリーの変数，x が連続量の変数の場合 y ~ w * x は，w の二つのカテゴリそれぞれに回帰式(intecept と slope) regression slope を推定しろという命令となる。 [[John Fox:http://socserv.mcmaster.ca/jfox/]] p.126 // S-plus * [S] EES の評価途中にスペースがある [#yb52040f] 式に改行やスペースが挟間っていると，そこで処理が止まってしまう。 * [S] ESS で操作中で画像の選択がある [#t7d79e80] 例えば plot を含めて，式を一括評価する時 par(ask=T) と設定した後だと，yes か no の答えを得ることができないので，そこで止まってします。 * データの並び替え [#q9c7bb32] [[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.417のデータならば > fishery xv yv 1 13.9885004 4.56836585 2 89.7028342 0.29710071 3 9.3151971 3.90699292 4 9.2132163 4.38109971 5 34.7017719 1.85936367 のようなデータフレームならば fishery[order(fishery$xv),] とすれば，xv の昇順にデータが整列する。 * [S] 内積とノルム 2005 08 23 [#f3b0090e] a b 内積は t(a) %*% b としても sum(a * b) としても良い。 ノルムは vecnorm(a); vecnorm(b) #で求めるか # mynorm if(p=="inf") return(max(abs(x))) sum(abs(x)^p)^(1/p)) } で求める。 * adj の設定 2005 08 25 [#ub65cb39] research/tex/p2005/2005/dunne.R より text(E.normalized.B.dunne2[1, E.normalized.B.dunne2[2,] > 0.0], E.normalized.B.dunne2[2, E.normalized.B.dunne2[2,] > 0.0], adj=c(1.5,-.5), labels=( colnames( E.normalized.B.dunne2[, E.normalized.B.dunne2[2,] adj=c(X,Y) X を正にすると，点は左に，Yを正にすると点は下に移動 * 負の2項分布 2005 08 30 [#w1336322] [[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.483 例えば 成功率 0.6, 成功までの試行数が 5 として(すなわち６回目に成功)， 100回の試行を行なう。 count table(count) 0 1 2 3 4 5 6 7 8 9 8 13 20 19 18 7 9 3 2 1 この場合 失敗無しで５回の成功が出た試行が 8 回 失敗 1 で５回の成功が出た試行が 13 回 あるのを意味する。 # 2006 07 24 追加 工学のためのデータサイエンス, p.49 - p.50 R の関数 _nbinom の引数 size, prob は 平均 mean が 'size * (1 - prob)/prob' を満たし，また 確率 prob が 'scale/(1 + scale)' を満たすように取られる d <- p(0, 1612), rep(1,164), rep(2,71), rep(3,47), rep(4,28), rep(5,17), rep(6,12), rep(7,12), rep(8,5), rep(9,7), rep(10,6), rep(11,3), rep(12,3), rep(13, 13)) mean(d) # = 0.6005 var(d) # = 3.262531 prob = 0.1155 * (1 - 0.1553) / 0.1553 # = 0.6282218 = size *(1 - prob)/prob scale = 0.1838523 / ( 1 + 0.1838523) # = 0.1553 = scale/ (scale + 0.1553) * lm, glm オブジェクトの summary における t 検定,z 検定 2005 09 02 [#t54cd3ab] S-Plus では t 値が，R では z 値がデフォルト [[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.534 [[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] S-plus.R, p.534 また S-Plus では，lm の出力は Pr(|t|) を出力するが，glm ではデフォルトでは出力されない。 * rgl.bbox の引数 2005 09 02 [#s091d17e] rgl.bbox(color="#999999", emission="#FFA500", specular="#FFF0F5", shinines=5,alpha=0.8) color="#999999" は座標の数値ラベル色 emission="#FFA500" はボックスの壁の色 specular="#FFF0F5" は明るさ，数値が大きいほど明るい shinines=5 は反射の具合 * 関数グラフの書き方 2005 09 19 [#q2e33709] xyplot(x,y,type="l") 同じことだが curve((x^2 - 6*x + 10), 0,6) * 2項分布を Poisson 分布で近似する 2005 09 28 [#h5577017] 例えば2項分布で，0.00001 で成功する賭けを 100000 回行なった時， ちょうど２０回成功する確率は pbinom(20, 100000, .0001) - pbinom(19, 100000, .0001) [1] 0.001865335 これをポアソン分布で近似するには ppois(20, 10) - ppois(19, 10) [1] 0.001866081 * データフレームのコピー 2005 10 08 [#s5f3e4ac] [[John Fox:http://socserv.mcmaster.ca/jfox/]] p.63 - 64 > search() [1] ".GlobalEnv" "package:car" "package:methods" [4] "package:stats" "package:graphics" "package:grDevices" [7] "package:utils" "package:datasets" "Autoloads" 10] "package:base" > cooperation エラー：オブジェクト "cooperation" は存在しません > attach(Guyer) # Guyerデータには変数 operationがある > objects() [1] "Guyer" "hump.data" "hump.model" "last.warning" > cooperation [1] 49 64 37 52 68 54 61 79 64 29 27 58 52 41 30 40 39 44 34 44 > cooperation > objects()新たに変数 cooperation が作成されている [1] "Guyer" "cooperation" "hump.data" "hump.model" > rm(cooperation)# 新たに作成された cooperation を消し > cooperation[1] Guyer変数のcooperation[1]を変更したつもりでも > objects()実際にはコピーができている [1] "Guyer" "cooperation" "hump.data" "hump.model" > rm(cooperation) > objects() [1] "Guyer" "hump.data" "hump.model" "last.warning" * png などの画像への保存 2005 11 09 [#b2a9096e] 始めに png(file="gazo.png") #として device を平き polt(y ~ x) # として画像を作成し dev.off() # として device を閉じる。 ここで始めて保存される。 * 標準誤差 2005 12 02 [#zdcdd38c] 標準誤差とは，平均値の推定値の範囲 標準偏差とは，データが平均を中心にどれだけばらついているかの目安 * Boxplot の notch の解釈 2005 12 23 [#mcf9c930] Boxplot の箱は 75 - 25 %タイルを示しているが，あるデータのこの範囲が，他のデータの平均値と重なっていない場合，二つのmedian は有意に異なる。Tukey による [[Crawley Statistics Using R:http://www.amazon.co.jp/Statistics-Introduction-Michael-J-Crawley/dp/0470022981/]] p.167 * [S] 関数内で式を表示させない 2005 12 24 [#e11ee55f] 関数内では，付置の式以外は表示される。 特に print 文では注意が必要である。 これを防ぐため，invisibleを使う。 [[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.35 渋谷，柴田 p.34 * [S] S-plus + ESS 環境でのヘルプ 2005 12 24 [#d40dc98b] キーワードにカーソルを合わせ，C-c C-v で表示させる。 ?キーワードでは表示されない // <<- はグローバル変数への代入 // // 舟尾 p.183