后验/预测害虫取样的数学模型

题目

假设你是一位在新的环境中进行害虫取样的生态学家。

你在测试地区安置了100个陷阱，第二天回去检查他们。你发现有37个陷阱被触发捕获害虫。一旦陷阱触发，它就不能再继续捕获其他害虫直到复位。如果你重设陷阱，两天内回来，你预期发现多少被触发的陷阱？计算这一后验/预测的分布。

解答：

1、每天触发陷阱害虫数的概率分布为泊松分布

X事件:假设lamda=X为每天平均触发陷阱数量

B事件：一天抓到k=37只虫

P(X)P(B|X)/P(B)=P(X|B)

P(B|X)=

P(X)=1/100

P(B)=

得到不同lamda下的后验概率P(X|B)

 

2、根据lamda的后验概率分布，用不同lamda和其相应的概率计算下次触发陷阱的分布并混合分布，得到下次触发陷阱的总体分布。

 

import thinkbayes
import thinkplot


class Euro(thinkbayes.Suite):
    """Represents hypotheses about the probability of heads."""

    def Likelihood(self, data, hypo):
        """Computes the likelihood of the data under the hypothesis.

        hypo: integer value of x, the probability of heads (0-100)
        data: tuple of (number of heads, number of tails)
        """
        lam = hypo
        k = data
        like = thinkbayes.EvalPoissonPmf(k, lam)
        return like


def MakeGoalPmf(suite, high=100):
        """Makes the distribution of goals scored, given distribution of lam.

        suite: distribution of goal-scoring rate
        high: upper bound

        returns: Pmf of goals per game
        """
        metapmf = thinkbayes.Pmf()

        for lam, prob in suite.Items():
            pmf = thinkbayes.MakePoissonPmf(lam, high)
            metapmf.Set(pmf, prob)

        mix = thinkbayes.MakeMixture(metapmf, name=suite.name)
        return mix


def main():
    suite = Euro(xrange(0, 100))
    heads = 37
    suite.Update(heads)
    print suite.Mean()
    thinkplot.Pmf(suite)
    thinkplot.Show()           #得到图一
    goal = MakeGoalPmf(suite)
    thinkplot.Pmf(goal)
    thinkplot.Show()           #得到图二
    print goal.Mean()


if __name__ == '__main__':
    main()

           图一：不同lamda下的后验概率P(X|B)，=37.999999999963755

       图二: 下次触发陷阱数量的总体分布，=37.999999225214054

参考书籍：《贝叶斯思维：统计建模的PYTHON学习法》