Rのメモ
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
-[[英語学研究2008-2]] | [[英語学研究2009-2]] | [[英語学研...
#contents
**ヘルプ参照 [#z5ba9ed0]
> help.start() # 全体のヘルプを出す場合
> help(rnorm) # function のヘルプを参照する場合、た...
**GUIの言語を変えたい [#pcb26f21]
-Rの起動オプション内で指定するのが楽な模様です。
-起動オプションを付与するには、Rを起動するショートカット...
> LANGUAGE=en #英語で起動したい場合
> "C:\Program Files\R\R-3.0.0\bin\x64\Rgui.exe" LANGUAGE...
**Function と argument [#g5e0231d]
-log(x) log() が function, x が argument
-plot(x,y) plot() が function, x, y が argument
**Vector ベクトル [#gffc451f]
-Numerical vector: 数値のベクトル
-Character vector: 文字列のベクトル
> wordlist <- c("able", "ability", "abuse")
> wordlist
[1] "able" "ability" "abuse"
**パッケージの読み込み [#kea71742]
-コマンドラインでライブラリを読み込む場合:library(xxxx) ...
> library(survival)
要求されたパッケージ splines をロード中です
次のパッケージを付け加えます: 'survival'
The following object(s) are masked from package:I...
lung
-library を閉じる場合:
> detach("package:survival")
-利用可能なlibrary を参照する場合:
> library()
**テキストファイルからの読み込み [#o74de443]
-(例)c:\temp にある lemma.num を bnclemma というデータ...
> bnclemma <- read.table("c:/temp/lemma.num", header=T)
header=T で1行目をヘッダとして読み込む
**テキストファイルへの書き出し [#k6a104af]
>write.table (データフレーム名, file="ファイル名")
**データエディタの起動 [#d9e3bbd5]
-作成したデータフレームをスプレッドシートで編集したい場合:
> lm <- edit(bnclemma)
--編集したデータフレームは lm という名前(何でもよい)に...
--bnelemma という元のフレームはそのまま
-上書きで編集するには fix() とする
> fix(bnclemma)
-まったく新しいファイルを書くには
>dd <- data.frame()
>fix(dd)
**ワークスペースの履歴・保存 [#i65b2a7d]
-ls() 関数
-今までの実行したワークスペースのコマンド一覧を出す
> ls()
[1] "bmi" "d" "energy" "exp.lean...
[7] "h" "height" "hh" "intake.p...
[13] "mylist" "oops" "pain" "thue2" ...
[19] "x" "y" "ylim"
-ワークスペースをファイルに保存する
>save.image()
**データフレームの部分加工 [#n225bbb0]
>icci2 <-icci[c(1,2)]
#icci のデータフレームの第1、2コラムを取り出して icci2 ...
**データフレームのサブセットの取り出し [#h0326953]
> icci_A1 <- subset(icci, CEFR == "A1")
#subset()関数で、データ icci の CEFRの変数の値が A1 のも...
**データフレームのマージ [#g73594ce]
>icci <- merge(icci1, icci2)
# これでお互いのファイルに共通する変数名を参照しながらマ...
# もしそれぞれにない個体があった場合に、それをすべて含め...
>icci <- merge(icci1, icci2, all=TRUE)
**データフレームのマージ(行の結合) [#aad41245]
> total <- rbind(A1_sample, A2_sample, B1_sample)
# rbind で2つ以上のデータフレームをくっつける。カラム変...
**データフレームからランダムに X 件を抽出 [#z177b7b2]
> A1_sample <- icci_A1[sample(1:nrow(icci_A1), 1000, rep...
# icci_A1 の行データ(1から個体数すべて)のなかから1000...
**分布 [#w2f350a6]
-Density or point probability 密度確率・
--the probability density function 確率密度関数
--the probability distribution function 確率分布関数
-確率分布関数=確率密度関数を積分したもの
-確率密度関数=確率分布を微分したもの
--統計データが十分にたくさんあれば、
--ヒストグラム = 確率密度関数になり
--累積頻度分布 = 確率分布関数になる。
***正規曲線の書き方1 [#a43c63aa]
> x <- seq(-4,4,0.1)
-4 から 4 まで 0.1ずつの感覚で sequence を作って x に格納
> plot(x, dnorm(x),type="l")
X軸に x を Y 軸に x にそって確率密度の正規曲線を生成
type="l" は line かな
> plot(x, dnorm(x),type="h")
同じプロットを h にすると histogram のようにそれぞれの x...
***正規曲線の書き方2 [#dbcc2dff]
> curve(dnorm(x), from=-4, to=4)
これでも同じカーブを描くことができる
***累積分布関数 [#wd8f1286]
-正規分布のどの位置にあるか、その位置までの cumulative di...
--平均132、標準偏差が13の正規分布の中で、160とい...
---正規分布上の値のとる確率は何パーセントか、を pnorm が...
---それを 1 から引くことで、160以上の値をとる確率を得ている
> 1-pnorm(160,mean=132,sd=13)
[1] 0.01562612
> 1-pbinom(15,size=20,prob=.5)
[1] 0.005908966
**記述統計 [#kaed79a2]
-サンプルデータとしてランダムに50個の正規分布を成す数値...
> x <- rnorm(50)
> x
[1] 1.11403983 -1.15165895 0.80461867 1.71557037 0....
[8] -0.10692887 -0.96407973 1.16309160 0.61545475 0....
[15] 0.12446743 0.77148385 1.19391658 -0.78081663 -0....
[22] -0.15209363 -0.49532705 -0.52829549 0.64414085 1....
[29] 1.55049255 0.21063906 0.22529564 0.34462737 0....
[36] 1.15710520 -0.69795906 -1.12325445 -2.08717916 0....
[43] 1.04982555 0.46580345 0.55212419 -0.53356749 -0....
[50] -0.06345214
-個別の記述統計の出し方
--平均
> mean(x)
[1] 0.02739456
-欠損値があって飛ばしたい場合の処理
> library(ISwR) #ISwR というライブラリを読み込む...
> data(juul) #data juul を読み込む
> juul
age menarche sex igf1 tanner testvol
1 NA NA NA 90 NA NA
2 NA NA NA 88 NA NA
3 NA NA NA 164 NA NA
4 NA NA NA 166 NA NA
5 NA NA NA 131 NA NA
--juul のデータの igf1 を参照したいが、いちいち juul$igf1...
> attach(juul) #こうすると juul$…を省略できる
> mean(igf1)
[1] NA #NA は欠損値があるので平均値を出せない、とエ...
> mean(igf1,na.rm=T) #"na.rm=T" not available, remove ...
[1] 340.168
--標準偏差
> sd(x)
[1] 0.876461
--分散
> var(x)
[1] 0.768184
--中央値(メディアン)
> median(x)
[1] 0.06215406
--分位数(quantile)
> quantile(x)
0% 25% 50% 75% 10...
-2.08717916 -0.68625583 0.06215406 0.63696932 1.71557...
---デフォルトではこのように正規分布の下から 25,50, 75,100...
---もし quantile を特定の観察したい割合で切りたい場合には...
> pvec <- seq(0,1,0.1)
pvec という変数に 0 から 1 までを 0.1 区切りに格納
> pvec
[1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
> quantile(x,pvec)
0% 10% 20% 30% 40% ...
-2.08717916 -0.97999720 -0.77189339 -0.52987709 -0.124...
70% 80% 90% 100%
0.49169967 0.77811082 1.16617410 1.71557037
---quantile に2つ目の argument で何%で切るかを 0.00-1.0...
**記述統計を一発で出す方法 [#fc2b1dfa]
-summary 関数は 最小値、第1分位点、メジアン、平均、第3...
> summary(igf1)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
25.0 202.2 313.5 340.2 462.8 915.0 321.0
-データフレーム全体の summary も可能
> summary(juul)
age menarche sex ...
Min. : 0.170 Min. : 1.000 Min. :1.000 Min. ...
1st Qu.: 9.053 1st Qu.: 1.000 1st Qu.:1.000 1st Q...
Median :12.560 Median : 1.000 Median :2.000 Media...
Mean :15.095 Mean : 1.476 Mean :1.534 Mean ...
3rd Qu.:16.855 3rd Qu.: 2.000 3rd Qu.:2.000 3rd Q...
Max. :83.000 Max. : 2.000 Max. :2.000 Max. ...
NA's : 5.000 NA's :635.000 NA's :5.000 NA's ...
testvol
Min. : 1.000
1st Qu.: 1.000
Median : 3.000
Mean : 7.896
3rd Qu.: 15.000
Max. : 30.000
NA's :859.000
**ヒストグラムを書く [#xb4b861d]
> x <- rnorm(100) # normal density of pseudo-ran...
> x
[1] -0.23852918 0.80007985 2.59488947 1.92574080 0....
[8] -0.07618785 -1.59706302 0.82031697 -1.20817887 1....
[15] 0.75213738 -0.06597683 0.23352029 -0.16270841 -1....
[22] -0.12607458 0.09170498 -0.10259290 -0.45916474 0....
[29] 1.45689028 0.64155579 -0.33764720 -1.25446597 0....
[36] -0.26736349 -0.97858879 1.45124404 -1.22135682 0....
[43] -0.45558029 1.61702489 -1.51530741 1.12869612 1....
[50] -0.73521276 -0.12948352 -0.09378984 0.93129898 0....
[57] 0.25044045 1.01388932 0.82988086 -1.59370935 -0....
[64] -0.38521575 1.30834830 1.20246317 -0.03461022 -0....
[71] 0.99960792 0.87056316 1.17321910 0.02667156 0....
[78] -0.85243182 -0.87548459 0.20321376 -1.46930772 1....
[85] -0.68914995 1.06938790 1.29973824 0.65570765 0....
[92] -0.25530870 0.46701640 1.27102185 1.73760250 1....
[99] -0.32908069 -0.49896537
> hist(x, freq=F) # write histog...
> curve(dnorm(x),add=T) # write the no...
-縦軸の目盛りが切れてしまった場合の処理(少し複雑)
> h <- hist(x, plot=F)
> ylim <- range(0, h$density, dnorm(0))
> hist(x, freq=F, ylim=ylim)
> curve(dnorm(x), add=T)
**変数:数値と文字列の操作 [#a703d170]
-たとえば energy というデータは expend, stature という2...
> data(energy)
> energy
expend stature
1 9.21 obese
2 7.53 lean
3 7.48 lean
4 8.08 lean
5 8.09 lean
6 10.15 lean
7 8.40 lean
8 10.88 lean
9 6.13 lean
10 7.90 lean
11 11.51 obese
12 12.79 obese
13 7.05 lean
14 11.85 obese
15 9.97 obese
16 7.48 lean
17 8.79 obese
18 9.69 obese
19 9.68 obese
20 7.58 lean
21 9.19 obese
22 8.11 lean
-lean の1行だけとりだす
> exp.lean <- energy$expend[energy$stature=="lean"]
-obese の行だけとりだす
> exp.obese <- energy$expend[energy$stature=="obese"]
> exp.lean
[1] 7.53 7.48 8.08 8.09 10.15 8.40 10.88 6.13 7.9...
> exp.obese
[1] 9.21 11.51 12.79 11.85 9.97 8.79 9.69 9.68 9.19
-同様のことが以下の split()関数でもできる
> l <-split(energy$expend, energy$stature)
>
> l
$lean
[1] 7.53 7.48 8.08 8.09 10.15 8.40 10.88 6.13 7.9...
$obese
[1] 9.21 11.51 12.79 11.85 9.97 8.79 9.69 9.68 9.19
**データ・フレームの一部取り出し [#h2e2b342]
(例)verbs というデータから AnimacyOfTheme が animate ...
> verbs[verbs$AnimacyOfTheme == "animate",]
RealizationOfRec Verb AnimacyOfRec AnimacyOfTheme L...
58 NP give animate animate ...
100 NP give animate animate ...
143 NP give inanimate animate ...
390 NP lend animate animate ...
506 NP give animate animate ...
736 PP trade animate animate ...
> subset(verbs, AnimacyOfTheme == "animate")
RealizationOfRec Verb AnimacyOfRec AnimacyOfTheme L...
58 NP give animate animate ...
100 NP give animate animate ...
143 NP give inanimate animate ...
390 NP lend animate animate ...
506 NP give animate animate ...
736 PP trade animate animate ...
**標準化得点を出す方法 [#rb5c6194]
> x=rnorm(10,3,0.1)
> x
[1] 2.987822 3.087673 2.972487 2.901589 3.056770 3.07972...
> x-mean(x)
[1] -0.01775045 0.08210025 -0.03308569 -0.10398376 0.0...
> y = x-mean(x)
> y/sd(x)
[1] -0.2184789 1.0105193 -0.4072305 -1.2798693 0.63015...
> scale(x)
[,1]
[1,] -0.2184789
[2,] 1.0105193
[3,] -0.4072305
[4,] -1.2798693
[5,] 0.6301557
[6,] 0.9126354
[7,] -0.1403719
[8,] -0.2493192
[9,] -1.6753157
[10,] 1.4172752
attr(,"scaled:center")
[1] 3.005573
attr(,"scaled:scale")
[1] 0.08124561
**さまざまなプロット [#o8ba5c7c]
-棒グラフ
--ratings 中の1つの変数 Length に関して crosstab を作り...
>barplot(xtabs( ~ ratings$Length), xlab="word length", c...
-ヒストグラム
--ratings 中の1つの変数 Length に関して、ヒストグラムを...
--この truehist というコマンドは MASS というパッケージ中...
>truehist(ratings$Length, xlab="word length", col="grey")
-一般的な描画関数 plot()
>plot(h)
>plotsort(lexdec$RT), ylab = "log RT") # RT という変数...
-箱ひげ図を書く
>boxplot(lexdec$RT) # ふつうに RT という変数を箱...
>boxplot(exp(lexdec$RT)) # RT は対数なのでそれをふつう...
-すべての変数間のプロットを一発で出す → scatterplot matri...
>pairs(rating[ , -c(1, 6:8, 10:14)])
-カテゴリーデータをモザイク図にする
>mosaicplot(verbs.xtabs, main = "dative")
--上位階層にしたいものを先に xtab で処理しておく
-上位階層を変数でくくって箱髭図を書く
**複数の表を並べて作成する場合 [#cb1be9cd]
-par()関数を使うと、以後作成する表を1つ1つ並べて表示し...
>par(mfrow = c(3,2))
--2列3行で表を並べてくれる
**表を画像にはき出す場合 [#rcdf67f6]
-jpeg で barplot.jpeg というファイルに書き出す
>jpeg("barplot.jpeg", width = 400, height=420)
>truehist ... # ここで作表を行う
>dev.off() # 画像モードを抜ける
**単語リストの操作 [#v2bfd808]
-テキスト・ファイルを読み込む
(例)c:\temp にある lemma.num を bnclemma というデータ...
> bnclemma <- read.table("c:/temp/lemma.num", header=T)
header=T で1行目をヘッダとして読み込む
-データフレームを編集する
--作成したデータフレームをスプレッドシートで編集したい場...
> lm <- edit(bnclemma)
---編集したデータフレームは lm という名前(何でもよい)に...
---lemmalist として lm というデータフレームを作成した
-単語リストの抽出
--データフレームを行列として参照する
---ランク10位の単語([,3]はwordのカラムを示す)を参照
> lm[10,3]
[1] to
---ランク10位の品詞([,4])を参照
> lm[10,4]
[1] prep
Levels: a adv conj det infinitive-marker interjection mo...
※このようにアウトプットの最後に、その変数に何レベルある...
---lm 中の word のカラムの最初の4個を取り出す
> lm$word[1:4]
[1] the be of and
5464 Levels: a abandon abbey ability able abnormal aboli...
---lm 中の pos のカラムの1番目を取り出す
> lm$pos[1]
[1] det
Levels: a adv conj det infinitive-marker interjection mo...
---lm 中のwordの項目数をカウントする → length()を使用
> length(lm$word)
[1] 6318
---rank が100位以内の word を取り出す
> lm$word[lm$rank < 100]
[1] the be of and a in to ...
[12] i that you he on with do ...
[23] but from they his that she or ...
[34] say will would can if their go ...
[45] her make who as out up see...
[56] some could so him year into its...
[67] than more about now last your me ...
[78] should these people also well any onl...
[89] way look like use her such how...
---rank が100位以内で、pos が名詞のものを取り出す
> lm$word[lm$rank < 100 & lm$pos == "n"]
[1] time year people way
※ [ ] で条件設定し、& で2つの条件を組み合わせている
※ logical operator としては、&, |(=or), !(=not) が使える
---rank が100位以内で、pos が名詞または形容詞のものを...
> lm$word[lm$rank < 100 & (lm$pos == "n"|lm$pos == "a")]
[1] time year last other people new way goo...
※ logical operator | が使われていることに注意
---10位以内のリストをすべてとりだす
> lm[lm$rank <10,]
rank freq word pos
1 1 6187267 the det
2 2 4239632 be v
3 3 3093444 of prep
4 4 2687863 and conj
5 5 2186369 a det
6 6 1924315 in prep
7 7 1620850 to infinitive-marker
8 8 1375636 have v
9 9 1090186 it pron
※10位まで取り出すならば、<11 とする
※[]の中で条件を付けるが、最後の ,] とカンマをつけるとこ...
---100位以内で品詞が動詞である行だけをとりだす
> lm[(lm$rank <101 & lm$pos == "v"),]
rank freq word pos
2 2 4239632 be v
8 8 1375636 have v
18 18 559596 do v
34 34 333518 say v
40 40 249540 go v
44 44 220940 get v
46 46 217268 make v
51 51 191661 see v
52 52 185534 know v
54 54 179220 take v
64 64 153881 think v
66 66 151871 come v
76 76 131417 give v
90 90 111058 look v
92 92 108820 use v
100 100 98899 find v
-100位以内の単語の頻度プロットを作る
> lm100 <- lm[(lm$rank <101),] # lm100 というデータフレ...
> lm100
rank freq word pos
1 1 6187267 the det
2 2 4239632 be v
3 3 3093444 of prep
4 4 2687863 and conj
5 5 2186369 a det
6 6 1924315 in prep
7 7 1620850 to infinitive-marker
8 8 1375636 have v
9 9 1090186 it pron
10 10 1039323 to prep
11 11 887877 for prep
12 12 884599 i pron
13 13 760399 that conj
14 14 695498 you pron
15 15 681255 he pron
16 16 680739 on prep
17 17 675027 with prep
18 18 559596 do v
19 19 534162 at prep
20 20 517171 by prep
※出力はこんな感じ
プロットの関数は plot()
> plot(lm100$rank, log(lm100$freq), type="n")
※ランキングを横軸、頻度の自然対数を縦軸にとる。
※type="n" はプロットの記号を消す、という意味。後で単語...
> text(lm100$rank, log(lm100$freq), as.character(lm100$w...
※この text() という関数が交差する点に as.character() ...
※cex=で文字サイズを指定する
--頻度リスト&プロット:応用 [#c0377bc4]
---100語以内の動詞のプロットを作る
> vlm <- lm100[lm100$pos == "v",] #動詞を100位以内か...
> vlm
rank freq word pos
2 2 4239632 be v
8 8 1375636 have v
18 18 559596 do v
34 34 333518 say v
40 40 249540 go v
44 44 220940 get v
46 46 217268 make v
51 51 191661 see v
52 52 185534 know v
54 54 179220 take v
64 64 153881 think v
66 66 151871 come v
76 76 131417 give v
90 90 111058 look v
92 92 108820 use v
100 100 98899 find v
> plot(vlm$rank, log(vlm$freq), type="n") #横軸ラン...
> text(vlm$rank, log(vlm$freq), as.character(vlm$word), ...
※vlm$rank と vlm$freq の交差する点に vlm$word から取って...
※cex=0.9 で少し大きめに
‐サブコーパス間のランク差の統計情報を得る
# 複数の頻度表を、見出し語・ランク・頻度・コーパスの4コ...
sum <- numSummary(brownclone[,"Rank"], groups=brownclone...
statistics=c("mean", "sd", "quantiles"), quantiles=c(0,...
write.csv(sum$table, file="numsummary_brownclone.csv")
終了行:
-[[英語学研究2008-2]] | [[英語学研究2009-2]] | [[英語学研...
#contents
**ヘルプ参照 [#z5ba9ed0]
> help.start() # 全体のヘルプを出す場合
> help(rnorm) # function のヘルプを参照する場合、た...
**GUIの言語を変えたい [#pcb26f21]
-Rの起動オプション内で指定するのが楽な模様です。
-起動オプションを付与するには、Rを起動するショートカット...
> LANGUAGE=en #英語で起動したい場合
> "C:\Program Files\R\R-3.0.0\bin\x64\Rgui.exe" LANGUAGE...
**Function と argument [#g5e0231d]
-log(x) log() が function, x が argument
-plot(x,y) plot() が function, x, y が argument
**Vector ベクトル [#gffc451f]
-Numerical vector: 数値のベクトル
-Character vector: 文字列のベクトル
> wordlist <- c("able", "ability", "abuse")
> wordlist
[1] "able" "ability" "abuse"
**パッケージの読み込み [#kea71742]
-コマンドラインでライブラリを読み込む場合:library(xxxx) ...
> library(survival)
要求されたパッケージ splines をロード中です
次のパッケージを付け加えます: 'survival'
The following object(s) are masked from package:I...
lung
-library を閉じる場合:
> detach("package:survival")
-利用可能なlibrary を参照する場合:
> library()
**テキストファイルからの読み込み [#o74de443]
-(例)c:\temp にある lemma.num を bnclemma というデータ...
> bnclemma <- read.table("c:/temp/lemma.num", header=T)
header=T で1行目をヘッダとして読み込む
**テキストファイルへの書き出し [#k6a104af]
>write.table (データフレーム名, file="ファイル名")
**データエディタの起動 [#d9e3bbd5]
-作成したデータフレームをスプレッドシートで編集したい場合:
> lm <- edit(bnclemma)
--編集したデータフレームは lm という名前(何でもよい)に...
--bnelemma という元のフレームはそのまま
-上書きで編集するには fix() とする
> fix(bnclemma)
-まったく新しいファイルを書くには
>dd <- data.frame()
>fix(dd)
**ワークスペースの履歴・保存 [#i65b2a7d]
-ls() 関数
-今までの実行したワークスペースのコマンド一覧を出す
> ls()
[1] "bmi" "d" "energy" "exp.lean...
[7] "h" "height" "hh" "intake.p...
[13] "mylist" "oops" "pain" "thue2" ...
[19] "x" "y" "ylim"
-ワークスペースをファイルに保存する
>save.image()
**データフレームの部分加工 [#n225bbb0]
>icci2 <-icci[c(1,2)]
#icci のデータフレームの第1、2コラムを取り出して icci2 ...
**データフレームのサブセットの取り出し [#h0326953]
> icci_A1 <- subset(icci, CEFR == "A1")
#subset()関数で、データ icci の CEFRの変数の値が A1 のも...
**データフレームのマージ [#g73594ce]
>icci <- merge(icci1, icci2)
# これでお互いのファイルに共通する変数名を参照しながらマ...
# もしそれぞれにない個体があった場合に、それをすべて含め...
>icci <- merge(icci1, icci2, all=TRUE)
**データフレームのマージ(行の結合) [#aad41245]
> total <- rbind(A1_sample, A2_sample, B1_sample)
# rbind で2つ以上のデータフレームをくっつける。カラム変...
**データフレームからランダムに X 件を抽出 [#z177b7b2]
> A1_sample <- icci_A1[sample(1:nrow(icci_A1), 1000, rep...
# icci_A1 の行データ(1から個体数すべて)のなかから1000...
**分布 [#w2f350a6]
-Density or point probability 密度確率・
--the probability density function 確率密度関数
--the probability distribution function 確率分布関数
-確率分布関数=確率密度関数を積分したもの
-確率密度関数=確率分布を微分したもの
--統計データが十分にたくさんあれば、
--ヒストグラム = 確率密度関数になり
--累積頻度分布 = 確率分布関数になる。
***正規曲線の書き方1 [#a43c63aa]
> x <- seq(-4,4,0.1)
-4 から 4 まで 0.1ずつの感覚で sequence を作って x に格納
> plot(x, dnorm(x),type="l")
X軸に x を Y 軸に x にそって確率密度の正規曲線を生成
type="l" は line かな
> plot(x, dnorm(x),type="h")
同じプロットを h にすると histogram のようにそれぞれの x...
***正規曲線の書き方2 [#dbcc2dff]
> curve(dnorm(x), from=-4, to=4)
これでも同じカーブを描くことができる
***累積分布関数 [#wd8f1286]
-正規分布のどの位置にあるか、その位置までの cumulative di...
--平均132、標準偏差が13の正規分布の中で、160とい...
---正規分布上の値のとる確率は何パーセントか、を pnorm が...
---それを 1 から引くことで、160以上の値をとる確率を得ている
> 1-pnorm(160,mean=132,sd=13)
[1] 0.01562612
> 1-pbinom(15,size=20,prob=.5)
[1] 0.005908966
**記述統計 [#kaed79a2]
-サンプルデータとしてランダムに50個の正規分布を成す数値...
> x <- rnorm(50)
> x
[1] 1.11403983 -1.15165895 0.80461867 1.71557037 0....
[8] -0.10692887 -0.96407973 1.16309160 0.61545475 0....
[15] 0.12446743 0.77148385 1.19391658 -0.78081663 -0....
[22] -0.15209363 -0.49532705 -0.52829549 0.64414085 1....
[29] 1.55049255 0.21063906 0.22529564 0.34462737 0....
[36] 1.15710520 -0.69795906 -1.12325445 -2.08717916 0....
[43] 1.04982555 0.46580345 0.55212419 -0.53356749 -0....
[50] -0.06345214
-個別の記述統計の出し方
--平均
> mean(x)
[1] 0.02739456
-欠損値があって飛ばしたい場合の処理
> library(ISwR) #ISwR というライブラリを読み込む...
> data(juul) #data juul を読み込む
> juul
age menarche sex igf1 tanner testvol
1 NA NA NA 90 NA NA
2 NA NA NA 88 NA NA
3 NA NA NA 164 NA NA
4 NA NA NA 166 NA NA
5 NA NA NA 131 NA NA
--juul のデータの igf1 を参照したいが、いちいち juul$igf1...
> attach(juul) #こうすると juul$…を省略できる
> mean(igf1)
[1] NA #NA は欠損値があるので平均値を出せない、とエ...
> mean(igf1,na.rm=T) #"na.rm=T" not available, remove ...
[1] 340.168
--標準偏差
> sd(x)
[1] 0.876461
--分散
> var(x)
[1] 0.768184
--中央値(メディアン)
> median(x)
[1] 0.06215406
--分位数(quantile)
> quantile(x)
0% 25% 50% 75% 10...
-2.08717916 -0.68625583 0.06215406 0.63696932 1.71557...
---デフォルトではこのように正規分布の下から 25,50, 75,100...
---もし quantile を特定の観察したい割合で切りたい場合には...
> pvec <- seq(0,1,0.1)
pvec という変数に 0 から 1 までを 0.1 区切りに格納
> pvec
[1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
> quantile(x,pvec)
0% 10% 20% 30% 40% ...
-2.08717916 -0.97999720 -0.77189339 -0.52987709 -0.124...
70% 80% 90% 100%
0.49169967 0.77811082 1.16617410 1.71557037
---quantile に2つ目の argument で何%で切るかを 0.00-1.0...
**記述統計を一発で出す方法 [#fc2b1dfa]
-summary 関数は 最小値、第1分位点、メジアン、平均、第3...
> summary(igf1)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
25.0 202.2 313.5 340.2 462.8 915.0 321.0
-データフレーム全体の summary も可能
> summary(juul)
age menarche sex ...
Min. : 0.170 Min. : 1.000 Min. :1.000 Min. ...
1st Qu.: 9.053 1st Qu.: 1.000 1st Qu.:1.000 1st Q...
Median :12.560 Median : 1.000 Median :2.000 Media...
Mean :15.095 Mean : 1.476 Mean :1.534 Mean ...
3rd Qu.:16.855 3rd Qu.: 2.000 3rd Qu.:2.000 3rd Q...
Max. :83.000 Max. : 2.000 Max. :2.000 Max. ...
NA's : 5.000 NA's :635.000 NA's :5.000 NA's ...
testvol
Min. : 1.000
1st Qu.: 1.000
Median : 3.000
Mean : 7.896
3rd Qu.: 15.000
Max. : 30.000
NA's :859.000
**ヒストグラムを書く [#xb4b861d]
> x <- rnorm(100) # normal density of pseudo-ran...
> x
[1] -0.23852918 0.80007985 2.59488947 1.92574080 0....
[8] -0.07618785 -1.59706302 0.82031697 -1.20817887 1....
[15] 0.75213738 -0.06597683 0.23352029 -0.16270841 -1....
[22] -0.12607458 0.09170498 -0.10259290 -0.45916474 0....
[29] 1.45689028 0.64155579 -0.33764720 -1.25446597 0....
[36] -0.26736349 -0.97858879 1.45124404 -1.22135682 0....
[43] -0.45558029 1.61702489 -1.51530741 1.12869612 1....
[50] -0.73521276 -0.12948352 -0.09378984 0.93129898 0....
[57] 0.25044045 1.01388932 0.82988086 -1.59370935 -0....
[64] -0.38521575 1.30834830 1.20246317 -0.03461022 -0....
[71] 0.99960792 0.87056316 1.17321910 0.02667156 0....
[78] -0.85243182 -0.87548459 0.20321376 -1.46930772 1....
[85] -0.68914995 1.06938790 1.29973824 0.65570765 0....
[92] -0.25530870 0.46701640 1.27102185 1.73760250 1....
[99] -0.32908069 -0.49896537
> hist(x, freq=F) # write histog...
> curve(dnorm(x),add=T) # write the no...
-縦軸の目盛りが切れてしまった場合の処理(少し複雑)
> h <- hist(x, plot=F)
> ylim <- range(0, h$density, dnorm(0))
> hist(x, freq=F, ylim=ylim)
> curve(dnorm(x), add=T)
**変数:数値と文字列の操作 [#a703d170]
-たとえば energy というデータは expend, stature という2...
> data(energy)
> energy
expend stature
1 9.21 obese
2 7.53 lean
3 7.48 lean
4 8.08 lean
5 8.09 lean
6 10.15 lean
7 8.40 lean
8 10.88 lean
9 6.13 lean
10 7.90 lean
11 11.51 obese
12 12.79 obese
13 7.05 lean
14 11.85 obese
15 9.97 obese
16 7.48 lean
17 8.79 obese
18 9.69 obese
19 9.68 obese
20 7.58 lean
21 9.19 obese
22 8.11 lean
-lean の1行だけとりだす
> exp.lean <- energy$expend[energy$stature=="lean"]
-obese の行だけとりだす
> exp.obese <- energy$expend[energy$stature=="obese"]
> exp.lean
[1] 7.53 7.48 8.08 8.09 10.15 8.40 10.88 6.13 7.9...
> exp.obese
[1] 9.21 11.51 12.79 11.85 9.97 8.79 9.69 9.68 9.19
-同様のことが以下の split()関数でもできる
> l <-split(energy$expend, energy$stature)
>
> l
$lean
[1] 7.53 7.48 8.08 8.09 10.15 8.40 10.88 6.13 7.9...
$obese
[1] 9.21 11.51 12.79 11.85 9.97 8.79 9.69 9.68 9.19
**データ・フレームの一部取り出し [#h2e2b342]
(例)verbs というデータから AnimacyOfTheme が animate ...
> verbs[verbs$AnimacyOfTheme == "animate",]
RealizationOfRec Verb AnimacyOfRec AnimacyOfTheme L...
58 NP give animate animate ...
100 NP give animate animate ...
143 NP give inanimate animate ...
390 NP lend animate animate ...
506 NP give animate animate ...
736 PP trade animate animate ...
> subset(verbs, AnimacyOfTheme == "animate")
RealizationOfRec Verb AnimacyOfRec AnimacyOfTheme L...
58 NP give animate animate ...
100 NP give animate animate ...
143 NP give inanimate animate ...
390 NP lend animate animate ...
506 NP give animate animate ...
736 PP trade animate animate ...
**標準化得点を出す方法 [#rb5c6194]
> x=rnorm(10,3,0.1)
> x
[1] 2.987822 3.087673 2.972487 2.901589 3.056770 3.07972...
> x-mean(x)
[1] -0.01775045 0.08210025 -0.03308569 -0.10398376 0.0...
> y = x-mean(x)
> y/sd(x)
[1] -0.2184789 1.0105193 -0.4072305 -1.2798693 0.63015...
> scale(x)
[,1]
[1,] -0.2184789
[2,] 1.0105193
[3,] -0.4072305
[4,] -1.2798693
[5,] 0.6301557
[6,] 0.9126354
[7,] -0.1403719
[8,] -0.2493192
[9,] -1.6753157
[10,] 1.4172752
attr(,"scaled:center")
[1] 3.005573
attr(,"scaled:scale")
[1] 0.08124561
**さまざまなプロット [#o8ba5c7c]
-棒グラフ
--ratings 中の1つの変数 Length に関して crosstab を作り...
>barplot(xtabs( ~ ratings$Length), xlab="word length", c...
-ヒストグラム
--ratings 中の1つの変数 Length に関して、ヒストグラムを...
--この truehist というコマンドは MASS というパッケージ中...
>truehist(ratings$Length, xlab="word length", col="grey")
-一般的な描画関数 plot()
>plot(h)
>plotsort(lexdec$RT), ylab = "log RT") # RT という変数...
-箱ひげ図を書く
>boxplot(lexdec$RT) # ふつうに RT という変数を箱...
>boxplot(exp(lexdec$RT)) # RT は対数なのでそれをふつう...
-すべての変数間のプロットを一発で出す → scatterplot matri...
>pairs(rating[ , -c(1, 6:8, 10:14)])
-カテゴリーデータをモザイク図にする
>mosaicplot(verbs.xtabs, main = "dative")
--上位階層にしたいものを先に xtab で処理しておく
-上位階層を変数でくくって箱髭図を書く
**複数の表を並べて作成する場合 [#cb1be9cd]
-par()関数を使うと、以後作成する表を1つ1つ並べて表示し...
>par(mfrow = c(3,2))
--2列3行で表を並べてくれる
**表を画像にはき出す場合 [#rcdf67f6]
-jpeg で barplot.jpeg というファイルに書き出す
>jpeg("barplot.jpeg", width = 400, height=420)
>truehist ... # ここで作表を行う
>dev.off() # 画像モードを抜ける
**単語リストの操作 [#v2bfd808]
-テキスト・ファイルを読み込む
(例)c:\temp にある lemma.num を bnclemma というデータ...
> bnclemma <- read.table("c:/temp/lemma.num", header=T)
header=T で1行目をヘッダとして読み込む
-データフレームを編集する
--作成したデータフレームをスプレッドシートで編集したい場...
> lm <- edit(bnclemma)
---編集したデータフレームは lm という名前(何でもよい)に...
---lemmalist として lm というデータフレームを作成した
-単語リストの抽出
--データフレームを行列として参照する
---ランク10位の単語([,3]はwordのカラムを示す)を参照
> lm[10,3]
[1] to
---ランク10位の品詞([,4])を参照
> lm[10,4]
[1] prep
Levels: a adv conj det infinitive-marker interjection mo...
※このようにアウトプットの最後に、その変数に何レベルある...
---lm 中の word のカラムの最初の4個を取り出す
> lm$word[1:4]
[1] the be of and
5464 Levels: a abandon abbey ability able abnormal aboli...
---lm 中の pos のカラムの1番目を取り出す
> lm$pos[1]
[1] det
Levels: a adv conj det infinitive-marker interjection mo...
---lm 中のwordの項目数をカウントする → length()を使用
> length(lm$word)
[1] 6318
---rank が100位以内の word を取り出す
> lm$word[lm$rank < 100]
[1] the be of and a in to ...
[12] i that you he on with do ...
[23] but from they his that she or ...
[34] say will would can if their go ...
[45] her make who as out up see...
[56] some could so him year into its...
[67] than more about now last your me ...
[78] should these people also well any onl...
[89] way look like use her such how...
---rank が100位以内で、pos が名詞のものを取り出す
> lm$word[lm$rank < 100 & lm$pos == "n"]
[1] time year people way
※ [ ] で条件設定し、& で2つの条件を組み合わせている
※ logical operator としては、&, |(=or), !(=not) が使える
---rank が100位以内で、pos が名詞または形容詞のものを...
> lm$word[lm$rank < 100 & (lm$pos == "n"|lm$pos == "a")]
[1] time year last other people new way goo...
※ logical operator | が使われていることに注意
---10位以内のリストをすべてとりだす
> lm[lm$rank <10,]
rank freq word pos
1 1 6187267 the det
2 2 4239632 be v
3 3 3093444 of prep
4 4 2687863 and conj
5 5 2186369 a det
6 6 1924315 in prep
7 7 1620850 to infinitive-marker
8 8 1375636 have v
9 9 1090186 it pron
※10位まで取り出すならば、<11 とする
※[]の中で条件を付けるが、最後の ,] とカンマをつけるとこ...
---100位以内で品詞が動詞である行だけをとりだす
> lm[(lm$rank <101 & lm$pos == "v"),]
rank freq word pos
2 2 4239632 be v
8 8 1375636 have v
18 18 559596 do v
34 34 333518 say v
40 40 249540 go v
44 44 220940 get v
46 46 217268 make v
51 51 191661 see v
52 52 185534 know v
54 54 179220 take v
64 64 153881 think v
66 66 151871 come v
76 76 131417 give v
90 90 111058 look v
92 92 108820 use v
100 100 98899 find v
-100位以内の単語の頻度プロットを作る
> lm100 <- lm[(lm$rank <101),] # lm100 というデータフレ...
> lm100
rank freq word pos
1 1 6187267 the det
2 2 4239632 be v
3 3 3093444 of prep
4 4 2687863 and conj
5 5 2186369 a det
6 6 1924315 in prep
7 7 1620850 to infinitive-marker
8 8 1375636 have v
9 9 1090186 it pron
10 10 1039323 to prep
11 11 887877 for prep
12 12 884599 i pron
13 13 760399 that conj
14 14 695498 you pron
15 15 681255 he pron
16 16 680739 on prep
17 17 675027 with prep
18 18 559596 do v
19 19 534162 at prep
20 20 517171 by prep
※出力はこんな感じ
プロットの関数は plot()
> plot(lm100$rank, log(lm100$freq), type="n")
※ランキングを横軸、頻度の自然対数を縦軸にとる。
※type="n" はプロットの記号を消す、という意味。後で単語...
> text(lm100$rank, log(lm100$freq), as.character(lm100$w...
※この text() という関数が交差する点に as.character() ...
※cex=で文字サイズを指定する
--頻度リスト&プロット:応用 [#c0377bc4]
---100語以内の動詞のプロットを作る
> vlm <- lm100[lm100$pos == "v",] #動詞を100位以内か...
> vlm
rank freq word pos
2 2 4239632 be v
8 8 1375636 have v
18 18 559596 do v
34 34 333518 say v
40 40 249540 go v
44 44 220940 get v
46 46 217268 make v
51 51 191661 see v
52 52 185534 know v
54 54 179220 take v
64 64 153881 think v
66 66 151871 come v
76 76 131417 give v
90 90 111058 look v
92 92 108820 use v
100 100 98899 find v
> plot(vlm$rank, log(vlm$freq), type="n") #横軸ラン...
> text(vlm$rank, log(vlm$freq), as.character(vlm$word), ...
※vlm$rank と vlm$freq の交差する点に vlm$word から取って...
※cex=0.9 で少し大きめに
‐サブコーパス間のランク差の統計情報を得る
# 複数の頻度表を、見出し語・ランク・頻度・コーパスの4コ...
sum <- numSummary(brownclone[,"Rank"], groups=brownclone...
statistics=c("mean", "sd", "quantiles"), quantiles=c(0,...
write.csv(sum$table, file="numsummary_brownclone.csv")
ページ名: