第二届东亚天文统计学研讨会笔记

我参加了Eric Feigelson老师主讲的第二届东亚天文统计学研讨会，这里留个笔记，以免我忘记
我之所以知道Eric老师，除了因为他写了天文统计教材外，也是因为看到过他在八十年代写的如何理解上下限数据点的文章，吃饭的时候跟Eric聊起他这个工作的时候，Eric说哦你还看过这个啊，我过几天会讲
课件都能在https://zenodo.org/record/1290840#.W1FQ3_kzbg8这里找到，教材是Eric老师和G. Jogesh Babu老师写的 Modern Statistical Methods for Astronomy: With R Applications 简称MSMA
Eric老师说，现代统计学和天文学的关系大概是这么回事儿：
Eric有只拖拉机他从来也不骑，有一天心血来潮骑着去赶集，他骑着拖拉机兔兔兔兔心里正得意，突然发现拖拉机停下来了，只好从工具箱找扳手改锥，拆开外壳仔细一看，链子断了，这好办，接上就好了，可是车链子没法用扳手改锥接上，倒是能用铁丝栓在一起，但是车开着不太舒服，报修后来了个修理工，带了一箱子工具，拿出来个他不认识的东西，断了的链子两头放一起，一夹，链子接上了，修好了，就这样
Eric能看明白当时的情况，但是没有合适的工具，因此一点也奇怪自己修不好拖拉机，Eric说，天文学研究经常能看到断了的链子，但是天文学家不经常有合适的工具，而统计学家早就有了非常合适的工具，但是他们不知道天文学家在哪里掉链子
从这个例子出发，Eric介绍了统计学研究的现状，比如统计和数学的关系，比如各种各样上世纪后期发展起来的分析方法，比如找系外行星这种事情，傅里叶变换当然可以找到光变周期认证掩食等等，ARIMA也可以找到周期，而且效率更高
种种迹象表明，统计学家研究出来的这么多的方法可以使我们越来越对得起数据
可是现代统计方法有很多不友好的地方，比如说很复杂很难学，为了解决这个问题，Eric隆重推荐了R，由于R有很多写好的库，可以方便的进行统计分析，R的库背后都有学术论文论证其合理性，因此R非常适合用于做现代统计分析
图片来自Eric老师的课件
Eric在第一天第一个上午介绍了这些现状后，很耐心的讲了讲R的一些语法作为热身，比如说怎么画散点儿图，怎么画轮廓等等，这就算认识了R，就可以操练了
这次讲座的形式是，老师先讲会儿课，之后带着学生用R体会一下现代统计方法，老师讲的东西大约是这些，当然这门课涉及到的话题远不止这些：
Local regression
这个部分主要是讲一些非参数统计，要有让数据为自己代言的感觉( to let the data speak for themselves — R. A. Fisher)
比如这种散的数据点，也没个规则说是什么分布，非参数估计就可以根据现在这些数据点的分布总结出一条红线和虚线来展示数据分布特征，做法就是上面那几行程序，载入了库和数据后，一个命令就分析好了
Statistical inference
这部分主要讲统计的需求，比如得到了一个样本后，要做哪些分析，要如何验证分析的可靠性，比如说KS检验，CM检验和AD检验的适用范围，比如bootstrap,一个好的估计应该具有的良好品质，最小开方拟合，极大似然估计和贝叶斯
老师隆重推荐了 Kelly 的
Some Aspects of Measurement Error in Linear Regression of Astronomical Data http://adsabs.harvard.edu/abs/2007ApJ…665.1489K
这篇文章
Regression
这部分主要解释如何找到合适的模型拟合数据，其过程要有物理分析，模型选择，拟合，拟合后各种检验，确认模型的确对得起数据
图片来自Eric老师的课件
最后一段最后那点儿是Eric的朋友，统计系的Babu老师补充的，说这样才是一场认真的统计分析
Data Mining
由于天文学家不停的在对观测到的数据做分类，比如星系/恒星分类，比如超新星分类，比如图像选源等等，这就有了各种各样的方法，比如依据SED形状发展起来的color-color图选源，比如从图像形态做区分的点源展源等等，现代统计学在这个课题上的进展一个是神经元网络，机器学习的应用，使得我们可以简单的对海量数据进行有效分类，一个是聚集程度的分析
这里推荐的是Everitt et al. Cluster Analysis
Bayesian
这个环节比较大，讲了好久好久，跟你们学的差不多，Eric老师评论说
很多天文学家都说自己在用贝叶斯估计，但是其实就是极大似然估计，如果你真的只是想得到个拟合结果和误差的话，最小开方拟合的方法都能给个误差，请不要滥用 贝叶斯方法 这个词，老老实实承认自己其实就是做了个极大似然估计，让 贝叶斯方法 这个词真正用在那些需要这个方法的课题上
推荐的是一本新书：
Bayesian Models for Astrophysical Data using R, Python, JAGS and Stan
Time Series Analysis
传统的方法当然就是自相关函数，傅里叶变换，Eric老师介绍了新的方法叫做
Autoregressive moving average model，以及一些别的方法，都有现成的R程序库可以直接调用，非常方便，Eric也讲了一下他组里用这种方法处理开普勒卫星光变数据找到很多原来没注意到的系外行星的事情，预计今年文章就会发表
Censoring and truncation in astronomical surveys
天文学研究中经常有上下限的数据，如何充分利用这类数据也是个不大不小的问题，我读博时的室友，高能所统计第一人，华北油田第一前锋 L博士 的第一个科研课题就是如何面对这类数据，于是我也是读博时就早早了解到这个话题，却久久没有得到个说法，所以一直比较在意
今天Eric先讲了个段子：他还在宾夕法尼亚大学念书的时候，遇到了有上限的数据，并且需要做个拟合，然而他不会
Eric 1985,Statistical methods for astronomical data with upper limits
有一天他在学校里闲逛，看到统计系，就进去问了一下说你们知道这种只有上限没有误差棒的数据点怎么参与拟合吗，统计系的人非常惊讶，从书架上抽出一本书说知道啊，而且我们刚写了本书，封面几乎就是这张图，然后建议他 Babu 老师聊聊，于是Eric老师和Babu老师这一聊就聊了三十年，聊的停不下来
Eric老师用统计系的老师提供的分析方法写了 Statistical methods for astronomical data with upper limits 这篇文章，他说其实他当时不懂统计，就是把统计专业的东西搬过来给天文学家看了一下，写完文章他又去统计系，说你们还会统计点儿别的什么东西吗 有功夫儿咱们掺和掺和，行嘛，掺和掺和就是热惹热惹
这种上下限的数据点，更有统计感的说法叫做Truncated Sample或者Censored Data，处理这种数据的统计工具，叫做survival analysis，中文叫生存分析
这个名字的来源据说是这种感觉：在人口普查这类行当里，要统计一些东西，比如一年吃多少大米，说起来简单，就看看这家人一年买多少米就好了，可是有时候年初这家人的个数和年底家人个数不一样，这时候就有了上下限，调查人口数量的时候会有存活率之类的概念，这种有上下限的分析也就叫做了生存分析
有了工具，当然就可以分析了，比如Eric的俩文章:
http://adsabs.harvard.edu/abs/1985ApJ…293..192F
http://adsabs.harvard.edu/full/1986ApJ…306..490I
由于探测总会有个极限，在光度函数这类统计工作中也需要考虑这种有截断的数据对结果的贡献，这里Eric吐槽说天文学家总是把描述光度函数的这个Gamma函数叫做Schechter函数
除了介绍了一些近代光度函数的处理方法，Eric老师特意表扬了
Lynden-Bell老师1971年 A method of allowing for known observational selection in small samples applied to 3CR quasars 文章的附录，说这个做法很好很难得
推荐的书是：
Lee, E. T. & Wang, J. W. (2013) Statistical Methods for Survival Data Analysis, 4th ed., Wiley
Helsel, D. R. (2012) Statistics for Censored Environmental Data Using Minitab and R, 2nd ed., Wiley
Moore, D. F. (2016) Applied Survival Analysis Using R, Springer
False-Discovery Rate
FDR控制失误率的方法，可以用在选源等鉴定相关的工作里
最初的文章是：Controlling the False-Discovery Rate in Astrophysical Data Analysishttp://adsabs.harvard.edu/abs/2001AJ….122.3492M
后续的工作有很多，比如 Vio, R.; Vergès, C.; Andreani, P.的
The correct estimate of the probability of false detection of the matched filter in weak-signal detection problems . II. Further results with application to a set of ALMA and ATCA datahttp://adsabs.harvard.edu/abs/2017A%26A…604A.115V
等等，我这里写这个是因为Vio被Eric好评了一下，说Vio出品，必属精品
还有一些看标题感觉很有意思的文章:
DUCHAMP: a 3D source finder for spectral-line data
Noise-based Detection and Segmentation of Nebulous Objects
Darth Fader: Using wavelets to obtain accurate redshifts of spectra at very low signal-to-noise
Darth Fader这个文章作者之一Starck, J被Eric反复点赞，Eric说他是个牛人，他写的东西都是完美的，我说对我知道他，他写了一本图像处理的书叫 Astronomical Image and Data Analysis ，我怎么都看不懂
最后一个报告里，Eric老师讲了一些天文学家常犯的错误，比如胡用KS检验，过分的resample，不得体的用回归分析，被遗忘的残差分析，不充分的模型选择，滥用贝叶斯分析，乱用误用friends of friends分析方法等等
回到第一个故事，Eric想说的是，统计学家已经研究出了很好的工具用于分析各种各样的数据，天文学家有大量的数据可以做统计分析，R语言就是统计学家的工具箱，方便天文学家使用，从这几天体验来看，R至少很适合体验统计方法
最后一个环节是领毕业证并且和Eric合影留念
主要内容就这些了，他还推荐了这个网站：https://asaip.psu.edu/ 我放在阅读原文了
这几天我问了Eric很多问题，还有印象的是：
1，误差棒不对称的数据怎么参与拟合
答：这个很困难，因为违反了一些拟合的基本需求，所以我要说我不知道
问：那我能不能用模拟的方法，假设个分布，然后沿着误差棒撒一堆点作为假的数据，假设这些点的误差都一样，拟合很多次，这样这样的
答：可以啊，反正没辙，这么着也是个办法
问：我是不是可以更进一步的说，所有误差分析之类的事情，用模拟误差的方法总能得到靠谱的误差，一方面误差传递公式的近似要在 误差比较小 的时候成立，另一方面有时候很难分析出来个误差，尤其是信息不是那么全的时候
答，对对，模拟的方法得到误差最靠谱了，好多量的误差来源太复杂了，想不清楚
2，如果数据点误差不独立怎么办
终于解释明白了在经历了傅里叶变换的干涉数据是一张照片，各个pixel是不独立的
答：还真是啊，这个我没想过，不知道啊，我想想统计学家会怎么办
(此时Eric突然沉默，仿佛变身成了宾大统计教授)
可能FDR会有帮助，你看一下这些这些 (就是上面FDR的那个文章)
3，数据需求量如何估计
比如说我有一个星系样本，就5个星系，研究出来个结论的话，是不是不太可信， 因为统计量太少了，我想画个轮廓图，然后就发现没啥好轮廓的，就5个点儿，干脆画五个圈儿吧
答：哈哈哈 是啊，5个听起来太少了
问：如果太少的话，再多观测一个星系，是不是对现状也没什么改变
答：对啊，也才六个，还是画圈儿
问：那我要看多少星系，得到个结论就显得比较可靠了，这有个判断准则没有？我写观测申请经常要论证这个样本的星系个数需求，但是总觉得说的不好，我经常说需要10个，因为10/根号10听起来大于3了 (我好像在统计课上学过一个判据说样本容量的，是陈希孺那本书，我有空去书店找找看)
答：没有，你要自己定这个准则，五个我也觉得少，十个我觉得还好，可是画contour还是不够，但是五个的确太不够了
问：另一方面比如说我有一百万个星系的样本，我想研究点性质，有没有个什么大道理能直接帮我从这个样本里选100个星系，代表一百万星系性质，这样只要研究100个星系就好了，远离大数据困扰
答：这好像是同一个问题的另一方面，而且是个好问题，我也确实不知道
4，神经网络的方法非常高效，可是会不会导致结果过于依赖training sets？
自己答：会的，这就是看你怎么用这个工具，至少我们能把和trainning sample差不多的东西找出来并且分类，剩下样本如果不多了，再用别的方法分析一下就可以了，如果剩下无法分类的太多，就说明training sets不是特别代表这个样本，需要做更多观测得到好的trainning sets
Eric答：会的，所以 trainning sample 很重要，当然配合别的分析方法也很重要
5，一个充满了Truncated Sample的样本的median怎么算，数据点很弥散的时候我见到有一些天文科技工作者会画一个巨大的点和误差棒，展示一些整体趋势
答：不知道啊，没想过这个
肯定不能是直接就拿着数字median了，是不是可以把这些上限的数据点合在一起，做一个测量，不过还是要确认没有胡叠加数据，别把不同的噪音来源的数据简单相加，然后你试试
另外你说的事情好像是在re-sample，这个很危险啊一定小心
6，关于R的黑盒子感觉
问：我自从做观测以来，被告知如果不了解细节的话，不要随便相信跑黑盒子程序的结果，其实我不排斥黑盒子，但是如果我似乎完全不知道这种回归是怎么做的，也无从通过一些简单的测试判断是不是准确，尽管结果看着挺顺眼的
答：这个建议非常old school，我更建议你试试看，至少做个参考，毕竟R是一堆统计专家写的，每个库都有同行评审的文章做背景，还是挺靠谱的
问：那python也有这种库的话，你是不是就会强烈推荐python了？
答：我不清楚有没有，我也不太会用python，对我来说R已经足够用了
7，意义在哪里
问：比如您今天讲的直线拟合，不同的方法拟合的结果其实差不多，而且天文数据一直在积累，所以就算现在拟合的特别好，来了新的数据也仍然是同样的分布但是还是可能会有甚至超出误差范围内的变化，这样以来如果不同的分析手段得到的结果互相都在误差范围内，或者改变真的不大的话
具体到我的经验来说，我发现即使我把开方拟合的公式写错了，比如我没有考虑误差，都能几乎得到合理考虑各种误差的拟合结果，这么高级的分析手段是不是显得大材小用没有必要了
答：首先，我承认的确没有改变很多，只是说这样分析了以后对数据的利用更充分，其次我们不用一下新的工具的话，就永远不知道会是什么样的
问：现代统计方法是不是只会把文章结果从好变到更好，而不是从不能变成能，比如说选源，我们用3sigma判据来选源的话，那么不论什么方法，都不能选来选去的选出个2sigma的源，愣说有3sigma的置信概率，听起来对吗
答：对啊
问：如果我们只研究好的数据，用安全的误差棒，即使用一百年前的分析方法，总不至于结果是错的对吧
答：对啊
问：现代统计方法再怎么研究，也不至于得到个完全相反的结论吧，更何况天文学研究中有时候不做回归分析，甚至所谓线性拟合都做的很不认真，也不怎么做假设检验，因为相关的趋势是一看就知道的，如果看着都不明显就算相关显著性够高，也还是没谁敢认真的用这种相关，而这相关背后的物理也完全不会因为拟合的斜率更准了一点点儿就理解的更深刻了一点
而对于统计分析来说，如果分析了信噪比不是很高的数据的话，即使方法很优秀，哪怕陈独秀，这个结论仍然会等待更多的观测来验证，而不是统计分析了一场，之后就相信了
那么现代统计分析工具的必要性在哪里呢？
答：这是个好问题，我想说首先天文学家都太聪明了，也的确是一眼就看出来相关，哈哈哈哈哈
我是个编辑，是个统计编辑，我平时的工作就是如果有谁的文章统计方法显得花里胡哨的，审稿人评审后还会发到我这里，我不讨论任何物理，就点评其统计方法。这是我的一个工作
就我看到的这些文章，坦率的说，就像你说的，极少工作是统计方法 错了 导致文章结论需要改变，大部分时候是统计方法不合理，按照我的建议改了以后文章结论也的确没怎么变，有时候文章作者中的老年人会写邮件跟编辑抱怨说怎么搞的又来了个统计的审稿人，我写了一辈子文章从没听说过，杂志社也讨论过我这份儿工作存在的意义，有一次开会，屋子里坐了三十来个编辑，我们讨论这个事儿，我就问各位编辑，我的存在是不是使得文章质量变高了，文章科学意义更大了，大部分编辑都认可我的说法，也认为我的确使文章的统计部分写的更规范了，所以我还有这个工作，哈哈哈哈哈哈哈