from scipy import optimize from scipy import integrate import scipy as sp import numpy as np from math import exp Non=15 Noff=15 Ton=10 Toff=50 q=0.95 def PoissonPosterior(s): value=0 s=np.maximum(s,np.zeros(np.shape(s))) Norm=0 for i in range(1,Non+1): coeff=(1.+Toff/Ton)**i*sp.factorial(Non+Noff-i-1)/sp.factorial(Non-i) Norm+=coeff #print coeff,s,Ton,i value=coeff*(s*Ton)**(i-1)*np.exp(-s*Ton)/sp.factorial(i-1) return value*Ton/Norm def solve_me(x,func): # x is an array of the values you are solving for a,b = x fa=func(a) fb=func(b) integral_error = integrate.quad(func, a , b)[0] - q prob_difference = fb - fa print a,b,integral_error, prob_difference return np.array([integral_error, prob_difference]) results=optimize.fsolve(solve_me, [1., 1.5],args=(PoissonPosterior)) print "%d%s credible interval: [%.2f,%.2f]"%(q*100,"%",results[0],results[1]) ################## def density(s): if s<0.5: return 4*s else : return 4-4*s #for testing : results=optimize.fsolve(solve_me, [0.5, 0.5],args=(density)) # test that the computation of the pdf is correct from pylab import * x=np.arange(0,3,0.001) y=PoissonPosterior(x) plot(x,y) show()