[基本パッケージ中の統計関連関数の全体をカテゴリー別に紹介すること，およびそのヘルプドキュメントの「ほぼ忠実な」和訳(ちょっと詳しくなった人のためのもの)]
|
|
|
|
● インストール：The Comprehensive R Archive Network(CRAN https://cran.r-project.org/)からダウンロード (Windows ではR-x.x.x-win.exe のファイル)。詳しくは Rの初歩(奥村晴彦先生)
● パッケージのインストール(とりあえずは不要)
|
|
● インストール
|
|
Rを
インストールの上、
Windows スタートボタン-> R x64 で検索してRGuiを起動
●
実行
● クリアの方法
rm(list = ls(all.names = TRUE)) # オブジェクトを全て削除する if( dev.cur() > 1 ) dev.off() # プロット欄の図をクリア● エディタ
|
|
● 数値の代入
a <- 256; a # a の中に256 を代入して a を表示 a <- a + 1; a # a をインクリメントして a を表示 (b <- 123) # 括弧をつけると代入したものを表示してくれる● データの基本
if (条件){
式1
} else{
式2
}
●
for 文
for(i in 1:5) cat(i,"\n")
|
|
● データの読み込みと確認
国語 = c(46, 30, 24, 49) # 国語のデータの入力 数学 = c(71, 72, 4, 49) # 数学のデータの入力 tokuten = data.frame(国語,数学) # データを一つのオブジェクトとして扱う tokuten # 各の成分表示 tokuten$国語 # 国語だけをとり出す tokuten$理科 = c(98, 35, 21, 87) # 理科のデータの追加 tokuten # オブジェクトの成分を確認(理科が追加)
x <- 1:10
y <- 1:10
## pdf 形式で表示 pictex(), postscript(), png() もある
pdf("c:/temp/xy.pdf", family="Japan1") # 日本語フォントの字化け防止
plot(x,y) # 実際の描画
dev.off() # 描画のデバイスを閉じる(これでファイルに保存される)
|
|
## 関数名<-function(引数) {return(返り値) } の形式
wa <- function(x, y) {
wa1 <- x + y # 足し算
return(wa1)
}
wa(10,4) # 使い方
tanjyobi <- function(n){
year<-365 # 平年
prob = 1 - prod(seq(from = year, by = -1, length.out = n) / year)
return(prob)
}
### 関数の応用例
ninzu <- 40; tanjyobi(ninzu) # 40人いる場合の誕生日が一致する人がいる確率
flag <- 0 # 初めて1/2を超えるための人数(停止時間)
for(i in 1:ninzu) {
cat(sprintf(" %2d人 一致確率 =%.2f ", i,tanjyobi(i)))
if(tanjyobi(i)>0.5 && flag == 0) {
cat(sprintf(" %2d人あつまれば一致する人がいる確率が1/2より大きい！",i))
flag <- 1
}
cat("\n")
}
n <-1:(ninzu+10) # 人数のデータ
p <-lapply(n, tanjyobi) # lapply 関数で一致確率をリストにする
plot(n,p, type="o",
main="誕生日の一致確率",
xlab="人数", ylab="一致確率")
|
|
データx, y, xy を定義 x <- c(67.2, 87.3, 77.2, 78.3) # xにデータが入っているとする y <- c(21.2, 98.6, 78.1, 91.6) # yにデータが入っているとする z <- c(31.8, 91.8, 54.1, 77.1) # zにデータが入っているとする xy = data.frame(x,y) # データをオブジェクトとして扱う データx, y, xy の情報を表示 sum(x) # 合計 mean(x) # 平均 length(x) # データのサイズ max(x) # 最大値 min(x) # 最小値 median(x) # 中央値 range(x) # 最小値と最大値 fivenum(x) # Tukey の五数要約 five number summary (最小値、下ヒンジ、中央値、上ヒンジ、最大値) quantile(x) # 四分位数 type で指定できる。詳しくは?quantileで調べる! IQR(x) # 四分位偏差 sd(x) # 標準偏差 var(x) # 不偏分散 (n-1で割る) var(x)*(length(x)-1)/length(x) # 標本分散(nで割る) cov(x, y) # x,yの共分散 cor(x, y) # x,yの相関係数 sort(x) sort(x, decreasing=T)
ヒストグラム #hist(x) #dev.off() # ウィンドウを閉じる hist(x, breaks=seq(60,90,5), main="xの度数分布", xlab="階級", ylab="度数") # タイトル、x軸, y軸の名前 ヒストグラムを縦に表記 par(mfrow=c(2,1)) # 画面2x1にヒストグラム hist(xy$x) # オブジェクトの要素としての出力 hist(xy$y) # オブジェクトの要素としての出力 ヒストグラムの軸を設定 par(mfrow=c(1,2)) # 画面1x2にヒストグラム hist(kokusuu$国語, xlim=c(0,100),ylim=c(0,12)) hist(kokusuu$数学, xlim=c(0,100),ylim=c(0,12))
散布図 plot(x,y)
3次元の散布図
install.packages("scatterplot3d") # install していないときにはインストールしておく
library(scatterplot3d) # scatterplot3dをよびだす。
x <- c(67.2, 87.3, 77.2, 78.3) # xにデータが入っているとする
y <- c(21.2, 98.6, 78.1, 91.6) # yにデータが入っているとする
z <- c(31.8, 91.8, 54.1, 77.1) # zにデータが入っているとする
scatterplot3d(x,y,z)
箱ひげ図 #boxplot(x) boxplot(y) boxplot(xy) boxplot(x,xlab="", ylab="", col="red", horizontal = T) # 箱を赤、横向き
単独のグラフ max<- 5; min<- -max; mu = 0; sigma = 1; curve( dnorm(x, mu, sigma), from=min, to=max ) # 標準正規分布 curve(dt(x, 2), add=TRUE, lty="dotted") # 自由度2 のt分布を追加
|
|
library(rmutil) rlevy(10)
library(Pareto) plot(dPareto(1:10,1,1),type="o")
|
|
関数の定義と表示1 f<-function(x) 2*sin(x) #関数の定義 curve(f(x),-pi,pi) # グラフの描画
関数の定義と表示2
entropy.func <- function(x){
-x*log2(x)-(1-x)*log2(1-x)
}
plot(entropy.func, 0, 1, type="l") # plot(関数名, 下限, 上限)
分布関数の表示(場合分けでの定義の例)
a<-3; xmax <- 5; ymax <- .1;
My.exp <- function(x){
if( x>0 )
F=1-exp( -a*x )
else
F=0
return(F)
}
h<-Vectorize(My.exp) # 場合分けの際はベクトルにする
curve(h,xlim=c(-xmax,xmax),
ylim=c(-ymax,1+ymax),
main="指数分布の分布関数",
ylab="",
lwd = 2, # 太さが2
)
abline(h = 0);abline(v = 0);
abline(h = 1, lty="dotted", col="blue")
abline(v = a, lty="dotted", col="blue")
# グンベル分布の分布関数、密度関数の微分での表示
range<-6
f1 <-deriv(~exp(-exp(-x)), "x", func=T) # グンベル分布の分布関数
f2 <-function(x) attr(f1(x),"gradient") # グンベル分布の密度関数
curve(f1(x),-range,range,
main="グンベル分布の分布関数と密度関数", # 全体のラベル
xlab="", ylab="", # 軸のラベルはつけない
lty = 1, # 実線を指定
col="blue")
curve(f2(x),-range,range,
lty=2, # 点線を指定
col="red",
add=T)
abline(v = 0) # y 軸を描く
abline(h = 0) # x 軸を描く
legend("topleft", # 左上に凡例を書く
legend=c("分布関数", "密度関数"), # 凡例のコメント
lty=c(1,2), # 線の種類を指定
col=c("blue", "red") # 色を指定
)
対数パレート分布
正規分布の密度関数の重ね書き
ran <- 8;
mean1 <- 0; sd1 <- 1;
mean2 <- 1; sd2 <- 2;
curve(dnorm(x,mean1,sd1), -ran, ran,
main="2つの正規分布の密度関数",
ylab="確率",
add = TRUE, # 重ね書きOK
col="blue",
type="l")
curve(dnorm(x,mean2,sd2), -ran, ran,
add = TRUE, # 重ね書きOK
col="red",
lty="dotted", # 線をドットにする
lwd = 2, # 太さが2
type="l")
ポアソン分布の頻度関数
lambda <- 2; range1 <- 10;
barplot(dpois(0:range1,lambda), names.arg=0:range1,xlab="x", ylab="prob",
main = "Poisson dist")
n <- 10; p<- 1/6;
plot(0:n, dbinom(0:n, n, p), # 二項分布
xlab="試行回数",
ylab="確率",
main="二項分布の頻度関数",
type="h", # 各点からx軸までの垂線プロットを行う．
lwd=7) #線分の幅 (line width) を番号で指定する．値が大きいほど太い線になる
# 二項分布の頻度関数 Bin(n,p)の表記
# Benjamin, Fluet (2000), Amer. Math. Monthly, 107, 6, 560-562.
n <- 9
p <- 2/3
x <- 0:n
# 確率関数
prob_f <- dbinom(x, size = n, prob = p)
# ヒストグラム
bp <- barplot(
prob_f,
names.arg = x,
# col = "skyblue",
xlab = "成功回数",
ylab = "確率",
main = paste("二項分布 Bin(", n, ", ", p, ")", sep = ""),
ylim = c(0, max(prob_f) * 1.5)
)
# 棒の上に確率を表記
text(
x = bp, # 棒の中央位置（barplot が返す値）
y = prob_f, # 棒の高さ
labels = round(prob_f, 4), # ラベルの内容（小数第3位まで表示）
pos = 3, # 棒の上に表示（位置指定）
cex = 0.7, # 文字サイズ
col = "black" # ラベルの色
)
f <- function(p){
p^10-20*p^19+19*p^20
}
curve(f(x), main="秀をもらう確率の差",
xlab="確率", ylab="確率の差",
xlim=c(0,1))
abline(h=0) # x軸
abline(v=0.8, col='red', lwd=1, lty=2) # 縦線, 赤, 太さ1, 点線
abline(v=0.9, col='blue', lwd=1, lty=2) # 縦線, 赤, 太さ1, 点線
legend("topleft", # 左上に凡例を書く
legend=c("A君", "B君"), # 凡例のコメント
lty=c(2,2), # 線の種類を指定
col=c("red", "blue") # 色を指定
)
# 中田、内藤、p.115 例5.6 (対数尤度関数の表示)
n<-5; k<-3
logL <- function(p){
p^k*(1-p)^(n-k)
}
curve(logL(x), main="対数尤度関数",
xlab="確率", ylab="対数尤度",
xlim=c(0,1))
abline(v=k/n, col='red', lwd=1, lty=2) # 縦線, 赤, 太さ1, 点線
密度関数と乱数によるヒストグラムを同時に描く(中田、内藤, p.101, 章末問題4.19)
N <- 10000; x <- runif(N); y <- -log(x)
hist(y,prob=TRUE,
main="一様分布の-log x による変換と指数分布",
xlab="", ylab="確率")
rng <- 10; curve(exp(-x), 0, rng, add=TRUE)
裾確率に色を塗る(多角形とみなす) xrange = 4; # xの範囲 min <- 2; #塗りつぶす軸を指定 i <- 10 #台形に分けるメッシュの細かさ eps <-0.01 # 文字をずらすときの幅 mu <- 0; sd <- 1; #標準正規分布 curve(dnorm(x, mu, sd),xlab="", ylab="", -xrange, xrange) xseq <- seq(min, xrange, length=i) # x の範囲を等分割 yseq <- dnorm(xseq) # 対応するy の値 # データを多角形にする(第1,4座標を同じにすることで実現) x1 <- c(min, xseq, xrange, min) y1 <- c( 0, yseq, dnorm(xrange, mu, sd), 0) # 多角形のデータに色を塗る polygon(x1, y1, col="gray") # 境界線とその説明 lines(c(min, min), c(0, dnorm(min, mu, sd)+eps), lwd = 2) text(min, dnorm(min, mu, sd)+eps, paste(min,"σ")) abline(h=0) #積分範囲の矢印 arrows( mid, 0.3, mid, dnorm(mid, mu, sd)/2, col="blue") #積分範囲の説明 mid <- (min+xrange)/2; text(mid+eps, 0.3+eps, paste(min,"σ以上の確率"))
|
|
# p.121 例5.9 (区間推定, 正規母集団, 母分散既知)
kichi_shinraikukan <- function(n, sigma2, xbar, shinraido) {
z <- qnorm(shinraido) # z<-qnorm(p) <=> P(N(0,1) <= z)=p
sgm <- sqrt(sigma2) # 母標準偏差
haba <- z*sgm/sqrt(n) # 標本偏差からの幅
shitagawa <- xbar-haba # 下側信頼限界
uegawa <- xbar+haba # 上側信頼限界
cat(sprintf("### 母平均の信頼区間(母分散既知)\n"))
cat(sprintf("信頼区間の幅/2\t: %.3f\n", haba))
cat(sprintf("信頼区間\t: [%.3f, %.3f]\n", shitagawa, uegawa))
} # kichi_shinraikukanの定義終了
#kichi_shinraikukan(10, 0.05^2, 1.021, 0.975)#例5.9
kichi_shinraikukan(10, 0.02^2, 5.001, 0.975) #問5.10
# michi_shinraikukan (標本の大きさ, 標本不偏分散, 標本平均, 信頼係数)
michi_shinraikukan <- function(n, uvar, xbar, shinraikeisu) {
p <- shinraikeisu+(1-shinraikeisu)/2
t <- qt(p, n-1) # t<-qt(p) <=> P(t(n-1) <= t)=p
haba <- t*sqrt(uvar)/sqrt(n) # 標本偏差からの幅
shitagawa <- xbar-haba # 下側信頼限界
uegawa <- xbar+haba # 上側信頼限界
cat(sprintf("### 母平均の信頼区間(母分散未知)\n"))
cat(sprintf("パーセント点\t: %.3f\n", t))
cat(sprintf("標本平均\t: %.3f\n", xbar))
cat(sprintf("頼区間の幅/2\t: %.3f\n", haba))
cat(sprintf("信頼区間\t: [%.2f, %.2f]\n", shitagawa, uegawa))
} # michi_shinraikukanの定義終了
#michi_shinraikukan(30, 145.36, 29.25, 0.95)#例5.11
##### 実際のデータでの母平均の信頼区間
toi5.12 = c(21,20.2,24,18.8,16.4,10.3,16.4,15.8,17.4,15.4,27.2,30.8,28.4,20,
22.2,22.6,20.6,26.4,25,24.8,34.8)
n <- length(toi5.12); uvar <- var(toi5.12); xbar <- mean(toi5.12) # 代入
michi_shinraikukan(n, uvar, xbar, 0.95)#問5.12
## R の t.test を使用すれば関数を定義する必要なし!
t.test(toi5.12,conf.level=0.95) # t検定を行うとともに95%信頼区間を出力
出力の一部
t = 17.254, df = 20, p-value = 1.773e-13 # t の実現値, 自由度, p値 alternative hypothesis: true mean is not equal to 0 # 両側対立仮説 95 percent confidence interval: 19.19378 24.47289 # 信頼区間 sample estimates: mean of x 21.83333 # 標本平均
# bobunsan_shinraikukan(標本の大きさ, 標本不偏分散, 信頼係数)
bobunsan_shinraikukan <- function(n, uvar, shinraikeisu) {
p <- shinraikeisu+(1-shinraikeisu)/2
chi2low <- qchisq(1-p, n-1) # c <-qchisq(p) <=> P(chisq(n-1) <= c)=p
chi2upp <- qchisq(p , n-1) # c <-qchisq(p) <=> P(chisq(n-1) <= c)=p
shitagawa <- (n-1)*uvar/chi2upp # 下側信頼限界
uegawa <- (n-1)*uvar/chi2low # 上側信頼限界
cat(sprintf("### 母分散の信頼区間\n"))
cat(sprintf("chi2下側\t: %.2f\nchi2上側\t: %.2f\n", chi2low, chi2upp))
cat(sprintf("信頼区間\t: [%.3f, %.3f]\n", shitagawa, uegawa))
} # bobunsan_shinraikukanの定義終了
bobunsan_shinraikukan(30,145.36,0.95)
toi5.12 = c(21,20.2,24,18.8,16.4,10.3,16.4,15.8,17.4,15.4,27.2,30.8,28.4,20,
22.2,22.6,20.6,26.4,25,24.8,34.8)
n <- length(toi5.12); uvar <- var(toi5.12)
bobunsan_shinraikukan(n,uvar,0.95)#問5.13
bobunsan_shinraikukan(n,uvar,0.90)#問5.13
############# R の関数として以下のみでOK
library(EnvStats) # パッケージを読み込む(初めての場合はインストールする)
varTest(toi5.12, alternative = "two.sided", conf.level=0.95) # カイ二乗分散検定
varTest(toi5.12, alternative = "two.sided", conf.level=0.90) # カイ二乗分散検定
出力の一部
Chi-Squared = 672.51, df = 20, p-value < 2.2e-16 # chi2 の実現値, 自由度, p値 alternative hypothesis: true variance is not equal to 1 # 両側対立仮説 95 percent confidence interval: 19.68143 70.12014 # 信頼区間 sample estimates: variance 33.62533 # 標本不偏分散
# hiritsu_shinraikukan (標本の大きさ, 成功回数, 信頼係数)
hiritsu_shinraikukan <- function(n, ncount, shinraikeisu) {
p <- shinraikeisu+(1-shinraikeisu)/2
z <- qnorm(p) # z<-qnorm(p) <=> P(N(0,1) <= z)=p
phat <- ncount/n # p の点推定
haba <- z*sqrt(phat*(1-phat)/n) # 標本偏差からの幅
shitagawa <- phat - haba # 下側信頼限界
uegawa <- phat + haba # 上側信頼限界
cat(sprintf("### 比率の信頼区間\n"))
cat(sprintf("信頼区間の幅/2\t: %.3f\n", haba))
cat(sprintf("信頼区間\t: [%.3f, %.3f]\n", shitagawa, uegawa))
} # hiritsu_shinraikukanの定義終了
hiritsu_shinraikukan(200,9,0.95) # 例5.13
hiritsu_shinraikukan(1000,482,0.95) # 問5.14
############# R の関数として以下のみでOK
prop.test(482,1000,p=0.5,alt="t",conf.level=0.95,correct=F)
# 482/1000, p=0.5 は比率が0.5の帰無仮説(信頼区間を求める際は無関係)
# correct=F は連続補正しないことを意味する
出力の一部
data: 482 out of 1000, null probability 0.5
X-squared = 1.296, df = 1, p-value = 0.2549
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
0.4511585 0.5129793 # 信頼区間
sample estimates:
p
0.482
# prop.test の出力は (0.4511585+0.5129793)/2=0.4820689なので左右対称の信頼区間ではなく、何らかの修正で計算されている(要チェック)
|
|
# 中田、内藤、p.138 例5.17, 問5.22 章末問題5.10 (正規母集団, 母分散既知 検定)
# kichi_test (標本の大きさ, 母分散, 標本平均, 有意水準, 帰無仮説, 対立仮説)
kichi_test <- function(n, sigma2, xbar, yuuisuijyun, mu0, alt) {
shinraikeisu<- 1-yuuisuijyun
p <- shinraikeisu # 片側検定
z <- qnorm(p) # z<-qnorm(p) <=> P(N(0,1) <= z)=p, z_alpha/2
p2 <- shinraikeisu+(1-shinraikeisu)/2 # 両側検定
z2 <- qnorm(p2)
sgm <- sqrt(sigma2) # 母標準偏差
ztest <- sqrt(n)*(xbar-mu0)/sgm # 検定統計量の実現値
cat(sprintf("### 母平均の検定: "))
switch(alt, # 棄却1, 採択0
"左側" = {
pchi <- pnorm(ztest); # p値の計算
if(ztest < -z) kikyaku<-1 else kikyaku<-0;
cat(sprintf("H0: mu=%.3f vs mu<%.3f \n", mu0, mu0));
cat(sprintf("左側検定：棄却域=(-inf, %.3f]\n", -z))
},
"右側" = {
pchi <- 1-pnorm(ztest); # p値の計算
if(ztest > z) kikyaku<-1 else kikyaku<-0;
cat(sprintf("H0: mu=%.3f vs mu>%.3f \n", mu0, mu0));
cat(sprintf("左側検定：棄却域=[%.3f, inf)\n", z))
},
"両側" = {
pchi <- 2*(1-pnorm(abs(ztest))); # p値の計算
if(abs(ztest) > z) kikyaku<-1 else kikyaku<-0;
cat(sprintf("H0: mu=%.3f vs mu!=%.3f \n", mu0, mu0));
cat(sprintf("両側検定：棄却域=
(-inf, %.3f]or [%.3f, inf)\n", -z2, z2))
}
)
cat(sprintf("検定統計量の実現値: %.3f, p値: %.3f\n", ztest, pchi))
if(kikyaku==1){
cat(sprintf("有意水準%.3fで棄却されます\n",yuuisuijyun))
}else{
cat(sprintf("有意水準%.3fでは棄却されません\n",yuuisuijyun))
}
} # kichi_testの定義終了
#kichi_test (標本の大きさ, 母分散, 標本平均, 有意水準, 帰無仮説, 対立仮説)
kichi_test(n=16, sigma2=5.31^2, xbar=997.79,
yuuisuijyun=0.05, mu0=1000, alt="左側") #p.138 例5.17
# kichi_test(n=10, sigma2=1.5^2, xbar=49,
# yuuisuijyun=0.01, mu0=50, alt="左側") #p.139 問5.22
kichi_test(n=50, sigma2=36, xbar=160,
yuuisuijyun=0.05, mu0=158, alt="両側") #p.145 章末問題5.10
t.test(toi5.12, alternative = "less", mu = 25, conf.level = 0.95)出力の一部
t = -2.5025, df = 20, p-value = 0.01056 # t0 の実現値, 自由度, p値
alternative hypothesis: true mean is less than 25 # 対立仮説
95 percent confidence interval:
-Inf 24.01577
sample estimates:
mean of x
21.83333 # 標本平均
# 中田、内藤、 p.143 例5.20, 問5.25 (正規母集団, 比率の検定；大標本)
# hiritsu_test (標本の大きさ, 成功回数, 有意水準, 帰無仮説, 対立仮説)
hiritsu_test <- function(n, ncount, yuuisuijyun, p0, alt) {
shinraikeisu<- 1-yuuisuijyun
p <- shinraikeisu # 片側検定
z <- qnorm(p) # z<-qnorm(p) <=> P(N(0,1) <= z)=p, z_alpha/2
p2 <- shinraikeisu+(1-shinraikeisu)/2 # 両側検定
z2 <- qnorm(p2)
xbar <- ncount/n # 比率
ztest <- sqrt(n)*(xbar-p0)/sqrt(p0*(1-p0)) # 検定統計量の実現値
cat(sprintf("### 母比率の検定: "))
switch(alt, # 棄却1, 採択0
"左側" = {
pchi <- pnorm(ztest); # p値の計算
if(ztest < -z) kikyaku<-1 else kikyaku<-0;
cat(sprintf("H0: mu=%.3f vs mu<%.3f \n", p0, p0));
cat(sprintf("左側検定：棄却域=(-inf, %.3f]\n", -z))
},
"右側" = {
pchi <- 1-pnorm(ztest); # p値の計算
if(ztest > z) kikyaku<-1 else kikyaku<-0;
cat(sprintf("H0: mu=%.3f vs mu>%.3f \n", p0, p0));
cat(sprintf("左側検定：棄却域=[%.3f, inf)\n", z))
},
"両側" = {
pchi <- 2*(1-pnorm(abs(ztest))); # p値の計算
if(abs(ztest) > z) kikyaku<-1 else kikyaku<-0;
cat(sprintf("H0: mu=%.3f vs mu!=%.3f \n", p0, p0));
cat(sprintf("両側検定：棄却域=
(-inf, %.3f]or [%.3f, inf)\n", -z2, z2))
}
)
cat(sprintf("検定統計量の実現値: %.3f, p値: %.3f\n", ztest, pchi))
if(kikyaku==1){
cat(sprintf("有意水準%.3fで棄却されます\n",yuuisuijyun))
}else{
cat(sprintf("有意水準%.3fでは棄却されません\n",yuuisuijyun))
}
} # hiritsu_testの定義終了
# hiritsu_test (標本の大きさ, 成功回数, 有意水準, 帰無仮説, 対立仮説)
hiritsu_test(n=500, ncount=30,
yuuisuijyun=0.01, p0=0.04, alt="両側") #p.143 例5.20
hiritsu_test(n=600, ncount=219,
yuuisuijyun=0.01, p0=0.4, alt="両側") #p.144 問5.25
|
|
############################################################
# 回帰直線の係数の出力
# 中田、内藤、p.150 例6.1, p.157 問6.7, p.158 章末問題6.11
############################################################
rm(list = ls(all.names = TRUE)) # オブジェクトを全て削除する
options(encoding = "shift-jis") # 漢字コードをshift-jisにする
#options(encoding = "utf-8") # 漢字コードをutf-8にする
## kaiki_chokusen_keisu(x,y): x, y の回帰直線だけを出力する
rm(list = ls(all.names = TRUE))
kaiki_chokusen_keisu <- function(x, y) {
n <- length(x)
xsum<- sum(x); ysum<- sum(y)
x2sum<- sum(x*x); y2sum<- sum(y*y); xysum<- sum(x*y)
bhat <- (n*xysum - xsum*ysum)/(n*x2sum-xsum^2)
ahat <- ysum/n-bhat*xsum/n
lm.xy <- lm(y ~ x) # R の関数としての回帰直線
cat(sprintf("例6.1の方法で求めた回帰直線 \n
y = %.3f + (%.4f) x\n", ahat, bhat))
#cat(sprintf("\nRの関数による回帰直線の詳細 \n"))
#summary(lm.xy)
} # kaiki_chokusen_keisuの定義終了
# #### p.150 例6.1
# jinkowariai<-read.csv("c:/R/hyo1_1.csv")
# #jinkowariai
# x <- jinkowariai$老年人口割合
# y <- jinkowariai$年少人口割合
# kaiki_chokusen_keisu(x,y)
# #### p.157 問6.7
# komeshuryo<-read.csv("c:/R/toi6_7.csv")
# #komeshuryo
# x <- komeshuryo$年
# y <- komeshuryo$収量
# kaiki_chokusen_keisu(x,y)
#### p.158 章末問題6.11
nenreizakohi<-read.csv("c:/R/shomatsu6_11.csv") # データの読み込み
#nenreizakohi
x <- nenreizakohi$年齢
y <- nenreizakohi$座高比
kaiki_chokusen_keisu(x,y)
############################################################
# 回帰直線における推定(定理6.4)
# 中田、内藤、p.150 例6.1, p.157 問6.7, p.158 章末問題6.11
############################################################
rm(list = ls(all.names = TRUE)) # オブジェクトを全て削除する
options(encoding = "shift-jis") # 漢字コードをshift-jisにする
#options(encoding = "utf-8") # 漢字コードをutf-8にする
kaiki_chokusen_shinraikukan <- function(x, y, vx, shinraikeisu) {
n <- length(x) # 標本の大きさ
xsum<- sum(x); ysum<- sum(y)
x2sum<- sum(x*x); y2sum<- sum(y*y); xysum<- sum(x*y)
xbar <- xsum/n
bhat <- (n*xysum - xsum*ysum)/(n*x2sum-xsum^2)
ahat <- ysum/n-bhat*xsum/n
sx2 <- x2sum/n-(xsum/n)^2 # xの標本分散
sy2 <- y2sum/n-(ysum/n)^2 # yの標本分散
r <- (n*xysum-xsum*ysum)/sqrt((n*x2sum-xsum^2)*(n*y2sum-ysum^2))
S <- n*sy2*(1-r^2) # p.150 (6.14) Sの定義どおり
sigmahat2 <- S/(n-2) # p.156 (6.22)
hnx <- sqrt((1+(vx-xbar)^2/sx2)* sigmahat2/n) # 定理6.4
#cat(sprintf("n = %d, S(a,b) = %.2f, sigmahat2= %.3f, t0test = %.3f\n",
# n, S, sigmahat2, t0test)) # 途中出力
abx <- ahat+ bhat*vx
p2 <- shinraikeisu+(1-shinraikeisu)/2 #
haba <- hnx * qt(p2, n-2)
shitagawa <- abx - haba # 下側信頼限界
uegawa <- abx + haba # 上側信頼限界
cat(sprintf("### 回帰直線における推定: "))
cat(sprintf("回帰直線 \n
y = %.3f + (%.4f) x\n", ahat, bhat))
cat(sprintf("x=%.1f における%.1f%%信頼区間\t: [%.3f, %.3f]\n",
vx, shinraikeisu*100, shitagawa, uegawa))
} # kaiki_chokusen_testの定義終了
# #### # 中田、内藤、p.156 例6.3
# jinkowariai<-read.csv("c:/R/hyo1_1.csv")
# x <- jinkowariai$老年人口割合
# y <- jinkowariai$年少人口割合
# yuuisuijyun <- 0.05
# kaiki_chokusen_shinraikukan(x,y, 40, 0.95)
# #### # 中田、内藤、p.157 問6.7
# komeshuryo<-read.csv("c:/R/toi6_7.csv")
# #komeshuryo
# x <- komeshuryo$年
# y <- komeshuryo$収量
# kaiki_chokusen_shinraikukan(x,y, 95, 0.95) # 教科書のt分布表とは少しずれる！
#### # 中田、内藤、p.158 章末問題6.11
nenreizakohi<-read.csv("c:/R/shomatsu6_11.csv")
#nenreizakohi
x <- nenreizakohi$年齢
y <- nenreizakohi$座高比
kaiki_chokusen_shinraikukan(x,y, 15, 0.99)
# 教科書の問題にはないが１５歳の信頼区間を出力
############################################################
# 回帰直線における検定(定理6.5)
# 中田、内藤、p.150 例6.1, p.157 問6.7, p.158 章末問題6.11
############################################################
rm(list = ls(all.names = TRUE)) # オブジェクトを全て削除する
options(encoding = "shift-jis") # 漢字コードをshift-jisにする
#options(encoding = "utf-8") # 漢字コードをutf-8にする
kaiki_chokusen_test <- function(x, y, yuuisuijyun) {
n <- length(x) # 標本の大きさ
shinraikeisu<- 1-yuuisuijyun
p2 <- shinraikeisu+(1-shinraikeisu)/2 # 両側検定
z <- qt(p2, n-2) # z<-qt(p,df) <=> P(t(df) <= z)=p
xsum<- sum(x); ysum<- sum(y)
x2sum<- sum(x*x); y2sum<- sum(y*y); xysum<- sum(x*y)
bhat <- (n*xysum - xsum*ysum)/(n*x2sum-xsum^2)
sx2 <- x2sum/n-(xsum/n)^2 # xの標本分散
sy2 <- y2sum/n-(ysum/n)^2 # yの標本分散
r <- (n*xysum-xsum*ysum)/sqrt((n*x2sum-xsum^2)*(n*y2sum-ysum^2))
S <- n*sy2*(1-r^2) # p.150 (6.14) Sの定義どおり
sigmahat2 <- S/(n-2) # p.156 (6.22)
t0test <- bhat/sqrt(sigmahat2/(n*sx2))# 検定統計量の実現値
#cat(sprintf("n = %d, S(a,b) = %.3f, sigmahat2= %.3f, t0test = %.3f\n",
# n, S, sigmahat2, t0test)) # 途中出力
pchi <- 2*(1-pt(abs(t0test),n-2)) # p値の計算
cat(sprintf("### 回帰直線における検定: "))
if(abs(t0test) > z) kikyaku<-1 else kikyaku<-0;
cat(sprintf("H0: mu=0 vs mu!=0 \n"));
cat(sprintf("両側検定：棄却域=(-inf, %.3f] or [%.3f, inf)\n", -z, z))
cat(sprintf("検定統計量の実現値: %.3f, p値: %.3f\n", t0test, pchi))
if(kikyaku==1){
cat(sprintf("有意水準%.3fで棄却されます\n",yuuisuijyun))
}else{
cat(sprintf("有意水準%.3fでは棄却されません\n",yuuisuijyun))
}
} # kaiki_chokusen_testの定義終了
# #### p.156 例6.3
# jinkowariai<-read.csv("c:/R/hyo1_1.csv")
# x <- jinkowariai$老年人口割合
# y <- jinkowariai$年少人口割合
# yuuisuijyun <- 0.05
# kaiki_chokusen_test(x, y, yuuisuijyun)
# #### p.157 問6.7
# komeshuryo<-read.csv("c:/R/toi6_7.csv")
# #komeshuryo
# x <- komeshuryo$年
# y <- komeshuryo$収量
# yuuisuijyun <- 0.05
# kaiki_chokusen_test(x, y, yuuisuijyun)
#### p.158 章末問題6.11
nenreizakohi<-read.csv("c:/R/shomatsu6_11.csv")
#nenreizakohi
x <- nenreizakohi$年齢
y <- nenreizakohi$座高比
kaiki_chokusen_test(x, y, 0.01)
############################################################
# 回帰直線のグラフ表示
# 中田、内藤、 p.150 例6.1, p.157 問6.7, p.158 章末問題6.11
############################################################
rm(list = ls(all.names = TRUE)) # オブジェクトを全て削除する
options(encoding = "shift-jis") # 漢字コードをshift-jisにする
#options(encoding = "utf-8") # 漢字コードをutf-8にする
if( dev.cur() > 1 ) dev.off()
kaiki_chokusen_hyoji <- function(x, y) {
m <- lm(y ~ x) # 回帰直線 p. 156
# newx <- data.frame(x = seq(20, 50, 0.5)) # x軸の最小、最大、ステップ
newx <- data.frame(x = seq(min(x), max(x), 0.5)) # x軸の最小、最大、ステップ
conf.interval <- predict(m, newdata = newx,
interval = 'confidence', level = 0.95) #信頼区間
plot(x, y,
#data = trees,
col = "red",
pch = 16, # プロットマーカーを塗りつぶす
# asp=1, # アスペクト比
# main="人口割合", xlab="老年人口割合ｘ100", ylab="年少人口割合ｘ100",
)
#abline(m, lwd=3)
lines(newx$x, conf.interval[, 1], col = 'green', lty=1, lwd=3) # 実線, 太さ3
lines(newx$x, conf.interval[, 2], col = 'blue', lty=2, lwd=3) # 点線, 太さ3
lines(newx$x, conf.interval[, 3], col = 'blue', lty=2, lwd=3) # 点線, 太さ3
#summary(m);coef(m)
} # kaiki_chokusen_hyojiの定義終了
#jinkowariai<-read.csv("c:/R/hyo1_1.csv")
#jinkowariai; summary(jinkowariai)
#####
### p.9 表1.1 人口割合の散布図
jinkowariai<-read.csv("c:/R/hyo1_1.csv")
x <- jinkowariai$老年人口割合
y <- jinkowariai$年少人口割合
kaiki_chokusen_hyoji(x,y)
#### p.157 問6.7
# komeshuryo<-read.csv("c:/R/toi6_7.csv")
# #komeshuryo
# x <- komeshuryo$年
# y <- komeshuryo$収量
# kaiki_chokusen_hyoji(x,y)
# # #### p.158 章末問題6.11
# nenreizakohi<-read.csv("c:/R/shomatsu6_11.csv")
# #nenreizakohi
# x <- nenreizakohi$年齢
# y <- nenreizakohi$座高比
# kaiki_chokusen_hyoji(x,y)
|
|
# 破産の問題 (中田、内藤p.70, 例3.19)
# 破産確率、破産するまでの回数の平均、そのヒストグラム
trial <- 1000 # 試行回数
a <- 8 # A君の所持金
N <- 20 # A君とB君の合計の初期金
p <- 0.5 # A君の勝つ確率
hasan <- function(a, N, p){
shojikin <- a # 初期のA君の所持金
hkaisu <- 0 # 破産するまでの回数
while (shojikin > 0 & shojikin < N) {
pm1 <- sample(c(-1,1),1,prob=c(1-p,p)) # +1, -1
shojikin <- shojikin + pm1
hkaisu <- hkaisu + 1
}
if (shojikin == 0)
flag <- 1 # A君破産なら1
else
flag <- 0 # B君破産なら0
hasan_jikan<-hkaisu # 破産するまでの回数
return(c(flag, hasan_jikan)) # 破産、回数を返り値にする
} # hasan 終了
########### trial回繰り返す!
hasan_data <- replicate(trial, hasan(a, N, p)) # hasan をtrial回繰り返す
########### 結果の出力のための変数準備
num_Awin <- sum(hasan_data[1,]) # A君の勝利数
mean_Awin <- mean(hasan_data[1,]) # A君の勝つ確率の標本平均
ideal_mean_Awin <- 1-a/N # 理論値(p=1/2)
mean_kaisu <- mean(hasan_data[2,]) # 回数の標本平均
ideal_mean_kaisu <- a*(N-a) # 理論値(p=1/2)
out<-sprintf(
"############### 破産の問題 #######################
# 試行\t\t\t\t= %d
# A君の初期の所持金\t\t= %d
# A君とB君の合計の所持金\t= %d
# A君の勝つ確率\t\t\t= %.2f
# A君の破産数\t\t\t= %d
# A君の破産率\t\t= %5.2f (理論値: %.2f)
# 破産までの平均回数\t= %5.2f (理論値: %5.2f)
##################################################\n",
trial, # 試行
a, # A君の初期の所持金
N, # A君とB君の合計の初期金
p, # A君の勝つ確率
num_Awin, #A君の破産数
mean_Awin, #A君の破産率
ideal_mean_Awin, #A君の破産率(理論値)
mean_kaisu, #破産までの平均回数
ideal_mean_kaisu #破産までの平均回数(理論値)
)
# ヒストグラム
hist(hasan_data[2,], # 破産までの回数のヒストグラム
main=paste("破産までの回数のヒストグラム (試行", trial, "回)\n",
"A君の所持金", a,"円, 所持金の合計", N,"円"
),
xlab="破産までの回数", ylab="度数"
)
# 理論値を書き込む
abline(v = ideal_mean_kaisu, col = "red", lwd = 2, lty = "dotted")
#### 結果の出力
cat(out)
|
|
N<-10 #コイン投げ coin <- c(0,1); coinN <- sample(coin, size = N, replace = TRUE); coinN N<-10 #サイコロ投げ dice <- c(1,2,3,4,5,6); diceN <- sample(dice, size = N, replace = TRUE); diceN
#サイコロ投げで1の割合が1/6に近づく様子(大数の法則)
### サイコロ投げで1の割合が1/6に近づく様子(大数の法則)
#set.seed(40) # 必要ならば乱数の種を固定する
###
n <- 1000 # サイコロを投げる回数
p <- 1/6 # サイコロで1が出る確率
dice1 <- sample(0:1, n, replace=TRUE, prob=c(1-p,p)) # replace=TRUE で復元抽出
succ1 <- cumsum(dice1) # 1が出た回数の合計
prop <- succ1/(1:n) # 1からn回のうちの1が出た回数の割合
plot(1:n,prop,type="l",xlab="サイコロ投げの試行回数",
ylab="1の割合",
main = "1000回のうちの1の割合")
abline(h = p, col="red", lwd=2, lty="dotted") # 理論値の横線、lwd:線の太さ
legend(800,0.07,legend = c(expression(bar(X)[n]), expression(mu)), #凡例
col = c("black", "red"), lty = c("solid", "dotted"), lwd = 2)
### サイコロ投げで1が出る標本平均が1/6に近づく様子の標本サイズによる比較
#rm(list = ls(all.names = TRUE)) # オブジェクトを全て削除する
p <- 1/6 # 1が出る確率
mu <- p # mu = E(Be(p))
sample_mean <- function(n){ # サイコロの目のn回の標本平均
mean(sample(0:1, n, replace=TRUE, prob=c(1-p,p))) # Be(p)の標本平均
}
trial <- 1000 # 試行回数
path_length <- c(10, 100) # パスの長さ
op <- par(fg = "orange", col = "black", mfrow = c(1, 2)) # 色付け、画面を分割
for(n in path_length){ # n:パスの長さ
xbar <- replicate(trial, sample_mean(n)) # 標本平均をtrial 回
hist(xbar, breaks = 20, xlim = c(0,1), # x の範囲を(0,1)にする
col = "blue", border = "blue",
main = paste0("標本サイズ",n), xlab = "標本平均", ylab = "頻度")
abline(h = 0, lwd = 2, lty = "solid")
abline(v = mu, col = "red", lwd = 2, lty = "dotted")
}
par(op) # 初期化
|
|
################
# X1,...,Xn sim Exp(1), (iid), Xbar=(X_1+...X_n)/n
# Xbar のn=1,2,..のヒストグラムのアニメーション（中心極限定理）
# パラメータ：
# n = 1,2,..,rpt回繰り返す
# 乱数を trl 個抽出
# アニメの時間間隔はkankaku秒
# 現在のフォルダにmovienm として保存
################
#rm(list = ls(all.names = TRUE)) # オブジェクトを全て削除する
#パッケージの読み込み(ない場合は[パッケージ]からインストールする)
library("animation")
### パラメータの設定
rpt <- 10 #アニメーションの枚数
trl <- 10000 # 繰り返しの回数
kankaku <- 1.0 # 切替間隔(秒)
movienm <-"unifmean.gif" #画像のファイル名
unif_mean <- function(size, trial){ # 関数の定義
#trial <- 1000 # 繰り返しの回数
#size <- i # 標本平均の大きさ（ここを変化させる）
width <- 1/(size*10) # 階級の幅
#x <- matrix(runif(trial*size), nrow=trial, ncol=size) # 一様乱数をtrial×size 個用意
###### 一様分布のを他の分布に変える
x <- matrix(rexp(trial*size, rate =1), nrow=trial, ncol=size) # 標準指数分布
##################################################################
xm <- apply(x, 1, mean) # 行ごとの平均
xma <- floor(min(xm)) # 最小値
xmb <- ceiling(max(xm)) # 最大値
xmc <- max(abs(xma),abs(xmb)) # 最大値と最小値の絶対値の大きい方
m <- mean(xm); # 平均
s <- sd(xm); # 標準偏差
hist(xm, breaks=seq(xma,xmb,by=width), # 最小から最大までwidthきざみ
xlim=c(0,xmc), # 横軸の範囲
ylim=c(0,1.2/(sqrt(2*pi)*s)), # 縦軸の範囲
freq=FALSE, xlab="", ylab="頻度", # 軸のラベル
#main = paste("[0,1]一様rvの", size, "個の標本平均","繰り返し回数:", trial) # メインのタイトル
main = paste("標準指数rvの", size, "個の標本平均","繰り返し回数:", trial) # メインのタイトル
)
curve(dnorm(x, m, s), add=T, col=2) # 標準正規分布の赤での重ね書き
}# unif_mean の定義の終了
### gif ファイルにアニメーションを保存
saveGIF({
for (i in 1:rpt){
#処理内容を記述する
unif_mean(i, trl)
}},
interval = kankaku,
movie.name = movienm
)
|
|
x<-rnorm(1000000) hist(x)
x<-rcauchy(1000000) hist(x)軸の範囲などで以下のように調整が必要である。
############################################
##### コーシー分布 のヒストグラム
# 標本平均がコーシー分布に従うため size を換えてもヒストグラムは変化しない!
############################################
trial <- 10000 # 繰り返しの回数
size <- 1 # 標本の大きさ(標本平均もコーシー分布)
width <- 0.1 # 階級の幅
half_range <- 10 # ヒストグラムの範囲の半分
x <- matrix(data=rcauchy(trial*size), ncol=size) # コーシー乱数をtrial×size 個用意
###### コーシー分布の性質(乱数のとり方を変える!)
#x <- matrix(data=1/rcauchy(trial*size), ncol=size) # 逆数もコーシー分布
#x <- matrix(data=rnorm(trial*size)/rnorm(trial*size), ncol=size) # N(0,1)/N(0,1)はコーシー分布
#x <- matrix(data=tan(runif(trial*size,-pi/2,pi/2)), ncol=size) #tan(一様(-π/2,π/2)はコーシー分布
##################################################################
xm <- apply(x, 1, mean) # 行ごとの平均
xma <- floor(min(xm)) # 最小値
xmb <- ceiling(max(xm)) # 最大値
hist(xm, breaks=seq(xma,xmb,by=width), # 最小から最大までwidthきざみ
xlim=c(-half_range,half_range), # 横軸の範囲
# ylim=c(0,1), # 縦軸の範囲
freq=FALSE, xlab="", ylab="頻度", # 軸のラベル
main="コーシー分布のヒストグラム") # メインのタイトル
curve(dcauchy(x, 0, 1), add=T, col=2) # 標準コーシー分布の赤での重ね書き
上記のコメントの部分の変更
f(x)= 8*(log(x))*((log(x))^2+pi^2)/(3*pi^4*(x^2-1)) (x>0, x!=1)
4/(3*pi^2) (x=1)
0 (x<=0)
であることのヒストグラムでのチェック
N <- 10000 ave <- numeric(N) # 0 をN個の配列で用意 x <- rcauchy(N) # コーシー分布に従う乱数を N 個 配列 x に格納 for (n in 1:N ) ave[n] <- mean(x[1:n]) # 1からnまでの平均をとる plot(1:N, ave) # データを表示 abline(h=0) # h = 0 の横線を入れる
##################
# Kolmogorov-Feller の定理(Gut-Kolmogorov-Feller の定理)の応用例
# (Gut, Probability, 2013, Springer, Th.4.4.2, p.281, b_n=n*log n)
# X_k:(iid) 2/(π(1+x^2)) (x>0) 右側標準コーシー分布
# ならば
# lim (X1+...+Xn)/(n*log n) = 2/π (確率収束)
##################
N <- 10000
exactlaw <- numeric(N) # 0 をN個の配列で用意
x <- abs(rcauchy(N)) # コーシー分布の絶対値に従う乱数を N 個 配列 x に格納
for (n in 1:N ) exactlaw[n] <- sum(x[1:n])/(n*log(n)) # 足してn*log(n)で割る!
plot(1:N, exactlaw,
xlab="回数n", ylab="X1+...+Xn/n*log(n)",
ylim=c(0,1), # 縦軸の範囲
) # データを表示
abline(h=2/pi) # 収束先の h = 2/pi の横線を入れる
### たいして良いシミュレーションでもない?
|
|
# コーシー分布の標本データの最大値の分布 gev.fit # https://search.r-project.org/CRAN/refmans/ismev/html/gev.fit.html rm(list = ls(all.names = TRUE)) # 変数クリア if( dev.cur() > 1 ) dev.off() # グラフクリア library(ismev) # 極値統計用パッケージ N<-10000 # 標本サイズ num_year <- 50 # まとめた年数 size = N/num_year # 各年のデータのサイズ x1<- rcauchy(N) # 標準コーシー分布からのN個の実現値 x<- matrix(x1, nrow=num_year, ncol=size) # num_year x size 行列にする xmax <- apply(x, 1, max) # 各行の最大値 x.test.ismev <- gev.fit(xmax) # GEV分布 の最尤フィッティング gev.diag(x.test.ismev) # グラフ表示(gev.fitの適用後)
|
|
# ビュホンの針
rm(list = ls(all.names = TRUE)) # 変数クリア
if( dev.cur() > 1 ) dev.off() # グラフクリア
n <- 300 # 針の本数
NL <- 1 # 針の長さ
W <- 1 # 平行線の縦の間隔
GS <- 5 # 画面サイズ
# 横線の y 座標：0, W, 2W, ..., [GS/W]
plot(0,
xlim=c(0,GS), ylim=c(0,GS),
type="n", #プロットしない
axes=F, # 軸を描かない
ann=F) # 軸ラベルなし
# 平行線をひく
for(i in 0:floor(GS/W)) abline(h=W*i,
lwd=2, # 太さ２
col="blue")
px <- runif(n, 0, GS) # x 座標n個
py <- runif(n, 0, GS) # y 座標n個
angle <- runif(n, 0, pi) # 角度 n個 (本来はpiを使うのは反則)
HL <- NL/2 # 針の長さの半分
pp<- cbind(px - HL*cos(angle), px + HL*cos(angle),
py- HL*sin(angle), py + HL*sin(angle)) # 端点の座標を配列にする
# y座標の始点、終点を求めて平行線と交わっているかどうかを調べる
ytab <- data.frame(pp[,3],pp[,4])
ymax <- apply(ytab, 1, max) # y座標の大きい方
ymin <- apply(ytab, 1, min) # y座標の小さい方
# y座標の始点、終点の間に平行線が存在すればcheck!
check <- (ymin < W*floor(ymax)/W)
#針を描く関数を定義
needle_line <- function(k){
lines(pp[k,1:2], pp[k,3:4],
col=rgb(check[k],0,0.5))# rgb の引数は0..1(色はテキトー)
}
#針を描く
void <-sapply(1:n, needle_line)
# piの近似
cat("pi =", 2*NL/(W*mean(check)), "\n")
|
|