2015年4月28日火曜日

Rのhistとggplot(geom_histogram)を比較する

ggplot2パッケージって使ったことありませんでした。記法がとっつきにくいというのもあったのですが、あの見た目を素直に「きれい」とは認めたくなかったというか。「ggplotの方がかっこいい」という人への反発みたいな部分も多分にあるのですが。

「白地に線画」みたいなビジュアルが好きなんです。下手に色をつけたらかえって趣味が悪くなったりするじゃないですか。自分自身センスがないことを自覚しているので、色やら塗りつぶしやら飾りっぽいものやらは極力使わないようにしていたんですよね。

でも今回、ggplot2パッケージを使ってみようと思ったのは、ベースとなるレイヤーを用意して、次に重ねたいレイヤーを(場合によっては複数)用意して、最後にバンっとplotする、っていうコンセプトがなんだか使いやすいのではないかという気がしたから。なので、徐々に試していこうかなと思っています。

前置きが長くなりました。

今回はヒストグラムですが、多くの人に馴染みのあるデフォルトのhistでの描画と、ggplot2パッケージを使った場合とで、どんな風に違うのか試してみました。

正規分布に従う乱数を発生させて、普通にヒストグラムを描いてみます。


data1 <- rnorm(1000)
 
hist(data1)

通常のhist関数で描画

ggplot2パッケージを使うと↓こんな描き方になります。

install.packages("ggplot2")
library(ggplot2)
 
df1 <- data.frame(data1) # データフレームにしておく必要がある
 
g <- ggplot(df1, aes(x=data1)) # 対象がdata1列であることを指定
g <- g + geom_histogram()      # ヒストグラムを描画することを指定
plot(g)                        # 描画

ggplot2パッケージを使ったもの

先ほどの、デフォルトのhistと階級幅を合わせる(=0.5)には、binwidthを指定します↓

g <- ggplot(df1, aes(x=data1))
g <- g + geom_histogram(binwidth=0.5)
plot(g)

ggplot2(階級幅を指定)

ちなみにまとめて書いてしまう方法もあって、これくらいレイヤー数だったら、この方がシンプルかもしれません。

ggplot(df1, aes(x=data1)) + geom_histogram(binwidth=0.5)


で、見た目どうでしょうか、黒のベタ塗りはけっこう好みが分かれそうです。

でも、もちろん、指定すればどうにでもできます。

g <- ggplot(df1, aes(x=data1))
g <- g + geom_histogram(colour="black", fill="white", binwidth=0.5)
plot(g)

ggplot2(ボーダーと塗りつぶしの色を指定)

colourで枠線の色を、fillで塗りつぶしの色を指定しました。私はこっちの方が好きですね。




2015年4月17日金曜日

Rでブートストラップ法

「ブートストラップ」という言葉は聞いたことはあるけど、具体的にどんなことなのかはよく知らない、というのが私のレベルです。

が、「統計学入門」(東京大学教養学部統計学教室 編)という教科書↓

統計学入門 (基礎統計学)
統計学入門 (基礎統計学)

の中の演習問題に出てきた(でも、教科書の本編には解説はない)ので、Rでちょっとやってみました。

第3章 練習問題
3.4 <ブートストラップ> 以下の手法をパソコンで行え。

i)  [1, 11]に属する、整数の乱数を発生させる手続きを考えよ。

ii)  i)を11回繰返して実行し、11個のランダムな番号を得たのち、表3.1、図3.1のデータからその番号のデータを取り出し(ただし、同一番号があれば重複して取り出し、そのデータは複数個あると考える)、相関係数rを計算せよ。

iii) ii)を200通り繰返し、相関係数の値 r1, r2, ..., r200 を得たとき、そのヒストグラムを作れ。以上の統計手法をブートストラップという。

「表3.1、図3.1」というのは、教科書の本編のところに出てきた表と図で、相関関係が認められる例として紹介されていたものです。

統計学入門 (基礎統計学) より

統計学入門 (基礎統計学) より

このデータを使う必要があるので、まずは入力しましょう。

> bro <- c(71, 68, 66, 67, 70, 71, 70, 73, 72, 65, 66)
> sis <- c(69, 64, 65, 63, 65, 62, 65, 64, 66, 59, 62)

書籍に載っていた図と同じような感じでプロットしてみました↓

> plot(bro, sis, xlim=c(68,70), ylim=c(59,70), asp=1, pch=16)

書籍に載っていたデータをRでプロットしたもの

書籍に出てくる図とほぼ同じものができました。ただ、点が10個しかなく1個少ないです。これは(70, 65)という全く同じデータが2つあって同じ座標にプロットされているためです。書籍の方では少しずらして、分かるようになっていました。

あと、書籍だと(64,62)あたりに点がありますが(一番左側のもの)、これは間違いでしょうね。点を数えてみると12個あるし・・・

使うデータは入力したので、演習問題をやっていきましょう。

i)  [1, 11]に属する、整数の乱数を発生させる手続きを考えよ。

乱数の従う分布について特に言及していないので、一様分布でしょうね。

まず思いつくのは、

> ceiling(runif(1, 0, 11))   # 1~11の整数の乱数

runif関数が返す値をceiling関数で整数化してやる感じですね。

でも、Rで整数の乱数を発生させるんだったら、sample関数を使った方がシンプルな気がします↓

> sample(1:11, 1, replace=T) # 1~11の整数の乱数

1, 2, ..., 11 の中から1個選ぶ、復元抽出(replace=T)で、という感じですね。

ii)  i)を11回繰返して実行し、11個のランダムな番号を得たのち、

これはsample関数の引数を変えるだけでOKですね。

> rnd <- sample(1:11, 11, replace=T)
> rnd
 [1]  8  1  8  2  4 11  5  3  9  8 10

表3.1、図3.1のデータからその番号のデータを取り出し(ただし、同一番号があれば重複して取り出し、そのデータは複数個あると考える)、

番号に応じたデータを取り出すには、先ほど発生させたrndをインデックスに指定ればいいですね。

> bro[rnd]
 [1] 73 71 73 68 67 66 70 66 72 73 65

> sis[rnd]
 [1] 64 69 64 64 63 62 65 65 66 64 59

相関係数rを計算せよ。

> cor(bro[rnd], sis[rnd])
[1] 0.5357521

iii) ii)を200通り繰返し、相関係数の値 r1, r2, ..., r200 を得たとき、そのヒストグラムを作れ。以上の統計手法をブートストラップ(太字)という。

> r <- numeric(200) # 200回分の格納用
>
> for ( i in 1:200 ) {
+   rnd <- sample(1:11, 11, replace=T)
+   r[i] <- cor(bro[rnd], sis[rnd])
+ }
>
> hist(r)

ブートストラップ法で作成したヒストグラム

これで宿題は終わりました。単位は取れそうです。

元のデータの相関係数をそのまま出すと、

> cor(bro, sis)
[1] 0.5580547

になるのですが、ヒストグラムのモードは違う場所(0.65あたり)に出てますね。分布も左に歪んでいるし。

以上の統計手法をブートストラップという。

とのことですが、ヒストグラムの意味(モードのずれや歪み)や使い方については、もうちょっと勉強しないといけなさそうです。




2015年4月14日火曜日

Rで y = x のグラフを描きたいとき

Rでグラフを描きたいときにはcurveを使いますよね。

y = x^2 のグラフなら↓こんな感じ

curve(x^2, -1, 1)


じゃあ、y = x のグラフのグラフは、↓これでいけるかと思いきや・・・

curve(x, -1, 1)

 以下にエラー eval(expr, envir, enclos) :
   関数 "x" を見つけることができませんでした


xだけじゃ、関数だと思ってくれなかったんでしょうか?

とりあえず、↓こうすれば大丈夫ではあります。

f <- function(x){x} # y=x の関数を定義しておく

curve(f, -1, 1)


でも、なんか面倒臭いですよね。

で、実は括弧で囲めばOKということに最近気付きました。

curve((x), -1, 1)


しかしまあ、「y = x」のグラフを描くことって、滅多にないでしょうけど。




Rの数式描画の中で変数を展開したいときはbquoteを使う

ループを使って複数のグラフを描くような状況。

par(mfrow=c(2,2)) # 2行2列で作図する
 
for(i in 1:4){
  curve( x^i, -1, 1, asp=1, ylab="y" )
}

複数のグラフを描いてみた(タイトルなし)

グラフを区別するためにmainオプションでタイトルを付けましょう↓

par(mfrow=c(2,2))
 
for(i in 1:4){
  curve( x^i, -1, 1, asp=1, ylab="y",
         main=paste("y=x^", i) )
}

タイトルをつけてみたが格好悪い

う~ん、伝わりはするけど、"y=x^2"みたいな見た目が格好悪い。ちゃんと数式(数字の肩に指数が乗っている感じ)で表現したいですよね。

そういうときは、expression関数で、いけるかな・・・?

par(mfrow=c(2,2))
 
for(i in 1:4){
  curve( x^i, -1, 1, asp=1, ylab="y",
         main=expression(y==x^i) )
}

expressionだと変数iが展開されない

ありゃ? 数式の形にはなりましたが、変数iが展開されずそのまま「i」と出ちゃいました。

やっと本題です。

このiを展開(i = 1, 2, 3, 4)したいときには、bquote関数が使えます↓

par(mfrow=c(2,2))
 
for(i in 1:4){
  curve( x^i, -1, 1, asp=1, ylab="y",
         main=bquote(y==x^.(i)) )
}

bquoteを使えば、変数を展開できる

はい。無事、所望の表現になりました。

bquoteはexpressionと同じように数式表現を解釈してくれるのですが、ピリオドと括弧「.( )」で囲まれた部分については変数を展開してくれます。(という仕様のようです。私もそんなに使い込んだわけではないので・・・)




2015年4月10日金曜日

ベンフォードの法則をRでシミュレーションする

電気料金の請求書の金額の1桁目をたくさん集めて、その分布を見たらどうなっているでしょう。

数字の出方に法則性があるようにも思えないから、1から9まで一様になっているのでは? という気もしますが、そうではないというのがベンフォードの法則です。

ベンフォードの法則 - Wikipedia

ベンフォードの法則は、自然界に出てくる多くの(全てのではない)数値の最初の桁の分布が一様ではない、ある特定のものになっているというものである。
(中略)
この直感に反するような結果は、電気料金の請求書、住所の番地、株価、人口の数値、死亡率、川の長さ、物理・数学定数、冪乗則で表現されるような過程(自然界ではとても一般的なものである)など、様々な種類の数値の集合に適用できることがわかっている。


どのように偏って分布しているかというと、↓こんな割合になるとのこと。


ベンフォードの法則 - Wikipedia より


本当にそうなるか、Rでシミュレーションしてみましょう。

数値が満たすべき分布について、ウィキペディアには、

論理的には、数値が対数的に分布しているときは常に最初の桁の数値がこのような分布で出現する。

とあります。いまいちピンとこなかったのですが、あわせて掲載されていた図でなんとなく理解しました。



ふむふむ、対数スケールのグラフ上に、見た目で一様になるように分布すればいいのね。

理解を助けるために、まず、対数の数直線を書いてみましょう。


x <- c(1, 10, 100, 1000)
y <- rep(0, length(x))
plot(x, y, axes=F, xlab="", ylab="", log="x")
abline(h=0)
text(x,y, x, pos=1)



runif関数で発生させた乱数をそのまま使っちゃだめですよね。それでは数値的な一様になってしまいますから。

数直線に指数表記を追加すると分かりやすいです。

a <- round(log10(x), 3) # 指数の肩の数字
for ( i in 1:length(x) ) {
  text(x[i],y[i], bquote(10^.(a[i])), pos=3)
}



指数の肩のところがリニアになっているわけですから、発生させた一様な乱数を指数の肩のところに入れて得られた数値をサンプルとすればよさそうです。

n <- 1000000        # サンプルサイズ(これくらいだとやや時間かかる)
u <- runif(n, 0, 3) # 0~3の範囲の一様乱数
x <- 10^u           # 乱数の指数の肩に
msd <- as.integer(substr(x, 1, 1)) # 最上位桁の取り出し
round(table(msd) / n * 100, 1)     # 割合をパーセントで表示


これの実行結果が↓これです。じゃん。

msd
   1    2    3    4    5    6    7    8    9
30.1 17.6 12.5  9.7  7.9  6.7  5.8  5.1  4.6

おお、前述のウィキペディアにあげてあった表の数値と全く同じ結果になりました。

最後にヒストグラムを描いてみましょう。

h <- hist(msd, breaks=seq(0.5, 9.5, 1), xaxt="n", yaxt="n",
          main="ベンフォードの法則", xlab="最上位桁の数字", ylab="出現割合")
percent <- paste(round(h$counts / n * 100, 1), "%", sep="")
text(h$mids, h$counts, percent, pos=1)
axis(side=1, at=1:9)



うん、ここまでやると納得できましたね。




2015年4月9日木曜日

Rで最上位桁を取り出す

なんでこんなことをしたいのかは、また別の機会に書くとして、やりたいのは数字の最上位桁の数字を取り出したいというものです。123だったら1、31415だったら「3」を取り出すという感じです。

整数商を使えばいけそうです。123から1を取り出すには100で割る、31415から3を取り出すには10000で割る。

  123 ÷  100 = 1 余り  23

 31415 ÷ 10000 = 3 余り 1415

割る数は、10×(桁数-1)ですね。↓こんな感じでしょうか。

n <- floor(log10(x)) # 桁数-1
msd <- x %/% 10^n    # 整数商


ちなみに、MSDというのはMost Significant Digit(最上位桁、最上位の数字)の略です。

正しく結果が出ているか、乱数で確認してみましょう。

set.seed(1) # 再現性のために
x <- sample(1:1000000, 20)

n <- floor(log10(x)) # 桁数-1
msd <- x %/% 10^n    # 整数商

data.frame(x, msd) # 並べて比較するためにデータフレームに

        x msd
1  265509   2
2  372124   3
3  572853   5
4  908206   9
5  201682   2
6  898386   8
7  944670   9
8  660794   6
9  629110   6
10  61786   6
11 205973   2
12 176555   1
13 687015   6
14 384099   3
15 769831   7
16 497692   4
17 717608   7
18 991890   9
19 380029   3
20 777431   7


うまくいっているようですね。

でも、最上位桁を求めるために、常用対数とfloorと整数商を動員するのか。なんだか仰々しいというか。

要は1文字目を取り出せばいいんでしょ、ってところは後から気づきました。

↓文字列型に変換して、1文字目を取り出して、それを整数型に変換する、という考え方です。このほうがシンプルですね。

msd2 <- as.integer(substr(as.character(x), 1, 1))

data.frame(x, msd, msd2)

        x msd msd2
1  265509   2    2
2  372124   3    3
3  572853   5    5
4  908206   9    9
5  201682   2    2
6  898386   8    8
7  944670   9    9
8  660794   6    6
9  629110   6    6
10  61786   6    6
11 205973   2    2
12 176555   1    1
13 687015   6    6
14 384099   3    3
15 769831   7    7
16 497692   4    4
17 717608   7    7
18 991890   9    9
19 380029   3    3
20 777431   7    7


うん、同じ結果になりました。

タイプ数は大して変わりませんが、ロジックがシンプルな分、後者の方がいいような気がしますね。

ちなみにas.characterは省略できる(勝手に文字列として解釈してsubstrを適用してくれる)みたいですが、私は書いておく派ですね。




2015年4月8日水曜日

Rの正規表現で最短マッチ(最短一致)させる

特にRに固有のことというわけではありませんし、正規表現に詳しい人にとっては当たり前のことだとは思いますが、私自身が結構つまづいたので・・・

↓こんな感じのベクトル(住所リストのイメージ)があったとします。

住所リスト <- c( "佐藤町中村1-2",
                 "高橋町新田2-3",
                 "佐々木町本町3-4",
                 "伊藤町原4-5",
                 "山本町町田5-6",
                 "長谷川町本郷6-7" )

住所リスト


[1] "佐藤町中村1-2"   "高橋町新田2-3"   "佐々木町本町3-4"
[4] "伊藤町原4-5"     "山本町町田5-6"   "長谷川町本郷6-7"


で、それぞれから冒頭の「○○町」の部分を取り除きたいとします。

町名の長さも一律ではないので、こういうときは正規表現が便利ですよね。

私が最初に思いついた正規表現が以下です。

sub("^.+町", "", 住所リスト, useBytes=F)

冒頭()から、任意の文字()が、1つ以上()あって、「」で終わるような文字列、というつもりです。

ですが、うまく行きませんでした。↓実行してみると・・・

sub("^.+町", "", 住所リスト, useBytes=F)

[1] "中村1-2" "新田2-3" "3-4"     "原4-5"   "田5-6"   "本郷6-7"

町名より下のレベルの住所に「町」という文字を含んでいると、そこまで余分にマッチしてしまったというわけです↓

"佐々木町本町3-4"
"山本町町田5-6"

これは正規表現のデフォルトの動きが最長マッチ(最長一致)になっているからなんですね。

これを最短マッチ(最短一致)で動くように指定してやれば、意図通りに動いてくれそうです。

最短マッチにするには1個以上の繰り返しの「+」に「?」を付け足してやればOKです。

sub("^.+?町", "", 住所リスト, useBytes=F)

[1] "中村1-2" "新田2-3" "本町3-4" "原4-5"   "町田5-6" "本郷6-7"

意図通りの結果になりました。同様に、0個以上の繰り返しの場合には「*?」を使えばOKです。

上記は、マッチした文字列を除去する例でした。では、マッチした文字列を抽出するにはどうすればいいでしょうか。

私はこの方法くらいしか思いつかなかったのですが、なんだかまどろっこしい。もっといいやり方があるのではないかという気もしていますが・・・

# マッチする開始位置を得る

位置 <- regexpr("^.+?町", 住所リスト, useBytes=F)
位置


[1] 1 1 1 1 1 1
attr(,"match.length")
[1] 3 3 4 3 3 4


# マッチする長さも、regexprの戻り値に入っているので、取り出す

長さ <- attr(位置, "match.length")
長さ

[1] 3 3 4 3 3 4

# マッチ位置とマッチ長さを使って、substrで抽出
substr(住所リスト, 位置, 長さ)

[1] "佐藤町"   "高橋町"   "佐々木町" "伊藤町"   "山本町"   "長谷川町"

無事、町名を取り出せたようです。

ちなみに、substrの第3引数には終了位置を指定する必要があります。上記の例では開始位置が1固定なので、終了位置=長さであるため、そのまま指定しましたが、開始位置が1固定でない場合は、

substr(住所リスト, 位置, 位置 + 長さ - 1)


とする必要があります。今回の例だと、実行結果は同じになります。

substr(住所リスト, 位置, 位置 + 長さ - 1)

[1] "佐藤町"   "高橋町"   "佐々木町" "伊藤町"   "山本町"   "長谷川町"




2015年4月7日火曜日

Rで最頻値(mode)を求める ~連続値データの場合~

Rで最頻値を求める方法として、よく見かけるのがtable関数を使った方法です↓

names(which.max(table(x)))

離散値データの場合は上記でうまくいくのですが、連続値データの場合はうまくいきません。

また離散値データであっても、取りうる値の種類に対してサンプルサイズが小さい場合もうまくいきません。

具体例をあげてみます。40人学級でのテストの点数をイメージしました。

# テストの点数っぽいサンプルデータを作る
set.seed(4) # (やや恣意的な)再現性のために
x <- round(rnorm(n=40, mean=50, sd=15))
x


 [1] 53 42 63 59 75 60 31 47 78 77 58 50 56 49 51 53 67 49 48 46
[21] 73 52 70 69 59 46 69 64 36 69 52 66 39 28 63 44 47 64 43 40


上記データに対して、tableを使って最頻値を求めると、

names(which.max(table(x)))

[1] "69"

想定している、"50"とかなりずれた値が出てしまいました。実現値がスパースなので、個別に最も現れた値を出すことに意味がないんですね。

↓階級を1点刻みにしたヒストグラムで見るとよく分かります。

# 1点刻みでヒストグラムを描く
hist(x, breaks=seq(20,80,1))


69点を取った人は3人で最も多いですが、たまたま感がありますよね。50点の周辺の方が頻度が高いように思えます。

すなおに階級に幅を持たせると、予想と近い最頻値が得られます。

# 階級幅はデフォルト(hist関数にまかせる)
hist(x)


上記の場合、40~50の階級が最も度数が多いので、この中間値を階級値として「最頻値は45である」という結論は直感に反していないと思われます。

階級値をきりのいい数字にしたい場合は、↓breaksを指定してやればOKです。

hist(x, breaks=seq(25, 85, 10))



最頻値は50の階級であるという、予想通りの結果となりました。

つまり、連続値(やそれに近い離散値)における最頻値を求めるとは、ヒストグラムを描いてその峰になっているところを求めるということであり、tableを使った方法だとうまくいかない。

というのが、長い前置きでした。

で、ヒストグラムの峰の階級値を知るには、hist関数の戻り値に入っている、counts(度数)とmids(階級値)を使えばOKです↓

h1 <- hist(x)
h1$mids[which.max(h1$counts)]


[1] 45

階級の取り方を指定する場合↓

h2 <- hist(x, breaks=seq(25, 85, 10))
h2$mids[which.max(h2$counts)]


[1] 50

階級の取り方が適切かどうかは目で見て判断しなきゃいけないので、hist関数を実行した際に表示される描画結果を見てチェックするという必要はありそうです。




2015年4月1日水曜日

所得の分布は対数正規分布に従っているのか?

日本の世帯ごとの所得をヒストグラムにすると↓こんな感じになります。

所得金額階級別にみた世帯数の相対度数分布(平成21年調査)
2 所得の分布状況|厚生労働省 より


右側に裾が長く延びているので、正規分布ではないようです。

よく、所得は対数正規分布に従う、なんて聞いたりします。所得額の対数をとってそれで度数分布表を作るような感じですね。

実際に試してみましょう。

グラフで示されていた値を、データ化します。

rate <- 0.01 * c(6.6, 12.7, 13.9, 13.3, 10.0, 8.9, 7.1, 6.2, 5.1, 3.9,
                 3.0,  2.1,  1.7,  1.3,  0.9, 0.9, 0.5, 0.4, 0.2, 0.2)
income <- seq(50, 1950, 100)
df <- data.frame(income, rate)
df
   income  rate
1      50 0.066
2     150 0.127
3     250 0.139
4     350 0.133
5     450 0.100
6     550 0.089
7     650 0.071
8     750 0.062
9     850 0.051
10    950 0.039
11   1050 0.030
12   1150 0.021
13   1250 0.017
14   1350 0.013
15   1450 0.009
16   1550 0.009
17   1650 0.005
18   1750 0.004
19   1850 0.002
20   1950 0.002

「2000万円以上」の階級はオープンエンドになっているので省きました。所得額の中央の値を階級値(単位:万円)として、rateはパーセントでなく割合で表しています。

そのままplotしてみると↓

plot(df$income, df$rate)
所得ごとの度数をプロット

当然、厚労省のページに載ってたような形が出てきます。

データフレームに、対数をとった所得の列を追加してみます↓

df$income_log <- log10(df$income)
df
   income  rate income_log
1      50 0.066   1.698970
2     150 0.127   2.176091
3     250 0.139   2.397940
4     350 0.133   2.544068
5     450 0.100   2.653213
6     550 0.089   2.740363
7     650 0.071   2.812913
8     750 0.062   2.875061
9     850 0.051   2.929419
10    950 0.039   2.977724
11   1050 0.030   3.021189
12   1150 0.021   3.060698
13   1250 0.017   3.096910
14   1350 0.013   3.130334
15   1450 0.009   3.161368
16   1550 0.009   3.190332
17   1650 0.005   3.217484
18   1750 0.004   3.243038
19   1850 0.002   3.267172
20   1950 0.002   3.290035

この対数所得の列でプロットしてみると↓

plot(df$income_log, df$rate)
所得の対数をとって、度数をプロット

なんだか正規分布っぽくなりました。

では、この対数所得に対して、平均と分散を求めてみましょう。度数分布表から平均や分散を求める方法は、高校で習いましたよね?

m <- sum( df$income_log * df$rate)          # 平均
v <- sum( (df$income_log - m)^2 * df$rate ) # 分散

m
[1] 2.560478
v
[1] 0.1332302

対数をとったものは、平均 2.560478、分散 0.1332302 になることが分かりました。

この値を使って100万世帯分の所得データのサンプルを作ってみましょう。

n <- 1000000 # とりあえず100万世帯
r <- rnorm(n, m, sqrt(v))
head(r, 20)
 [1] 2.448054 3.062347 2.342119 2.527496 2.299730
 [6] 3.189821 2.096890 3.187842 2.440127 2.915923
[11] 2.553295 2.497428 2.379083 2.774537 2.668117
[16] 2.475265 2.716102 2.395672 2.960818 2.258748

上記は対数をとった所得になっているので、元に戻すと↓

head(10^r, 20)
 [1]  280.5781 1154.3754  219.8460  336.8965  199.4023
 [6] 1548.1791  124.9941 1541.1390  275.5034  823.9926
[11]  357.5157  314.3603  239.3772  595.0269  465.7116
[16]  298.7205  520.1181  248.6978  913.7295  181.4461

1人目から、280万円、1154万円、219万円、・・・という感じです。

できれば2人目くらいになりたいものです。

ヒストグラムを描いてみましょう↓

hist(10^r, breaks=seq(0, 100000, 100), xlim=c(0,2000), ylim=c(0,200000), col="lightgray")
対数正規分布に従う乱数のプロット

なんだか、っぽい感じですね。

比較できるように割合の数値を書き込んでみます↓

h <- hist(10^r, breaks=seq(0, 100000, 100), xlim=c(0,2000), ylim=c(0,200000), col="lightgray")
text(h$mids, h$counts, round(h$counts / n * 100, 1), pos=3, cex=0.8)
対数正規分布に従う乱数のプロット(割合を表示)

う~ん、割合の数字がけっこう違いますね。ピークの位置も微妙に違うし、尖り具合も違うし。

「まあ、そんなうまくはいかないよな」な例でした。

発生させた乱数でも分かるように、ちょっと低所得者層が多めに出ているようですね。

でも、当たらずといえども遠からずな感じだと思いません?