Skip to content
Snippets Groups Projects
Commit 2cc12ed0 authored by Julian Lenz's avatar Julian Lenz
Browse files

Finally fixed plotting.py

parents 19271a23 3d5b776e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 100 additions and 298 deletions
plotting.py
+ 100
− 298
View file @ 2cc12ed0
......@@ -10,15 +10,16 @@ import numpy as np
import matplotlib.cm as cm
from matplotlib.pyplot import figure
from matplotlib import rcParams
import matplotlib.pyplot as plt
from numbers import Number
from .styling import Fig, tight, newFig, ax, colors, names
from .aux import iterTest, gaussian
from .styling import Fig, tight, ax, colors, names
from .aux import iterTest
def errorbar3D(ax, x, y, z, xerr=None, yerr=None, zerr=None, **kwargs):
"""
Plot 3D errorbars.
This is a simple function to plot errorbars in 3D which is not natively
supported by matplotlib. It just plots the data points via the scatter
method and adds a short line corresponding to the error in every direction
......@@ -61,7 +62,6 @@ def errorbar3D(ax, x, y, z, xerr=None, yerr=None, zerr=None, **kwargs):
im:
Return value of ax.scatter(x,y,z,...)
"""
cmap = kwargs.pop("cmap", None)
if cmap is not None and type(cmap) == str:
cmap = cm.get_cmap(cmap)
......@@ -148,17 +148,69 @@ def errorbar3D(ax, x, y, z, xerr=None, yerr=None, zerr=None, **kwargs):
class ErrorSleeve:
"""
Compute the error sleeve, i.e. confidence band for a function.
This class takes a function, its estimated (fitted) parameters and their covariance
matrix and computes a confidence band from it. This is done via bootstrapping under
the assumption of normally distributed parameters (according to their covariance
matrix). Bounds can be given that the parameters should obey. However, due to the
assumption of normal distribution this just redraws parameters that violate these
bounds (which is far from fool-proof). Use this feature with care.
--------------------------------------------------
Parameters:
func - function:
The function that was fitted. It is assumed be of signature:
func(x,*p).
p - iterable of parameters:
Provides the parameters that were found in the fit (usually) via
p,pe=scipy.optimize.curve_fit(func,x,y,...).
pe - np.array (1D or 2D):
If 2D, it is interpreted as the covariance matrix as returned from
the call:
p,pe=scipy.optimize.curve_fit(func,x,y,...).
If 1D, it is interpreted as the standard deviations of the
parameters 'p' and is converted into a covariance matrix via:
pe=np.diag(pe**2).
N - int:
Number of bootstrap samples at every point.
Default:1000
part - int:
Percentage of the confidence interval.
Default: 65, i.e. 65% confidence interval
bounds - array-like:
Bounds for the parameters.
"""
safetyMax: int = 100
def __init__(self, func, p, pe, part=0.65, N=1000, bounds=None):
def __init__(self, func, p, pe, part=65, N=1000, bounds=None):
self.func = func
self.p = p
self.pe = pe
if len(np.shape(self.pe)) == 1:
self.pe = np.diag(self.pe ** 2)
self.N = N
self.part = part
self.bounds = np.inf * np.ones((2,) + np.shape(self.p))
self.bounds[0] *= -1.0
self.errorFunction = self.makeErrorFunction()
if bounds is None:
self.bounds = np.inf * np.ones((2,) + np.shape(self.p))
self.bounds[0] *= -1.0
else:
self.bounds = np.array(bounds)
self.errorFunction = self._makeErrorFunction()
def _normalizeInput(self, X):
if not iterTest(X):
......@@ -184,319 +236,69 @@ class ErrorSleeve:
safetyCount += 1
if safetyCount == self.safetyMax:
raise Exception("Something went wrong here.")
print(ptmp.shape)
return ptmp
def makeErrorFunction(self):
def _makeErrorFunction(self):
def errorsBootstrap(x):
x = self._normalizeInput(x)
ySample = np.array([self.func(x, *p) for p in self._drawParameters()])
print(ySample.shape)
return np.vstack(
[
np.percentile(ySample, (1 - self.part) / 2, axis=1),
np.percentile(ySample, self.part + (1 - self.part) / 2, axis=1),
]
).T
return (
np.vstack(
[
np.percentile(ySample, (100 - self.part) / 2, axis=0),
np.percentile(
ySample, self.part + (100 - self.part) / 2, axis=0
),
]
).T
- self.func(x, *self.p)[:, np.newaxis]
)
return errorsBootstrap
def plot(self, ax, x, color=None, alpha=0.3, label=None):
"""
Plot the error sleeve.
This method plots the error sleeve onto the given axes `ax`. If `ax` is an
instance of FancyAxes, it takes the current plotting color as a default.
Parameters
----------
ax : matplotlib axes
Axes to plot onto.
x : array-like
Values on which to evaluate the errors.
color : matplotlib color, optional
Color of the plot. The default is the current style's color if `ax` is a
FancyAxes and else 'r'.
alpha : float between 0 and 1, optional
Same as for other matplotlib plot functions. The default is 0.3.
label : str, optional
Label of the error sleeve. The default is None.
Returns
-------
matplotlib.collections.PolyCollection
Return value of ax.fill_between(...).
"""
if color is None:
try:
color = ax.getCurrentStyle("plot")["color"]
except AttributeError:
# ax is not a FancyAxes
color = "r"
e = self.errorFunction(x)
e = self.errorFunction(x) + self.func(x, *self.p)[:, np.newaxis]
return ax.fill_between(
x, e[:, 0], e[:, 1], alpha=alpha, color=color, linewidth=0, label=label
)
def values(self, x):
"""Return the errors at the positions x."""
return self.errorFunction(x)
def errorSleeve(func, p, pe, **kwargs):
"""
NOT YET FULLY DOCUMENTED!
Computes a confidence interval of a function given estimated parameters and
their covariance matrix (resp. variances).
The method 'bootstrap' draws just brute force for every parameter in the
range of the given variances and samples this. It can take a lot of time
to get smooth confidence curves.
The method 'analytic' calculates the standard deviation from the assumption
of multi-variately normal distributed parameters by an integration.
--------------------------------------------------
Parameters:
func - function:
The function that was fitted. It is assumed be of signature:
func(x,*p).
p - iterable of parameters:
Provides the parameters that were found in the fit (usually) via
p,pe=scipy.optimize.curve_fit(func,x,y,...).
pe - np.array (1D or 2D):
If 2D, it is interpreted as the covariance matrix as returned from
the call:
p,pe=scipy.optimize.curve_fit(func,x,y,...).
If 1D, it is interpreted as the standard deviations of the
parameters 'p' and is converted into a covariance matrix via:
pe=np.diag(pe**2).
method - str 'bootstrap':
Chooses the method to use. 'bootstrap' randomly samples the
distribution at every point with parameters drawn from a multivariate
normal distribution. It then returns the left and right edge such
that all values out of this region are below 'rel' * maximum. At
the moment, only 'bootstrap' is implemented.
Default: 'bootstrap'
rel - float (0 <= rel <= 1):
Relative height at which the width is calculated, i.e. the
confidence interval returned ensures that away from it the
distribution is below 'rel' * maximum.
N - int:
Number of bootstrap samples at every point.
Default:1000
--------------------------------------------------
Return Values:
errors - function:
A function that calculates the lower and upper bound of the
confidence interval:
[[lower(x1),upper(x1)],[lower(x2),upper(x2)],...]=errors([x1,x2,...]).
more - list:
In specific scenarios additional output is helpful. This is returned
inside of this list. In the following scenarios it has the following
elements:
- more==[] if nothing was further specialized.
- If a grid was used, the second last element is always the
grid and the last element is always errors(grid).
- If anything was plotted, the first element is always the axis
in use followed by output from the plotting methods, i.e. with
plot==2 one would get [ax,im,line,e], where 'im' is the
return value of ax.fill_between(...) and line is the return
value of ax.plot(...).
--------------------------------------------------
"""
method = kwargs.pop("method", "bootstrap")
label = kwargs.pop("label", None)
grid = kwargs.pop("grid", None)
plot = kwargs.pop(
"plot", 0 if grid is None else 1
) # 0 no plotting, 1 plot just sleeve, 2 also function
if plot not in (0, 1, 2):
raise Exception("plot parameter not understood.")
ax = kwargs.pop("ax", None)
xlim = kwargs.pop("xlim", None)
ngrid = kwargs.pop("ngrid", 300)
alpha = kwargs.pop("alpha", 0.3)
color = kwargs.pop("color", "r")
bounds = kwargs.pop("bounds", None)
rel = kwargs.pop("rel", 0.5)
if method == "bootstrap":
N = kwargs.pop("N", 1000)
Bins = kwargs.pop("bins", int(np.max([np.floor(N / 100), 10])))
smear = kwargs.pop("smear", 0)
if smear % 2 == 1:
smear = smear + 1
if smear > 0:
weight = kwargs.pop("smearWeight", "gaussian")
if type(weight) == str:
if weight == "gaussian":
sigma = np.sqrt(-(smear ** 2) / (4 * np.log(0.01)))
weights = gaussian(
np.linspace(-smear / 2, smear / 2, smear), 1, sigma, 0
)
weights *= 1 / np.sum(weights)
weight = lambda x: gaussian(x, 1, sigma, 0)
elif weight == "uniform":
weights = np.ones(smear) / smear
weight = lambda x: 1
else:
raise Exception("No valid smearing weight.")
else:
try:
weights = weight(np.linspace(-smear / 2, smear / 2, smear))
except TypeError:
print("The smearing weight you have given is not a string")
print("and probably also not a function.")
raise
else:
weight = lambda x: 1
weights = np.array([1])
safetyMax = kwargs.pop("safetyMax", 100)
if bounds is None:
bounds = np.array([(-np.inf, np.inf) for i in range(len(p))]).T
else:
bounds = np.array(bounds)
if plot > 0:
if ax is None:
f, ax = newFig()
if grid is None:
if xlim is None:
raise Exception("You must either provide grid or xlim to plot.")
else:
grid = np.linspace(xlim[0], xlim[1], ngrid)
if np.shape(pe) != (len(p), len(p)):
# pe are just standard deviations; make them a covariance matrix.
pe = np.diag(pe ** 2)
def errorsBootstrap(X, showHist=False):
"""
showHist - bool:
For debugging only! Shows the histogram at every point.
"""
if not iterTest(X):
X = np.array([X])
ret = np.empty((len(X), 2))
tmp = np.empty(N)
for i, x in enumerate(X):
ptmp = np.random.multivariate_normal(mean=p, cov=pe, size=N)
safetyCount = 0
while (
not np.all(ptmp > bounds[0, :])
or not np.all(ptmp < bounds[1, :])
and safetyCount < safetyMax
):
mask = np.invert(np.all(ptmp > bounds[0, :], axis=1))
ptmp[mask] = np.random.multivariate_normal(
mean=p, cov=pe, size=np.sum(mask)
)
mask = np.invert(np.all(ptmp < bounds[1, :], axis=1))
ptmp[mask] = np.random.multivariate_normal(
mean=p, cov=pe, size=np.sum(mask)
)
safetyCount += 1
if safetyCount == safetyMax:
raise Exception("Something went wrong here.")
for n in range(N):
tmp[n] = func(x, *ptmp[n, :])
if showHist:
plt.figure()
plt.hist(tmp, bins=Bins)
plt.vlines(func(x, *p), *plt.gca().get_ylim())
tmpHist, bins = np.histogram(tmp, bins=Bins)
if not np.any(np.array(np.shape(np.nonzero(tmpHist))) > 1):
argLeft = np.nonzero(tmpHist)[0]
argRight = np.nonzero(tmpHist)[0]
else:
arg = np.where(bins >= func(x, *p))[0][0]
if arg == 0:
argLeft = 0
else:
try:
mask = tmpHist[:arg] > rel * tmpHist[arg]
except IndexError:
mask = tmpHist[:arg] > rel * tmpHist[arg - 1]
if len(mask) == 0 or mask[0]:
argLeft = 0
elif not np.any(mask):
argLeft = arg
else:
argLeft = np.where(mask)[0][0]
if arg == len(tmpHist) - 1:
argRight = len(tmpHist)
else:
try:
mask = tmpHist[arg:][::-1] > rel * tmpHist[arg]
except IndexError:
mask = tmpHist[arg:][::-1] > rel * tmpHist[arg - 1]
if len(mask) == 0 or mask[0]:
argRight = len(tmpHist)
elif not np.any(mask):
argRight = arg + 1
else:
argRight = len(tmpHist) - np.where(mask)[0][0]
ret[i] = [bins[argLeft], bins[argRight]]
if len(X) > smear and smear > 0:
ret[int(smear / 2) : -int(smear / 2), :] = np.sum(
[
(weights * ret[j : j + smear, :].T).T
for j in range(len(ret) - smear)
],
axis=1,
)
for i in range(int(smear / 2)):
ret[i, :] = np.sum(
(
weight(np.arange(-i, int(smear / 2)))
/ np.sum(weight(np.arange(-i, int(smear / 2))))
* ret[: int(smear / 2) + i, :].T
).T,
axis=0,
)
ret[-i - 1, :] = np.sum(
(
weight(np.arange(-i, int(smear / 2)))
/ weight(np.arange(-i, int(smear / 2)))
* ret[-int(smear / 2) - i :, :].T
).T,
axis=0,
)
return ret
errors = errorsBootstrap
else:
raise Exception('Method not understood. Choose from "bootstrap".')
if plot == 0:
if grid is None:
return errors, []
else:
return errors, [grid, errors(grid)]
e = errors(grid)
if plot == 1:
im = ax.fill_between(
grid, e[:, 0], e[:, 1], alpha=alpha, color=color, linewidth=0, label=label
)
return errorsBootstrap, [ax, im, grid, e]
elif plot == 2:
im = ax.fill_between(
grid, e[:, 0], e[:, 1], alpha=alpha, color=color, linewidth=0
)
line = ax.plot(grid, func(grid, *p), color=color, label=label)
return errorsBootstrap, [ax, im, line, grid, e]
else:
raise Exception("plot had an unexpected value.")
return None
class BrokenAxis:
"""
The BrokenAxis class is designed to yield plots that look as if one (or both)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment