贝叶斯统计:仿真

         按照以往例子，还是从最简单的例子讲起。

         话说衡量肿瘤生长率是按照多长时间肿瘤长大一倍来算，当然没有直接用这个数，而是用它的倒数，学名叫做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的相关系数。从最后计算结果上来看，使用了相关系数后，所得到的肿瘤生长记录使长得慢的肿瘤长得更慢，长得快的长得更快，从而使其和不适用相关系数所得到的成长记录相比较起来，差别蛮大。