logsum.py

# Illustrates why log-sums are better than direct sums for probability integrals
# when probabilities are initially computed on the log scale

import pymc
import numpy 

lp_big = numpy.random.normal(loc=0,scale=2,size=1000)
print numpy.log(numpy.mean(numpy.exp(lp_big)))
print pymc.flib.logsum(lp_big)-numpy.log(len(lp_big))

lp_small = numpy.random.normal(loc=-3000,scale=2,size=1000)
print numpy.log(numpy.mean(numpy.exp(lp_small)))
print pymc.flib.logsum(lp_small)-numpy.log(len(lp_small))