ラベル R入門 の投稿を表示しています。 すべての投稿を表示
ラベル R入門 の投稿を表示しています。 すべての投稿を表示

2012年6月20日水曜日

R入門02

関数の作成

> varp <- function(x){
+ 標準分散 <- var(x)*(length(x)-1)/length(x)
+ 標準分散
+ }
> x <- c(10,13,8,15,8)
> x
[1] 10 13  8 15  8
> var(x)
[1] 9.7
> varp(x)
[1] 7.76

 

関数名 <- function(x){

処理

}

が大まかな流れらしい。“+”は勝手に補完される。

だったらいちいち標準分散なんて変数作って、そいつを表示させるようなまどろこしい・・・てか馬鹿丸出しな関数つくんなと思っていたら

> varpp <- function(x){
+ var(x)*(length(x)-1)/length(x)
+ }
> varpp(x)
[1] 7.76

これでもきちんと動いた。大丈夫か本当に?

“標準分散”ていう変数に代入して残しておきたいのなら

標準分散<-varpp(x)で大丈夫なはずだし、どうしてこんなにめんどくさいことをするんだろうか?

もしかしたら出来ないのではないのかという妙な疑問が湧いたので試してみた

> 標準<-varpp(x)
> 標準
[1] 7.76

普通にできるじゃないか・・・大丈夫か本当にw

 

というか、function(x)の文字列をこのwindows live writerで入力しようとしたら

(x)の部分で

Image113

こんな画像になった。これは仕様なのか?俺としてはこっちのほうが戦慄が走り、背筋が凍りそうになった。



2011年4月23日土曜日

R入門08

2つの変数の記述統計

国語の点数が高いほど、英語の点数が高い傾向がある

・・・量的変数どうしの関係=相関

質的変数どうしの関係=連関

 

散布図の描写

Image125

a<-c(・・・・・)

b<-c(・・・・・)

plot(a,b)

 

共分散

・・・偏差の積の平均(此処で言う偏差とは平均からの偏差)

> sum((a-mean(a))*(b-mean(b)))/length(a)
[1] 7.55

または

> mean((a-mean(a))*(b-mean(b)))
[1] 7.55

一発出る関数はcov

> cov(a,b)
[1] 7.947368

しかしながら、これは先ほど求めた共分散と微妙に値が異なっている。

これは先程のようにnで割るか、n-1で割るかで差が出てくるから。

不偏共分散

 

共分散=不偏共分散*n-1/n

> cov(a,b)*(length(a)-1)/length(a)
[1] 7.55

 

今度は身長と体重で共分散その様子を見てみる。

> m<-c(1.5,1.6,1.7,1.8,1.9)
> cm<-m*100
> k<-c(50,70,60,80,90)
> cov(cm,k)
[1] 225
> cov(m,k)
[1] 2.25

 

相関係数

rxy=sxy/sxsy

> a
[1]  6 10  6 10  5  3  5  9  3  3 11  6 11  9  7  5  8  7  7  9
> b
[1] 10 13  8 15  8  6  9 10  7  3 18 14 18 11 12  5  7 12  7  7
> cov(a,b)/(sd(a)*sd(b))
[1] 0.749659

こいつを一発で出す関数はcor()

> cor(a,b)
[1] 0.749659

無相関に近づくと相関係数は0に近づく。

相関係数の値に明確な基準はないけど目安として

相関係数の絶対値が

0.2以下・・・ほとんど相関なし

0.4以下・・・弱い相関

0.7以下・・・中程度の相関あり

それ以上・・・強い相関あり

 

故に、先程算出した相関係数を見る限りa,bは強い相関があると考えられる。

しかし、極少数の外れ値によって相関係数が大きく左右されることがあるので、相関係数も大事だが、必ず散布図を描写すること。

 

クロス集計表

次は連関を見ていこう

質的変数どうしの関係。

aに関して

> a<-c("S","S","A","A","S","S","S","S","S","A","A","S","A","S","S","A","S","S","S","S")
> a
[1] "S" "S" "A" "A" "S" "S" "S" "S" "S" "A" "A" "S" "A" "S" "S" "A" "S"
[18] "S" "S" "S"
> table(a)
a
A  S
6 14

 

bに関して

> b<-c("A","A","A","A","S","S","S","S","S","S","A","A","A","S","A","S","S","S","S","S")
> b
[1] "A" "A" "A" "A" "S" "S" "S" "S" "S" "S" "A" "A" "A" "S" "A" "S" "S"
[18] "S" "S" "S"
> table(b)
b
A  S
8 12

 

abの連関

> table(a,b)
   b
a    A  S
  A  4  2
  S  4 10

 

ファイ係数

相関係数の特別な場合で、1と0の二つの値からなる二値変数に対して計算冴えれる相関係数のこと。

> aa<-ifelse(a=="A",1,0)
> aa
[1] 0 0 1 1 0 0 0 0 0 1 1 0 1 0 0 1 0 0 0 0

> bb<-ifelse(b=="A",1,0)
> bb
[1] 1 1 1 1 0 0 0 0 0 0 1 1 1 0 1 0 0 0 0 0

このように1、0に変換した後の相関係数がファイ係数。

> cor(aa,bb)
[1] 0.3563483

 

練習問題

Image129

> (sum((a-mean(a))*(b-mean(b)))/length(a))/(sqrt(sum((a-mean(a))^2))*sqrt(sum((b-mean(b))^2)))
[1] 0.09092974

 

> aa<-c("A","S","S","A","S","A","A","S","A","A","S","A","S","A","S","S","A","A","S","S")
> table(a)
a
1  3  4  5  6  8 10 12
2  2  1  1  1  1  1  1
> table(aa)
aa
A  S
10 10
> bb<-c("A","S","A","A","S","S","S","S","A","A","A","A","S","S","A","S","S","A","S","S")
> table(bb)
bb
A  S
9 11
> table(aa,bb)
   bb
aa  A S
  A 6 4
  S 3 7

 

> aaa<-ifelse(aa=="S",1,0)
> bbb<-ifelse(bb=="S",1,0)
> cor(aaa,bbb)
[1] 0.3015113



R入門07

練習問題

Image119

Image120

余計な物が写りこんでるが気にするな。

> mean(a)
[1] 112

> sqrt(sum((a-mean(a))^2)/length(a))
[1] 77.82031

> mean(b)
[1] 65
> sqrt(sum((b-mean(b))^2)/length(b))
[1] 31.70173

 

> (a-mean(a))/sqrt(sum((a-mean(a))^2)/length(a))
[1] -0.6682061 -0.1542014 -0.7967072 -0.9252084 -0.7967072  1.5163138
[7]  0.1028009  1.6448149  1.1308103 -1.0537096
> (b-mean(b))/sqrt(sum((b-mean(b))^2)/length(b))
[1] -0.4731602 -0.1577201 -0.7886004 -0.4731602  1.1040405  0.4731602
[7] -1.1040405 -1.4194807  1.1040405  1.7349208



R入門06

標準化

z得点=(データの値-平均)/標準偏差

> a<-c(13,14,7,12,10,6,8,15,4,14,9,6,10,12,5,12,8,8,12,15)
> a
[1] 13 14  7 12 10  6  8 15  4 14  9  6 10 12  5 12  8  8 12 15
> mean(a)
[1] 10
> sqrt(mean((a-mean(a))^2))
[1] 3.331666

ここまでは標準偏差の式。

> (a-mean(a))/sqrt(mean((a-mean(a))^2))
[1]  0.9004503  1.2006005 -0.9004503  0.6003002  0.0000000 -1.2006005
[7] -0.6003002  1.5007506 -1.8009007  1.2006005 -0.3001501 -1.2006005
[13]  0.0000000  0.6003002 -1.5007506  0.6003002 -0.6003002 -0.6003002
[19]  0.6003002  1.5007506

 

> sqrt(mean(((a-mean(a))/sqrt(mean((a-mean(a))^2))-mean((a-mean(a))/sqrt(mean((a-mean(a))^2))))^2))
[1] 1

 

偏差値=z得点*10+50
> 10*(a-mean(a))/sqrt(mean((a-mean(a))^2))+50
[1] 59.00450 62.00600 40.99550 56.00300 50.00000 37.99400 43.99700 65.00751 31.99099 62.00600 46.99850 37.99400 50.00000 56.00300
[15] 34.99249 56.00300 43.99700 43.99700 56.00300 65.00751

 

偏差値の平均は50、標準偏差が10になる事を確認
> mean(10*(a-mean(a))/sqrt(mean((a-mean(a))^2))+50)
[1] 50
> sqrt(mean(((10*(a-mean(a))/sqrt(mean((a-mean(a))^2))+50)-mean(10*(a-mean(a))/sqrt(mean((a-mean(a))^2))+50))^2))
[1] 10



R入門05

散布度

データを見る際には中心的位置の他にどれだけ散らばっていれば分かればなおよい。

それを散布度という。

> b<-mean(a)
> b
[1] 10.8
> a-b
[1] -0.8  2.2 -2.8  4.2 -2.8
> c<-a-b
> d<-c*c
> d
[1]  0.64  4.84  7.84 17.64  7.84
> e<-sum(d)
> e
[1] 38.8
> f<-length(a)
> f
[1] 5
> g<-e/f
> g
[1] 7.76
> sqrt(g)
[1] 2.785678

なんでこの本ではこんなにまどろっこしい事をするんだろうとかなり疑問に思った。

> sqrt(sum((a-mean(a))^2)/length(a))
[1] 2.785678

これで良くないか?

> var(a)
[1] 9.7

なぜ違うの?

>これは不偏分散という。ついさっき求めたのは標本分散という。

>標本分散・・・偏差の2乗/データ数

>不偏分散・・・偏差の2乗/(データ数-1)

∴不偏分散から標本分散を求めるには

> var(a)*(length(a)-1)/length(a)
[1] 7.76

こうすればよい

 

> sd(a)
[1] 3.114482

これは不偏分散の平方根

ここから標本分散を求めるには

> ((sd(a)^2)*(length(a)-1)/length(a))
[1] 7.76

こうすればよい。

 

平均偏差

個別のデータ - 平均の絶対値の平均

> mean(abs(a-mean(a)))
[1] 2.56

 

レンジ(範囲)

> max(a)-min(a)
[1] 7



R入門04

1つの変数の記述統計

データを構成するモノ・・・質的変数

CF.二値変数(とりうる値が2通りのみ) ⇔ 量的変数

別に同値じゃないけど分かるよね、対比て意味ね。

 

質的変数と量的変数の区別は統計解析上重要らしい。

 

データの視覚的表現

> a<-c("C","B","B","A","B","C","A","D","D","D","A","A","A","C","B","D","C","C","B","D")
> a
[1] "C" "B" "B" "A" "B" "C" "A" "D" "D" "D" "A" "A" "A" "C" "B" "D" "C"
[18] "C" "B" "D"
> table(a)
a
A B C D
5 5 5 5

 

> b<-c(13,14,7,12,10,6,8,15,4,14,9,6,10,12,5,12,8,8,12,15)
> b
[1] 13 14  7 12 10  6  8 15  4 14  9  6 10 12  5 12  8  8 12 15
> hist(b)

Image118

代表値の求め方

・平均

> 10+13+8+15+8
[1] 54
> sum(10,13,8,15,8)
[1] 54
> a<-c(10,13,8,15,8)
> sum(a)
[1] 54
> sum(a)/5
[1] 10.8
> sum(a)/length(a)
[1] 10.8
> mean(a)
[1] 10.8

 

・中央値

> median(a)
[1] 10

データを大きさの順に並べた場合真ん中に位置する値の事

奇数の場合は真ん中が存在しうるが、データ数が偶数の場合には真ん中の二つの値の平均を持って中央値することが一般的

 

・最頻値

> table(a)
a
8 10 13 15
2  1  1  1

から8だと分かる。

最も頻繁に観測される値。

最頻値が2つあるような分布(二峰型分布)の場合、最頻値を代表値として用いるのは適切でないと思われます。

複数の最頻値がある場合、例えばSASではそれらの最小値を最頻値として表示しますが、いずれにせよ、こういう場合は最頻値以外を代表値として用いるほうが無難です。

データに外れ値なのがある場合など、中央値が用いられる場合はありますが、一般的には、代表値としては平均を用いると考えて差し支えありません。

 

 

だとよ。



R入門03

関数をパッケージする。

まずはファイル>新しいスクリプトを開く

すると以下のようなエディット窓が開く

Image114

こいつに次のような標準分散を求める関数をカキカキ。

varp(x)<-function(x){
    var(x)*(length(x)-1)/length(x)
}

で保存するんだけど。

Image115

普通のソフトはこれだけ入力すれば、後は拡張子は自動的に付けてくれるものなんだけどこいつは付けてくれない。

だから、自分で.rと拡張子を入力する必要が出てくる。

このまま保存しても、基本的には問題ないんだけど、このままだと拡張子の全くつかないファイルとして保存されることになる。

だから次回このファイルを読み込む際に、すべてのファイルという指定をつけないと認識してくれないので、やはりめんどくさいが.rは付けることにしよう。

で、Rgui.exeを再起動

Image117

> x<-c(11,464,12,55,123)
> varp(x)

と入力すると以下の文字列が返される。

>エラー:  関数 "varp" を見つけることができませんでした

つまり、この段階ではvarpという関数は定義されていないと言うことだ。

> source("標準分散.r")
> varp(x)
[1] 29054

ここでsource(“ファイル名”)を実行させることで、そのファイルの関数を利用出来るようになる。

あっそwすごいね~しびれるぅ。



R入門01

> 2+4
[1] 6

 

> 20/5
[1] 4

 

> sqrt(200)
[1] 14.14214

 

なんか便利な関数があるアルヨてな話

> summary(c(173,178,180,183,182,174,179,179,174,192))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  173.0   175.0   179.0   179.4   181.5   192.0

 

変数に代入するのにいちいち<-ていう演算子ぶち込まないといけないよねてな話

> 身長 <- c(173,178,180,183,182,174,179,179,174,192)
> 身長
[1] 173 178 180 183 182 174 179 179 174 192

 

でこいつは本には書いてなかったのだが、一回代入すると、多分ソフトを再起動させるか、同じ変数名に別の数値を代入しない限り、その値は保持されるみたい。

> 身長 <- c(173,178,180,183,182,174,179,179,174,192)
> 身長
[1] 173 178 180 183 182 174 179 179 174 192
> sqrt(身長)
[1] 13.15295 13.34166 13.41641 13.52775 13.49074 13.19091 13.37909 13.37909
[9] 13.19091 13.85641
> summary(身長)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  173.0   175.0   179.0   179.4   181.5   192.0

 

ここで発見したこと。どうやらこの言語には予約語という概念は存在しないらしい。

要するにsqrt(sqrt)ということも、sqrtという変数に正のsqrt関数有効範囲内の数値が代入されていれば有効であるということだ。

なんとも中途半端な言語だなあ。本当に使っていて大丈夫なのか?これ?

> matrix <- matrix(c(173,178,180,183,182,174,179,179,174,192,17000,14000,9000,50000,30000,12000,900,2100,1000,25000),10,2)
> matrix
      [,1]  [,2]
[1,]  173 17000
[2,]  178 14000
[3,]  180  9000
[4,]  183 50000
[5,]  182 30000
[6,]  174 12000
[7,]  179   900
[8,]  179  2100
[9,]  174  1000
[10,]  192 25000

で、このmatrix(c(x,y,・・・・・)n,m)

これでnm行列の作成と言うところだね。

まあ、統計言語としちゃあ仕方ないんだけど、俺的には配列変数を応用してやったほうが直感的で使い良いと思うんだけどね。

Image111

test.csvの中身

> hoks <- read.csv("test.csv")
> hoks
    X height salary
1   1    173  17000
2   2    178  14000
3   3    180   9000
4   4    193  50000
5   5    182  30000
6   6    174  12000
7   7    179    900
8   8    179   2100
9   9    174   1000
10 10    192  25000

csvファイルから直接読み込めるのはちょっと便利かもね。



R入門

なにこれ?

>フリーのソフトウェア。自由に改変、再配布可能。

 

どんなやつ?具体的に。

>こんなの、貧相なUI。

>Image107

>しかし、この貧相さがグッと来る。貧相といえば・・・

>ああここはそういう事をいう場所じゃなかったね~

 

インスコ

>ここから http://www.r-project.org/:公式

>Image104

 

文字化け?

>インストーラを走らせると、Comodoオオトカゲさんを使っているとこう言うのが出る。

>Image109

>こんなモノをリメンバーしても意味ないから、チェック外してインストローラとして実行させる。

>で、こんなんが出る。

>ffdshowでもあったけど、なんでこうなるのだ?●ね。

Image108

 

話変わるけど、今シュタゲでも見てるけどなんぞあれ?

ストーリはまあまあ、椎名へきる=花澤香菜で視聴確定したけど・・・

正直見ていてしんどい。あんなコッテコテのWEBランゲージ使うやついないだろしかし。

あんな典型的なてか、あんな言葉リアルで言っている「香具師」なんていないでしかしお寿司だお(^^)。

 

視聴中にため息が耐えないアニメだが、3話の終わりごろからようやくオッと思う展開が。

あのWEBランゲージが耳につかないほどの超展開、面白い展開を期待したいところだな。

原作をやっていないから、所見だからとても楽しみだなあ。アニメ終わったら原作でもプレイしてみようか。