按照以往例子,还是从最简单的例子讲起。
话说衡量肿瘤生长率是按照多长时间肿瘤长大一倍来算,当然没有直接用这个数,而是用它的倒数,学名叫做reciprocal doubling time(RDT),换句话说,如果RDT=1,说明肿瘤在一年内长大了一倍,RDT=2说明,翻了两番,也就长大4倍,如果-1,表明缩小了一半。
现在,一个美国老兵来找我,说他现在肿瘤的尺寸是多大,他现在已经退役3291天了,请计算下他当初退役时候肿瘤已经生长的可能性有多大,背景是如果肿瘤在服役时候生长,那么政府给的补偿会多。
现在开始处理这个问题了。
首先还是建模,我假设肿瘤以一个恒定的RDT生长,那么如果肿瘤是以一般速度生长,那么在他退役那天,这个肿瘤有多大?
然后作者统计了52个病人的RDT,然后取了一个中值,之后针对当前肿瘤大小、RDT、以及从退役到当前时间长短,计算出他退役时候肿瘤直径大概是6cm。
根据这个结果,显然在退役时候这个老兵已经长出肿瘤了。
让我们认真回想下这个问题,在建模过程中我们可能存在一些问题,导致最后分析中出错。
其中一个问题在于,我们假设RDT是恒定不变的,即便我们改了另外一个模型出来,划分时间段,每个时间段的RDT是随机抽取的,但这个时候每个时间段的RDT仍然是独立的,但实际情况是,连续时间段的RDT可能是相关的,举个例子来说,前一个时间段中RDT如果比较大的话,下一个时间段RDT很大可能性数值也很大,所以:我们需要构建一个序列相关数值。
既然想清楚了,那就让我们开始做吧。
过程是这样的:
(1) 从高斯分布中产生一个相关值,这很容易办到,我们以上一个数值为条件,来计算下一个数值的分布情况
(2) 使用高斯CDF将每个数值转化成它的积累概率
(3) 根据给定的cdf,将每个积累概率转化为相关数值。
当然,光看这个过程似乎不太好懂,少数高手不在此列。相关代码如下。
def CorrelatedGenerator(cdf, rho):
"""Generates a sequence of values from cdf with correlation.
Generates a correlated standard Gaussian series, then transforms to
values from cdf
cdf: distribution to choose from
rho: target coefficient of correlation
"""
def Transform(x):
"""Maps from a Gaussian variate to a variate with the given CDF."""
p = thinkbayes.GaussianCdf(x)
y = cdf.Value(p)
return y
# for the first value, choose from a Gaussian and transform it
x = random.gauss(0, 1)
yield Transform(x)
# for subsequent values, choose from the conditional distribution
# based on the previous value
sigma = math.sqrt(1 - rho**2)
while True:
x = random.gauss(x * rho, sigma)
yield Transform(x)
我觉得这段代码中,挺神奇的是rho,也就是前一时间内,RDT和后一时间段RDT的相关系数。从最后计算结果上来看,使用了相关系数后,所得到的肿瘤生长记录使长得慢的肿瘤长得更慢,长得快的长得更快,从而使其和不适用相关系数所得到的成长记录相比较起来,差别蛮大。