Python scipy.stats.norm 模块，fit() 实例源码

def test_megkde_2d_basic():
    # Draw from normal, fit KDE, see if sampling from kde's pdf recovers norm
    np.random.seed(1)
    data = np.random.multivariate_normal([0, 1], [[1.0, 0.], [0., 0.75 ** 2]], size=10000)
    xs, ys = np.linspace(-4, 4, 50), np.linspace(-4, 4, 50)
    xx, yy = np.meshgrid(xs, ys, indexing='ij')
    samps = np.vstack((xx.flatten(), yy.flatten())).T
    zs = MegKDE(data).evaluate(samps).reshape(xx.shape)
    zs_x = zs.sum(axis=1)
    zs_y = zs.sum(axis=0)
    cs_x = zs_x.cumsum()
    cs_x /= cs_x[-1]
    cs_x[0] = 0
    cs_y = zs_y.cumsum()
    cs_y /= cs_y[-1]
    cs_y[0] = 0
    samps_x = interp1d(cs_x, xs)(np.random.uniform(size=10000))
    samps_y = interp1d(cs_y, ys)(np.random.uniform(size=10000))
    mu_x, std_x = norm.fit(samps_x)
    mu_y, std_y = norm.fit(samps_y)
    assert np.isclose(mu_x, 0, atol=0.1)
    assert np.isclose(std_x, 1.0, atol=0.1)
    assert np.isclose(mu_y, 1, atol=0.1)
    assert np.isclose(std_y, 0.75, atol=0.1)

def fit(self, X, y):
        self.base_models_ = [list() for x in self.base_models]
        self.meta_model_ = clone(self.meta_model)
        kfold = KFold(n_splits=self.n_folds, shuffle=True, random_state=15)

        # train cloned base models then create out-of-fold predictions that are needed to train the cloned meta-model
        out_of_fold_predictions = np.zeros((X.shape[0], len(self.base_models)))
        for i, model in enumerate(self.base_models):
            for train_index, holdout_index in kfold.split(X, y):
                instance = clone(model)
                self.base_models_[i].append(instance)
                instance.fit(X[train_index], y[train_index])
                y_pred = instance.predict(X[holdout_index])
                out_of_fold_predictions[holdout_index, i] = y_pred

        # now train the cloned  meta-model using the out-of-fold predictions as new feature
        self.meta_model_.fit(out_of_fold_predictions, y)
        return self

    # do the predictions of all base models on the test data and use the averaged predictions as 
    #meta-features for the final prediction which is done by the meta-model

def test_megkde_1d_basic():
    # Draw from normal, fit KDE, see if sampling from kde's pdf recovers norm
    np.random.seed(0)
    data = np.random.normal(loc=0, scale=1.0, size=2000)
    xs = np.linspace(-3, 3, 100)
    ys = MegKDE(data).evaluate(xs)
    cs = ys.cumsum()
    cs /= cs[-1]
    cs[0] = 0
    samps = interp1d(cs, xs)(np.random.uniform(size=10000))
    mu, std = norm.fit(samps)
    assert np.isclose(mu, 0, atol=0.1)
    assert np.isclose(std, 1.0, atol=0.1)

def test_megkde_1d_uniform_weight():
    # Draw from normal, fit KDE, see if sampling from kde's pdf recovers norm
    np.random.seed(0)
    data = np.random.normal(loc=0, scale=1.0, size=2000)
    xs = np.linspace(-3, 3, 100)
    ys = MegKDE(data, weights=np.ones(2000)).evaluate(xs)
    cs = ys.cumsum()
    cs /= cs[-1]
    cs[0] = 0
    samps = interp1d(cs, xs)(np.random.uniform(size=10000))
    mu, std = norm.fit(samps)
    assert np.isclose(mu, 0, atol=0.1)
    assert np.isclose(std, 1.0, atol=0.1)

def test_megkde_1d_changing_weights():
    # Draw from normal, fit KDE, see if sampling from kde's pdf recovers norm
    np.random.seed(0)
    xs = np.linspace(-3, 3, 1000)
    weights = norm.pdf(xs)
    ys = MegKDE(xs, weights=weights).evaluate(xs)
    cs = ys.cumsum()
    cs /= cs[-1]
    cs[0] = 0
    samps = interp1d(cs, xs)(np.random.uniform(size=10000))
    mu, std = norm.fit(samps)
    assert np.isclose(mu, 0, atol=0.1)
    assert np.isclose(std, 1.0, atol=0.1)

def logistic_fidelity(self):
        #group data and assign state labels
        gnd_features = np.hstack([np.real(self.ground_data.T),
                                np.imag(self.ground_data.T)])
        ex_features = np.hstack([np.real(self.excited_data.T),
                                np.imag(self.excited_data.T)])
        #liblinear wants arrays in C order
        features = np.ascontiguousarray(np.vstack([gnd_features, ex_features]))
        state = np.ascontiguousarray(np.hstack([np.zeros(self.ground_data.shape[1]),
                                                np.ones(self.excited_data.shape[1])]))
        #Set up logistic regression with cross-validation using liblinear.
        #Cs sets the inverse of the regularization strength, which will be optimized
        #through cross-validation. Uses the default Stratified K-Folds
        #CV generator, with 3 folds.
        #This is set up to be as consistent with the MATLAB implementation
        #as I can make it. --GJR
        Cs = np.logspace(-1,2,5)
        logreg = LogisticRegressionCV(Cs, cv=3, solver='liblinear')
        logreg.fit(features, state) #fit the model
        predictions = logreg.predict(features) #in-place classification
        score = logreg.score(features,state) #mean accuracy of classification
        N = len(predictions)
        S = np.sum(predictions == state) #how many we got right
        #now calculate confidence intervals
        c = 0.95
        flo = betaincinv(S+1, N-S+1, (1-c)/2., )
        fhi = betaincinv(S+1, N-S+1, (1+c)/2., )
        logger.info(("In-place logistic regression fidelity: " +
                "{:.2f}% ({:.2f}, {:.2f})".format(100*score, 100*flo, 100*fhi)))

def _rerender(self):
        nmr_maps = len(self.maps_to_show)
        if self._show_trace:
            nmr_maps *= 2

        grid = GridSpec(nmr_maps, 1, left=0.04, right=0.96, top=0.94, bottom=0.06, hspace=0.2)

        i = 0
        for map_name in self.maps_to_show:
            samples = self._voxels[map_name]

            if self._sample_indices is not None:
                samples = samples[:, self._sample_indices]

            title = map_name
            if map_name in self.names:
                title = self.names[map_name]

            if isinstance(self._nmr_bins, dict) and map_name in self._nmr_bins:
                nmr_bins = self._nmr_bins[map_name]
            else:
                nmr_bins = self._nmr_bins

            hist_plot = plt.subplot(grid[i])
            n, bins, patches = hist_plot.hist(np.nan_to_num(samples[self.voxel_ind, :]), nmr_bins, normed=True)
            plt.title(title)
            i += 1

            if self._fit_gaussian:
                mu, sigma = norm.fit(samples[self.voxel_ind, :])
                bincenters = 0.5*(bins[1:] + bins[:-1])
                y = mlab.normpdf(bincenters, mu, sigma)
                hist_plot.plot(bincenters, y, 'r', linewidth=1)

            if self._show_trace:
                trace_plot = plt.subplot(grid[i])
                trace_plot.plot(samples[self.voxel_ind, :])
                i += 1

def display_distrib(pd, feature):
    plt.figure()
    sns.distplot(pd[feature].dropna() , fit=norm);
    (mu, sigma) = norm.fit(pd[feature].dropna())    

    plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)], loc='best')
    plt.ylabel('Frequency')
    plt.title('SalePrice distribution')
    plt.show()

def __init__(self, models):
        self.models = models

    # we define clones of the original models to fit the data in

def fit(self, X, y):
        self.models_ = [clone(x) for x in self.models]

        # Train cloned base models
        for model in self.models_:
            model.fit(X, y)

        return self

    # now we do the predictions for cloned models and average them

def __init__(self, base_models, meta_model, n_folds=5):
        self.base_models = base_models
        self.meta_model = meta_model
        self.n_folds = n_folds

    # we again fit the data on clones of the original models

def central_limit(rvs, n, length):
    rv_mean = np.zeros(length)
    for i in range(n):
        rv = rvs(size=length)
        rv_mean = rv_mean + rv
    rv_mean = rv_mean / n
    gaussian_params = norm_dist.fit(rv_mean)
    gaussian = norm_dist(gaussian_params[0], gaussian_params[1])
    return rv_mean, gaussian

def graphing(what, data, who):

    kek = data[who][what]    
    fig = figure()
    ax = fig.add_subplot(111)
    if '\xa0' in str(kek):
       kek = (int(kek.replace(u'\xa0', '')))

    data = stuff[what]
    mean, sigma = norm.fit(data)
    plt.hist(data, bins = 1000, normed=True, alpha=0.4)

    xmin, xmax = plt.xlim()
    x = np.linspace(xmin, xmax, 1000)
    p = norm.pdf(x, mean, sigma)
    plt.plot(x, p, 'k', linewidth=2)
    title = f'{what}\nFit results: mean = {round(mean,2)}, sigma = {round(sigma,2)}'
    plt.title(title)

    two, eight = norm.interval(.999, loc = mean, scale = sigma)
    plt.xlim(0, eight)

    new = norm.pdf(kek, mean, sigma)
    plt.plot([kek, kek], [0, new], color='r', linestyle='--')
    if what == 'Death/Game':
        at = round((1-norm.cdf(kek, mean, sigma))*100, 2)
    else:
        at = round(norm.cdf(kek, mean, sigma)*100, 2)
    text = f'You are better than \n{at}% players'
    ax.annotate(text, xy=(kek, new), xytext=(0.7, 0.8), textcoords='axes fraction',
                arrowprops=dict(facecolor='g', color='g'))

    what = what.replace('/', '-')
    #plt.show()
    plt.savefig(f'Y:/python/graphs/data/{what}.png')
    #plt.close()

def scales(stuff):
    limits = {}
    for key in stuff.keys():
        try:
            data = stuff[key]
            mean, sigma = norm.fit(data)

            one, nine = norm.interval(.99, loc = mean, scale = sigma)
            two, eight = norm.interval(.9, loc = mean, scale = sigma)
            three, seven = norm.interval(.7, loc = mean, scale = sigma)
            four, six = norm.interval(.5, loc = mean, scale = sigma)
            five = mean

            if key != 'Karma':
                if four < 0:
                    four = five/2
                if three < 0:
                    three = four/2
                if two < 0:
                    two = three/2
                if one < 0:
                    one = two/2

            if key == 'Fleet Strength':
                if nine > 200:
                    nine = 200
            if key == 'Eff rating':
                if nine > 9000:
                    nine = 9000

            if key == 'Death/Game' or key == 'Total deaths':
                limits[key] = {one: 9, two: 8, three: 7, four: 6, five: 5,
                               six: 4, seven: 3, eight: 2, nine: 1}      
            else:
                limits[key] = {one: 1, two: 2, three: 3, four: 4, five: 5,
                               six: 6, seven: 7, eight: 8, nine: 9}                

        except: continue
    return limits

def fit_distribution(self):
    y = np.zeros(self.num_sent)
    iter = 0
    for key in self.stats:
      iter_end = self.stats[key].num_sents + iter
      y[iter:iter_end] = key
      iter = iter_end
    mu, std = norm.fit(y)
    self.distr = norm(mu, std)

def __get_est_price_mode_pe(realtime, code, years):

    q_str = 'code==' + '\"' + code + '\"'

    pe_obj = {}
    eps = 0
    for y in range(datetime.now().year - years, datetime.now().year + 1):
        for q in range(1, 5):
        quarter = str(y)+'q'+str(q)
        if (os.path.exists(PREFIX+'/'+quarter+'.csv')):
        r = __pd_read_report(quarter)
        if (len(r.query(q_str)) > 0):
            # save all pe history and latest eps
                    tmp_pe, tmp_eps = __get_pe_and_eps(code, quarter)
                    if (isinstance(tmp_pe, float) & isinstance(tmp_eps, float)):
                        pe_obj[quarter] = tmp_pe
                        eps = tmp_eps
                        log.debug('%s pe: %.2f, eps: %.2f', quarter, pe_obj[quarter], eps)
                    else:
                        log.warn('skip %s', quarter)
                        continue
    #sorted(pe_obj)
    #log.debug(pe_obj)
    arr = pe_obj.values()
    mu, std = norm.fit(arr)

    if (realtime):
    d = datetime.now()
    today = __pd_read_today_all()
    close = round(today.query(q_str).trade.values[0], 2)
    else:
    k, d = __get_k_data_of_last_trade_day(code)
    close = round(k.close.values[0], 2)
    log.info('%s price: %.2f @ pe %.2f, eps %.2f', d.strftime("%Y-%m-%d"), close, close/eps, eps)
    log.info('mu, std: %.2f, %.2f', mu, std)

    growth = __get_growth(code, years)

    left = __estimation_formula_bg_dynamic(growth, eps, mu - std)
    centrum = __estimation_formula_bg_dynamic(growth, eps, mu)
    right = __estimation_formula_bg_dynamic(growth, eps, mu + std)
    value = __estimation_formula_bg_dynamic(growth, eps, 8.5)

    log.info('est dynamic: %.2f~%.2f~%.2f', left, centrum, right)
    log.info('est value: %.2f', value)
    log.info('range from left: %.2f%%', (close-left)/left*100.0)
    log.info('position: %.2f%%', (close-left)/(right-left)*100.0)

    return left, centrum, right, value

def __get_est_price_mode_pb(realtime, code, years):

    q_str = 'code==' + '\"' + code + '\"'

    pb_obj = {}
    bvps = 0
    for y in range(datetime.now().year - years, datetime.now().year + 1):
        for q in range(1, 5):
        quarter = str(y)+'q'+str(q)
        if (os.path.exists(PREFIX+'/'+quarter+'.csv')):
        r = __pd_read_report(quarter)
        if (len(r.query(q_str)) > 0):
            # save all pb history and latest bvps
                    tmp_pb, tmp_bvps = __get_pb_and_bvps(code, quarter)
                    if (isinstance(tmp_pb, float) & isinstance(tmp_bvps, float)):
                        pb_obj[quarter] = tmp_pb
                        bvps = tmp_bvps
                        log.debug('%s pb: %.2f, bvps: %.2f', quarter, pb_obj[quarter], bvps)
                    else:
                        log.warn('skip %s', quarter)
                        continue
    #sorted(pb_obj)
    #log.debug(pb_obj)
    arr = pb_obj.values()
    mu, std = norm.fit(arr)

    if (realtime):
    d = datetime.now()
    today = __pd_read_today_all()
    close = round(today.query(q_str).trade.values[0], 2)
    else:
    k, d = __get_k_data_of_last_trade_day(code)
    close = round(k.close.values[0], 2)
    log.info('%s price: %.2f @ pb %.2f, bvps %.2f', d.strftime("%Y-%m-%d"), close, close/bvps, bvps)
    log.info('mu, std: %.2f, %.2f', mu, std)

    left = __estimation_formula_pb(bvps, mu - std)
    centrum = __estimation_formula_pb(bvps, mu)
    right = __estimation_formula_pb(bvps, mu + std)
    value = __estimation_formula_pb(bvps, 1.0)

    log.info('est dynamic: %.2f~%.2f~%.2f', left, centrum, right)
    log.info('est value: %.2f', value)
    log.info('range from left: %.2f%%', (close-left)/left*100.0)
    log.info('position: %.2f%%', (close-left)/(right-left)*100.0)

    return left, centrum, right, value

def show(self, voxel_ind=0, names=None, maps_to_show=None, to_file=None, block=True, maximize=False,
             show_trace=True, nmr_bins=20, window_title=None, show_sliders=True, fit_gaussian=True,
             figure_options=None, sample_indices=None):
        """Show the samples per voxel.

        Args:
            voxel_ind (int): the voxel to show the samples from.
            names (dict): A list of names for the different maps. Use as ``{map_name: display_name}`` that is,
                 the key is the name of the map in the volumes dictionary and the display name is the string that will
                 be used as title for that map.
            maps_to_show (:class:`list`): A list of maps to show.
                The items in this list must correspond to the keys in the volumes dictionary.
            to_file (string, optional, default None): If to_file is not None it is supposed
                to be a filename where the image will be saved.
                If not set to None, nothing will be displayed, the results will directly be saved.
                Already existing items will be overwritten.
            block (boolean): If we want to block after calling the plots or not. Set this to False if you
                do not want the routine to block after drawing. In doing so you manually need to block.
            maximize (boolean): if we want to display the window maximized or not
            show_trace (boolean): if we show the trace of each map or not
            nmr_bins (dict or int): either a single value or one per map name
            show_sliders (boolean): if we show the slider or not
            fit_gaussian (boolean): if we fit and show a normal distribution (Gaussian) to the histogram or not
            window_title (str): the title of the window. If None, the default title is used
            figure_options (dict) options for the figure
            sample_indices (list): the list of sample indices to use
        """
        figure_options = figure_options or {'figsize': (18, 16)}
        self._figure = plt.figure(**figure_options)

        if names:
            self.names = names
        if maps_to_show:
            self.maps_to_show = maps_to_show
        self.voxel_ind = voxel_ind
        self._nmr_bins = nmr_bins or self._nmr_bins
        self._show_trace = show_trace
        self.show_sliders = show_sliders
        self._fit_gaussian = fit_gaussian
        self._sample_indices = sample_indices
        self._setup()

        if maximize:
            mng = plt.get_current_fig_manager()
            mng.window.showMaximized()

        if window_title:
            mng = plt.get_current_fig_manager()
            mng.canvas.set_window_title(window_title)

        if to_file:
            plt.savefig(to_file)
            plt.close()
        else:
            plt.draw()
            if block:
                plt.show(True)

def calibrateNoise(self):

        bg, bgstd, las, time, conc, ok = CalibrationDialog.setExt()

        _np.asarray(bg)
        _np.asarray(bgstd)
        _np.asarray(las)
        _np.asarray(time)
        _np.asarray(conc)

        x_3d = _np.array([conc, las, time])
        p0 = [1, 1]
        fitParamsBg, fitCovariances = curve_fit(fitFuncBg, x_3d, bg, p0)
        print(' fit coefficients :\n', fitParamsBg)

        # SET VALUES TO PARAMETER
        self.lasercEdit.setValue(fitParamsBg[0])
        self.imagercEdit.setValue(fitParamsBg[1])

        x_3dStd = _np.array([las, time, bg])
        p0S = [1, 1, 1]
        fitParamsStd, fitCovariances = curve_fit(fitFuncStd, x_3dStd, bgstd, p0S)

        print(' fit coefficients2:\n', fitParamsStd)

        self.EquationAEdit.setValue(fitParamsStd[0])
        self.EquationBEdit.setValue(fitParamsStd[1])
        self.EquationCEdit.setValue(fitParamsStd[2])

        # Noise model working point

        figure4 = plt.figure()

        # Background
        bgmodel = fitFuncBg(x_3d, fitParamsBg[0], fitParamsBg[1])
        ax1 = figure4.add_subplot(121)
        ax1.cla()
        ax1.plot(bg, bgmodel, 'o')
        x = _np.linspace(*ax1.get_xlim())
        ax1.plot(x, x)
        title = "Background Model:"
        ax1.set_title(title)

        # Std
        bgmodelstd = fitFuncStd(x_3dStd, fitParamsStd[0], fitParamsStd[1], fitParamsStd[2])
        ax2 = figure4.add_subplot(122)
        ax2.cla()
        ax2.plot(bgstd, bgmodelstd, 'o')
        x = _np.linspace(*ax2.get_xlim())
        ax2.plot(x, x)
        title = "Background Model Std:"
        ax2.set_title(title)

        figure4.show()

def run(self):
        """
        NOTE: Gibbs sampling is a way of selecting variables one at a time instead of all at once. This is beneficial in
        high dimensional variable space because it will increase the probability of accepting a new sample. It isn't
        implemented here, but it might be worth keeping in mind.
        """
        counter = 0
        while not self._check_convergence() and (counter <= self.steps or self.steps == 0):
            tmp_dump = self.dumpfile + ".temp"
            with open(tmp_dump, "wb") as ofile:
                dump_obj = [chain._dump_obj() for chain in self.chains]
                dill.dump(dump_obj, ofile, protocol=-1)
            shutil.move(tmp_dump, self.dumpfile)
            counter += 1
            child_list = OrderedDict()
            for chain in self.chains:  # Note that this will spin off (c * w) new processes, where c=chains, w=walkers
                for walker in chain.walkers:
                    # Start new process
                    func_args = []
                    for variable in walker.variables:
                        if walker.lava:
                            variable.draw_random()
                        else:
                            variable.draw_new_value(walker.heat)
                        func_args.append(variable.draw_value)

                    # Always add a new seed for the target function
                    func_args.append(self.rand_gen.randint(1, 999999999999999))
                    p = Process(target=self.mc_step_run, args=(walker, [func_args]))
                    p.start()
                    child_list[walker.name] = p

            # wait for remaining processes to complete
            while len(child_list) > 0:
                for _name, child in child_list.items():
                    if child.is_alive():
                        continue
                    else:
                        del child_list[_name]
                        break

            for chain in self.chains:
                # Get the normalized standard deviation among all historical walker scores for this chain
                history_series = pd.Series([score for walker in chain.walkers for score in walker.score_history])
                mu, std = norm.fit(history_series)
                for walker in chain.walkers:
                    self.step_parse(walker, std)
                    if self.best["score"] is None or walker.current_score > self.best["score"]:
                        self.best["score"] = walker.current_score
                        self.best["variables"] = OrderedDict([(x.name, x.current_value) for x in walker.variables])

            for chain in self.chains:
                chain.swap_hot_cold()
                # Send output to files
                if counter % self.sample_rate == 0:
                    chain.step_counter += 1
                    chain.write_sample()
        return