『ポアソン分布・ポアソン回帰・ポアソン過程』をRで再現してみる(6)

非定常ポアソン過程です。

図5.2(b)の、密度が変化する点配置を再現してみます。

library(ggplot2)
library(dplyr)
library(pracma)
library(MASS)

set.seed(123)

### 5.3
r <- function(x, y) {
  (-0.00002 * (x - 40)^2 + 0.05) * 0.04 * y
}

x.min <- 0
x.max <- 100
y.min <- 0
y.max <- 100
opt <- optim(par = c(1, 1),
             fn = function(par) -r(par[1], par[2]),
             method = "L-BFGS-B",
             lower = c(x.min, y.min),
             upper = c(x.max, y.max))
r.S <- -opt$value
N <- rpois(1, r.S * (x.max - x.min) * (y.max - y.min))
df <- data.frame(x = runif(N, x.min, x.max),
                 y = runif(N, y.min, y.max),
                 r = runif(N, 0, 1)) %>%
  filter(r(x, y) / r.S > r)

## 強度
integral2(r, x.min, x.max, y.min, y.max)

$Q
[1] 626.6667

$error
[1] 2.131628e-14

## 点の数
nrow(df)

[1] 654

##
ggplot(data = df, mapping = aes(x = x, y = y)) +
  geom_point() +
  xlim(x.min, x.max) +
  ylim(y.min, y.max) +
  coord_equal() +
  theme(aspect.ratio = 1)

2次元カーネル密度推定値と重ねてみます。2次元カーネル密度推定には、MASSパッケージのkde2d関数を使用しました。

### 5.6
## Kernel Density
n <- c(100, 100)
lx <- (x.max - x.min) / n[1]
ly <- (y.max - y.min) / n[2]
z <- kde2d(x = df$x, y = df$y, n = n,
           lims = c(x.min + lx / 2, x.max - lx / 2,
                    y.min + ly / 2, y.max - ly / 2))

df2 <- data.frame(x = rep(z$x, n[2]), y = rep(z$y, each = n[1]), z = c(z$z))
ggplot(data = df2, mapping = aes(x = x, y = y, fill = z), alpha = 0.5) +
  geom_tile() +
  geom_point(data = df, mapping = aes(x = x, y = y),
             colour = "black", fill = "black") +
  coord_equal() +
  theme_bw() +
  theme(aspect.ratio = 1)

タグ：R