Python uncertainties.correlated_values函数代码示例

    def test_correlated_values():
        "Correlated variables."

        u = uncertainties.ufloat((1, 0.1))
        cov = uncertainties.covariance_matrix([u])
        # "1" is used instead of u.nominal_value because
        # u.nominal_value might return a float.  The idea is to force
        # the new variable u2 to be defined through an integer nominal
        # value:
        u2, = uncertainties.correlated_values([1], cov)
        expr = 2 * u2  # Calculations with u2 should be possible, like with u

        ####################

        # Covariances between output and input variables:

        x = ufloat((1, 0.1))
        y = ufloat((2, 0.3))
        z = -3 * x + y

        covs = uncertainties.covariance_matrix([x, y, z])

        # "Inversion" of the covariance matrix: creation of new
        # variables:
        (x_new, y_new, z_new) = uncertainties.correlated_values(
            [x.nominal_value, y.nominal_value, z.nominal_value], covs, tags=["x", "y", "z"]
        )

        # Even the uncertainties should be correctly reconstructed:
        assert matrices_close(numpy.array((x, y, z)), numpy.array((x_new, y_new, z_new)))

        # ... and the covariances too:
        assert matrices_close(numpy.array(covs), numpy.array(uncertainties.covariance_matrix([x_new, y_new, z_new])))

        assert matrices_close(numpy.array([z_new]), numpy.array([-3 * x_new + y_new]))

        ####################

        # ... as well as functional relations:

        u = ufloat((1, 0.05))
        v = ufloat((10, 0.1))
        sum_value = u + 2 * v

        # Covariance matrices:
        cov_matrix = uncertainties.covariance_matrix([u, v, sum_value])

        # Correlated variables can be constructed from a covariance matrix, if
        # NumPy is available:
        (u2, v2, sum2) = uncertainties.correlated_values([x.nominal_value for x in [u, v, sum_value]], cov_matrix)

        # matrices_close() is used instead of _numbers_close() because
        # it compares uncertainties too:
        assert matrices_close(numpy.array([0]), numpy.array([sum2 - (u2 + 2 * v2)]))

def auswerten(name, d, n, t, z, V_mol, eps, raw):
    d *= 1e-3
    N = unp.uarray(n/t, np.sqrt(n)/t) - N_u

    if name=="Cu":
        tools.table((raw[0], raw[1], N), ("D/mm", "n", "(N-N_U)/\per\second"), "build/{}.tex".format(name), "Messdaten von {}.".format(name), "tab:daten{}".format(name), split=2, footer=r"$\Delta t = \SI{60}{s}$")#"(N-N_U)/\per\second"
    else:
        tools.table((raw[0], raw[1], raw[2], N), ("D/mm", "n", "\Delta t/s", "(N-N_U)/\per\second"), "build/{}.tex".format(name), "Messdaten von {}.".format(name), "tab:daten{}".format(name), split=2)
    mu = z * const.N_A / V_mol * 2 * np.pi * (const.e**2 / (4 * np.pi * const.epsilon_0 * const.m_e * const.c**2))**2 * ((1+eps)/eps**2 * ((2 * (1+eps))/(1+2*eps) - 1/eps * np.log(1+2*eps)) + 1/(2*eps) * np.log(1+ 2*eps) - (1+ 3*eps)/(1+2*eps)**2)

    params, pcov = curve_fit(fit, d, unp.nominal_values(N), sigma=unp.std_devs(N))
    params_ = unc.correlated_values(params, pcov)
    print("{}: N(0) = {}, µ = {}, µ_com = {}".format(name, params_[0], -params_[1], mu))

    sd = np.linspace(0, .07, 1000)

    valuesp = (fit(sd, *(unp.nominal_values(params_) + 10*unp.std_devs(params_)))).astype(float)
    valuesm = (fit(sd, *(unp.nominal_values(params_) - 10*unp.std_devs(params_)))).astype(float)

    #plt.xlim(0,7)
    plt.xlabel(r"$D/\si{mm}$")
    plt.ylabel(r"$(N-N_U)/\si{\per\second}$")
    plt.plot(1e3*sd, fit(sd, *params), 'b-', label="Fit")
    plt.fill_between(1e3*sd, valuesm, valuesp, facecolor='blue', alpha=0.125, edgecolor='none', label=r'$1\sigma$-Umgebung ($\times 10$)')
    plt.errorbar(1e3*d, unp.nominal_values(N), yerr=unp.std_devs(N), fmt='rx', label="Messdaten")
    plt.legend(loc='best')
    plt.yscale('linear')
    plt.tight_layout(pad=0)
    plt.savefig("build/{}.pdf".format(name))
    plt.yscale('log')
    plt.savefig("build/{}_log.pdf".format(name))
    plt.clf()

    def model_with_uncertainties(self):
        """Best fit model with uncertainties

        The parameters on the model will have the units of the model attribute.
        The covariance matrix passed on initialization must also have these
        units.

        TODO: Add to gammapy.spectrum.models

        This function uses the uncertainties packages as explained here
        https://pythonhosted.org/uncertainties/user_guide.html#use-of-a-covariance-matrix

        Examples
        --------
        TODO
        """
        import uncertainties
        pars = self.model.parameters

        # convert existing parameters to ufloats
        values = [pars[_].value for _ in self.covar_axis]
        ufloats = uncertainties.correlated_values(values, self.covariance)
        upars = dict(zip(self.covar_axis, ufloats))

        # add parameters missing in covariance
        for name in pars:
            upars.setdefault(name, pars[name].value)

        return self.model.__class__(**upars)

def salt2mu_aberr(x1=None,x1err=None,
                  c=None,cerr=None,
                  mb=None,mberr=None,
                  cov_x1_c=None,cov_x1_x0=None,cov_c_x0=None,
                  alpha=None,beta=None,
                  alphaerr=None,betaerr=None,
                  M=None,x0=None,sigint=None,z=None,peczerr=0.0005):
    from uncertainties import ufloat, correlated_values, correlated_values_norm
    alphatmp,betatmp = alpha,beta
    alpha,beta = ufloat(alpha,alphaerr),ufloat(beta,betaerr)

    sf = -2.5/(x0*np.log(10.0))
    cov_mb_c = cov_c_x0*sf
    cov_mb_x1 = cov_x1_x0*sf

    mu_out,muerr_out = np.array([]),np.array([])
    for i in range(len(x1)):

        covmat = np.array([[mberr[i]**2.,cov_mb_x1[i],cov_mb_c[i]],
                           [cov_mb_x1[i],x1err[i]**2.,cov_x1_c[i]],
                           [cov_mb_c[i],cov_x1_c[i],cerr[i]**2.]])
        mb_single,x1_single,c_single = correlated_values([mb[i],x1[i],c[i]],covmat)

        mu = mb_single + x1_single*alpha - beta*c_single + 19.36
        if sigint: mu = mu + ufloat(0,sigint)
        zerr = peczerr*5.0/np.log(10)*(1.0+z[i])/(z[i]*(1.0+z[i]/2.0))

        mu = mu + ufloat(0,np.sqrt(zerr**2. + 0.055**2.*z[i]**2.))
        mu_out,muerr_out = np.append(mu_out,mu.n),np.append(muerr_out,mu.std_dev)

    return(mu_out,muerr_out)

    def get_covariance_matrix(self):
        '''
        Get the covariance matrix from all contributions
        https://root.cern.ch/doc/master/classTUnfoldDensity.html#a7f9335973b3c520e2a4311d2dd6f5579
        '''
        import uncertainties as u
        from numpy import array, matrix, zeros
        from numpy import sqrt as np_sqrt
        if self.unfolded_data is not None:
            # Calculate the covariance from TUnfold
            covariance = asrootpy( 
                self.unfoldObject.GetEmatrixInput('Covariance'))

            # Reformat into a numpy matrix
            zs = list(covariance.z())
            cov_matrix = matrix(zs)

            # Just the unfolded number of events         
            inputs = hist_to_value_error_tuplelist(self.unfolded_data)
            nominal_values = [i[0] for i in inputs]         
            # # Unfolded number of events in each bin
            # # With correlation between uncertainties
            values_correlated = u.correlated_values( nominal_values, cov_matrix.tolist() )
            corr_matrix = matrix(u.correlation_matrix(values_correlated) )

            return cov_matrix, corr_matrix
        else:
            print("Data has not been unfolded. Cannot return unfolding covariance matrix")
        return

    def model_with_uncertainties(self):
        """Best fit model with uncertainties

        The parameters on the model will be in units ``keV``, ``cm``, and
        ``s``. Thus, when evaluating the model energies have to be passed in
        keV and the resulting flux will be in ``cm-2 s-1 keV-1``. The
        covariance matrix passed on initialization must also have these units.

        TODO: This is due to sherpa units, make more flexible
        TODO: Add to gammapy.spectrum.models

        This function uses the uncertainties packages as explained here
        https://pythonhosted.org/uncertainties/user_guide.html#use-of-a-covariance-matrix

        Examples
        --------
        TODO
        """
        import uncertainties
        pars = self.model.parameters

        # convert existing parameters to ufloats
        values = [pars[_].value for _ in self.covar_axis]
        ufloats = uncertainties.correlated_values(values, self.covariance)
        upars = dict(zip(self.covar_axis, ufloats))

        # add parameters missing in covariance
        for name in pars:
            upars.setdefault(name, pars[name].value)

        return self.model.__class__(**upars)

def integral():
    x1 = uc.ufloat( 20 /2 *10**(-3)  ,   1*10**(-3))
    x2 = uc.ufloat( 150/2 *10**(-3)  ,   1*10**(-3))
    L =  uc.ufloat( 175   *10**(-3)  ,   1*10**(-3))
    l =  uc.ufloat( 150   *10**(-3)  ,   1*10**(-3))
    I =  uc.ufloat(10 , 0.1)
    N =  3600

    A = N  / (2 *L  * (x2-x1)) 
    factor =  (x2 * (sqrt( x2**2 + ((L+l)/2)**2)- sqrt( x2**2 + ((L-l)/2)**2) ) \
              -x1 * (sqrt( x1**2 + ((L+l)/2)**2)- sqrt( x1**2 + ((L-l)/2)**2) )  \
    + (((L+l)/2)**2* log((x2 + sqrt(((L+l)/2)**2  + x2**2))/ (x1 + sqrt(((L+l)/2)**2  + x1**2)) ) \
    -((L-l)/(4))**2* log((x2 + sqrt(((L-l)/2)**2  + x2**2))/ (x1 + sqrt(((L-l)/2)**2  + x1**2)) )))

    print(factor)
    print(A * factor)

    a = np.load(input_dir +"a_3.npy")
    I = np.load(input_dir + "i_3.npy")
    I_fit = np.linspace(min(I),max(I),100) 

    S_a = 0.2
    S_I = 0.05

    weights = a*0 + 1/S_a 
    p, cov= np.polyfit(I,a, 1, full=False, cov=True, w=weights)

    (p1, p2) = uc.correlated_values([p[0],p[1]], cov)

    print(p1 / (A*factor))
    print(p1/2556/100 *oe *60)

def main():
    data = pd.read_csv(
        'data/messwerte_2.csv',
        skiprows=(1, 2, 3, 4, 5),
    )
    data['T'] = data['T'].apply(const.C2K)

    func = partial(linear, x0=data['t'].mean())

    params, cov = curve_fit(func, data['t'], data['T'])
    a, b = unc.correlated_values(params, cov)
    print('Rate = {} Kelvin per minute'.format(a))

    with open('build/rate.tex', 'w') as f:
        f.write(r'b = \SI{{{a.n:1.2f} +- {a.s:1.2f}}}{{\kelvin\per\minute}}'.format(a=a))

    T_max = 289.95 * u.kelvin
    relaxation_time = (
        (u.boltzmann_constant * T_max**2) /
        (0.526 * u.eV * a * u.kelvin / u.minute) *
        np.exp(-0.526 * u.eV / (u.boltzmann_constant * T_max))
    )
    print(relaxation_time.to(u.second))

    t = np.linspace(0, 60, 2)
    plt.plot(t, func(t, *params), label='Ausgleichsgerade', color='gray')
    plt.plot(data['t'], data['T'], 'x', ms=3, label='Messwerte')
    plt.xlabel(r'$t \mathbin{/} \si{\minute}$')
    plt.ylabel(r'$T \mathbin{/} \si{\kelvin}$')
    plt.legend(loc='lower right')

    plt.tight_layout(pad=0)
    plt.savefig('build/rate.pdf')

def ucurve_fit(f, x, y, **kwargs):
    if np.any(unp.std_devs(y) == 0):
        sigma = None
    else:
        sigma = unp.std_devs(y)

    popt, pcov = scipy.optimize.curve_fit(f, x, unp.nominal_values(y), sigma=sigma, **kwargs)

    return unc.correlated_values(popt, pcov)

def CdTe_Am():
    detector_name = "CdTe"
    sample_name = "Am"
    suff = "_" + sample_name + "_" + detector_name 
    npy_file = npy_dir + "spectra" + suff + ".npy"
    n = np.load(npy_file)[:500]
    s_n = np.sqrt(n)
    x = np.arange(len(n))
    
    # only fit
    p0 = [250000, 250, 25]
    red = (x > 235) * (x < 270)
    n_red = n[red]
    x_red = x[red]
    s_n_red = s_n[red]
    fit, cov_fit = curve_fit(gauss, x_red, n_red, p0=p0, sigma=s_n_red)
    fit  = np.abs(fit)
    fit_corr = uc.correlated_values(fit, cov_fit)
    s_fit = un.std_devs(fit_corr)
    # chi square
    n_d = 4
    chi2 = np.sum(((gauss(x_red, *fit) - n_red) / s_n_red) ** 2 )
    chi2_test = chi2/n_d
    fit_r, s_fit_r = err_round(fit, cov_fit)
    fit_both = np.array([fit_r, s_fit_r])
    fit_both = np.reshape(fit_both.T, np.size(fit_both))
    As[2], mus[2], sigmas[2] = fit
    s_As[2], s_mus[2], s_sigmas[2] = un.std_devs(fit_corr)

    all = np.concatenate((fit_both, [chi2_test]), axis=0)
    def plot_two_gauss():
        fig1, ax1 = plt.subplots(1, 1)
        if not save_fig:
            fig1.suptitle("Detector: " + detector_name + "; sample: " + sample_name)
        plot1, = ax1.plot(x, n, '.', alpha=0.3)   # histo plot
        ax1.errorbar(x, n, yerr=s_n, fmt=',', alpha=0.99, c=plot1.get_color(), errorevery=10) # errors of t are not changed!
        ax1.plot(x_red, gauss(x_red, *fit))
        ax1.set_xlabel("channel")
        ax1.set_ylabel("counts")
        textstr = 'Results of fit:\n\
                \\begin{eqnarray*}\
                A     &=& (%.0f \pm %.0f) \\\\ \
                \mu   &=& (%.1f \pm %.1f) \\\\ \
                \sigma&=& (%.1f \pm %.1f) \\\\ \
                \chi^2 / n_d &=& %.1f\
                \end{eqnarray*}'%tuple(all)
        ax1.text(0.65, 0.95, textstr, transform=ax1.transAxes, va='top', bbox=props)
        if show_fig:
            fig1.show()
        if save_fig:
            file_name = "detector" + suff
            fig1.savefig(fig_dir + file_name + ".pdf")
            fig1.savefig(fig_dir + file_name + ".png")
        return 0
    plot_two_gauss()
    return fit, s_fit

def poly_fit(x, y_e, x_range, p0):
    x_min, x_max = x_range
    mask = (x > x_min) * (x < x_max)
    x_fit = x[mask]
    y_fit = un.nominal_values(y_e[mask])
    y_sigma = np.sqrt(un.std_devs(y_e[mask]))
    coeff, cov = curve_fit(poly, x_fit, y_fit, p0=p0, 
                           sigma=y_sigma, absolute_sigma=True)
    c = uc.correlated_values(coeff, cov)
    return c

def parabola_fit(points):
    dims = points[0][0].shape[0]
    
    x = np.array([p[0] for p in points])
    f = np.array([p[1] for p in points])

    A = build_design_matrix(x, f)
    B = build_design_vector(f)[:,np.newaxis] # make column vector

    # Compute best-fit parabola coefficients using a singular value
    # decomposition.
    U, w, V = np.linalg.svd(A, full_matrices=False)
    V = V.T # Flip to convention used by Numerical Recipies
    inv_w = 1.0/w
    inv_w[np.abs(w) < 1e-6] = 0.0
    # Numpy version of Eq 15.4.17 from Numerical Recipies (C edition)
    coeffs = np.zeros(A.shape[1])
    for i in xrange(len(coeffs)):
        coeffs += (np.dot(U[:,i], B[:,0]) * inv_w[i]) * V[:,i]

    # Chi2 and probability for best fit and quadratic coefficents
    chi2_terms = np.dot(A, coeffs[:,np.newaxis]) - B
    chi2 = (chi2_terms**2).sum()
    ndf = len(points) - (1 + dims + dims * (dims + 1) / 2)
    prob = ROOT.TMath.Prob(chi2, ndf)

    # Covariance is from Eq 15.4.20
    covariance = np.dot(V*inv_w**2, V.T)

    # Pack the coefficients into ufloats
    ufloat_coeffs = uncertainties.correlated_values(coeffs, covariance.tolist())

    # Separate coefficients into a, b, and c
    a = ufloat_coeffs[0]
    b = ufloat_coeffs[1:dims+1]
    c = np.zeros(shape=(dims,dims), dtype=object)
    index = dims + 1
    for i in xrange(dims):
        for j in xrange(i, dims):
            c[i,j] = ufloat_coeffs[index]
            c[j,i] = ufloat_coeffs[index]
            if j != i:
                # We combined the redundant off-diagonal parts of c
                # matrix, but now we have to divide by two to 
                # avoid double counting when c is used later
                c[i,j] /= 2.0
                c[j,i] /= 2.0
            index += 1

    return a, np.array(b), c, chi2, prob    

def two_bw_fit(x, y_e, x_range, p0, fit=True):
    x_min, x_max = x_range
    mask = (x > x_min) * (x < x_max)
    x_fit = x[mask]
    y_fit = un.nominal_values(y_e[mask])
    y_sigma = np.sqrt(un.std_devs(y_e[mask]))
    if fit:
        coeff, cov = curve_fit(two_breit_wigner, x_fit, y_fit, p0=p0, 
                                  sigma=y_sigma, absolute_sigma=True)
        c = uc.correlated_values(coeff, cov)
        fit_peak = two_breit_wigner(x_fit, *coeff)
    else:
        fit_peak = two_breit_wigner(x_fit, *p0)
        c = un.uarray(p0, [0, 0, 0, 0])
    return x_fit, fit_peak, c

def verdet_constant():
        a = np.load(input_dir +"a_3.npy")
        I = np.load(input_dir + "i_3.npy")
        I_fit = np.linspace(min(I),max(I),100) 

        S_a = 0.2
        S_I = 0.05

        weights = a*0 + 1/S_a 
        p, cov= np.polyfit(I,a, 1, full=False, cov=True, w=weights)

        (p1, p2) = uc.correlated_values([p[0],p[1]], cov)

        print(p1/2556)
        print(p1/2556/100 *oe *60)

    def _ufloats(self):
        """Return dict of ufloats with covariance."""
        from uncertainties import correlated_values

        values = [_.value for _ in self.parameters]

        try:
            # convert existing parameters to ufloats
            uarray = correlated_values(values, self.covariance)
        except np.linalg.LinAlgError:
            raise ValueError("Covariance matrix not set.")

        upars = {}
        for par, upar in zip(self.parameters, uarray):
            upars[par.name] = upar

        return upars