Python scipy.stats.norm.fit() Examples

def doc_thres(model, data, model_type, scale = 2.):
    train_y = data['train_set_Y']
    if "mlp" in model_type:
        train_pred = model.predict(data['train_set_X'])
    else:
        train_pred = model.predict(data['train_set_idx_X'])
    #print train_y.shape, train_pred.shape 
    def fit(prob_pos_X):
        prob_pos = [p for p in prob_pos_X]+[2-p for p in prob_pos_X]
        pos_mu, pos_std = dist_model.fit(prob_pos)
        return pos_std
    stds = []
    for c in range(train_pred.shape[-1]):
        idx = [train_y == c]
        c_pred = train_pred[idx]
        c_prob = c_pred[:,c]
        std = fit(c_prob)
        stds.append(std)
    thres = [max(0.5, 1. - scale * std) for std in stds]
    return thres 

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_predict(self, X, y=None):
        """Fit the model according to the given training data and predict if a
        particular training sample is an outlier or not.

        Parameters
        ----------
        X : array-like of shape (n_samples, n_features)
            Training Data.

        y : ignored

        Returns
        -------
        y_pred : array-like of shape (n_samples,)
            Return -1 for outliers and +1 for inliers.
        """

        if hasattr(self, 'novelty'):
            check_novelty(self.novelty, 'fit_predict')

        return self.fit(X).predict() 

def plot_t_value_hist(
	img_path='~/ni_data/ofM.dr/l1/as_composite/sub-5703/ses-ofM/sub-5703_ses-ofM_task-EPI_CBV_chr_longSOA_tstat.nii.gz',
	roi_path='~/ni_data/templates/roi/DSURQEc_ctx.nii.gz',
	mask_path='/usr/share/mouse-brain-atlases/dsurqec_200micron_mask.nii',
	save_as='~/qc_tvalues.pdf',
	):
	"""Make t-value histogram plot"""

	f, axarr = plt.subplots(1, sharex=True)

	roi = nib.load(path.expanduser(roi_path))
	roi_data = roi.get_data()
	mask = nib.load(path.expanduser(mask_path))
	mask_data = mask.get_data()
	idx = np.nonzero(np.multiply(roi_data,mask_data))
	img = nib.load(path.expanduser(img_path))
	data = img.get_data()[idx]
	(mu, sigma) = norm.fit(data)
	n, bins, patches = axarr.hist(data,'auto',normed=1, facecolor='green', alpha=0.75)
	y = mlab.normpdf(bins, mu, sigma)

	axarr.plot(bins, y, 'r--', linewidth=2)
	axarr.set_title('Histogram of t-values $\mathrm{(\mu=%.3f,\ \sigma=%.3f}$)' %(mu, sigma))
	axarr.set_xlabel('t-values')
	plt.savefig(path.expanduser(save_as)) 

def draw_pdf(var1, var2=None, label1='var1', label2='var2', xlab='variable', \
        ylab='Normalized frequency', linewidth=2, framealpha=0.3, \
        markersize=5, fontsize=18, xmax=18, binfac=1): 
    
    bins = range(xmax+1)
    (mu, sigma) = norm.fit(var1)
    y1 = mlab.normpdf(bins, mu, sigma)
    figure()
    grid('on')
    plot(bins, y1, 'b-o', linewidth=linewidth, markersize=markersize, label=label1)
    if var2 != None:
        (mu, sigma) = norm.fit(var2)
        y2 = mlab.normpdf(bins, mu, sigma)
        plot(bins, y2, 'r-o', linewidth=linewidth, markersize=markersize, label=label2)
    xlabel(xlab, fontsize=fontsize)
    ylabel(ylab, fontsize=fontsize)
    xlim(0,xmax)
    tick_params(axis='both', which='major', labelsize=fontsize)
    if var2 != None:
        legend(fancybox=True, framealpha=framealpha, loc='upper right', fontsize=fontsize)
    
    return 

def visualize_stochastic_histogram(
        self, log_evidences, max_log_evidence, during_analysis
    ):

        if log_evidences is None:
            return

        plotter = self.plotter.plotter_with_new_output(
            path=self.plotter.output.path + "other/"
        )

        if self.plot_stochastic_histogram and not during_analysis:

            bins = conf.instance.general.get(
                "inversion", "stochastic_histogram_bins", int
            )

            (mu, sigma) = norm.fit(log_evidences)
            n, bins, patches = plt.hist(x=log_evidences, bins=bins, density=1)
            y = norm.pdf(bins, mu, sigma)
            plt.plot(bins, y, "--")
            plt.xlabel("log evidence")
            plt.title("Stochastic Log Evidence Histogram")
            plt.axvline(max_log_evidence, color="r")
            plt.savefig(
                f"{plotter.output.path}stochastic_histogram.png", bbox_inches="tight"
            ) 

def run(opts):

	scriptwd = os.path.dirname(os.path.dirname(os.path.realpath(__file__)))
	img, data = loadmgh(opts.input[0])
	outname = str(opts.output[0])
	hemi = str(opts.hemi[0])
	fwhm = float(opts.geodesicfwhm[0])
	sigma = fwhm / np.sqrt(8 * np.log(2))

	# load distances and indices 
	dist = np.load(opts.distanceslist[0])
	basenamefwhm = opts.distanceslist[0].split('_distances.npy',1)[0]
	indices = np.load('%s_indices.npy' % basenamefwhm)

	cortex_index = np.load('%s/adjacency_sets/%s_cortex_mask_index.npy' % (scriptwd,hemi))
	distmask_index_start = np.in1d(indices[:,0], cortex_index)
	distmask_index_end = np.in1d(indices[:,1], cortex_index)
	distmask_index = distmask_index_start & distmask_index_end

	# mask out non-cortex values
	indices_masked = indices[distmask_index]
	dist_masked = dist[distmask_index]

	if opts.correct_surface:
		multipler = opts.correct_surface[0]
		y = data[cortex_index].flatten()
		(mu, sigma) = norm.fit(y)
		cthresh = multipler*sigma + mu
		print("The upper threshold is: %1.4f" % cthresh)
		data[data>cthresh] = cthresh

	smoothed = calc_gd_fwhm(indices_masked.astype(np.int32, order='c'), dist_masked.astype(np.float32, order='c'), data[:,0,0].astype(np.float32, order = "c"), sigma)

	outdata = np.zeros_like(data)
	outdata[cortex_index,0,0] = smoothed[cortex_index]
	nib.save(nib.freesurfer.mghformat.MGHImage(outdata,img.affine),outname) 

def test_draw_samples_non_mock(self, plot=False):
        # Also make sure the non-mock sampler works
        dtype = np.float32
        num_samples = 100000

        mean = np.array([0.5, 0])
        covariance = np.array([[2, 0.5], [0.5, 2]])

        rv_shape = (2,)

        mean_mx = add_sample_dimension(mx.nd, mx.nd.array(mean, dtype=dtype))
        covariance_mx = add_sample_dimension(mx.nd, mx.nd.array(covariance, dtype=dtype))

        rand_gen = None
        var = MultivariateNormal.define_variable(shape=rv_shape, rand_gen=rand_gen, dtype=dtype).factor
        variables = {var.mean.uuid: mean_mx, var.covariance.uuid: covariance_mx}
        rv_samples_rt = var.draw_samples(F=mx.nd, variables=variables, num_samples=num_samples)

        assert array_has_samples(mx.nd, rv_samples_rt)
        assert get_num_samples(mx.nd, rv_samples_rt) == num_samples
        assert rv_samples_rt.dtype == dtype

        from scipy.stats import multivariate_normal
        if plot:
            plot_bivariate(samples=rv_samples_rt, dist=multivariate_normal, mean=mean, cov=covariance)

        # mean_est, scale_est = multivariate_normal.fit(rv_samples_rt.asnumpy().ravel())
        mean_est = np.mean(rv_samples_rt.asnumpy(), axis=0)
        cov_est = np.cov(rv_samples_rt.asnumpy(), rowvar=False)
        mean_tol = 1e-1
        covariance_tol = 1e-1
        assert np.allclose(mean, mean_est, atol=mean_tol)
        assert np.allclose(covariance, cov_est, covariance_tol) 

def test_draw_samples_non_mock(self, plot=False):
        # Also make sure the non-mock sampler works
        dtype = np.float32
        num_samples = 100000

        mean = np.array([0.5])
        variance = np.array([2])

        rv_shape = (1,)

        mean_mx = add_sample_dimension(mx.nd, mx.nd.array(mean, dtype=dtype))
        variance_mx = add_sample_dimension(mx.nd, mx.nd.array(variance, dtype=dtype))

        rand_gen = None
        var = Normal.define_variable(shape=rv_shape, rand_gen=rand_gen, dtype=dtype).factor
        variables = {var.mean.uuid: mean_mx, var.variance.uuid: variance_mx}
        rv_samples_rt = var.draw_samples(F=mx.nd, variables=variables, num_samples=num_samples)

        assert array_has_samples(mx.nd, rv_samples_rt)
        assert get_num_samples(mx.nd, rv_samples_rt) == num_samples
        assert rv_samples_rt.dtype == dtype

        from scipy.stats import norm
        if plot:
            plot_univariate(samples=rv_samples_rt, dist=norm, loc=mean[0], scale=np.sqrt(variance[0]))

        mean_est, scale_est = norm.fit(rv_samples_rt.asnumpy().ravel())
        mean_tol = 1e-1
        variance_tol = 1e-1
        assert np.abs(mean[0] - mean_est) < mean_tol
        assert np.abs(variance[0] - scale_est ** 2) < variance_tol 

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 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_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 visualize_hyper_galaxy(self, fit, hyper_fit, galaxy_image, contribution_map_in):
        hyper_plots.subplot_fit_hyper_galaxy(
            fit=fit,
            hyper_fit=hyper_fit,
            galaxy_image=galaxy_image,
            contribution_map_in=contribution_map_in,
            include=self.include,
            sub_plotter=self.sub_plotter,
        ) 

def visualize_fit_in_fits(self, fit):

        fits_plotter = self.plotter.plotter_with_new_output(
            path=self.plotter.output.path + "fit_imaging/fits/", format="fits"
        )

        fit_imaging_plots.individuals(
            fit=fit,
            plot_image=True,
            plot_noise_map=True,
            plot_signal_to_noise_map=True,
            plot_model_image=True,
            plot_residual_map=True,
            plot_normalized_residual_map=True,
            plot_chi_squared_map=True,
            plot_subtracted_images_of_planes=True,
            plot_model_images_of_planes=True,
            plot_plane_images_of_planes=True,
            include=self.include,
            plotter=fits_plotter,
        )

        if fit.inversion is not None:

            fits_plotter = self.plotter.plotter_with_new_output(
                path=self.plotter.output.path + "inversion/fits/", format="fits"
            )

            inversion_plots.individuals(
                inversion=fit.inversion,
                plot_reconstructed_image=True,
                plot_interpolated_reconstruction=True,
                plot_interpolated_errors=True,
                include=self.include,
                plotter=fits_plotter,
            ) 

def plot_fit_individuals(self, fit):

        fit_galaxy_plots.individuals(
            fit=fit,
            plot_image=self.plot_galaxy_fit_image,
            plot_noise_map=self.plot_galaxy_fit_noise_map,
            plot_model_image=self.plot_galaxy_fit_model_image,
            plot_residual_map=self.plot_galaxy_fit_residual_map,
            plot_chi_squared_map=self.plot_galaxy_fit_chi_squared_map,
            include=self.include,
            plotter=self.plotter,
        ) 

def plot_galaxy_fit_subplot(self, fit):
        if self.plot_subplot_galaxy_fit:
            fit_galaxy_plots.subplot_fit_galaxy(
                fit=fit, include=self.include, sub_plotter=self.sub_plotter
            ) 

def fit(self, X, y=None):
        """Fit the model according to the given training data.

        Parameters
        ----------
        X : array-like of shape (n_samples, n_features)
            Training data.

        y : ignored

        Returns
        -------
        self : object
            Return self.
        """

        self._check_params()

        X                     = self._check_array(X, estimator=self)

        self._fit(X)

        self.classes_         = np.array([NEG_LABEL, POS_LABEL])
        _, self.n_features_   = X.shape
        self.anomaly_score_   = self._anomaly_score(X)
        self.threshold_       = self._get_threshold()
        self.contamination_   = self._get_contamination()
        self.random_variable_ = self._get_random_variable()

        return self 

def _get_random_variable(self):
        """Get the RV object according to the derived anomaly scores."""

        loc, scale = norm.fit(self.anomaly_score_)

        return norm(loc=loc, scale=scale) 

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])
            try:
                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
            except IndexError:
                pass 

def get_stroke_properties(self, stroke_widths):
        if len(stroke_widths) == 0:
            return (0, 0, 0, 0, 0, 0)
        try:
            most_probable_stroke_width = mode(stroke_widths, axis=None)[0][0]
            most_probable_stroke_width_count = mode(stroke_widths, axis=None)[1][0]
        except IndexError:
            most_probable_stroke_width = 0
            most_probable_stroke_width_count = 0
        try:
            mean, std = norm.fit(stroke_widths)
            x_min, x_max = int(min(stroke_widths)), int(max(stroke_widths))
        except ValueError:
            mean, std, x_min, x_max = 0, 0, 0, 0
        return most_probable_stroke_width, most_probable_stroke_width_count, mean, std, x_min, x_max 

def fit_t(ret):
    parm = t.fit(ret)
    nu, mu_t, sig_t = parm
    nu = np.round(nu)
    return mu_t, sig_t, nu 

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 = int(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 fit_normal(ret):
    mu_norm, sig_norm = norm.fit(ret)
    return mu_norm, sig_norm