# Licensed under a 3-clause BSD style license - see LICENSE.rstimportnumpyasnpimportscipy.optimizefromgammapy.utils.interpolationimportinterpolate_profilefromgammapy.utils.rootsimportfind_rootsfrom.likelihoodimportLikelihood__all__=["confidence_scipy","covariance_scipy","optimize_scipy","stat_profile_ul_scipy",]defoptimize_scipy(parameters,function,store_trace=False,**kwargs):method=kwargs.pop("method","Nelder-Mead")pars=[par.factorforparinparameters.free_parameters]bounds=[]forparinparameters.free_parameters:parmin=par.factor_minifnotnp.isnan(par.factor_min)elseNoneparmax=par.factor_maxifnotnp.isnan(par.factor_max)elseNonebounds.append((parmin,parmax))likelihood=Likelihood(function,parameters,store_trace)result=scipy.optimize.minimize(likelihood.fcn,pars,bounds=bounds,method=method,**kwargs)factors=result.xinfo={"success":result.success,"message":result.message,"nfev":result.nfev,"trace":likelihood.trace,}optimizer=Nonereturnfactors,info,optimizerclassTSDifference:"""Fit statistic function wrapper to compute TS differences."""def__init__(self,function,parameters,parameter,reoptimize,ts_diff):self.stat_null=function()self.parameters=parametersself.function=functionself.parameter=parameterself.ts_diff=ts_diffself.reoptimize=reoptimizedeffcn(self,factor):self.parameter.factor=factorifself.reoptimize:self.parameter.frozen=Trueoptimize_scipy(self.parameters,self.function,method="L-BFGS-B")value=self.function()-self.stat_null-self.ts_diffreturnvaluedef_confidence_scipy_brentq(parameters,parameter,function,sigma,reoptimize,upper=True,**kwargs):ts_diff=TSDifference(function,parameters,parameter,reoptimize,ts_diff=sigma**2)lower_bound=parameter.factorifupper:upper_bound=parameter.conf_max/parameter.scaleelse:upper_bound=parameter.conf_min/parameter.scalemessage,success="Confidence terminated successfully.",Truekwargs.setdefault("nbin",1)roots,res=find_roots(ts_diff.fcn,lower_bound=lower_bound,upper_bound=upper_bound,**kwargs)result=(roots[0],res[0])ifnp.isnan(roots[0]):message=("Confidence estimation failed. Try to set the parameter.min/max by hand.")success=Falsesuffix="errp"ifupperelse"errn"return{"nfev_"+suffix:result[1].iterations,suffix:np.abs(result[0]-lower_bound),"success_"+suffix:success,"message_"+suffix:message,"stat_null":ts_diff.stat_null,}defconfidence_scipy(parameters,parameter,function,sigma,reoptimize=True,**kwargs):iflen(parameters.free_parameters)<=1:reoptimize=Falsewithparameters.restore_status():result=_confidence_scipy_brentq(parameters=parameters,parameter=parameter,function=function,sigma=sigma,upper=False,reoptimize=reoptimize,**kwargs,)withparameters.restore_status():result_errp=_confidence_scipy_brentq(parameters=parameters,parameter=parameter,function=function,sigma=sigma,upper=True,reoptimize=reoptimize,**kwargs,)result.update(result_errp)returnresult# TODO: implement, e.g. with numdifftools.Hessiandefcovariance_scipy(parameters,function):raiseNotImplementedError
[docs]defstat_profile_ul_scipy(value_scan,stat_scan,delta_ts=4,interp_scale="sqrt",**kwargs):"""Compute upper limit of a parameter from a likelihood profile. Parameters ---------- value_scan : `~numpy.ndarray` Array of parameter values. stat_scan : `~numpy.ndarray` Array of delta fit statistic values, with respect to the minimum. delta_ts : float, optional Difference in test statistics for the upper limit. Default is 4. interp_scale : {"sqrt", "lin"}, optional Interpolation scale applied to the fit statistic profile. If the profile is of parabolic shape, a "sqrt" scaling is recommended. In other cases or for fine sampled profiles a "lin" can also be used. Default is "sqrt". **kwargs : dict Keyword arguments passed to `~scipy.optimize.brentq`. Returns ------- ul : float Upper limit value. """interp=interpolate_profile(value_scan,stat_scan,interp_scale=interp_scale)deff(x):returninterp((x,))-delta_tsidx=np.argmin(stat_scan)norm_best_fit=value_scan[idx]roots,res=find_roots(f,lower_bound=norm_best_fit,upper_bound=value_scan[-1],nbin=1,**kwargs)returnroots[0]