ヒストグラムと密度の推定
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
// 舟尾
COLOR(red){SIZE(25){ヒストグラムと密度の推定}}
([[グラフィックス参考実例集]]に戻る。[[Rのグラフィックス...
#contents
~
-混合分布からのサンプリングコードは鮮やかですね。data3 <-...
-修正ありがとうございます。ついでにコメントも入れました。...
-多(2)変量密度推定に関するパッケージや関数についての情報...
-library(MASS) にある関数 kde2d() を参照してみてください...
-kde2d と bkde2D について例を挙げてみました． -- &new{20...
-あらら…，２変量のヒストグラムについて書くのを忘れていま...
-関数 generator() を最適化して頂いたので，早速掲載しまし...
-[[Rコードの最適化:http://www.okada.jp.org/RWiki/?R%A5%B3...
#comment
~
*真の密度とデータの準備 [#yb7a0eb5]
**密度 f(x) = 0.6φ(x)+0.4ψ(x) [#s5aac18f]
φを平均-1，分散1の正規分布に従う確率変数，ψを平均2，分散1...
f(x) = 0.6φ(x)+0.4ψ(x)
となる密度関数 truedensity() を定義する．
truedensity <- function (x) {
0.6/sqrt(2*pi)*exp(-(x+1)^2/2) + 0.4/sqrt(2*pi)*exp(-(...
}
> curve(truedensity, xlim=c(-6,6), ylim=c(0,0.3), col=2)
#ref(ヒストグラムと密度の推定/hist-00.png, center)
**f(x) に従う乱数を生成する関数 [#xc7aa680]
次に， f(x) に従う乱数を生成する関数を定義する．
generator_tmp <- function(n) {
data1 <- rnorm(n)-1 # φ(x) に従う乱数
data2 <- rnorm(n)+2 # ψ(x) に従う乱数
data3 <- (runif(n) <= 0.6) # 確率0.6でdata1を採択し，
data1*data3+data2*(1-data3) # 確率0.4でdata2を採択する...
}
上と同じ動作をする関数で，最適化された関数を以下に挙げる...
generator <- function(n) {
(rnorm(n)-1)+3*(runif(n) > 0.6)
}
> data <- generator(1000)
【[[参考:http://www.okada.jp.org/RWiki/?R%A5%B3%A1%BC%A5%...
f(x) = p*φ(x) + (1-p)*ψ(x)
に従う乱数を生成する関数 rmixnorm4() を定義すると以下のよ...
rmixnorm4 <- function(n, p1, N1, N2) {
n1 <- sum(runif(n) < p1)
c(rnorm(n1, mean=N1[1], sd=N1[2]), rnorm(n-n1, mean=N...
}
*ヒストグラムの作成 [#x4565714]
**単純なヒストグラム：Sturges (1926) の方法 [#n3c9cc0e]
関数 hist() でヒストグラムを表示することも出来る．区切り...
> data <- generator(1000)
> hist(data)
#ref(ヒストグラムと密度の推定/hist-01.png, center)
『適当に選択』するが，適切な選択とは云えない．上のグラフ...
> library(MASS)
> x <- generator(1000)
> truehist(x)
#ref(ヒストグラムと密度の推定/hist-03.png, center)
----
追加：ユーザが自分で階級数を適切に（！）指定してやれば tr...
The default for breaks is "Sturges": see nclass.Sturges. ...
ということになっているのだが。
> set.seed(777)
> x <- generator(1000)
> layout(matrix(1:2, 2))
> library(MASS)
> truehist(x)
> hist(x, nclass=20, freq=FALSE)
> layout(1)
#ref(falsehist.png)
----
**プラグイン法によるヒストグラム：Wandの方法 [#x16f8b96]
次に，パッケージ ''KernSmooth'' にある関数 dpih() を用い...
> library(KernSmooth)
KernSmooth 2.22 installed
Copyright M. P. Wand 1997
> x <- generator(1000)
> h <- dpih(x)
> bins <- seq(min(x)-0.1, max(x)+0.1+h, by=h)
> hist(x, breaks=bins)
#ref(ヒストグラムと密度の推定/hist-02.png, center)
ところで，データに対応する x 軸上の点に縦線で印を付ける場...
> data(faithful)
> attach(faithful)
> plot(density(eruptions, bw=0.15))
> rug(eruptions)
> rug(jitter(eruptions, amount = .01), side = 3, col = "...
> detach("faithful")
#ref(ヒストグラムと密度の推定/hist04.png, center)
**Average Shifted Histogram [#u14adf44]
ヒストグラムには区切り幅選択の他にもう一つの重大な欠陥が...
> library(KernSmooth)
> x <- generator(100)
> h <- dpih(x)
> bins <- seq(min(x)-2*h, max(x)+h, by=h)
> hist(x, breaks=bins)
#ref(ヒストグラムと密度の推定/hist05.png, center)
上の例において，端点を (区切り幅)/2 だけ右にずらしてみる．
> bins <- seq(min(x)-3*h/2, max(x)+3*h/2, by=h)
> hist(x, breaks=bins)
#ref(ヒストグラムと密度の推定/hist06.png, center)
一つの改善案として，いくつかヒストグラムを描いたものの平...
# my Average Shifted Hi...
myASH <- function(data,breaks) { # Scott (1992) ( ≒ Haed...
library(KernSmooth) # パッケージの呼び出し
h <- dpih(data) # 区切り幅の推定
delta <- h/breaks # ヒストグラムをズラす幅
binnum <- length(seq(min(data)-h-delta/2, max(data)+h+...
by=delta)) # ヒストグラムの棒の数
pts <- rep(0, binnum) # 生データを選別して入...
binnum <- c() # メモリを開放(されてい...
for (i in 1:breaks) { # breaks 回くり返し
tmpbins <- seq(min(data)-i*delta, max(data)+(breaks-...
by=h) # データをそれぞれの区...
tmppts <- hist(data, breaks=tmpbins, prob=T, plot=F...
for (j in 0:(length(tmppts)-1)) { # 選別は合計 break...
pts[breaks*j+i] <- tmppts[j+1] # それぞれの結果は...
} # ズラされて配列に...
}
bins <- seq(min(data)-3*h/2-delta, max(data)+3*h/2+del...
pts <- append(rep(0,breaks), pts) # bin : plot の x ...
pts <- append(pts, rep(0,breaks)) # pts : 両端を 0 ...
meanpts <- rep(0, length(pts)-breaks+1)
for (i in 1:(length(pts)-breaks+1)) {
for (j in 1:breaks) {
meanpts[i] <- meanpts[i]+ pts[i+j-1]
} # meanpts : plot ...
meanpts[i] <- meanpts[i]/breaks # 周辺 breaks 個の...
} # の y 座標となる
plot(bins, meanpts, type="l") # 結果を plot する
}
> myASH(x,5) # 5 つの histogram...
#ref(ヒストグラムと密度の推定/hist07.png, center)
*密度の推定 [#zcc59f1e]
前述のデータについてヒストグラムを描いたが，デコボコすぎ...
**カーネル推定 [#r04a6f52]
実は，myASH() の引数 breaks を無限大に飛ばすとカーネル推...
> data <- generator(1000)
> plot(density(data), xlim=c(-6,6), ylim=c(0,0.3))
> lines(density(data), xlim=c(-6,6), ylim=c(0,0.3)) # ...
>par(new=T)
curve(truedensity, xlim=c(-6,6), ylim=c(0,0.3), col=2)
#ref(ヒストグラムと密度の推定/kern-00.png, center)
関数 density() の書式は以下の通りである．
density(x, bw = "nrd0", adjust = 1,
kernel = c("gaussian", "epanechnikov", "rectangu...
"biweight", "cosine", "optcosine"),
window = kernel, width, give.Rkern = FALSE,
n = 512, from, to, cut = 3, na.rm = FALSE)
引数の説明は以下の通り．ただし R のバージョンによって異な...
で確認することをお勧めする．
-x ： 観測値のベクトル．この観測値の従う確率分布の密度を...
-bw ： 用いられる平滑化バンド幅を指定する．指定された(も...
--nrd0 ： 平滑化カーネル（ガウスカーネル）の標準偏差にな...
--nrd ： nrd0 を Scott (1992) の方法により一般的したもの...
--ucv ： バイアス無しのクロスバリデーション規準によりバン...
--bcv ： バイアス付きのクロスバリデーション規準によりバン...
--SJ-ste ： パイロット評価を使用して（方程式を解くことに...
--SJ-dpi ： パイロット評価を使用して（プラグイン法により...
-adjust ： "bw" で指定されたバンド幅は，実際には adjust*b...
-kernel ，window ： 推定に用いるウインドウを指定する文字...
（注）"cosine"は S の命令であるが，これは "optcosine" よ...
-width ： S 言語との互換性の為に存在する引数で，云うなれ...
-give.Rkern ： これを TRUE にすれば密度推定が行われず，代...
--x ： 密度関数が推定される n 個の点の座標値
--y ： 推定された密度関数値
--bw ： 用いられるバンド幅
--N ： 欠損値を除いた後の標本の大きさ（サンプルサイズ）
--call ： 結果を生み出した呼出し
--data.name ： 引数 x の deparsed name
--has.na ： 互換性のための論理値（常に FALSE ）
-n ： 密度関数の値を求める等間隔点の数．
-from ，to ： 密度を from から to までの区間を n-1 等分し...
-cut ： 推定された密度を，両端において 0 に下がらせるため...
-na.rm ： これを TRUE にすれば欠損値が ｘ から取り除かれ...
前で定義した generator() で f(x) に従う乱数を 1000 個生成...
（バンド幅は "SJ" により選択）．真の密度を赤，density() ...
過ぎる嫌いがある．滑らかになり過ぎないようにするにはbcvの...
> data <- generator(1000)
> plot(density(data,bw="SJ"), xlab="x", ylab="y", xlim=c...
> truedensity <- function (x) {
+ 0.6/sqrt(2*pi)*exp(-(x+1)^2/2)+ 0.4/sqrt(2*pi)*exp(-...
+ }
> curve(truedensity, xlim=c(-6,6), ylim=c(0,0.3), col=2)
> par(new=T)
> plot(density(data,bw="SJ"), xlim=c(-6,6), ylim=c(0,0.3...
#ref(ヒストグラムと密度の推定/kern-01.png, center)
**プラグイン法による密度推定 [#q0caa129]
パッケージ ''KernSmooth'' にある関数 dpik() を用いること...
> library(KernSmooth)
> h <- dpik(data)
> est <- bkde(data, bandwidth=h)
> plot(est, type="l", xlim=c(-6,6), ylim=c(0,0.3))
> par(new=T)
> curve(truedensity, xlim=c(-6,6), ylim=c(0,0.3), col=2)
#ref(ヒストグラムと密度の推定/kern-02.png, center)
*二変量データのヒストグラムと密度推定 [#ta85f364]
**二変量データのヒストグラム：イメージ風 [#hc6b12fd]
二変量データをイメージ的なヒストグラムで表すには以下のよ...
> library(fields)
> look<- image.count( precip$x, nrow=32, ncol=32)
> image.plot( look)
#ref(ヒストグラムと密度の推定/hist-image.png, center)
**二変量データのヒストグラム：marginal [#w52316be]
二変量データの marginal をヒストグラムで表すには以下のよ...
> library(ade4)
> data(rpjdl)
> coa1 <- dudi.coa(rpjdl$fau, scannf = FALSE, nf = 4)
> s.hist(coa1$li)
> s.hist(coa1$li, cgrid = 2, cbr = 3, adj = 0.5, clab = 0)
#ref(ヒストグラムと密度の推定/hist-marginal.png, center)
**二変量データのヒストグラム：頻度ポリゴン [#xa09f59d]
二変量データのヒストグラム(頻度ポリゴン)を描く場合は以下...
> library(gregmisc)
> # データ：無相関な2変量正規乱数
> x <- rnorm(2000, sd=4)
> y <- rnorm(2000, sd=1)
> # 遠近法プロット (persp) のためのデータを hist2d() を...
> h2d <- hist2d(x,y,show=FALSE, same.scale=TRUE, nbins=c...
> persp( h2d$x, h2d$y, h2d$counts,
+ ticktype="detailed", theta=60, phi=30,
+ expand=0.5, shade=0.5, col="cyan", ltheta=-30)
#ref(ヒストグラムと密度の推定/hist-polygon.png, center)
**二変量データのカーネル密度推定 [#bbb37576]
二変量データのカーネル密度推定を行う場合は関数 kde2d() や...
> library(MASS)
> data(geyser)
> attach(geyser)
> f1 <- kde2d(duration, waiting, n = 50, lims = c(0.5, 6...
> f2 <- kde2d(duration, waiting, n = 50, lims = c(0.5, 6...
+ h = c(width.SJ(duration), width.SJ(waiting...
> persp(f2, phi = 30, theta = 20, d = 5)
#ref(ヒストグラムと密度の推定/kern-04.png, center)
> library(KernSmooth)
KernSmooth 2.22 installed
Copyright M. P. Wand 1997
> data(geyser, package="MASS")
> x <- cbind(geyser$duration, geyser$waiting)
> est <- bkde2D(x, bandwidth=c(0.7,7))
> persp(est$fhat)
#ref(ヒストグラムと密度の推定/kern-03.png, center)
*カーネル密度推定のわかりやすい文献(なんでも掲示板2003-10...
R言語は多くの種類のカーネル密度推定法をサポートしている．
この手法をわかりやすく解説した書物は以下のようなものがあ...
-Haerdle, W. (1991) : Smoothing Techniques with Implement...
-A. Bowman and A. Azzalini(1997): Applied smoothing techn...
-W.N. Venables and B.D. Ripley(2002):Modern applied stati...
-B. W. Silverman(1986): Density estimation for statistics...
-ジェフリー S. シモノフ 著，竹澤邦夫・大森宏 訳(1999) 「...
-竹澤邦夫 著(2003) みんなのためのノンパラメトリック回帰(...
*lowess() による平滑化 [#g4e374c0]
関数 lowess() によってデータの平滑化を行うことも出来る．
lowess() は平滑結果の座標を与える成分 x と y のリストを返...
平滑化は関数 lines() を用いて元の散布図に追加することが出...
書式は以下の通り．
lowess(x, y, f=2/3, iter=3, delta=.01*diff(range(x)))
引数の説明は以下の通り．
-x ，y ： 散布図中のプロットの座標を与えるベクトル．単一...
-f ： 平滑幅．これは各位置での平滑に影響を及ぼすプロット...
-iter ： 実行されるべき頑健化繰り返しの数．小さな iter の...
-delta ： 互いに距離 delta 以内に位置する x の値は lowess...
以下に使用例を示す．
> data(cars)
> plot(cars, main = "lowess(cars)")
> lines(lowess(cars), col = 2)
> lines(lowess(cars, f = 0.2), col = 3)
> legend(5, 120, c(paste("f = ", c("2/3", ".2"))), lty =...
#ref(ヒストグラムと密度の推定/smooth-00.png, center)
終了行:
// 舟尾
COLOR(red){SIZE(25){ヒストグラムと密度の推定}}
([[グラフィックス参考実例集]]に戻る。[[Rのグラフィックス...
#contents
~
-混合分布からのサンプリングコードは鮮やかですね。data3 <-...
-修正ありがとうございます。ついでにコメントも入れました。...
-多(2)変量密度推定に関するパッケージや関数についての情報...
-library(MASS) にある関数 kde2d() を参照してみてください...
-kde2d と bkde2D について例を挙げてみました． -- &new{20...
-あらら…，２変量のヒストグラムについて書くのを忘れていま...
-関数 generator() を最適化して頂いたので，早速掲載しまし...
-[[Rコードの最適化:http://www.okada.jp.org/RWiki/?R%A5%B3...
#comment
~
*真の密度とデータの準備 [#yb7a0eb5]
**密度 f(x) = 0.6φ(x)+0.4ψ(x) [#s5aac18f]
φを平均-1，分散1の正規分布に従う確率変数，ψを平均2，分散1...
f(x) = 0.6φ(x)+0.4ψ(x)
となる密度関数 truedensity() を定義する．
truedensity <- function (x) {
0.6/sqrt(2*pi)*exp(-(x+1)^2/2) + 0.4/sqrt(2*pi)*exp(-(...
}
> curve(truedensity, xlim=c(-6,6), ylim=c(0,0.3), col=2)
#ref(ヒストグラムと密度の推定/hist-00.png, center)
**f(x) に従う乱数を生成する関数 [#xc7aa680]
次に， f(x) に従う乱数を生成する関数を定義する．
generator_tmp <- function(n) {
data1 <- rnorm(n)-1 # φ(x) に従う乱数
data2 <- rnorm(n)+2 # ψ(x) に従う乱数
data3 <- (runif(n) <= 0.6) # 確率0.6でdata1を採択し，
data1*data3+data2*(1-data3) # 確率0.4でdata2を採択する...
}
上と同じ動作をする関数で，最適化された関数を以下に挙げる...
generator <- function(n) {
(rnorm(n)-1)+3*(runif(n) > 0.6)
}
> data <- generator(1000)
【[[参考:http://www.okada.jp.org/RWiki/?R%A5%B3%A1%BC%A5%...
f(x) = p*φ(x) + (1-p)*ψ(x)
に従う乱数を生成する関数 rmixnorm4() を定義すると以下のよ...
rmixnorm4 <- function(n, p1, N1, N2) {
n1 <- sum(runif(n) < p1)
c(rnorm(n1, mean=N1[1], sd=N1[2]), rnorm(n-n1, mean=N...
}
*ヒストグラムの作成 [#x4565714]
**単純なヒストグラム：Sturges (1926) の方法 [#n3c9cc0e]
関数 hist() でヒストグラムを表示することも出来る．区切り...
> data <- generator(1000)
> hist(data)
#ref(ヒストグラムと密度の推定/hist-01.png, center)
『適当に選択』するが，適切な選択とは云えない．上のグラフ...
> library(MASS)
> x <- generator(1000)
> truehist(x)
#ref(ヒストグラムと密度の推定/hist-03.png, center)
----
追加：ユーザが自分で階級数を適切に（！）指定してやれば tr...
The default for breaks is "Sturges": see nclass.Sturges. ...
ということになっているのだが。
> set.seed(777)
> x <- generator(1000)
> layout(matrix(1:2, 2))
> library(MASS)
> truehist(x)
> hist(x, nclass=20, freq=FALSE)
> layout(1)
#ref(falsehist.png)
----
**プラグイン法によるヒストグラム：Wandの方法 [#x16f8b96]
次に，パッケージ ''KernSmooth'' にある関数 dpih() を用い...
> library(KernSmooth)
KernSmooth 2.22 installed
Copyright M. P. Wand 1997
> x <- generator(1000)
> h <- dpih(x)
> bins <- seq(min(x)-0.1, max(x)+0.1+h, by=h)
> hist(x, breaks=bins)
#ref(ヒストグラムと密度の推定/hist-02.png, center)
ところで，データに対応する x 軸上の点に縦線で印を付ける場...
> data(faithful)
> attach(faithful)
> plot(density(eruptions, bw=0.15))
> rug(eruptions)
> rug(jitter(eruptions, amount = .01), side = 3, col = "...
> detach("faithful")
#ref(ヒストグラムと密度の推定/hist04.png, center)
**Average Shifted Histogram [#u14adf44]
ヒストグラムには区切り幅選択の他にもう一つの重大な欠陥が...
> library(KernSmooth)
> x <- generator(100)
> h <- dpih(x)
> bins <- seq(min(x)-2*h, max(x)+h, by=h)
> hist(x, breaks=bins)
#ref(ヒストグラムと密度の推定/hist05.png, center)
上の例において，端点を (区切り幅)/2 だけ右にずらしてみる．
> bins <- seq(min(x)-3*h/2, max(x)+3*h/2, by=h)
> hist(x, breaks=bins)
#ref(ヒストグラムと密度の推定/hist06.png, center)
一つの改善案として，いくつかヒストグラムを描いたものの平...
# my Average Shifted Hi...
myASH <- function(data,breaks) { # Scott (1992) ( ≒ Haed...
library(KernSmooth) # パッケージの呼び出し
h <- dpih(data) # 区切り幅の推定
delta <- h/breaks # ヒストグラムをズラす幅
binnum <- length(seq(min(data)-h-delta/2, max(data)+h+...
by=delta)) # ヒストグラムの棒の数
pts <- rep(0, binnum) # 生データを選別して入...
binnum <- c() # メモリを開放(されてい...
for (i in 1:breaks) { # breaks 回くり返し
tmpbins <- seq(min(data)-i*delta, max(data)+(breaks-...
by=h) # データをそれぞれの区...
tmppts <- hist(data, breaks=tmpbins, prob=T, plot=F...
for (j in 0:(length(tmppts)-1)) { # 選別は合計 break...
pts[breaks*j+i] <- tmppts[j+1] # それぞれの結果は...
} # ズラされて配列に...
}
bins <- seq(min(data)-3*h/2-delta, max(data)+3*h/2+del...
pts <- append(rep(0,breaks), pts) # bin : plot の x ...
pts <- append(pts, rep(0,breaks)) # pts : 両端を 0 ...
meanpts <- rep(0, length(pts)-breaks+1)
for (i in 1:(length(pts)-breaks+1)) {
for (j in 1:breaks) {
meanpts[i] <- meanpts[i]+ pts[i+j-1]
} # meanpts : plot ...
meanpts[i] <- meanpts[i]/breaks # 周辺 breaks 個の...
} # の y 座標となる
plot(bins, meanpts, type="l") # 結果を plot する
}
> myASH(x,5) # 5 つの histogram...
#ref(ヒストグラムと密度の推定/hist07.png, center)
*密度の推定 [#zcc59f1e]
前述のデータについてヒストグラムを描いたが，デコボコすぎ...
**カーネル推定 [#r04a6f52]
実は，myASH() の引数 breaks を無限大に飛ばすとカーネル推...
> data <- generator(1000)
> plot(density(data), xlim=c(-6,6), ylim=c(0,0.3))
> lines(density(data), xlim=c(-6,6), ylim=c(0,0.3)) # ...
>par(new=T)
curve(truedensity, xlim=c(-6,6), ylim=c(0,0.3), col=2)
#ref(ヒストグラムと密度の推定/kern-00.png, center)
関数 density() の書式は以下の通りである．
density(x, bw = "nrd0", adjust = 1,
kernel = c("gaussian", "epanechnikov", "rectangu...
"biweight", "cosine", "optcosine"),
window = kernel, width, give.Rkern = FALSE,
n = 512, from, to, cut = 3, na.rm = FALSE)
引数の説明は以下の通り．ただし R のバージョンによって異な...
で確認することをお勧めする．
-x ： 観測値のベクトル．この観測値の従う確率分布の密度を...
-bw ： 用いられる平滑化バンド幅を指定する．指定された(も...
--nrd0 ： 平滑化カーネル（ガウスカーネル）の標準偏差にな...
--nrd ： nrd0 を Scott (1992) の方法により一般的したもの...
--ucv ： バイアス無しのクロスバリデーション規準によりバン...
--bcv ： バイアス付きのクロスバリデーション規準によりバン...
--SJ-ste ： パイロット評価を使用して（方程式を解くことに...
--SJ-dpi ： パイロット評価を使用して（プラグイン法により...
-adjust ： "bw" で指定されたバンド幅は，実際には adjust*b...
-kernel ，window ： 推定に用いるウインドウを指定する文字...
（注）"cosine"は S の命令であるが，これは "optcosine" よ...
-width ： S 言語との互換性の為に存在する引数で，云うなれ...
-give.Rkern ： これを TRUE にすれば密度推定が行われず，代...
--x ： 密度関数が推定される n 個の点の座標値
--y ： 推定された密度関数値
--bw ： 用いられるバンド幅
--N ： 欠損値を除いた後の標本の大きさ（サンプルサイズ）
--call ： 結果を生み出した呼出し
--data.name ： 引数 x の deparsed name
--has.na ： 互換性のための論理値（常に FALSE ）
-n ： 密度関数の値を求める等間隔点の数．
-from ，to ： 密度を from から to までの区間を n-1 等分し...
-cut ： 推定された密度を，両端において 0 に下がらせるため...
-na.rm ： これを TRUE にすれば欠損値が ｘ から取り除かれ...
前で定義した generator() で f(x) に従う乱数を 1000 個生成...
（バンド幅は "SJ" により選択）．真の密度を赤，density() ...
過ぎる嫌いがある．滑らかになり過ぎないようにするにはbcvの...
> data <- generator(1000)
> plot(density(data,bw="SJ"), xlab="x", ylab="y", xlim=c...
> truedensity <- function (x) {
+ 0.6/sqrt(2*pi)*exp(-(x+1)^2/2)+ 0.4/sqrt(2*pi)*exp(-...
+ }
> curve(truedensity, xlim=c(-6,6), ylim=c(0,0.3), col=2)
> par(new=T)
> plot(density(data,bw="SJ"), xlim=c(-6,6), ylim=c(0,0.3...
#ref(ヒストグラムと密度の推定/kern-01.png, center)
**プラグイン法による密度推定 [#q0caa129]
パッケージ ''KernSmooth'' にある関数 dpik() を用いること...
> library(KernSmooth)
> h <- dpik(data)
> est <- bkde(data, bandwidth=h)
> plot(est, type="l", xlim=c(-6,6), ylim=c(0,0.3))
> par(new=T)
> curve(truedensity, xlim=c(-6,6), ylim=c(0,0.3), col=2)
#ref(ヒストグラムと密度の推定/kern-02.png, center)
*二変量データのヒストグラムと密度推定 [#ta85f364]
**二変量データのヒストグラム：イメージ風 [#hc6b12fd]
二変量データをイメージ的なヒストグラムで表すには以下のよ...
> library(fields)
> look<- image.count( precip$x, nrow=32, ncol=32)
> image.plot( look)
#ref(ヒストグラムと密度の推定/hist-image.png, center)
**二変量データのヒストグラム：marginal [#w52316be]
二変量データの marginal をヒストグラムで表すには以下のよ...
> library(ade4)
> data(rpjdl)
> coa1 <- dudi.coa(rpjdl$fau, scannf = FALSE, nf = 4)
> s.hist(coa1$li)
> s.hist(coa1$li, cgrid = 2, cbr = 3, adj = 0.5, clab = 0)
#ref(ヒストグラムと密度の推定/hist-marginal.png, center)
**二変量データのヒストグラム：頻度ポリゴン [#xa09f59d]
二変量データのヒストグラム(頻度ポリゴン)を描く場合は以下...
> library(gregmisc)
> # データ：無相関な2変量正規乱数
> x <- rnorm(2000, sd=4)
> y <- rnorm(2000, sd=1)
> # 遠近法プロット (persp) のためのデータを hist2d() を...
> h2d <- hist2d(x,y,show=FALSE, same.scale=TRUE, nbins=c...
> persp( h2d$x, h2d$y, h2d$counts,
+ ticktype="detailed", theta=60, phi=30,
+ expand=0.5, shade=0.5, col="cyan", ltheta=-30)
#ref(ヒストグラムと密度の推定/hist-polygon.png, center)
**二変量データのカーネル密度推定 [#bbb37576]
二変量データのカーネル密度推定を行う場合は関数 kde2d() や...
> library(MASS)
> data(geyser)
> attach(geyser)
> f1 <- kde2d(duration, waiting, n = 50, lims = c(0.5, 6...
> f2 <- kde2d(duration, waiting, n = 50, lims = c(0.5, 6...
+ h = c(width.SJ(duration), width.SJ(waiting...
> persp(f2, phi = 30, theta = 20, d = 5)
#ref(ヒストグラムと密度の推定/kern-04.png, center)
> library(KernSmooth)
KernSmooth 2.22 installed
Copyright M. P. Wand 1997
> data(geyser, package="MASS")
> x <- cbind(geyser$duration, geyser$waiting)
> est <- bkde2D(x, bandwidth=c(0.7,7))
> persp(est$fhat)
#ref(ヒストグラムと密度の推定/kern-03.png, center)
*カーネル密度推定のわかりやすい文献(なんでも掲示板2003-10...
R言語は多くの種類のカーネル密度推定法をサポートしている．
この手法をわかりやすく解説した書物は以下のようなものがあ...
-Haerdle, W. (1991) : Smoothing Techniques with Implement...
-A. Bowman and A. Azzalini(1997): Applied smoothing techn...
-W.N. Venables and B.D. Ripley(2002):Modern applied stati...
-B. W. Silverman(1986): Density estimation for statistics...
-ジェフリー S. シモノフ 著，竹澤邦夫・大森宏 訳(1999) 「...
-竹澤邦夫 著(2003) みんなのためのノンパラメトリック回帰(...
*lowess() による平滑化 [#g4e374c0]
関数 lowess() によってデータの平滑化を行うことも出来る．
lowess() は平滑結果の座標を与える成分 x と y のリストを返...
平滑化は関数 lines() を用いて元の散布図に追加することが出...
書式は以下の通り．
lowess(x, y, f=2/3, iter=3, delta=.01*diff(range(x)))
引数の説明は以下の通り．
-x ，y ： 散布図中のプロットの座標を与えるベクトル．単一...
-f ： 平滑幅．これは各位置での平滑に影響を及ぼすプロット...
-iter ： 実行されるべき頑健化繰り返しの数．小さな iter の...
-delta ： 互いに距離 delta 以内に位置する x の値は lowess...
以下に使用例を示す．
> data(cars)
> plot(cars, main = "lowess(cars)")
> lines(lowess(cars), col = 2)
> lines(lowess(cars, f = 0.2), col = 3)
> legend(5, 120, c(paste("f = ", c("2/3", ".2"))), lty =...
#ref(ヒストグラムと密度の推定/smooth-00.png, center)
ページ名:
![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
