# -*- coding: utf-8 -*-
import logging
import numpy as np
from scipy.optimize import curve_fit
logger = logging.getLogger(__name__)
[docs]def uncertainties(pder, resid, unc, freep_name, plot=False):
binsize = np.median(np.abs(resid))
nbins = (np.max(resid) - np.min(resid)) / binsize + 1
nbins = int(nbins)
h, x = np.histogram(resid, bins=nbins)
x = x[:-1] + np.diff(x) / 2
_, a = gaussfit(x, h, nterms="constant")
sigma = a[2]
# 6 sigma threshold
resid_threshold = 6 * sigma
nfree = len(freep_name)
freep_unc = np.zeros(nfree)
freep_med = np.zeros(nfree)
freep_msig = np.zeros(nfree)
freep_psig = np.zeros(nfree)
for ifree in range(nfree):
gradlim = np.median(np.abs(pder[:, ifree]))
i = np.where(
(np.abs(pder[:, ifree]) > gradlim) & (np.abs(resid) < resid_threshold)
)[0]
ni = len(i)
if ni > 21:
ii = np.argsort(resid[i] / pder[i, ifree])
# Sort pixels according to the change of the ifree
# parameter needed to match the observations
ch_x = resid[i[ii]] / pder[i[ii], ifree]
# Weights of the individual pixels also sorted
ch_y = np.abs(pder[i[ii], ifree]) / unc[i[ii]]
# Cumulative weights
ch_y = np.cumsum(ch_y)
# Normalized cumulative weights
ch_y = ch_y / ch_y[-1]
# value of distribution at -sigma
freep_msig[ifree] = np.interp(0.5 - 0.6827 / 2, ch_x, ch_y)
# value of distribution at +sigma
freep_psig[ifree] = np.interp(0.5 + 0.6827 / 2, ch_x, ch_y)
# Mean sigma
freep_unc[ifree] = (freep_psig[ifree] - freep_msig[ifree]) * 0.5
# Median change
freep_med[ifree] = np.median(ch_y)
print(
"%s : %.3f %.3f %.3f %.3f"
% (
freep_name[ifree],
freep_med[ifree],
freep_msig[ifree],
freep_psig[ifree],
freep_unc[ifree],
)
)
if plot:
# TODO plot
pass
return freep_med, freep_msig, freep_psig, freep_unc
[docs]def gaussfit(x, y, nterms="none"):
"""
Fit a simple gaussian to data
gauss(x, a, mu, sigma) = a * exp(-z**2/2)
with z = (x - mu) / sigma
Parameters
----------
x : array(float)
x values
y : array(float)
y values
Returns
-------
gauss(x), parameters
fitted values for x, fit paramters (a, mu, sigma)
"""
z = lambda x, A1, A2: (x - A1) / A2
if nterms in ["none", 3]:
gauss = lambda x, A0, A1, A2: A0 * np.exp(-z(x, A1, A2) ** 2 / 2)
p0 = [max(y), 1, 1]
elif nterms in ["constant", 4]:
gauss = lambda x, A0, A1, A2, A3: A0 * np.exp(-z(x, A1, A2) ** 2 / 2) + A3
p0 = [max(y), 1, 1, 0]
elif nterms in ["linear", 5]:
gauss = (
lambda x, A0, A1, A2, A3, A4: A0 * np.exp(-z(x, A1, A2) ** 2 / 2)
+ A3
+ A4 * x
)
p0 = [max(y), 1, 1, 0, 0]
elif nterms in ["quadratic", 6]:
gauss = (
lambda x, A0, A1, A2, A3, A4, A5: A0 * np.exp(-z(x, A1, A2) ** 2 / 2)
+ A3
+ A4 * x
+ A5 * x ** 2
)
p0 = [max(y), 1, 1, 0, 0, 0]
popt, _ = curve_fit(gauss, x, y, p0=p0)
return gauss(x, *popt), popt