スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

Parsley!!

この記事を閲覧するにはパスワードが必要です
パスワード入力

Creatures

この記事を閲覧するにはパスワードが必要です
パスワード入力

標本分散と不偏分散

標本数$N$の標本分散が不偏分散の$\frac{N-1}{N}$倍となることの、elephantな「力技」による確認を行う。


標本分散は次の関数として定義される:
variance <- function(x){
mean(x^2 - mean(x)^2)
}

Rでvar()が返す値は不偏分散であり、var()=$\frac{N}{N-1}$variance()であることに注意する。


正規分布の場合:
pdf("gauss17Mar2013.pdf", height=5.5, width=7, pointsize=14)
x <- rnorm(100000,0,1)
hist(x,breaks=100,probability=T,main="Norm",col="#CDBEE9")
dev.off()
windows()

gauss17Mar2013_01-mini.png

test <- matrix(nrow=100000, ncol=200)
test2 <- matrix(nrow=100000, ncol=200)
plot <- matrix(nrow=200,ncol=1)
plot2 <- matrix(nrow=200,ncol=1)

for(j in c(5:200)){
for(i in c(1:100000)){
x <- rnorm(j,0,2)
test[i,j] <- variance(x)
test2[i,j] <- var(x)
}
plot[j,1] <- mean(test[,j])
plot2[j,1] <- mean(test2[,j])
}

pdf("normplot17Mar2013.pdf", height=5.5, width=7, pointsize=14)
plot(plot2[,1],xlim=c(4,200),ylim=c(3.1,4.1),xlab="N",ylab="variance",
main="Norm",col="#E82202",pch=20,ps=0.01)
par(new=T)
plot(plot[,1],xlim=c(4,200),ylim=c(3.1,4.1),xlab="N",ylab="variance",
col="#6013FB",pch=20,ps=0.01)
dev.off()
windows()

normplot17Mar2013_01-mini.png



Weibull分布の場合(※1)
pdf("weibull17Mar2013.pdf", height=5.5, width=7, pointsize=14)
x <- rweibull(100000,1,7)
hist(x,breaks=100,xlim=c(0,60),probability=T,main="Weibull",col="#CDBEE9")
dev.off()
windows()

weibull17Mar2013_01-mini.png

※1 なぜここでWeibull分布? と思われた方も多いだろう。特に理由は無いのだが、私はWeibull分布とtangledな関係にあるらしく、例に示すなら正規分布の次にWeibull分布、という約束に、今日からした。

for(j in c(5:200)){
for(i in c(1:100000)){
x <- rnorm(j,0,2)
test[i,j] <- variance(x)
test2[i,j] <- var(x)
}
plot[j,1] <- mean(test[,j])
plot2[j,1] <- mean(test2[,j])
}

pdf("normplot17Mar2013.pdf", height=5.5, width=7, pointsize=14)
plot(plot2[,1],xlim=c(4,200),ylim=c(3.1,4.1),xlab="N",ylab="variance",
main="Norm",col="#E82202",pch=20,ps=0.01)
par(new=T)
plot(plot[,1],xlim=c(4,200),ylim=c(3.1,4.1),xlab="N",ylab="variance",
col="#6013FB",pch=20,ps=0.01)
dev.off()
windows()

weibullplot17Mar2013_02-mini.png

多項式近似

sin(x)からランダムにズレた点群への多項式fitting(6次多項式)

num50_01-mini.png

過学習の例:

num7_01-mini.png



num <- 50
dim <- 7 # 6次多項式で近似する

x <- runif(num, min = 0, max = 5)
y <- sin(x)+runif(num, min = -0.2, max = 0.2)
plot(x,y,xlim=c(0,5),ylim=c(-1,1))

A <- matrix(nrow=dim, ncol=dim)
B <- matrix(nrow=dim, ncol=1)

for(i in c(1:dim)){
for(j in c(1:dim)){
A[i,j] <- sum(x**(i+j-2))
}
}

for(i in c(1:dim)){
B[i,1] <- sum(y*(x**(i-1)))
}

iA <- solve(A)

ans <- iA %*% B

polyapp <- function(x,ans){
y <- 0
for(i in c(1:length(ans))){
y <- y+ans[i]*x^(i-1)
}
return(y)
}

col <- c("#1B9E77", "#D95F02", "#7570B3")
main <- ""
xlab <- "x"
ylab <- "y"

pdf("test2.pdf", height=5.5, width=7, pointsize=14)
curve(sin(x),xlim=c(0,5),ylim=c(-1.4,1.4),xlab=xlab,ylab=ylab,
col="red",lty=2,lwd=1.5)
par(new=T)
curve(polyapp(x,ans),xlim=c(0,5),ylim=c(-1.4,1.4),xlab=xlab,ylab=ylab,
col=col[1],lwd=1.5)
par(new=T)
plot(x,y,xlim=c(0,5),ylim=c(-1.4,1.4),col=col[3])
dev.off()

Rデー

多重の条件に当てはまる行の抽出
今のところ下記以上にスマートに実行する方法は知らない。

> swiss
Fertility Agriculture Examination Education Catholic Infant.Mortality
Courtelary 80.2 17.0 15 12 9.96 22.2
Delemont 83.1 45.1 6 9 84.84 22.2
Franches-Mnt 92.5 39.7 5 5 93.40 20.2
  ・
  ・
  ・


> x <-swiss[swiss["Fertility"]<=72.0,]
> x
Fertility Agriculture Examination Education Catholic Infant.Mortality
Aigle 64.1 62.0 21 12 8.52 16.5
Aubonne 66.9 67.5 14 7 2.27 19.1
Avenches 68.9 60.7 19 12 4.43 22.7
  ・
  ・
  ・


> y <-x[x["Agriculture"]<=60.0,]
> y
Fertility Agriculture Examination Education Catholic Infant.Mortality
Grandson 71.7 34.0 17 8 3.30 20.0
Lausanne 55.7 19.4 26 28 12.11 20.2
La Vallee 54.3 15.2 31 20 2.15 10.8
  ・
  ・
  ・



逆行列の計算と行列の積

> A <- matrix(nrow=3, ncol=3)
> A

[,1] [,2] [,3]
[1,] NA NA NA
[2,] NA NA NA
[3,] NA NA NA

> A[,1] <- round(runif(3, min = 0, max = 9))
> A[,2] <- round(runif(3, min = 0, max = 9))
> A[,3] <- round(runif(3, min = 0, max = 9))
> A

[,1] [,2] [,3]
[1,] 6 4 7
[2,] 2 1 7
[3,] 7 1 2

> B <- solve(A)
> B
[,1] [,2] [,3]
[1,] -0.04347826 -0.008695652 0.1826087
[2,] 0.39130435 -0.321739130 -0.2434783
[3,] -0.04347826 0.191304348 -0.0173913

> round(A %*% B) # 行列の積

[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1

平均待ち時間の公式の導出

この記事を閲覧するにはパスワードが必要です
パスワード入力
Profile

zoa

Author:zoa
Metropolis (1927)

Calendar
02 | 2013/03 | 04
- - - - - 1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31 - - - - - -
Labels
Admin
Previous Posts
Recent Comments
Recent Trackback
Archive
このページのトップへ
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。