L-функция и ggplot в R. Я не вижу, где я ошибся - PullRequest
0 голосов
/ 13 апреля 2020

Я пытаюсь следовать этому туториалу, чтобы сделать симпатичный ggplot всего распределения (всего диапазона доверительного интервала), как они сделали.

Однако мой теоретический L-функция, кажется, сильно расходится, когда я ее строю. P в туториале - это мой poisson1, а L в туториале - это L_CSR1 в моем коде ниже. Когда я запускаю код ниже, я получаю следующий рисунок:

enter image description here

, который немного отличается от приведенного в турториале. Кажется, что значения моей черной линии неверны. Если я удаляю только кодовую строку (которая удаляет только значения L_CSR1)

+ geom_line(aes(r, iso), data = (L_CSR1), size = 1, col = "#707070")

снизу, по крайней мере, я получаю красивую фигуру, показывающую только доверительные полосы:

enter image description here

Я провел весь день, пытаясь выяснить, где проблема, но независимо от того, что я изменяю, я не могу понять это правильно. Я очень ценю, если кто-то может взглянуть на это.

Вот мой код:

library(spatstat)
library(raster)
library(reshape)
library(ggplot2)
library(maptools)

nclust = function(x0, y0, radius, n) {
  return(runifdisc(n, radius, centre=c(x0, y0)))
}

poisson1 = rPoissonCluster(40, 0.2, nclust, radius = 0.2, n = 5, win = owin(c(0,1.31),c(0,1)), lmax = NULL, nsim = 1, drop = TRUE, saveparents=TRUE)

L_CSR1 = Lest(poisson1, correction = "Ripley")

n = 9999
Lm = matrix(nrow = n, ncol = length(L_CSR1$r))        # Define an empty matrix
for(i in 1:n){
  Pmc = runifpoint(poisson1$n, win = owin(c(0,1.31),c(0,1)))
  Lmc = Lest(Pmc, nsim = n, border = "Ripley")      # Compute L
  Lm[i,] = Lmc$iso - Lmc$r                          # Set Ltheo to 0 by rotating plot
}


# Create a ppp object from Lm
Lmp <- ppp(rep(L_CSR1$r, each = n), as.vector(Lm), window = owin( c(min(L_CSR1$r), max(L_CSR1$r)), c(min(Lm), max(Lm))))

# Define an empty raster (grid). The extent is defined by the distance bands r and the
# minimum and maximum L-values.
r <- raster(nrows = 150, ncols = length(L_CSR1$r), xmn = min(L_CSR1$r), xmx = max(L_CSR1$r), ymn = min(Lm), ymx = max(Lm))

# Rasterize the ppp object. This creates a raster 'density' object where each cell (pixel) is assigned 
# the number of points.
poisson1.r <- rasterize(as.SpatialPoints.ppp(Lmp), r, field = 1, fun = sum)

# Next, convert the raster object to a matrix object. This step is needed if we want to normalize
# the cell values at their respective band distances.
poisson1.rm <- as.matrix(poisson1.r)

# Convert column values to normalized values
N <- apply(poisson1.rm, 2, function(x) {Fn = ecdf(x); Fn(x)})

# Convert matrix back to a raster
N.r <- raster(N, xmn = min(L_CSR1$r), xmx = max(L_CSR1$r), ymn = min(Lm), ymx = max(Lm))

# Use ggplot to plot the results

r.df <- as.data.frame(cbind(xyFromCell(N.r,1:ncell(N.r)), getValues(N.r)))
r.df <- na.exclude(r.df)
plotData <- melt(r.df, id = c("x","y"))
cols <- c('#AA0000', '#F03B20', '#FD8D3C', '#FECC5C', '#FFFFCC') 
brks <- c(0.05, .2, 0.5, 0.8, 0.95) 
ggplot(aes(x = x, y = y), data = plotData) + 
  geom_tile(aes(fill = value)) + 
  scale_fill_gradientn(colours = cols, breaks=brks) +
  theme(legend.position = "right") + xlab("r") + ylab(expression(hat("L"))) +
  geom_line(aes(r, iso), data = (L_CSR1), size = 1, col = "#707070")
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...