gns.nested_sampler module

exception gns.nested_sampler.ZTheorException

Bases: Exception

gns.nested_sampler.Eoft2OverEoftPowi(nLive, i)

gns.nested_sampler.Eoft2Powi(nLive, i)

gns.nested_sampler.EoftPowi(nLive, i)

gns.nested_sampler.LLhood(LhoodObj)

Returns .logpdf method of LhoodObj

gns.nested_sampler.Lhood(LhoodObj)

Returns .pdf method of LhoodObj

gns.nested_sampler.NestedRun(priorParams, LLhoodParams, paramNames, setupDict)

function which completes a NS run. parameters of priors and likelihood need to be specified, as well as a flag indication type of prior for each dimension and the pdf for the lhood. setupDict contains other setup parameters such as termination type & factor, method of finding new livepoint, details of how weights are calculated, how final Z contribution is added, and directory/file prefix for saved files.

gns.nested_sampler.NestedRunLinear(priorParams, LhoodParams, paramNames, setupDict)

function which completes a NS run. parameters of priors and likelihood need to be specified, as well as a flag indication type of prior for each dimension and the pdf for the lhood. setupDict contains other setup parameters such as termination type & factor, method of finding new livepoint, details of how weights are calculated, how final Z contribution is added, and directory/file prefix for saved files.

gns.nested_sampler.applyBoundary(point, lower, upper, differenceCorrection, pointCorrection)

give a point in or outside the domain, in case of point being in the domain it does nothing. When it is outside, it calculates a ‘distance’ from the domain according to differenceCorrection type, and then uses this ‘distance’ to transform the point into the domain by using either reflective or wrapping methods according to the value of pointCorrection. It is recommended that differenceCorrection and pointCorrection take the same values (makes most intuitive sense to me)

gns.nested_sampler.calc1mt(nLive, expectation='1mt', sampling=False, maxPoints=False)

calc value of 1-t from its pdf, from (supposedely equivalent) way of deriving form of pdf, or from E[.] or E[l(.)]

gns.nested_sampler.calc1mt2(nLive, expectation='1mt2', sampling=False, maxPoints=False)

calc value of (1-t)^2 from its pdf, from (supposedely equivalent) way of deriving form of pdf, or from E[.] or E[l(.)]

gns.nested_sampler.calcEofLogZ(EofZ, varZ)

Calc E[log(Z)] from E[Z] and Var[Z]. Assumes Z is log-normally distributed

gns.nested_sampler.calcEofZ(EofLogZ, varLogZ)

calc E[Z] from E[logZ] and var[logZ]. Assumes Z is log-normal

gns.nested_sampler.calcEofZ2FinalKeeton(finalLhoods, nLive, nest)

Averages over Lhood, which I don’t think is the correct thing to do as it doesn’t correspond to a unique parameter vector value. TODO: consider other ways of getting final contribution from livepoints with Keeton’s method

gns.nested_sampler.calcEofZ2FinalKeetonLog(finalLLhoods, nLive, nest)

TODO

gns.nested_sampler.calcEofZ2Keeton(Lhoods, nLive, nest)

Calculate second (raw) moment of Z from main NS loop (equation 22 Keeton)

gns.nested_sampler.calcEofZ2KeetonLog(LLhoods, nLive, nest)

TODO

gns.nested_sampler.calcEofZFinalKeeton(finalLhoods, nLive, nest)

Averages over Lhood, which I don’t think is the correct thing to do as it doesn’t correspond to a unique parameter vector value. TODO: consider other ways of getting final contribution from livepoints with Keeton’s method

gns.nested_sampler.calcEofZFinalKeetonLog(finalLLhoods, nLive, nest)

TODO

gns.nested_sampler.calcEofZKeeton(Lhoods, nLive, nest)

Calculate first moment of Z from main NS loop. According to paper, this is just E[Z] = 1. / nLive * sum_i^nest L_i * E[t]^i

gns.nested_sampler.calcEofZKeetonLog(LLhoods, nLive, nest)

TODO

gns.nested_sampler.calcEofZZFinalKeeton(Lhoods, finalLhoods, nLive, nest)

Averages over Lhood for contribution from final points, which I don’t think is the correct thing to do as it doesn’t correspond to a unique parameter vector value. TODO: consider other ways of getting final contribution from livepoints with Keeton’s method

gns.nested_sampler.calcEofZZFinalKeetonLog(deadPointsLLhood, livePointsLLhood, nLive, nest)

TODO

gns.nested_sampler.calcEofftArr(Eofft, nLive, n)

Calculates E[f(t)]^i then returns this with yield. Yield means next time function is called, it picks off from where it last returned, with same variable values as before returning. Note the function isn’t executed until the generator return by yield is iterated over Putting for loop here is faster than filling in blank array

gns.nested_sampler.calcEofts(nLive)

calculate expected values of t related variables to update Z and X moments

gns.nested_sampler.calcHErr(Z, ZErr, LhoodPiLogLIntegral, LhoodPiLogLErr)

Calculates error on H due to uncertainty of Z, HIntegrand and ln(Z). ignores possible correlation between Z and LhoodPiLogLIntegral

gns.nested_sampler.calcHKeeton(EofZ, Lhoods, nLive, nest)

Calculate H from KL divergence equation transformed to LX space as given in Keeton.

gns.nested_sampler.calcHKeetonLog(logEofZK, deadPointsLLhood, nLive, nest)

TODO

gns.nested_sampler.calcHSkillingK(EofZ, varZ, nLive)

Uses definition of error given in Skilling’s NS paper, ACCORDING to Keeton. Only valid in limit that Skilling’s approximation of var[log(Z)] = H / nLive being correct, and E[Z]^2 >> var[Z] so that log(1+x)~x approximation is valid. Also requires that Z is log-normally distributed I think this is only valid for NS loop contributions, not final part or total

gns.nested_sampler.calcHTheor(priorParams, priorFuncsPdf, LLhoodFunc, nDims, Z, ZErr, LhoodFunc=None)

Calculates HTheor from the KL divergence equation: H = int[P(theta) * ln(P(theta) / pi(theta))] = 1/Z * int[L(theta) * pi(theta) * ln(L(theta))] - ln(Z). For uniform priors, calculates volume and skips that part of integral (over pi(theta)). Uses same trick as ZTheor in that it composes a dictionary of functions to integrate, mapped to dimension(s) of theta vector to integrate along for given function. Passes this dictionary to function which nquad actually evaluates.

gns.nested_sampler.calcHTheorApprox(Z, nDims, priorVolume)

Only valid in limit that prior is hyperrectangle, and majority of lhood is contained in prior hypervolume such that limits of integration can be extended to +- infinity

gns.nested_sampler.calcHTotalKeeton(EofZ, Lhoods, nLive, nest, finalLhoods)

Calculates total value of H based on KL divergence equation transformed to LX space as given in Keeton. Uses H function used to calculate loop H value (but with total Z), and adapts final result to give HTotal

gns.nested_sampler.calcHTotalKeetonLog(logEofZFinalK, deadPointsLLhood, nLive, nest, livePointsLLhood)

TODO

gns.nested_sampler.calcInitTrialSig(livePoints)

calculate initial standard deviation based on width of domain defined by max and min parameter values of livepoints in each dimension

gns.nested_sampler.calcInnerSums(Lhoods, nLive, nest)

Second generator (yielding) function, which returns inner sum for outer index k

gns.nested_sampler.calcLogEofts(nLive)

Calculate log(E[t] - E[t^2]) as it is much easier to do so here than later having on log(E[t]) and log(E[t^2])

gns.nested_sampler.calcSums(Lhoods, nLive, nest)

Calculate double summation in equation (22) of Keeton using two generator (yielding) functions. First one creates array associated with index of inner sum (which is subsequently summed). Second one creates array of summed inner sums, which is then multiplied by array of L_k * E[t]^k terms to give outer summation terms. Outer summation terms are added together to give total of double sum.

gns.nested_sampler.calcSumsLoop(Lhoods, nLive, nest)

Calculate double summation in equation (22) of Keeton using double for loop (one for each summation). Inefficient (I think) but easy

gns.nested_sampler.calcVarLogZ(EofZ, varZ, method)

Uses propagation of uncertainty formula or relationship between log-normal r.v.s and the normally distributed log of the log-normal r.v.s to calculate var[logZ] from EofZ and varZ (taken from Wikipedia)

gns.nested_sampler.calcVarLogZSkilling(H, nLive)

Skilling works in log space throughout, including calculating the moments of log(*) i.e. E[f(log(*))]. Thus he derives a value for the variance of log(Z), through his discussions of Poisson fluctuations whilst exploring the posterior.

gns.nested_sampler.calcVarZ(varLogZ, method, EofZ=None, EofLogZ=None)

Uses propagation of uncertainty formula or relationship between log-normal r.v.s and the normally distributed log of the log-normal r.v.s to calculate var[Z] from EofZ and varLogZ (taken from Wikipedia)

gns.nested_sampler.calcVarZSkillingK(EofZ, nLive, H)

Uses definition of error given in Skilling’s NS paper, ACCORDING to Keeton. Only valid in limit that Skilling’s approximation of var[log(Z)] = H / nLive being correct, and E[Z]^2 >> var[Z] so that log(1+x)~x approximation is valid. Also requires that Z is log-normally distributed I think this is only valid for NS loop contributions, not final part or total

gns.nested_sampler.calcVariance(EofX, EofX2)

Calculate second moment of X from first moment and raw second moment

gns.nested_sampler.calcVarianceLog(logEofX, logEofX2)

Calc log(var(X)) from log(E[X]) and log(E[X^2]) Does logsubtractexp manually, so doesn’t account for possible underflow issues with exponentiating like np.logaddexp() does, but this shouldn’t be an issue for the numbers involved here.

gns.nested_sampler.calcZMomentsFinalKeeton(finalLhoods, nLive, nest)

calculate Z moments a-posteri with list of final Lhood points (ones remaining at termination of main loop), using equations given in Keeton

gns.nested_sampler.calcZMomentsFinalKeetonLog(livePointsLLhood, nLive, nest)

TODO

gns.nested_sampler.calcZMomentsKeeton(Lhoods, nLive, nest)

calculate Z moments a-posteri with full list of Lhoods used in NS loop, using equations given in Keeton

gns.nested_sampler.calcZMomentsKeetonLog(deadPointsLLhood, nLive, nest)

TODO

gns.nested_sampler.calcZTheor(priorParams, priorFuncsPdf, LhoodFunc, nDims)

numerically integrates L(theta) * pi(theta) over theta priorFuncs must be in same order as dimensions of LhoodFunc when it was fitted (and in same order as priorParams). Have to temporarily set np.seterr to warnings only as underflow always seems to occur when evaluating integral (make sense I guess). A bit slow, but not sure how I can make it faster tbh, as it will get exponentially slower with # of dimensions LhoodFunc and priorFuncsPdf can be .pdf() or .logpdf() methods

gns.nested_sampler.calcZTheorApprox(priorParams)

Only valid in limit that prior is hyperrectangle, and majority of lhood is contained in prior hypervolume such that limits of integration (domain of the sampling space defined by the prior) can be extended close enough +- infinity such that the Lhood integrates to 1 over this domain

gns.nested_sampler.calct(nLive, expectation='t', sampling=False, maxPoints=False)

calc value of t from its pdf, from (supposedely equivalent) way of deriving form of pdf, or from E[.] or E[l(.)]

gns.nested_sampler.calct2(nLive, expectation='t2', sampling=False, maxPoints=False)

calc value of t^2 from its pdf, from (supposedely equivalent) way of deriving form of pdf, or from E[.] or E[l(.)]

gns.nested_sampler.callGetDist(chainsFilePrefix, plotName, nParams)

produces triangular posterior plots using getDist for first nParams parameters from chains file as labelled in that file and in .paramnames

gns.nested_sampler.checkBoundaries(livePoint, priorParams)

For all parameters with fixed boundaries, (just uniform for NOW), ensures trial point in that dimension has value in allowed domain

gns.nested_sampler.checkIfAveragedLhood(nFinal, livePointsLhood, avLhood)

Checks if Lhood was averaged over in getLogWeightsLive or not. If it was, need to work with average LLhood value for remainder of calculations, if not then carry on working with nLive size array of LLhoods. For ZLiveType == ‘max’, average is just taken over array size one with max Lhood value in it, so it is still just max value

gns.nested_sampler.checkInputParamsShape(priorParams, LhoodParams, nDims)

checks prior and Lhood input arrays are correct shape

gns.nested_sampler.evalExpLogIntegrand(*args)

evaluates .logpdf() of Lhood/ prior functions, adds them together along with an arbitrary ‘large’ value given in paramsNIntLogFuncsNIntLogVal[-1] then exponentiates this value to avoid underflow when evaluating the pdfs/ multiplying them together. This should avoid any underflow at all, provided the given value of integrandLogVal is large enough There is some underflow you cannot avoid, regardless of the value of integrandLogVal. This is because the LLhood values range from e.g. O(-10) to ~ O(-10^3). Subtracting approximation of max value of integrand is also not helpful for the same reason (spans too many orders of magnitude) Using an integrandLogVal suitable for the former will still cause the latter to underflow, but using a value to prevent the latter from underflow will result in the former overflowing! Dynamically calculating integrandLogVal wouldn’t help, as you would need to include this number again after exponentiating (which would cause it to underflow), but before the value is added to integral (if factor isn’t constant, doesn’t commute with integral) n.b. underflow is usually better than overflow, as you don’t want to miss your most likely values args consists of 1) the vector of theta values for given call (from nquad) 2) list of functions which make up the integrand 3) ‘large’ number to be added to logarithm before exponentiated to give integrand

gns.nested_sampler.evalIntegrand(*args)

evaluates a-priori parameter fitted pdfs with given data. data has to be reshaped because scipy functions are annoying and require last axis to be number of dimensions. Last element of parametersAndIntegrandFuncs is list of pdf objects for Lhood and priors, which are evaluated and multiplied together to give the value of the integrand. Note in general, Lhood will have dimensionality nDims whereas each prior will have dimensionality 1. The value (of the key corresponding to the function) is the relevant slice of x to get the correct dimensions for each function. May suffer from underflow either when evaluating the .pdf() calls, or when multiplying them together. args consists of 1) the vector of theta values for given call (from nquad) 2) list of functions which make up the integrand

gns.nested_sampler.evalLogLExpLogIntegrand(*args)

Same as above evalExpLogIntegrand but multiplies by log(L) for calculating H using exp(log(integrand)) method. LLhoodFunc has to be passed in separately as well as in the dictionary of functions, so it can be used to evaluate LLhoodFunc(theta)*np.exp(integrandLogVal) args consists of 1) the vector of theta values for given call (from nquad) 2) list of functions which make up the integrand 3) ‘large’ number to be added to logarithm before exponentiated to give integrand 4) LLhood function required for log(L(theta)) part of integrand

gns.nested_sampler.fitLhood(LLhoodParams)

fit lhood (without data) for parameters to make future evaluations much faster.

gns.nested_sampler.fitPriors(priorParams)

Only currently handles one dimensional (independent priors). Scipy.stats multivariate functions do not have built in inverse CDF methods, so if I want to consider multivariate priors I may have to write my own code. Note scipy.stats.uniform takes parameters loc and scale where the boundaries are defined to be loc and loc + scale, so scale = upper - lower bound. Returns list of fitted prior objects length of nDims (one function for each parameter)

gns.nested_sampler.getEofWeightsLive(Eofw, LhoodStar, liveLhoods, trapezoidalFlag, ZLiveType)

as above but non-log space version

gns.nested_sampler.getEofZ2TotalKeeton(EofZ2, EofZ2Final)

get total from NS loop and final contributions

gns.nested_sampler.getEofZ2TotalKeetonLog(logEofZ2, logEofZ2Final)

TODO

gns.nested_sampler.getEofZLive(nLive, EofX, livePointsLhood, LhoodStar, ZLiveType, trapezoidalFlag)

as above but in linear space

gns.nested_sampler.getEofZTotalKeeton(EofZ, EofZFinal)

get total from NS loop and final contributions

gns.nested_sampler.getEofZTotalKeetonLog(logEofZ, logEofZFinal)

TODO

gns.nested_sampler.getEofftArr(Eofft, nLive, nest)

faster than creating array of zeroes and looping over

gns.nested_sampler.getEofwLive(nLive, EofX, ZLiveType)

as above but in non-log space

gns.nested_sampler.getFinalAverage(livePointsPhys, livePointsLLhood, X, nFinal, avLLhood, space)

gets final livepoint values and X value per remaining livepoint for average Z criteria NOTE Xfinal is a list not a numpy array space says whether you are working in linear or log space (X or logX)

gns.nested_sampler.getFinalContribution(verbose, ZLiveType, trapezoidalFlag, nFinal, EofZ, EofZ2, EofX, EofWeights, H, livePointsPhys, livePointsLhood, avLhood, liveLhoodMax, liveMaxIndex, LhoodStar, errorEval='recursive')

as above but in linear space

gns.nested_sampler.getFinalContributionLog(verbose, ZLiveType, trapezoidalFlag, nFinal, logEofZ, logEofZ2, logEofX, logEofWeights, H, livePointsPhys, livePointsLLhood, avLLhood, liveLLhoodMax, liveMaxIndex, LLhoodStar, errorEval='recursive')

Get final contribution from livepoints after NS loop has ended. Way of estimating final contribution is dictated by ZLiveType. Also updates H value and gets final weights (and physical values) for posterior this function could be quite taxing on memory as it has to copy all arrays/ lists across NOTE: for standard quadrature summation, average Lhood and average X give same values of Z (averaging over X is equivalent to averaging over L). However, correct posterior weights are given by latter method, and Z errors are different in both cases

gns.nested_sampler.getFinalMax(liveMaxIndex, livePointsPhys, liveLhoodMax, X)

get livepoint and physical livepoint values corresponding to maximum likelihood point in remaining livepoints. Note Xfinal is a list not a numpy array or a scalar Function works for log or linear space

gns.nested_sampler.getLivePointsPhysFinal(livePointsPhys, avLhood)

Get physical values associated with remaining contribution of livepoints. If LLhood isn’t averaged over (X is) this is just the input livepoint values, but if LLhood is averaged it is non-trivial, I.E. THE PHYSICAL VALUES ASSOCIATED WITH THIS POINT ARE MEANININGLESS Only relevant for ZLiveType == ‘average’ as for ‘max’ case, physical values are just that corresponding to max(L) These are needed for posterior samples of remaining contribution of livepoints

gns.nested_sampler.getLogEofWeight(logw, LLhood_im1, LLhood_i, trapezoidalFlag)

calculates logw + log(f(L_im1, L_i)) where f(L_im1, L_i) = L_i for standard quadrature and f(L_im1, L_i) = (L_im1 + L_i) / 2. for the trapezium rule

gns.nested_sampler.getLogEofWeightsLive(logEofw, LLhoodStar, liveLLhoods, trapezoidalFlag, ZLiveType)

From Will’s implementation, Z = sum (X_im1 - X_i) * 0.5 * (L_i + L_im1) Unsure whether you should treat final contribution using trapezium rule (when it is used for rest of sum). I think you should and in case of ZLiveType == ‘average ‘, the L values used are L + {L_live} and in the case of ZLiveType == ‘max’, the L values used are L* + {max(L_live)}. When trapezium rule isn’t used (for rest of sum), L values used are {L_live} in case of ZLiveType == ‘average ‘ and {max(L_live)} in case of ZLiveType == ‘max’. When ZLiveType == ‘average *’ there is an added complication of what the average is ‘taken over’ (for both trapezium rule and standard quadrature) i.e. over the prior volume or the likelihood. If ZLiveType == ‘average X’ the average is taken over X, meaning there are still nLive live log weights (equally spaced in X with values X / nLive) which for standard quadrature have values: {log(X / nLive) + log(L_1), …, log(X / nLive) + log(L_nLive)} and for trapezium rule: {log(X / nLive) + log((L + L_1) / 2. ), …, log(X / nLive) + log((L_nLive-1 + L_nLive) / 2. )} If ZLiveType == ‘average Lhood’ the average is taken over the remaining L values, meaning there is 1 live log weight with X value X (i.e. the L_average value is assumed to be at X = 0). For the standard quadrature method the live log weight thus has a value log(X) + log(sum_i^nLive[L_i] / nLive) and for the trapezoidal rule log(X) + log((L* + sum_i^nLive[L_i] / nLive) / 2.). When ZLiveType == ‘max’, the maximum is obviously taken over the remaining Lhoods. Thus there is only one live log weight. For standard quadrature this is log(X) + log(max(L_i) and for the trapezium rule it is log(X) + log((L* + max(L_i)) / 2.) If averaging over L, final livepoint needs to be attributed this L, so it is stored here under the variable avLLhood

gns.nested_sampler.getLogEofXLogEofw(nLive, X)

get increment (part of weight for posterior and evidence calculations) based on previous value of X, calculates latest X using t calculated from either expected value or sampling. Expected value can be of t (E[t]) or log(t) E[log(t)]. These are roughly the same for large nLive Sampling can take two forms: sampling from the pdf or taking the highest of U[0,1]^Nlive values (from which the pdf form is derived from), so they should in theory be the same.

gns.nested_sampler.getLogEofZLive(nLive, logEofX, livePointsLLhood, LLhoodStar, ZLiveType, trapezoidalFlag)

NOTE logWeightsLive here is an np array newLiveLLhoods has same shape as logWeightsLive (i.e. account for averageLhoodOrX value). If ZLiveType == ‘max’ avLLhood will just be the maximum LLhood value. there is no averaging to consider if ZLiveType == ‘max Lhood’. Could return live weights, but these need to be calculated again in final contribution function so don’t bother

gns.nested_sampler.getLogEofwLive(nLive, logEofX, ZLiveType)

Determines final logw based on ZLiveType and averageLhoodOrX, i.e. it determines whether final contribution is averaged/ maximised over L or averaged over X.

gns.nested_sampler.getMaxLhood(ZLiveType, livePointsLhood)

For ZLiveType == ‘max’ returns a 1 element array with maximum LLhood value, its value as a scalar, and the index of the max LLhood in the given array. For ZLiveType == ‘average’ it essentially does nothing

gns.nested_sampler.getNewLiveBlind(priorFuncsPpf, LhoodFunc, LhoodStar)

Blindly picks points isin U[0, 1]^D, converts these to physical values according to physical prior and uses as candidates for new livepoint until L > L* is found LhoodFunc can be Lhood or LLhood func, either works fine

gns.nested_sampler.getNewLiveMH(livePointsPhys, deadIndex, priorFuncsPdf, priorParams, LhoodFunc, LhoodStar)

gets new livepoint using variant of MCMC MH algorithm. From current livepoints (not including one to be excluded in current NS iteration), first picks a point at random as starting point. Next the standard deviation of the trial distribution is calculated from the width of the current livepoints (including the one to be excluded) as 0.1 * [max(param_value) - min(param_value)] in each dimension. A trial distribution (CURRENTLY GAUSSIAN) is centred on the selected livepoint with the calculated variance, and the output is used as a trial point. This point is kept with probability prior(trial) / prior(previous) if L_trial > L* and rejected otherwise. Each time a trial point is proposed, nTrials is incremented. If the trial point was accepted nAccept is incremented, if not nReject is. The trial distribution is updated based on nAccept, nReject such that the acceptance rate should be roughyl 50%. Once nTrials nTrials = maxTrials, if nAccept = 0 the whole process is repeated as failure to do so will mean the returned live point is a copy of another point. If not, the new livepoint is returned. One would hope that nAccept is > 1 to ensure that the sampling space is explored uniformly. LhoodFunc can be Lhood or LLhood func, either works fine

gns.nested_sampler.getPriorIntegrandAndLimits(priorParams, priorFuncsPdf, integrandFuncs)

Adds prior functions to integrandFuncs dict for ‘non-rectangular’ (uniform) dimensions, as long as a mapping to dimensions of data required to integrate over for that function. For uniform priors, calculates the hyperrectangular volume. Also creates list of limits in each dimension of integral. For Gaussian priors, sets limits to +- infinity

gns.nested_sampler.getPriorLogPdfs(priorObjs)

Takes list of fitted prior objects, returns list of objects’ .logpdf() methods

gns.nested_sampler.getPriorPdfs(priorObjs)

Takes list of fitted prior objects, returns list of objects’ .pdf() methods

gns.nested_sampler.getPriorPpfs(priorObjs)

Takes list of fitted prior objects, returns list of objects’ .ppf() methods

gns.nested_sampler.getTotal(deadPointsPhys, livePointsPhysFinal, deadPointsLhood, livePointsLhoodFinal, XArr, XFinalArr, weights)

gets final arrays of physical, llhood and X values for all accepted points in algorithm. This function mutates deadPointsPhys by appending numpy array livePointsPhysFinal. This is at the end of the program however, so it shouldn’t be an issue. Concatenate works on a list of numpy arrays (those corresponding to deadPoints should have shape (1, nDims) and there should be nest of them, the single numpy arrays corresponding to the final live points should have shape (nLive, nDims) if average of Z was used for final contribution or (1, nDims) if max of Z was used. Concatenating list of numpy arrays is much more efficient than using np.append() at each iteration.

gns.nested_sampler.getVarTotalKeeton(varZ, varZFinal, EofZ, EofZFinal, EofZZFinal)

Get total variance from NS loop and final contributions. For recursive method, since E[ZLive] = E[ZTot] etc., and assuming that the recurrence relations account for the covariance between Z and ZFinal, this is just varZFinal. For Keeton’s method, have to explicitly account for correlation as expectations for Z and ZLive are essentially calculated independently TODO: check if recurrence relations of Z and ZFinal properly account for correlation between two

gns.nested_sampler.getVarTotalKeetonLog(logVarZ, logVarZFinal, logEofZ, logEofZFinal, logEofZZFinal)

TODO

gns.nested_sampler.invPrior(livePoints, priorFuncsPpf)

take in array of livepoints each which has value isin[0,1], and has nDim dimensions. Output physical array of values corresponding to priors for each parameter dimension.

gns.nested_sampler.logAddArr(x, y, axis=None)

logaddexp where x is a scalar and y is an array. Returns a scalar in case that axis = None (exponentiates elements of array then adds them together). np implementation on its own returns an array i.e. doesn’t sum exponentiated elements of array. Axis specifies axis to do summation of exponentials over, and if specified returns an array with shape of the remaining dimensions.

gns.nested_sampler.logAddArr2(x, y, indexes=(None, ))

Alternative version of logAddArr that avoids over/ underflow errors of exponentiating the array y, to the same extent that np.logaddexp() does Note however that it is slower than logAddArr, so in cases where over/ underflow isn’t an issue, use that Loops over each specified element of y (using rowcol values) and adds to log of sums. By default loops over entire array, but to specify a certain row for e.g. a 2d array set indexes to (row_index, slice(None)) or for a certain column (slice(None), col_index)

gns.nested_sampler.logsubexp(a, b)

Only currently works for scalars. Calculates log(exp(a) - exp(b)). Subtracts max(a,b) from a and b before exponentiating in attempt to avoid underflow when a and b are small

gns.nested_sampler.main()

gns.nested_sampler.modToDomainReflect(point, lower, upper)

following ensures point isin [lower, upper]. This reflects according to (positive) difference between point and nearest bound and returns a ‘distance’ which can actually be used to get the correct value of the point within the domain. It makes most intuitive sense to me to use this when you want to reflect the points in the domain. Operation done to ensure reflecting is different based on whether the difference between the point and the nearest part of the domain is an odd or even (incl. 0) multiple of the boundary.

gns.nested_sampler.modToDomainWrap(point, lower, upper)

gns.nested_sampler.nDIntegratorHTheor(integrandFuncs, limitsList, integrandLogVal=200.0, LLhoodFunc=None)

As above but has to call nquad with a slightly different function representing integrand when using exp(log(integrand)) method for calculating, due to LLhood(theta) part of integrand

gns.nested_sampler.nDIntegratorZTheor(integrandFuncs, limitsList, integrandLogVal=200.0)

integrator used for calculating theoretical value of Z. integrand is function which evaluates to value of integrand at given parameter values (which are determined by nquad function). integrandFuncs is list of functions (Lhood & non-rectangular priors) which are multiplied together to give value of integrand. If integrandLogVal evaluates to true, does integration method which takes exp of log of integrand * some large number given by integrandLogVal to avoid underflow. The final integral result (and error) is then divided by exp(integrandLogVal) to get the final value.

gns.nested_sampler.plotPhysPosteriorIW(x, unnormalisedSamples, Z, space)

Plots posterior in physical space according to importance weights w(theta)L(theta) / Z. Doesn’t use KDE so isn’t true shape of posterior. If inputting logWeights/ logZ then set space == ‘log’

gns.nested_sampler.plotXPosterior(X, L, Z, space)

Plots X*L(X)/Z in log X space, not including KDE methods

gns.nested_sampler.printBreak()

tell user final contribution to sampling is being calculated

gns.nested_sampler.printFinalLivePoints(i, physValue, Lhood, ZLiveType, space)

print information about final livepoints used to calculate final contribution to Z/ posterior samples.

gns.nested_sampler.printSampleNum(numSamples)

Print number of samples used in sampling (including final livepoints used for posterior weights)

gns.nested_sampler.printTerminationUpdateInfo(nest, terminator)

Print update on termination status when evaluating by H value

gns.nested_sampler.printTerminationUpdateZ(EofZLive, endValue, terminationFactor, space)

Print update on termination status when evaluating by Z ratio

gns.nested_sampler.printTheoretical(ZTheor, ZTheorErr, HTheor, HTheorErr)

Outputs values for theoretical values of Z and H (and their errors)

gns.nested_sampler.printUpdate(nest, deadPointPhys, deadPointLhood, EofZ, livePointPhys, livePointLhood, space)

gives update on latest deadpoint and newpoint found to replace it

gns.nested_sampler.printZHValues(EofZ, EofZ2, varZ, H, space, stage, method)

print values of Z (including varios moments, variance) and H in either log or linear space, at a given stage and calculated by a given method

gns.nested_sampler.priorFuncsProd(livePoint, priorFuncsPdf)

calculates pdf of prior for each parameter dimension, then multiplies these together to get the pdf of the prior (i.e. the prior pdf assuming the parameters are independent) Works in linear space, but can easily be adapted if this ever leads to underflow errors (consider sum of log(pi(theta))).

gns.nested_sampler.testTrial(trialPoint, startPoint, trialLhood, startLhood, LhoodStar, priorFuncsPdf)

Check if trial point has L > L* and accept with probability prior(trial) / prior(previous)

gns.nested_sampler.tryTermination(verbose, terminationType, terminationFactor, nest, nLive, EofX, livePointsLhood, LhoodStar, ZLiveType, trapezoidalFlag, EofZ, H)

as above but in linear space

gns.nested_sampler.tryTerminationLog(verbose, terminationType, terminationFactor, nest, nLive, logEofX, livePointsLLhood, LLhoodStar, ZLiveType, trapezoidalFlag, logEofZ, H)

See if termination condition for main loop of NS has been met. Can be related to information value H or whether estimated remaining evidence is below a given fraction of the Z value calculated up to that iteration

gns.nested_sampler.updateEofX(Eoft, EofX)

Update value of raw first moment of X

gns.nested_sampler.updateEofX2(Eoft2, EofX2)

Update value of raw 2nd momement of X

gns.nested_sampler.updateEofX2Final(EofX, nFinal)

can’t be proved mathematically, just derived from recurrence relations

gns.nested_sampler.updateEofXFinal(EofX, nFinal)

can’t be proved mathematically, X is treated deterministically to be X / nLive

gns.nested_sampler.updateEofZ(EofZ, Eof1mt, EofX, L)

Update mean estimate of Z.

gns.nested_sampler.updateEofZ2(EofZ2, Eof1mt, EofZX, Eof1mt2, EofX2, L)

Update value of raw 2nd moment of Z based on Lhood value obtained in that NS iteration

gns.nested_sampler.updateEofZ2Final(EofZ2, EofX, EofZ, EofX2, L)

TODO: rewrite docstring

gns.nested_sampler.updateEofZFinal(EofZ, EofX, L)

TODO: rewrite docstring

gns.nested_sampler.updateEofZX(Eoft, EofZX, Eoft2, EofX2, L)

Updates raw ‘mixed’ moment of Z and X. Required to calculate E(Z)^2.

gns.nested_sampler.updateH(H, weight, ZNew, Lhood, Z)

Same as Skilling’s implementation but in linear space Handles FloatingPointErrors associated with taking np.log(0) (0 * log(0) = 0)

gns.nested_sampler.updateHLog(H, logWeight, logZNew, LLhood, logZ)

update H using previous value, previous and new log(Z) and latest weight Isn’t a non-log version as H propto log(L). As given in Skilling’s paper TODO: consider if trapezium rule should lead to different implementation

gns.nested_sampler.updateLogEofX(logEoft, logEofX)

as above but for log space

gns.nested_sampler.updateLogEofX2(logEoft2, logEofX2)

as above but for log space

gns.nested_sampler.updateLogEofX2Final(logEofX, nFinal)

can’t be proved mathematically, just derived from recurrence relations

gns.nested_sampler.updateLogEofXFinal(logEofX, nFinal)

can’t be proved mathematically, just derived from recurrence relations

gns.nested_sampler.updateLogEofZ(logEofZ, logEof1mt, logEofX, LL)

as above but for log space

gns.nested_sampler.updateLogEofZ2(logEofZ2, logEof1mt, logEofZX, logEof1mt2, logEofX2, LL)

as above but for log space

gns.nested_sampler.updateLogEofZ2Final(logEofZ2, logEofX, logEofZ, logEofX2, LL)

TODO: rewrite docstring

gns.nested_sampler.updateLogEofZFinal(logEofZ, logEofX, LL)

TODO: rewrite docstring

gns.nested_sampler.updateLogEofZX(logEoft, logEofZX, logEoftmEoft2, logEofX2, LL)

as above but for log space

gns.nested_sampler.updateLogZnXM(nLive, logEofZ, logEofZ2, logEofZX, logEofX, logEofX2, LL)

as above but for log space

gns.nested_sampler.updateLogZnXMoments(nLive, logEofZ, logEofZ2, logEofZX, logEofX, logEofX2, LLhoodStarOld, LLhoodStar, trapezoidalFlag)

as above but for log space

gns.nested_sampler.updateLogZnXMomentsF(nFinal, logEofZ, logEofZ2, logEofX, LL, errorEval)

TODO: rewrite docstring

gns.nested_sampler.updateLogZnXMomentsFinal(nFinal, logEofZ, logEofZ2, logEofX, LLhood_im1, LLhood_i, trapezoidalFlag, errorEval)

Wrapper around updateZnXMomentsF taking into account whether trapezium rule is used or not

gns.nested_sampler.updateTrialSigma(trialSigma, nAccept, nReject)

update standard deviation as in Sivia 2006. Apparently this ensures that ~50% of the points are accepted, but I’m sceptical.

gns.nested_sampler.updateZnXM(nLive, EofZ, EofZ2, EofZX, EofX, EofX2, L)

Update moments of Z and X based on their previous values, expected value of random variable t and Lhood value ((L_i + L_i-1) / 2. in case of trapezium rule). Used to calculate the mean and standard deviation of Z, and thus of log(Z) as well TODO: CONSIDER KEETON NON-RECURSIVE METHOD

gns.nested_sampler.updateZnXMoments(nLive, EofZ, EofZ2, EofZX, EofX, EofX2, LhoodStarOld, LhoodStar, trapezoidalFlag)

Wrapper around updateZnXM taking into account whether trapezium rule is used or not

gns.nested_sampler.updateZnXMomentsF(nFinal, EofZ, EofZ2, EofX, L, errorEval)

TODO: rewrite docstring TODO: CONSIDER KEETON NON-RECURSIVE METHOD WHICH EXPLICITLY ACCOUNTS FOR CORRELATION BETWEEN EOFZ AND EOFZLIVE

gns.nested_sampler.updateZnXMomentsFinal(nFinal, EofZ, EofZ2, EofX, Lhood_im1, Lhood_i, trapezoidalFlag, errorEval)

Wrapper around updateZnXMomentsF taking into account whether trapezium rule is used or not

gns.nested_sampler.writeOutput(outputFile, totalPointsPhys, totalPointsLhood, weights, XArr, paramNames, space)

writes a summary file which contains values for all sampled points. Also writes files needed for getDist.