From 4b60116a0da65614070953367372c15876b677f0 Mon Sep 17 00:00:00 2001 From: grayson rich Date: Tue, 13 Dec 2016 15:48:57 -0500 Subject: [PATCH] simple TOF model test actually works - o, glorious day.. --- tests/mcModelIntegration.py | 20 +++-- tests/simpleTOFmodel.py | 146 +++++++++++++++++++++++++++++++----- 2 files changed, 142 insertions(+), 24 deletions(-) diff --git a/tests/mcModelIntegration.py b/tests/mcModelIntegration.py index 1f6e48c..1719c67 100644 --- a/tests/mcModelIntegration.py +++ b/tests/mcModelIntegration.py @@ -5,6 +5,8 @@ # # intended to demo and develop MC-based integration of models # needed to numerically evaluate likelihoods / probabilities +# +# all the 'magic' here is in the functions lnlike and generateModelData import numpy as np from numpy import inf @@ -32,6 +34,13 @@ distance_cellLength = 2.86 # cm, length of gas cell distance_zeroDegLength = 3.81 # cm, length of 0deg detector +############## +# vars for binning of TOF +tof_nBins = 25 +tof_minRange = 175.0 +tof_maxRange = 200.0 +tof_range = (tof_minRange,tof_maxRange) + def getDDneutronEnergy(deuteronEnergy, labAngle = 0): """Get the energy of neutrons produced by DDN reaction @@ -100,7 +109,7 @@ def lnlike(params, observables, nDraws=2000000): """ #print('checking type ({}) and length ({}) of params in lnlikefxn'.format(type(params),len(params))) evalData=generateModelData(params, nDraws) - evalHist, evalBinEdges = np.histogram(evalData[:,3], 25, (175.0,200.0), + evalHist, evalBinEdges = np.histogram(evalData[:,3], tof_nBins, tof_range, density=True) logEvalHist = np.log(evalHist) #print(logEvalHist) @@ -179,7 +188,8 @@ def lnlike(params, observables, nDraws=2000000): plot.xlabel('Location in cell (cm)') plot.show() -histTOFdata, tof_bin_edges = np.histogram(fakeDataSet[:,3], 25, (175.0,200.0), +histTOFdata, tof_bin_edges = np.histogram(fakeDataSet[:,3], tof_nBins, + tof_range, density=True) plot.figure(20) plot.scatter(tof_bin_edges[:-1], histTOFdata, color='k') @@ -201,8 +211,8 @@ def lnlike(params, observables, nDraws=2000000): # make our vector of 'n' -observedVectorN, observed_bin_edges = np.histogram(fakeObsData[:,3], 25, - (175.0,200.0)) +observedVectorN, observed_bin_edges = np.histogram(fakeObsData[:,3], tof_nBins, + tof_range) loghist = np.log(histTOFdata) @@ -255,7 +265,7 @@ def lnlike(params, observables, nDraws=2000000): plot.figure() plot.subplot(221) -plot.hist(fakeObsData[:,3], 25, (175.0,200.0)) +plot.hist(fakeObsData[:,3], tof_nBins, tof_range) plot.subplot(223) plot.scatter(multiplier,nll_sigmas) plot.xlabel('fraction sigma TRUE') diff --git a/tests/simpleTOFmodel.py b/tests/simpleTOFmodel.py index 5fdf344..875dd03 100644 --- a/tests/simpleTOFmodel.py +++ b/tests/simpleTOFmodel.py @@ -7,6 +7,7 @@ # doesn't include all the realistic features needed to do actual analysis import numpy as np +from numpy import inf from scipy.integrate import quad import scipy.optimize as optimize import matplotlib.pyplot as plot @@ -31,6 +32,13 @@ distance_cellLength = 2.86 # cm, length of gas cell distance_zeroDegLength = 3.81 # cm, length of 0deg detector +############## +# vars for binning of TOF +tof_nBins = 25 +tof_minRange = 175.0 +tof_maxRange = 200.0 +tof_range = (tof_minRange,tof_maxRange) + def getDDneutronEnergy(deuteronEnergy, labAngle = 0): """Get the energy of neutrons produced by DDN reaction @@ -58,7 +66,69 @@ def getTOF(mass, energy, distance): return tof -def evaluateModel() +def generateModelData(params, nSamples): + """Generate some fake data from our mode with num samples = nSamples + params is an array of the parameters, [e0, e1, sigma] + Returns a tuple of len(nSamples), [x, ed, en, tof] + """ + initialEnergy, eLoss, sigma = params + data_x=np.random.uniform(low=0.0, high=distance_cellLength, size=nSamples) + data_ed= np.random.normal(loc=initialEnergy + eLoss*data_x, + scale=sigma) + data_en = getDDneutronEnergy(data_ed) + + neutronDistance = distance_cellToZero + (distance_cellLength - data_x) + neutronTOF = getTOF(mass_neutron, data_en, neutronDistance) + effectiveDenergy = (initialEnergy + data_ed)/2 + deuteronTOF = getTOF( mass_deuteron, effectiveDenergy, data_x ) + data_tof = neutronTOF + deuteronTOF + + data = np.column_stack((data_x,data_ed,data_en,data_tof)) + return data + +def lnlike(params, observables, nDraws=1000000): + """Evaluate the log likelihood given a set of params and observables + Observables is a vector; a histogrammed time distribution + Params is a list of [initial D energy, E_D loss, sigma] + nDraws is the number of points drawn from (energy,location) distribution\ + which are used to produce the PDF evaluations at different TOFs + """ + #print('checking type ({}) and length ({}) of params in lnlikefxn'.format(type(params),len(params))) + evalData=generateModelData(params, nDraws) + evalHist, evalBinEdges = np.histogram(evalData[:,3], tof_nBins, tof_range, + density=True) + logEvalHist = np.log(evalHist) + #print(logEvalHist) + # find what TOFs have zero observed data + # we'll use this to handle cases where we might wind up with -inf*0 + # likelihood is fine if PDF is 0 somewhere where no data is found + # without checks though, ln(PDF=0)=-inf, -inf*0 = nan + # however, if PDF is 0 (lnPDF=-inf) where there IS data, lnL should be -inf + zeroObservedIndices = np.where(observables == 0)[0] + for idx in zeroObservedIndices: + if logEvalHist[idx] == -inf: + logEvalHist[zeroObservedIndices] = 0 + + loglike = np.dot(logEvalHist,observables) + return loglike + + + +def lnprior(theta): + e_0, e_1, sigma = theta + if 800 < e_0 < 1200 and -200