在R中基于dwp包拟合的xep02分布抽样的技术求助

通过dwp包拟合得到xep02概率密度函数，拟合参数如下：

> Kbatmod$xep02
Distribution: xep02
Formula: ncarc ~ log(r) + I(r^2) + offset(log(exposure))

Parameters:
           b0            b2 
-0.8396640654 -0.0004653923 

Coefficients:
  (Intercept)        log(r)        I(r^2) 
-2.1154557406 -0.8396640654 -0.0004653923 

Variance:
              (Intercept)        log(r)        I(r^2)
(Intercept)  3.034712e-01 -1.111308e-01  4.550014e-05
log(r)      -1.111308e-01  4.608998e-02 -2.451855e-05
I(r^2)       4.550014e-05 -2.451855e-05  2.531856e-08

需要围绕(0,0)生成点云：方位角取0-360随机值，距离符合xep02分布，但dwp包无对应抽样函数，使用RVCompare的sampleFromDensity时报错：

> dwpPDF <- Kbatmod$xep02
> PDFsamples <- sampleFromDensity(dwpPDF, 100, c(0,100))
Error in density.default(X[[i]], ...) : 
  need at least 2 points to select a bandwidth automatically

1. 手动定义xep02的概率密度函数（PDF） #

xep02分布的强度函数为λ(r) = exp(β₀ + β₁*log(r) + β₂*r²)，概率密度函数是强度函数除以0到最大距离（这里取100）的积分归一化常数。从拟合结果提取系数后定义PDF：

# 提取模型系数
beta0 <- -2.1154557406
beta1 <- -0.8396640654
beta2 <- -0.0004653923

# 定义xep02概率密度函数
xep02_pdf <- function(r) {
  lambda <- exp(beta0 + beta1 * log(r) + beta2 * r^2)
  # 计算归一化常数（0到100的积分，下限取0.001避免log(0)报错）
  norm_const <- integrate(function(x) exp(beta0 + beta1 * log(x) + beta2 * x^2), lower = 0.001, upper = 100)$value
  lambda / norm_const
}

2. 用逆变换抽样生成距离样本 #

通过累积分布函数（CDF）的逆函数实现抽样，步骤如下：

# 定义CDF函数
xep02_cdf <- function(r) {
  integrate(xep02_pdf, lower = 0.001, upper = r)$value
}

# 定义逆CDF函数（用uniroot求解）
xep02_inv_cdf <- function(p) {
  uniroot(function(r) xep02_cdf(r) - p, interval = c(0.001, 100))$root
}

# 生成100个符合xep02分布的距离样本
set.seed(123) # 可选，保证结果可复现
n_samples <- 100
u_samples <- runif(n_samples)
distance_samples <- sapply(u_samples, xep02_inv_cdf)

3. 生成随机方位角并转换为笛卡尔坐标 #

# 生成0-360度的随机方位角（转换为弧度）
azimuth_rad <- runif(n_samples, 0, 2 * pi)

# 转换为以(0,0)为中心的笛卡尔坐标
x_coords <- distance_samples * cos(azimuth_rad)
y_coords <- distance_samples * sin(azimuth_rad)

# 整理为数据框
point_cloud <- data.frame(X = x_coords, Y = y_coords)

4. 导出到ArcGIS并统计点数 #

• 将点云数据导出为CSV：

  write.csv(point_cloud, "xep02_point_cloud.csv", row.names = FALSE)
  

• 在ArcGIS中导入CSV为点图层，使用Count Points in Polygons工具或空间连接功能，统计落入预定义多边形内的点数。

报错原因说明 #

RVCompare::sampleFromDensity需要输入数值型的密度数据（如向量形式的密度值），而非dwp包的模型对象，因此直接传入Kbatmod$xep02会导致函数无法识别数据结构，触发带宽选择的错误。手动实现抽样逻辑更可靠且无需依赖第三方包的兼容格式。

内容的提问来源于stack exchange，提问作者BugaB