Python Fit.est_errors方法代码示例

# 需要导入模块: from sherpa.fit import Fit [as 别名]
# 或者: from sherpa.fit.Fit import est_errors [as 别名]
    def fit(self, current):
        self._fit.model.thawedpars = current

        # nm = NelderMead()
        # nm.config['iquad'] = 0
        # nm.config['finalsimplex'] = 1

        #lm = LevMar()
        #lm.config['maxfev'] = 5

        cv = Covariance()

        # Use the fit method defined before called get_draws().  This way the user
        # does not have to pass in the fitting method and method options.
        fit = Fit(self._fit.data, self._fit.model, self._fit.stat, self._fit.method,
                  cv)

        fit_result = fit.fit()
        covar_result = fit.est_errors()
        sigma = np.array(covar_result.extra_output)

        if np.isnan(sigma).any():
            raise CovarError("NaNs found in covariance matrix")

        # cache the fitting scales
        self._sigma = sigma

# 需要导入模块: from sherpa.fit import Fit [as 别名]
# 或者: from sherpa.fit.Fit import est_errors [as 别名]
    def fit(self):
        """Fit spectrum"""
        from sherpa.fit import Fit
        from sherpa.models import ArithmeticModel, SimulFitModel
        from sherpa.astro.instrument import Response1D
        from sherpa.data import DataSimulFit

        # Translate model to sherpa model if necessary
        if isinstance(self.model, models.SpectralModel):
            model = self.model.to_sherpa()
        else:
            model = self.model

        if not isinstance(model, ArithmeticModel):
            raise ValueError('Model not understood: {}'.format(model))

        # Make model amplitude O(1e0)
        val = model.ampl.val * self.FLUX_FACTOR ** (-1)
        model.ampl = val

        if self.fit_range is not None:
            log.info('Restricting fit range to {}'.format(self.fit_range))
            fitmin = self.fit_range[0].to('keV').value
            fitmax = self.fit_range[1].to('keV').value

        # Loop over observations
        pha = list()
        folded_model = list()
        nobs = len(self.obs_list)
        for ii in range(nobs):
            temp = self.obs_list[ii].to_sherpa()
            if self.fit_range is not None:
                temp.notice(fitmin, fitmax)
                if temp.get_background() is not None:
                    temp.get_background().notice(fitmin, fitmax)
            temp.ignore_bad()
            if temp.get_background() is not None:
                temp.get_background().ignore_bad()
            pha.append(temp)
            # Forward folding
            resp = Response1D(pha[ii])
            folded_model.append(resp(model) * self.FLUX_FACTOR)

        data = DataSimulFit('simul fit data', pha)
        fitmodel = SimulFitModel('simul fit model', folded_model)

        log.debug(fitmodel)
        fit = Fit(data, fitmodel, self.statistic)
        fitresult = fit.fit()
        log.debug(fitresult)
        # The model instance passed to the Fit now holds the best fit values
        covar = fit.est_errors()
        log.debug(covar)

        for ii in range(nobs):
            efilter = pha[ii].get_filter()
            shmodel = fitmodel.parts[ii]
            self.result[ii].fit = _sherpa_to_fitresult(shmodel, covar, efilter, fitresult)

# 需要导入模块: from sherpa.fit import Fit [as 别名]
# 或者: from sherpa.fit.Fit import est_errors [as 别名]
class test_sim(SherpaTestCase):

    @requires_fits
    @requires_xspec
    def setUp(self):
        from sherpa.astro.io import read_pha
        from sherpa.astro.xspec import XSwabs, XSpowerlaw
        # self.startdir = os.getcwd()
        self.old_level = logger.getEffectiveLevel()
        logger.setLevel(logging.CRITICAL)

        pha = self.make_path("refake_0934_1_21_1e4.fak")
        # rmf = self.make_path("ccdid7_default.rmf")
        # arf = self.make_path("quiet_0934.arf")

        self.simarf = self.make_path("aref_sample.fits")
        self.pcaarf = self.make_path("aref_Cedge.fits")

        data = read_pha(pha)
        data.ignore(None, 0.3)
        data.ignore(7.0, None)

        rsp = Response1D(data)
        self.abs1 = XSwabs('abs1')
        self.p1 = XSpowerlaw('p1')
        model = rsp(self.abs1 * self.p1)

        self.fit = Fit(data, model, CStat(), NelderMead(), Covariance())

    def tearDown(self):
        # os.chdir(self.startdir)
        if hasattr(self, 'old_level'):
            logger.setLevel(self.old_level)

    @requires_xspec
    @requires_data
    def test_pragbayes_simarf(self):
        mcmc = sim.MCMC()

        self.abs1.nh = 0.092886
        self.p1.phoindex = 0.994544
        self.p1.norm = 9.26369

        mcmc.set_sampler("PragBayes")
        mcmc.set_sampler_opt("simarf", self.simarf)
        mcmc.set_sampler_opt("p_M", 0.5)
        mcmc.set_sampler_opt("nsubiter", 7)

        covar_results = self.fit.est_errors()
        cov = covar_results.extra_output

        niter = 10
        stats, accept, params = mcmc.get_draws(self.fit, cov, niter=niter)
        # try:
        #     assert (covar_results.parmaxes < params.std(1)).all()
        # except AssertionError:
        #     print 'covar: ', str(covar_results.parmaxes)
        #     print 'param: ', str(params.std(1))
        #     raise

    @requires_xspec
    @requires_data
    def test_pragbayes_pcaarf(self):
        mcmc = sim.MCMC()

        self.abs1.nh = 0.092886
        self.p1.phoindex = 0.994544
        self.p1.norm = 9.26369

        mcmc.set_sampler("pragBayes")
        mcmc.set_sampler_opt("simarf", self.pcaarf)
        mcmc.set_sampler_opt("p_M", 0.5)
        mcmc.set_sampler_opt("nsubiter", 5)

        covar_results = self.fit.est_errors()
        cov = covar_results.extra_output

        niter = 10
        stats, accept, params = mcmc.get_draws(self.fit, cov, niter=niter)

# 需要导入模块: from sherpa.fit import Fit [as 别名]
# 或者: from sherpa.fit.Fit import est_errors [as 别名]
class test_sim(SherpaTestCase):

    @unittest.skipIf(not has_fits_support(),
                     'need pycrates, pyfits')
    @unittest.skipIf(not has_package_from_list('sherpa.astro.xspec'),
                     "required sherpa.astro.xspec module missing")
    def setUp(self):
        try:
            from sherpa.astro.io import read_pha
            from sherpa.astro.xspec import XSwabs, XSpowerlaw
        except:
            return
        #self.startdir = os.getcwd()
        self.old_level = logger.getEffectiveLevel()
        logger.setLevel(logging.CRITICAL)

        datadir = SherpaTestCase.datadir
        if datadir is None:
            return

        pha = os.path.join(datadir, "refake_0934_1_21_1e4.fak")
        rmf = os.path.join(datadir, "ccdid7_default.rmf")
        arf = os.path.join(datadir, "quiet_0934.arf")
        
        self.simarf = os.path.join(datadir, "aref_sample.fits")
        self.pcaarf = os.path.join(datadir, "aref_Cedge.fits")

        data = read_pha(pha)
        data.ignore(None,0.3)
        data.ignore(7.0,None)

        rsp = Response1D(data)
        self.abs1 = XSwabs('abs1')
        self.p1 = XSpowerlaw('p1')
        model = rsp(self.abs1*self.p1)
        
        self.fit = Fit(data, model, CStat(), NelderMead(), Covariance())


    def tearDown(self):
        #os.chdir(self.startdir)
        if hasattr(self,'old_level'):
            logger.setLevel(self.old_level)

    @unittest.skipIf(not has_package_from_list('sherpa.astro.xspec'),
                     "required sherpa.astro.xspec module missing")
    @unittest.skipIf(test_data_missing(), "required test data missing")
    def test_pragbayes_simarf(self):
        datadir = SherpaTestCase.datadir
        if datadir is None:
            return

        mcmc = sim.MCMC()

        self.abs1.nh = 0.092886
        self.p1.phoindex = 0.994544
        self.p1.norm = 9.26369

        mcmc.set_sampler("PragBayes")
        mcmc.set_sampler_opt("simarf", self.simarf)
        mcmc.set_sampler_opt("p_M", 0.5)
        mcmc.set_sampler_opt("nsubiter", 7)

        covar_results = self.fit.est_errors()
        cov = covar_results.extra_output

        niter = 10
        stats, accept, params = mcmc.get_draws(self.fit, cov, niter=niter)
        # try:
        #     assert (covar_results.parmaxes < params.std(1)).all()
        # except AssertionError:
        #     print 'covar: ', str(covar_results.parmaxes)
        #     print 'param: ', str(params.std(1))
        #     raise

    @unittest.skipIf(not has_package_from_list('sherpa.astro.xspec'),
                     "required sherpa.astro.xspec module missing")
    @unittest.skipIf(test_data_missing(), "required test data missing")
    def test_pragbayes_pcaarf(self):
        datadir = SherpaTestCase.datadir
        if datadir is None:
            return

        mcmc = sim.MCMC()

        self.abs1.nh = 0.092886
        self.p1.phoindex = 0.994544
        self.p1.norm = 9.26369

        mcmc.set_sampler("pragBayes")
        mcmc.set_sampler_opt("simarf", self.pcaarf)
        mcmc.set_sampler_opt("p_M", 0.5)
        mcmc.set_sampler_opt("nsubiter", 5)

        covar_results = self.fit.est_errors()
        cov = covar_results.extra_output

        niter = 10
        stats, accept, params = mcmc.get_draws(self.fit, cov, niter=niter)

# 需要导入模块: from sherpa.fit import Fit [as 别名]
# 或者: from sherpa.fit.Fit import est_errors [as 别名]
    def fit(self):
        """Fit spectrum"""
        from sherpa.fit import Fit
        from sherpa.models import ArithmeticModel, SimulFitModel
        from sherpa.astro.instrument import Response1D
        from sherpa.data import DataSimulFit

        # Reset results
        self._result = list()

        # Translate model to sherpa model if necessary
        if isinstance(self.model, models.SpectralModel):
            model = self.model.to_sherpa()
        else:
            model = self.model

        if not isinstance(model, ArithmeticModel):
            raise ValueError('Model not understood: {}'.format(model))

        # Make model amplitude O(1e0)
        val = model.ampl.val * self.FLUX_FACTOR ** (-1)
        model.ampl = val

        if self.fit_range is not None:
            log.info('Restricting fit range to {}'.format(self.fit_range))
            fitmin = self.fit_range[0].to('keV').value
            fitmax = self.fit_range[1].to('keV').value

        # Loop over observations
        pha = list()
        folded_model = list()
        nobs = len(self.obs_list)
        for ii in range(nobs):
            temp = self.obs_list[ii].to_sherpa()
            if self.fit_range is not None:
                temp.notice(fitmin, fitmax)
                if temp.get_background() is not None:
                    temp.get_background().notice(fitmin, fitmax)
            temp.ignore_bad()
            if temp.get_background() is not None:
                temp.get_background().ignore_bad()
            pha.append(temp)
            log.debug('Noticed channels obs {}: {}'.format(
                ii, temp.get_noticed_channels()))
            # Forward folding
            resp = Response1D(pha[ii])
            folded_model.append(resp(model) * self.FLUX_FACTOR)

        if (len(pha) == 1 and len(pha[0].get_noticed_channels()) == 1):
            raise ValueError('You are trying to fit one observation in only '
                             'one bin, error estimation will fail')

        data = DataSimulFit('simul fit data', pha)
        log.debug(data)
        fitmodel = SimulFitModel('simul fit model', folded_model)
        log.debug(fitmodel)

        fit = Fit(data, fitmodel, self.statistic, method=self.method_fit)

        fitresult = fit.fit()
        log.debug(fitresult)
        # The model instance passed to the Fit now holds the best fit values
        covar = fit.est_errors()
        log.debug(covar)

        for ii in range(nobs):
            efilter = pha[ii].get_filter()
            # Skip observations not participating in the fit
            if efilter != '':
                shmodel = fitmodel.parts[ii]
                result = _sherpa_to_fitresult(shmodel, covar, efilter, fitresult)
                result.obs = self.obs_list[ii]
            else:
                result = None
            self._result.append(result)

        valid_result = np.nonzero(self.result)[0][0]
        global_result = copy.deepcopy(self.result[valid_result])
        global_result.npred = None
        global_result.obs = None
        all_fitranges = [_.fit_range for _ in self._result if _ is not None]
        fit_range_min = min([_[0] for _ in all_fitranges])
        fit_range_max = max([_[1] for _ in all_fitranges])
        global_result.fit_range = u.Quantity((fit_range_min, fit_range_max))
        self._global_result = global_result

# 需要导入模块: from sherpa.fit import Fit [as 别名]
# 或者: from sherpa.fit.Fit import est_errors [as 别名]

#.........这里部分代码省略.........
        # model.ref = 0.1
        # fit = Fit(data, model, Cash(), NelderMead())

        # NOTE: We cannot use the Levenbergr-Marquart optimizer in Sherpa
        # because it relies on the fvec return value of the fit statistic (we
        # return None). The computation of fvec is not straightforwad, not just
        # stats per bin. E.g. for a cash fit the sherpa stat computes it
        # according to cstat
        # see https://github.com/sherpa/sherpa/blob/master/sherpa/include/sherpa/stats.hh#L122

        self._sherpa_fit = Fit(data,
                               SherpaModel(self),
                               SherpaStat(self),
                               NelderMead())
        fitresult = self._sherpa_fit.fit()
        log.debug(fitresult)
        self._make_fit_result()

    def _make_fit_result(self):
        """Bundle fit results into `~gammapy.spectrum.SpectrumFitResult`.

        It is important to copy best fit values, because the error estimation
        will change the model parameters and statval again
        """
        from . import SpectrumFitResult

        model = self.model.copy()
        if self.background_model is not None:

            bkg_model = self.background_model.copy()
        else:
            bkg_model = None

        statname = self.stat

        results = []
        for idx, obs in enumerate(self.obs_list):
            fit_range = self.true_fit_range[idx]
            statval = np.sum(self.statval[idx])
            stat_per_bin = self.statval[idx]
            npred_src = copy.deepcopy(self.predicted_counts[idx][0])
            npred_bkg = copy.deepcopy(self.predicted_counts[idx][1])

            results.append(SpectrumFitResult(
                model=model,
                fit_range=fit_range,
                statname=statname,
                statval=statval,
                stat_per_bin=stat_per_bin,
                npred_src=npred_src,
                npred_bkg=npred_bkg,
                background_model=bkg_model,
                obs=obs
            ))

        self.result = results

    def est_errors(self):
        """Estimate parameter errors."""
        if self.err_method == 'sherpa':
            self._est_errors_sherpa()
        else:
            raise NotImplementedError('{}'.format(self.err_method))

        for res in self.result:
            res.covar_axis = self.covar_axis
            res.covariance = self.covariance
            res.model.parameters.set_parameter_covariance(self.covariance, self.covar_axis)

    def _est_errors_sherpa(self):
        """Wrapper around Sherpa error estimator."""
        covar = self._sherpa_fit.est_errors()
        self.covar_axis = [par.split('.')[-1] for par in covar.parnames]
        self.covariance = copy.deepcopy(covar.extra_output)

    def run(self, outdir=None):
        """Run all steps and write result to disk.

        Parameters
        ----------
        outdir : Path, str
            directory to write results files to (if given)
        """
        log.info('Running {}'.format(self))

        self.fit()
        self.est_errors()

        if outdir is not None:
            self._write_result(outdir)

    def _write_result(self, outdir):
        outdir = make_path(outdir)
        outdir.mkdir(exist_ok=True, parents=True)

        # Assume only one model is fit to all data
        modelname = self.result[0].model.__class__.__name__
        filename = outdir / 'fit_result_{}.yaml'.format(modelname)
        log.info('Writing {}'.format(filename))
        self.result[0].to_yaml(filename)

# 需要导入模块: from sherpa.fit import Fit [as 别名]
# 或者: from sherpa.fit.Fit import est_errors [as 别名]
class test_sim(SherpaTestCase):

    _fit_results_bench = {
        'rstat': 89.29503933428586,
        'qval': 0.0,
        'succeeded': 1,
        'numpoints': 100,
        'dof': 95,
        'nfev': 93,
        'statval': 8483.0287367571564,
        'parnames': ['p1.gamma', 'p1.ampl', 'g1.fwhm',
                     'g1.pos', 'g1.ampl'], 
        'parvals': numpy.array(
            [1.0701938169914813,
             9.1826254677279469,
             2.5862083052721028,
             2.601619746022207,
             47.262657692418749])
        }

    _x = numpy.arange(0.1, 10.1, 0.1)
    _y = numpy.array(
        [ 114, 47, 35, 30, 40, 27, 30, 26, 24, 20, 26, 35,
          29, 28, 34, 36, 43, 39, 33, 47, 44, 46, 53, 56,
          52, 53, 49, 57, 49, 36, 33, 42, 49, 45, 42, 32,
          31, 34, 18, 24, 25, 11, 17, 17, 11,  9,  8,  5,
           4, 10,  3,  4,  6,  3,  0,  2,  4,  4,  0,  1,
           2,  0,  3,  3,  0,  2,  1,  2,  3,  0,  1,  0,
           1,  0,  0,  1,  3,  3,  0,  2,  0,  0,  1,  2,
           0,  1,  0,  1,  1,  0,  1,  1,  1,  1,  1,  1,
           1,  0,  1,  0
          ]
        )
    _err = numpy.ones(100)*0.4


    def setUp(self):
        data = Data1D('fake', self._x, self._y, self._err)

        g1 = Gauss1D('g1')
        g1.fwhm.set(1.0, _tiny, _max, frozen=False)
        g1.pos.set(1.0, -_max, _max, frozen=False)
        g1.ampl.set(1.0, -_max, _max, frozen=False)
        p1 = PowLaw1D('p1')
        p1.gamma.set(1.0, -10, 10, frozen=False)
        p1.ampl.set(1.0, 0.0, _max, frozen=False)
        p1.ref.set(1.0, -_max, _max, frozen=True)      
        model = p1 + g1

        method = LevMar()
        method.config['maxfev'] = 10000
        method.config['ftol'] = float(_eps)
        method.config['epsfcn'] = float(_eps)
        method.config['gtol'] = float(_eps)
        method.config['xtol'] = float(_eps)
        method.config['factor'] = float(100)

        self.fit = Fit(data, model, Chi2DataVar(), method, Covariance())
        results = self.fit.fit()

        for key in ["succeeded", "numpoints", "nfev"]:
            assert self._fit_results_bench[key] == int(getattr(results, key))

        for key in ["rstat", "qval", "statval", "dof"]:
            assert numpy.allclose(float(self._fit_results_bench[key]),
                                  float(getattr(results, key)),
                                  1.e-7, 1.e-7)

        for key in ["parvals"]:
            try:
                assert numpy.allclose(self._fit_results_bench[key],
                                      getattr(results, key),
                                      1.e-4, 1.e-4)
            except AssertionError:
                print 'parvals bench: ', self._fit_results_bench[key]
                print 'parvals fit:   ', getattr(results, key)
                print 'results', results
                raise

        covresults = self.fit.est_errors()
        self.dof = results.dof
        self.mu = numpy.array(results.parvals)
        self.cov = numpy.array(covresults.extra_output)
        self.num = 10


    def test_student_t(self):
        multivariate_t(self.mu, self.cov, self.dof, self.num)
        
    def test_cauchy(self):
        multivariate_cauchy(self.mu, self.cov, self.num)       

    def test_parameter_scale_vector(self):
        ps = ParameterScaleVector()
        ps.get_scales(self.fit)

    def test_parameter_scale_matrix(self):
        ps = ParameterScaleMatrix()
        ps.get_scales(self.fit)

#.........这里部分代码省略.........

# 需要导入模块: from sherpa.fit import Fit [as 别名]
# 或者: from sherpa.fit.Fit import est_errors [as 别名]

#.........这里部分代码省略.........
            input coordinates (dependent for 1D fits or independent for 2D fits)
        z : array or list of arrays (optional)
            input coordinates (dependent for 2D fits)
        xbinsize : array or list of arrays (optional)
            an array of xbinsizes in x  - this will be x -/+ (binsize  / 2.0)
        ybinsize : array or list of arrays (optional)
            an array of xbinsizes in y  - this will be y -/+ (ybinsize / 2.0)
        err : array or list of arrays (optional)
            an array of errors in dependent variable
        bkg : array or list of arrays (optional)
            this will act as background data
        bkg_sale : float or list of floats (optional)
            the scaling factor for the dataset if a single value
            is supplied it will be copied for each dataset
        **kwargs :
            keyword arguments will be passed on to sherpa fit routine

        Returns
        -------
        model_copy : `astropy.modeling.FittableModel` or a list of models.
            a copy of the input model with parameters set by the fitter
        """

        tie_list = []
        try:
            n_inputs = models[0].n_inputs
        except TypeError:
            n_inputs = models.n_inputs

        self._data = Dataset(n_inputs, x, y, z, xbinsize, ybinsize, err, bkg, bkg_scale)

        if self._data.ndata > 1:

            if len(models) == 1:
                self._fitmodel = ConvertedModel([models.copy() for _ in xrange(self._data.ndata)], tie_list)
                # Copy the model so each data set has the same model!
            elif len(models) == self._data.ndata:
                self._fitmodel = ConvertedModel(models, tie_list)
            else:
                raise Exception("Don't know how to handle multiple models "
                                "unless there is one foreach dataset")
        else:
            if len(models) > 1:
                self._data.make_simfit(len(models))
                self._fitmodel = ConvertedModel(models, tie_list)
            else:
                self._fitmodel = ConvertedModel(models)

        self._fitter = Fit(self._data.data, self._fitmodel.sherpa_model, self._stat_method, self._opt_method, self._est_method, **kwargs)
        self.fit_info = self._fitter.fit()

        return self._fitmodel.get_astropy_model()

    def est_errors(self, sigma=None, maxiters=None, numcores=1, methoddict=None, parlist=None):
        """
        Use sherpa error estimators based on the last fit.

        Parameters
        ----------
        sigma: float
            this will be set as the confidance interval for which the errors are found too.
        maxiters: int
            the maximum number of iterations the error estimator will run before giving up.
        methoddict: dict
            !not quite sure couldn't figure this one out yet!
        parlist: list
            a list of parameters to find the confidance interval of if none are provided all free
            parameters will be estimated.
        """
        if self._fitter is None:
            ValueError("Must complete a valid fit before errors can be calculated")
        if sigma is not None:
            self._fitter.estmethod.config['sigma'] = sigma
        if maxiters is not None:
            self._fitter.estmethod.config['maxiters'] = maxiters
        if 'numcores' in self._fitter.estmethod.config:
            if not numcores == self._fitter.estmethod.config['numcores']:
                self._fitter.estmethod.config['numcores'] = numcores

        self.error_info = self._fitter.est_errors(methoddict=methoddict, parlist=parlist)
        pnames = [p.split(".", 1)[-1] for p in self.error_info.parnames]  # this is to remove the model name
        return pnames, self.error_info.parvals, self.error_info.parmins, self.error_info.parmaxes

    @property
    def _est_config(self):
        """This returns the est.config property"""
        return self._est_method.config

    @property
    def _opt_config(self):
        """This returns the opt.config property"""
        return self._opt_method.config

    # Here is the MCMC wrapper!
    @format_doc("{__doc__}\n{mcmc_doc}",mcmc_doc=SherpaMCMC.__doc__)
    def get_sampler(self, *args, **kwargs):
        """
        This returns and instance of `SherpaMCMC` with it's self as the fitter
        """
        return SherpaMCMC(self, *args, **kwargs)

# 需要导入模块: from sherpa.fit import Fit [as 别名]
# 或者: from sherpa.fit.Fit import est_errors [as 别名]
class test_sim(SherpaTestCase):
    def setUp(self):
        # self.startdir = os.getcwd()
        self.old_level = logger.getEffectiveLevel()
        logger.setLevel(logging.CRITICAL)

        datadir = SherpaTestCase.datadir
        if datadir is None:
            return

        pha = os.path.join(datadir, "refake_0934_1_21_1e4.fak")
        rmf = os.path.join(datadir, "ccdid7_default.rmf")
        arf = os.path.join(datadir, "quiet_0934.arf")

        self.simarf = os.path.join(datadir, "aref_sample.fits")
        self.pcaarf = os.path.join(datadir, "aref_Cedge.fits")

        data = read_pha(pha)
        data.ignore(None, 0.3)
        data.ignore(7.0, None)

        rsp = Response1D(data)
        self.abs1 = XSwabs("abs1")
        self.p1 = XSpowerlaw("p1")
        model = rsp(self.abs1 * self.p1)

        self.fit = Fit(data, model, CStat(), NelderMead(), Covariance())

    def tearDown(self):
        # os.chdir(self.startdir)
        logger.setLevel(self.old_level)

    @needs_data
    def test_pragbayes_simarf(self):
        datadir = SherpaTestCase.datadir
        if datadir is None:
            return

        mcmc = sim.MCMC()

        self.abs1.nh = 0.092886
        self.p1.phoindex = 0.994544
        self.p1.norm = 9.26369

        mcmc.set_sampler("PragBayes")
        mcmc.set_sampler_opt("simarf", self.simarf)
        mcmc.set_sampler_opt("p_M", 0.5)
        mcmc.set_sampler_opt("nsubiter", 7)

        covar_results = self.fit.est_errors()
        cov = covar_results.extra_output

        niter = 10
        stats, accept, params = mcmc.get_draws(self.fit, cov, niter=niter)
        # try:
        #     assert (covar_results.parmaxes < params.std(1)).all()
        # except AssertionError:
        #     print 'covar: ', str(covar_results.parmaxes)
        #     print 'param: ', str(params.std(1))
        #     raise

    @needs_data
    def test_pragbayes_pcaarf(self):
        datadir = SherpaTestCase.datadir
        if datadir is None:
            return

        mcmc = sim.MCMC()

        self.abs1.nh = 0.092886
        self.p1.phoindex = 0.994544
        self.p1.norm = 9.26369

        mcmc.set_sampler("pragBayes")
        mcmc.set_sampler_opt("simarf", self.pcaarf)
        mcmc.set_sampler_opt("p_M", 0.5)
        mcmc.set_sampler_opt("nsubiter", 5)

        covar_results = self.fit.est_errors()
        cov = covar_results.extra_output

        niter = 10
        stats, accept, params = mcmc.get_draws(self.fit, cov, niter=niter)