• 新鲜事
  • 《预测核辐射强度的分布》求批

目前,我找到了一个卫星影像数据,精度挺高的 10米乘以10米的,感觉可以用来介绍背景。

Cloud2016 我在 Google Maps 上从左到右随便取了 3 个点,把坐标输进 Open-Elevation API,分别返回了13米、6米和11米,可以把你的网格位置都请求一遍,看看精度是否足够画图:

https://api.open-elevation.com/api/v1/lookup?locations=11.152198,166.842824
https://api.open-elevation.com/api/v1/lookup?locations=11.156872,166.870147
https://api.open-elevation.com/api/v1/lookup?locations=11.173753,166.897657

    nan.xiao 我看了下它的代码,数据源应该是来自 https://srtm.csi.cgiar.org/ 。90 米精度的,意味着 90×9090 \times 90 的区域内是一个海拔值,而朗格拉普岛本来就很小,所以,这个数据源在这里应该不可用。能不能找到更高精度的,比如 10 米的,肯定就够用了。

    6 天 后

    现在遇到几个问题,求助各位。目前,这几个问题以注释/注意的形式放在正文中

    1. nlme 包的函数 Variogram() 拟合结果中,n.pairs 总对数不是 157*156/2 = 12246 而是 156*155/2 = 12090,不知道是哪个地方有问题?函数 Variogram() 帮助文档显示 n*(n-1)/2 对的。
    2. 我将之前的简单克里金插值(无权重约束)换成普通克里金插值(带权重约束),公式是从《Statistics for Spatial Data》(Cressiee 1993) 第 122 页公式 3.2.16 仿写的。全书的数学公式符号与 Cressiee 的书不同,我看了书籍材料之后,手推了一遍,手推仿写的可能有问题,也可能是我的实现代码有问题,理由是普通克里金插值预测方差图看着有问题,见本书图 19.17 的上子图。不知道能否帮我更正?

    我先自报家门,回答两个大家可能关心的问题。也想和大家探讨,我这写作的思路是否合理?

    1. 其实,有很多 R 包可以用来介绍空间数据分析的内容,为什么选用 nlme 包和 spaMM 包?

    空间随机效应是随机效应的一种,本章想将所有内容纳入广义线性混合效应模型框架下,参数的估计和预测都是基于频率派的。nlme 包可以求解带空间随机效应的线性混合效应模型,nlme 包高亮突出广义最小二乘和限制极大似然法,它们是很基础和核心的内容,因此,本章特意选择 nlme 包,但是对空间数据,nlme 包不支持预测,这部分就要手推公式、写代码实现普通克里金插值预测的过程。最后,与 spaMM 包实现的空间广义线性混合效应模型比较。spaMM 包基于似然函数推断的(属于频率派的),限制极大似然法估计方差-协方差矩阵的参数,nlme 包也是。nlme 包是 R 内置的,spaMM 包计算很高效,依赖都相对比较少且轻,安装方便。

    1. 从变差函数拟合、网格划分、克里金插值预测、画图,10-20 行代码就可以了,举例如下,为什么书中写了那么多?

    因不想做调包侠,想搞通一条道路。计算过程都被 R 包函数封装了,而且不同的 R 包大多是不同的人开发的,各有自己的一套理论视角,以 nlme 包作为参照,有的参数估计结果与之完全一样,有的需要做一些转化才能获得一样的结果,有的方法不同只能获得近似的结果,更别提还有贝叶斯派方法。

    library(gstat)
    # 变差模型结构
    vario <- variogram(object = log(counts / time) ~ 1, data = rongelap_sf)
    # 变差模型的协方差结构和参数初值
    fit_vgm <- vgm(psill = 1, model = "Exp", range = 200, nugget = 1)
    fit_vgm
    
    # 拟合变差
    fit_vario <- fit.variogram(object = vario, model = fit_vgm)
    fit_vario
    
    library(abind)
    library(stars)
    rongelap_coastline_sf <- st_as_sf(rongelap_coastline, coords = c("cX", "cY"), dim = "XY")
    rongelap_coastline_sfp <- st_cast(st_combine(st_geometry(rongelap_coastline_sf)), "POLYGON")
    
    rongelap_coastline_buffer <- st_buffer(rongelap_coastline_sfp, dist = 50)
    
    rongelap_grid_stars <- st_bbox(rongelap_coastline_buffer) |>
      st_as_stars(nx = 150, ny = 75) |>
      st_crop(rongelap_coastline_sfp)
    rongelap_grid_stars
    
    
    # 预测
    fit_krige_pred <- krige(
      formula = log(counts / time) ~ 1, locations = rongelap_sf,
      newdata = rongelap_grid_stars, model = fit_vario
    )
    
    fit_krige_pred
    
    library(ggplot2)
    
    # 预测值
    ggplot() +
      geom_stars(
        data = fit_krige_pred, aes(fill = exp(var1.pred)),
        na.action = na.omit
      ) +
      scale_fill_viridis_c(
        option = "C", breaks = 0:12,
        guide = guide_colourbar(
          barwidth = 1, barheight = 15
        )
      ) +
      theme_bw() +
      labs(x = "横坐标", y = "纵坐标", fill = "辐射强度")
    
    # 预测方差
    ggplot() +
      geom_stars(
        data = fit_krige_pred, aes(fill = var1.var),
        na.action = na.omit
      ) +
      scale_fill_viridis_c(
        option = "C", breaks = 0.1 * 0:12 / 4,
        guide = guide_colourbar(
          barwidth = 1, barheight = 15
        )
      ) +
      theme_bw() +
      labs(x = "横坐标", y = "纵坐标", fill = "预测方差")

      如果问题 2 不是问题,就有一套解释:从现在的预测方差图来看,克里金插值欠拟合,采样的地方方差比较大,未采样的地方方差小。spaMM 正好相反,过拟合,采样的地方方差小,未采样的地方方差大。

      空间混合效应模型对我这种数学水平的人来说有点太难了,所以对以上具体的技术问题并没有特别好的思路。不过作为参加过几个正式/非正式书籍项目的人,我可以对如何组织材料提一个大方向的建议:

      • 楼主可以在书中正文里保持简洁,把每个方法的基本假设,如何正确调包,如何解读结果讲清楚即可。
      • 然后,可以把所有需要遵循第一性原理风格的手推公式和底层实现问题在章节末尾给读者留成习题 (exercise),另起一个新的 Quarto book 项目,提供习题解答。
      • 这样,既改善了正文的长度,保持了读者的阅读兴趣,又可以让对底层感兴趣的读者了解原理从头实现,还可以把所有完整细节都放到习题解中,而不必受到书中章节篇幅的限制。

      以上仅供参考,其实目前的长度和风格还可以,可能后半部分的数学公式有点略多 😂

        nan.xiao 这个建议非常好哈 😄 。目前,我故意把大段的数学公式放在后面的,不想深入了解的读者,这部分是可以跳过的。拓展阅读部分,应该画点图表,详细地对比 R 包的。REML 的部分弄完,我就按上面的建议重新捋一下组织结构。我其实最头疼的就是组织结构,因为材料越挖越多,越堆越多,把原来的组织结构弄变形了。

        Cloud2016 对上面第一个问题有初步答案了,函数 Variogram() 有一个选项 collapse 它默认取值是 "quantiles" ,如果设置为 "none" 它会返回每个点对的变差函数值,是 12246 个,符合预期。这个选项 collapse 顾名思义,就是要把每个点对的变差,按照不同的方式卷/分组起来计算。quantiles 意味着根据点对的距离分布,按分位点划分距离区间。至于为什么按照这个方式划分之后,每个区间的点对数加起来少了 156 个,可能要去抠源代码了。

        下面写点体会

        1. 我怀疑我可能对一些细节抠得过细了,比如上面这个例子。写书的过程中,常常会陷入一些细节,因为想把每一个步骤,每一个逻辑,每一段文字保证是对的,至少自己觉得都符合逻辑而且对。写书而不是看书或听课,导致会浮现出来超多的细节,这也有可能是自己本身对问题了解得太少。
        2. 如果把遇到的所有细节抠明白肯定能成为这方面的「专家」,但是写书的成稿日期会无限延长,所以,是要交稿还是要成为专家,需要权衡一下。
        3. 抠细节,很多都涉及统计计算、数值优化和代码实现了,跟数据、分析和统计关系越来越远,写书的内容有走偏的可能。书中,空间线性混合效应模型的似然函数是非凸的、多模态的。为了更深的了解这一点,经过一番查找,我又去复现了 B. D. Ripley (就是那个R 语言核心团队的核心,也是空间统计领域的名人)1987 的一篇论文的结果(极高的引用量),作为数值优化章节的示例。又发现一个有意思的现象,普通优化求解器都无法找到全局最优解,但是 nlme 包的 gls() 函数好像可以,可见其内部除了初值选得好,还干了点别的初始化的事。可见,统计理论与代码实现之间又隔了几条街。

        感觉目前写得有点乱,像是在记笔记不是写书,所以,我把核污染分析同时挂到博客上了,博客是自留地,随便搞,书稿待调整。

        因为感觉写得有点乱,所以,接下来,需要思考和调整一下,关于书籍内容的,写哪些,写到什么程度。欢迎大家讨论、多给意见。

        1. 侧重分析思维:把数据分析的思路和过程突出出来。涉及到的函数和参数讲解一下,背后理论一概不讲,或一言以蔽之,给个参考文献。
        2. 侧重统计理论:参数估计的方法,数学公式推导。比如 nan.xiao 建议放在习题里留待读者自己去研究,另起一个 Quarto 项目。
        3. 侧重软件工具:尽量把用到的每个函数的每个参数都讲清楚作用。

        我感觉在文中一旦突出一个,其他的篇幅必然要压缩,都突出等于没有突出,像现在的稿子。

          Cloud2016 以上是我作为写作者的视角,特别想听听其他读者或写作者的视角。

          Cloud2016 我很赞同“侧重分析思维”这一项。据我收集的小样本,近年来群众反响较好的“实战”类建模书籍有 ISLR 和 François Chollet 的 Deep Learning with Python,可以参考一下他们使用的结构。

            目前的书稿放了 7 个东西,混在一起没有做好拆分,结构不清晰,进而给人混乱的感觉。其实,拆分好之后,就清晰了。

            1. 空间数据分析的思路和流程。
            2. 高斯模型。比如空间线性混合效应模型、经验变差函数拟合、克里金插值预测。
            3. 非高斯模型。广义线性模型、广义线性混合效应模型、空间广义线性混合效应模型。
            4. 理论部分。关于参数估计和模型预测,不同方法的理论。
            5. 不同模型之间的效果比较。高斯模型之间,非高斯模型之间。同一模型的不同参数估计方法比较。
            6. 模型拓展。比如添加块金效应,这些成分导致模型复杂化,参数估计需要一些技巧。
            7. 穿插其中的代码说明。

            有些部分要么删减要么变成习题。总算捋清楚了。感觉是把别人的一本书弄成自己的一章,可不就乱了吗?

            4 个月 后

            克里斯托弗·诺兰的新电影《奥本海默》上映了,奥本海默被称为原子弹之父,与之合作的泰勒被称为氢弹之父,二战后,前苏联、美国、英国、法国、中国,陆续紧锣密鼓地、争先恐后地研制核武器。科学家大都是爱国的,却都不幸沦为政治家、野心家的工具,奥本海默、爱因斯坦、波尔等一批从事理论和实验的著名物理科学家感到深深的内疚、有罪。
            1954 年美国在太平洋比基尼环礁上的那次氢弹核试验是人类历史上规模最大的一次实验,这颗氢弹的当量相当于广岛原子弹的 1000 倍。由于实际当量超出预估的 2.5 倍,气象条件也不利,产生大面积核辐射,一些基地、居民都没来得及撤离,几百公里外的朗格拉普环礁也受到核辐射,半个世纪后,岛上的核辐射强度仍然很高,不适合人类居住。