統計学2022 Lecture 7: 正規分布におけるパラメータの推定

Rによる演習

正規分布の確率密度関数の概形

#正規分布の確率密度関数を描く関数を作成
NormalPDF1 <- function(mu, sigma, From, To){ 
   curve(dnorm(x,mu,sigma), xlim=c(From,To), ylim=c(0,0.2),  ylab="Density", 
        main=paste("N(",mu,",",sigma^2,")")) 
   lines(c(mu,mu),c(0,dnorm(mu,mu,sigma)),lty=2)
}

#上記の関数を実行
par(mar=c(2,2,2,2), mfrow=c(2,2)) #mfrow=c(2,2)で画面を2×2に分割し，行から順に埋めていく
NormalPDF1(15,2,0,35)
NormalPDF1(20,2,0,35)
NormalPDF1(15,4,0,35)
NormalPDF1(20,4,0,35)

ライントランセクト法による鯨類資源量の推定値の確率分布

発見関数の定義

#detection function
g.HN <- function(PD, sigma, maxPD=999){ exp(-0.5*(PD/sigma)^2 )*(PD <= maxPD) }

#probability density function for perpendicular distance
f.HN <- function(PD, sigma, mu, maxPD=999){
 g <- exp(-0.5*(PD/sigma)^2 )*(PD <= maxPD)
 f <- g/mu; f
}

maxPD<-3
curve(g.HN(PD=x, sigma=0.5), lty=1, xlim=c(0, maxPD), xlab="横距離", ylab="発見確率")
curve(g.HN(PD=x, sigma=1) , lty=2, add=T)
curve(g.HN(PD=x, sigma=2) , lty=3, add=T)
legend(2, 1, legend=c("sigma=0.5", "sigma=1", "sigma=2"), lty=c(1,2,3), bty= "n", cex=0.8)

データの生成

#Distribution of animals and systematic allocation of transect lines
Nanimals <- 500
Area <- 100^2  
maxPD <- 5

set.seed(2013)
x <- runif(n=Nanimals, min=0, max=100)
y <- runif(n=Nanimals, min=0, max=100)

plot(0, 0, xlim=c(0,100), ylim=c(0,100), xlab="", ylab="", type="n")
points(x, y, pch=20, cex=1, col="gray")

Ntransects <- 3
tx <- c(20, 50, 80)
rectangles <- cbind(rep(2*maxPD, Ntransects), rep(100, Ntransects))
symbols(tx, rep(50,Ntransects), rectangles=rectangles, inches=F, add=T)

#発見プロセス

sigma <- 1
nvec <- numeric(Ntransects)
PDobs.list <- as.list(NULL)
PDobs.list[[1]] <- NULL

Index.Detected <- numeric(Nanimals)
for(i in 1:Ntransects){
 #Perpendicular distance and detection
 PD <- abs(x-tx[i])
 Detected <- as.logical( rbinom(Nanimals, 1, g.HN(PD, sigma, maxPD)) )
 nvec[i] <- sum(Detected)
 PDobs.list[[i]] <- PD[Detected]
 #draw lines and points detected
 lines( rep(tx[i], 2), c(0,100), lty=2)
 points(x[Detected], y[Detected], pch=4, cex=0.4, col="blue")
}

#nvec
PDobs.list

## [[1]]
##  [1] 0.13616423 0.45743668 0.50169476 0.29958731 0.53606140 0.01122308
##  [7] 0.08266618 0.61389015 0.72777497 0.97268858
## 
## [[2]]
##  [1] 2.26330995 0.68141106 0.35452985 0.20079068 1.84240532 0.03090585
##  [7] 0.52523906 0.63804407 1.30995212 0.08796044 0.60562722 0.08136069
## [13] 1.42675790 0.32980486 0.15148339 0.78227120
## 
## [[3]]
##  [1] 1.5050417 0.1410921 0.9914358 0.4067632 0.9260823 0.1051426 1.1646966
##  [8] 1.2144121 0.8109604 0.7468023 1.0707869

PDobs <- unlist(PDobs.list)

発見横距離の分布

hist(PDobs, xlab="観測された横距離")

sigma.est <- sqrt(mean(PDobs^2)); sigma.est

## [1] 0.852965

mu.est <- sqrt(pi/2)*sigma.est; mu.est 

## [1] 1.069033

A <- 100^2
L <- 3*100

N.est <- sum(nvec)*A/(2*mu.est*L)
N.est

## [1] 576.8453

資源量推定値の分布シミュレーション

par(mfrow=c(1,2))

set.seed(2022)

Nanimals <- 500

Nit <- 5000

sigma <- 0.5
mu <- sqrt(pi/2)*sigma

sigma.est <- numeric(Nit)
mu.est <- numeric(Nit)
N.est <- numeric(Nit)

for(h in 1:Nit){

x <- runif(n=Nanimals, min=0, max=100)
y <- runif(n=Nanimals, min=0, max=100)

tx <- c(20, 50, 80)
Ntransects <- length(tx)

nvec <- numeric(Ntransects)
PDobs.list <- as.list(NULL)
PDobs.list[[1]] <- NULL

Index.Detected <- numeric(Nanimals)

if(h <= 5) plot(0, 0, xlim=c(0,100), ylim=c(0,100), xlab="", ylab="", type="n")
for(i in 1:Ntransects){
 PD <- abs(x-tx[i])
 Detected <- as.logical( rbinom(Nanimals, 1, g.HN(PD, sigma, maxPD)) )
 nvec[i] <- sum(Detected)
 PDobs.list[[i]] <- PD[Detected]
 if(h <= 5){
    points(x, y, pch=20, cex=0.4, col="gray")
    rectangles <- cbind(rep(2*maxPD, Ntransects), rep(100, Ntransects))
    symbols(tx, rep(50,Ntransects), rectangles=rectangles, inches=F, add=T)
    lines( rep(tx[i], 2), c(0,100), lty=2)
    points(x[Detected], y[Detected], pch=4, cex=0.4, col="red") 
 }
}
   
nvec
PDobs.list
PDobs <- unlist(PDobs.list)

A <- 100^2
L <- Ntransects*100

sigma.est[h] <- sqrt(mean(PDobs^2))
mu.est[h] <- sqrt(pi/2)*sigma.est[h]
N.est[h] <- sum(nvec)*A/(2*mu.est[h]*L)

if(h <= 5){
 hist(PDobs, col="orange", freq=F, main=paste("N-hat = ",round(N.est[h],1)))
 curve(f.HN(x,sigma.est[h],mu.est[h]), add=T, lwd=2)
}


} #end of for loop in "h"

N.mean <- mean(N.est)
N.sd <- sd(N.est)

logN.mean <- mean(log(N.est))
logN.sd <- sd(log(N.est))

hist(N.est, col="skyblue", freq=F, breaks=50); 
points(Nanimals, 0, pch=19, cex=1.5, col="red")
curve(dnorm(x, N.mean, N.sd), col="red", add=T)

hist(log(N.est), col="skyblue", freq=F, breaks=50); 
points(log(Nanimals), 0, pch=19, cex=1.5, col="red")
curve(dnorm(x, logN.mean, logN.sd), col="red", add=T)