使用Stan估计多边形面积的均值与标准差：模型结果存疑求助

问题描述 #

我拥有某多边形顶点的采样(x,y)坐标（未知真实顶点坐标），想要计算该多边形面积的均值和标准差。自行编写的Stan模型可编译运行，但结果不符合预期，请求排查问题。

原模型代码 #

data {
  int<lower=0> N;            // Number of data points（应为：单个采样的顶点数）
  int<lower=0> C;            // 采样的多边形数量
  matrix[C, N] x;               // x coordinates（每行对应一个采样的N个顶点x坐标）
  matrix[C, N] y;               // y coordinates（每行对应一个采样的N个顶点y坐标）
}

parameters {
  real<lower=0> mean_a;
  real<lower=0> sd_a;
}

model {
  real area = 0.5 * abs(sum(x[1:N-1,] .* y[2:N,] - y[1:N-1,] .* x[2:N,]));

  // Prior for mean_a and sd_a
  mean_a ~ normal(0, 50);
  sd_a ~ cauchy(0, 20);
  
  // Likelihood
  target += normal_lpdf(area | mean_a, sd_a);
}

示例数据集 #

{
  "N": 3,
  "C": 4,
  "x": [
    [1, 2, 1],
    [4, 4.5, 4.7],
    [2, 1.5, 1.5],
    [3.8, 4.2, 4]
  ],
  "y": [
    [1, 2, 1],
    [4, 4.5, 4.7],
    [2, 1.5, 1.5],
    [3.8, 4.2, 4]
  ]
}

问题分析与修正方案 #

核心错误1：面积计算逻辑完全错误 #

原模型混淆了数据维度：x和y是[C, N]矩阵，每行对应一个完整的多边形采样（N个顶点），但原代码直接对矩阵的前N-1行与后N行运算，把所有顶点混排计算出一个“伪面积”，没有分别计算C个采样多边形的真实面积，导致仅生成1个面积值，无法拟合分布的均值和标准差。

正确的面积计算需遍历每个采样，单独计算单个多边形的面积：

vector[C] areas;
for (c in 1:C) {
  areas[c] = 0.5 * abs(sum(x[c, 1:N-1] .* y[c, 2:N] - y[c, 1:N-1] .* x[c, 2:N]));
}

核心错误2：似然函数使用错误 #

原模型仅用1个错误的面积值拟合正态分布，正确做法是用C个采样面积的整体来拟合，即对每个面积值应用正态似然：

target += normal_lpdf(areas | mean_a, sd_a);

次要问题：先验与参数约束不匹配 #

• mean_a加了<lower=0>约束，却使用normal(0,50)先验，会自动截断为半正态分布，直接使用half_normal(0, 50)更清晰。
• sd_a是标准差，必须非负，原代码用cauchy(0,20)需截断，应改为half_cauchy(0, 20)；同时过宽的先验可能导致收敛问题，可根据预期面积规模调整（比如half_cauchy(0, 5)）。

修正后的完整模型 #

data {
  int<lower=3> N;            // 单个多边形的顶点数（至少3个）
  int<lower=1> C;            // 多边形采样数量
  matrix[C, N] x;               // 每行对应一个采样的N个顶点x坐标
  matrix[C, N] y;               // 每行对应一个采样的N个顶点y坐标
}

parameters {
  real<lower=0> mean_a;
  real<lower=0> sd_a;
}

model {
  vector[C] areas;
  
  // 计算每个采样多边形的面积
  for (c in 1:C) {
    areas[c] = 0.5 * abs(sum(x[c, 1:N-1] .* y[c, 2:N] - y[c, 1:N-1] .* x[c, 2:N]));
  }
  
  // 先验分布
  mean_a ~ half_normal(0, 50);
  sd_a ~ half_cauchy(0, 5);
  
  // 似然函数：所有采样面积服从正态分布
  areas ~ normal(mean_a, sd_a);
}

generated quantities {
  // 可选：生成预测的面积值，用于验证
  real pred_area = normal_rng(mean_a, sd_a);
}

额外说明 #

• 需确保每个采样的顶点按顺时针或逆时针顺序排列，否则即使取绝对值，也可能计算出错误的面积。
• 示例数据集中第一个采样的顶点是(1,1),(2,2),(1,1)，属于退化三角形（面积为0），真实采样数据需检查顶点有效性。

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