"Definition of a univariate parameter."
import equadratures.plot as plot
import numpy as np
import scipy as sc
[docs]
class ParentParameter(object):
""" This class defines a univariate parameter.
Parameters
----------
lower : float, optional
Lower bound for the parameter.
upper : float, optional
Upper bound for the parameter.
order : int, optional
Order of the parameter.
param_type : str, optional
The type of distribution that characterizes the parameter (see [1, 2]). Options include `chebyshev (arcsine) <https://en.wikipedia.org/wiki/Arcsine_distribution>`_, `gaussian <https://en.wikipedia.org/wiki/Normal_distribution>`_,
`truncated-gaussian <https://en.wikipedia.org/wiki/Truncated_normal_distribution>`_, `beta <https://en.wikipedia.org/wiki/Beta_distribution>`_,
`cauchy <https://en.wikipedia.org/wiki/Cauchy_distribution>`_, `exponential <https://en.wikipedia.org/wiki/Exponential_distribution>`_,
`uniform <https://en.wikipedia.org/wiki/Uniform_distribution_(continuous)>`_, `triangular <https://en.wikipedia.org/wiki/Triangular_distribution>`_, `gamma <https://en.wikipedia.org/wiki/Gamma_distribution>`_,
`weibull <https://en.wikipedia.org/wiki/Weibull_distribution>`_, `rayleigh <https://en.wikipedia.org/wiki/Rayleigh_distribution>`_,
`pareto <https://en.wikipedia.org/wiki/Pareto_distribution>`_, `lognormal <https://en.wikipedia.org/wiki/Log-normal_distribution>`_,
`students-t <https://en.wikipedia.org/wiki/Student%27s_t-distribution>`_, `logistic <https://en.wikipedia.org/wiki/Log-normal_distribution>`_,
`gumbel <https://en.wikipedia.org/wiki/Gumbel_distribution>`_, `chi <https://en.wikipedia.org/wiki/Chi_distribution>`_ and `chi-squared <https://en.wikipedia.org/wiki/Chi-squared_distribution>`_.
If no string is provided, a ``uniform`` distribution is assumed. Data-driven and custom analytical parameters can also be constructed by setting this option to ``data`` and ``analytical`` and providing a **weight_function** (see examples).
shape_parameter_A : float, optional
Most of the aforementioned distributions are characterized by two shape parameters. For instance, in the case of a ``gaussian`` (or ``truncated-gaussian``), this represents the mean. In the case of a beta distribution this represents the alpha value. For a ``uniform`` distribution this input is not required.
shape_parameter_B : float, optional
This is the second shape parameter that characterizes the distribution selected. In the case of a ``gaussian`` or ``truncated-gaussian``, this is the variance.
data : numpy.ndarray, optional
A data-set with shape (number_of_data_points, 2), where the first column comprises of parameter values, while the second column corresponds to the data observations. This input should only be used with the ``Analytical`` distribution.
endpoints : str, optional
If set to ``both``, then the quadrature points and weights will have end-points, based on Gauss-Lobatto quadrature rules. If set to ``upper`` or ``lower`` a Gauss-Radau rule is used to compute one end-point at either the upper or lower bound.
weight_function: Weight, optional
An instance of Weight, which contains a bespoke analytical or data-driven weight (probability density) function.
Examples
--------
A uniform parameter
>>> param = eq.Parameter(distribution='uniform', lower=-2, upper=2., order=3)
A beta parameter
>>> param = eq.Parameter(distribution='beta', lower=-2., upper=15., order=4,
>>> shape_parameter_A=3.2, shape_parameter_B=1.7)
A data-driven parameter
>>> pdf = eq.Weight( stats.gaussian_kde(data, bw_method='silverman'),
>>> support=[-3, 3.2])
>>> param = eq.Parameter(distribution='analytical',
>>> weight_function=pdf, order=2)
References
----------
1. Xiu, D., Karniadakis, G. E., (2002) The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations. SIAM Journal on Scientific Computing, 24(2), `Paper <https://epubs.siam.org/doi/abs/10.1137/S1064827501387826?journalCode=sjoce3>`__
2. Gautschi, W., (1985) Orthogonal Polynomials-Constructive Theory and Applications. Journal of Computational and Applied Mathematics 12 (1985), pp. 61-76. `Paper <https://www.sciencedirect.com/science/article/pii/037704278590007X>`__
"""
def __init__(self, distribution, order, variable, endpoints):
self.distribution = distribution
self.variable = variable
self.order = order
self.endpoints = endpoints
self.weight_function = weight_function
self._set_distribution()
self._set_bounds()
self._set_moments()
if self.endpoints is not None:
if (self.distribution.bounds[0] == -np.inf) and (self.distribution.bounds[1] == np.inf) and (self.endpoints.lower() == 'both'):
raise(ValueError, 'Parameter: The lower bound for your distribution is -infinity and the upper bound is infinity. Furthermore, you have selected the to have both endpoints. These options are incompatible!')
if (self.distribution.bounds[0] == -np.inf) and (self.endpoints.lower() == 'lower'):
raise(ValueError, 'Parameter: The lower bound for your distribution is -infinity and you have selected the lower bound option in the endpoints. These options are incompatible!')
if (self.distribution.bounds[1] == np.inf) and (self.endpoints.lower() == 'upper'):
raise(ValueError, 'Parameter: The upper bound for your distribution is infinity and you have selected the upper bound option in the endpoints. These options are incompatible!')
def _set_distribution(self):
""" Private function that sets the distribution. """
if self.name.lower() == 'gaussian' or self.name.lower() == 'normal':
self.distribution = Gaussian(self.shape_parameter_A, self.shape_parameter_B)
elif self.name.lower() == 'uniform':
self.distribution = Uniform(self.lower, self.upper)
elif self.name.lower() == 'triangular':
self.distribution = Triangular(self.lower, self.upper, self.shape_parameter_A)
elif self.name.lower() == 'analytical' or self.name.lower() == 'data':
self.distribution = Analytical(self.weight_function)
elif self.name.lower() == 'beta':
self.distribution = Beta(self.lower, self.upper, self.shape_parameter_A, self.shape_parameter_B)
elif self.name.lower() == 'truncated-gaussian':
self.distribution = TruncatedGaussian(self.shape_parameter_A, self.shape_parameter_B, self.lower, self.upper)
elif self.name.lower() == 'cauchy':
self.distribution = Cauchy(self.shape_parameter_A, self.shape_parameter_B)
elif self.name.lower() == 'exponential':
self.distribution = Exponential(self.shape_parameter_A)
elif self.name.lower() == 'gamma':
self.distribution = Gamma(self.shape_parameter_A, self.shape_parameter_B)
elif self.name.lower() == 'weibull':
self.distribution = Weibull(self.shape_parameter_A, self.shape_parameter_B)
elif self.name.lower() == 'arcsine' or self.name.lower() == 'chebyshev':
self.distribution = Chebyshev(self.lower, self.upper)
elif self.name.lower() == 'rayleigh':
self.distribution = Rayleigh(self.shape_parameter_A)
elif self.name.lower() == 'chi-squared':
self.distribution = Chisquared(self.shape_parameter_A)
elif self.name.lower() == 'chi':
self.distribution = Chi(self.shape_parameter_A)
elif self.name.lower() == 'pareto':
self.distribution = Pareto(self.shape_parameter_A)
elif self.name.lower() == 'gumbel':
self.distribution = Gumbel(self.shape_parameter_A, self.shape_parameter_B)
elif self.name.lower() == 'logistic':
self.distribution = Logistic(self.shape_parameter_A, self.shape_parameter_B)
elif self.name.lower() == 'students-t' or self.name.lower() == 't' or self.name.lower() == 'studentt':
self.distribution = Studentst(self.shape_parameter_A)
elif self.name.lower() == 'lognormal' or self.name.lower() == 'log-normal':
self.distribution = Lognormal(self.shape_parameter_A)
else:
distribution_error()
self.mean = self.distribution.mean
self.variance = self.distribution.variance
[docs]
def plot_orthogonal_polynomials(self, ax=None, order_limit=None, number_of_points=200, show=True):
""" Plots the first few orthogonal polynomials. See :meth:`~equadratures.plot.plot_orthogonal_polynomials` for full description. """
return plot.plot_orthogonal_polynomials(self,ax,order_limit,number_of_points,show)
[docs]
def plot_pdf(self, ax=None, data=None, show=True, lim_range=True):
""" Plots the probability density function for a Parameter. See :meth:`~equadratures.plot.plot_pdf` for full description. """
return plot.plot_pdf(self,ax, data, show, lim_range)
[docs]
def plot_cdf(self, ax=None, show=True, lim_range=True):
""" Plots the cumulative density function for a Parameter. See :meth:`~equadratures.plot.plot_cdf` for full description. """
return plot.plot_cdf(self,ax, show, lim_range)
[docs]
def get_pdf(self, points=None):
""" Computes the probability density function associated with the Parameter.
Parameters
----------
points : numpy.ndarray, optional
Values of the parameter at which the PDF must be evaluated.
Returns
-------
numpy.ndarray
If ``points!=None``. ndarray containing the probability density function evaluated at the points in ``points``.
tuple
If ``points=None``. A tuple (`x`, `pdf`), where `pdf` is the probability density function evaluated at the points in `x`.
"""
if points is None:
x = self.x_range_for_pdf
return x, self.get_pdf(x)
else:
return self.get_pdf(points)
[docs]
def get_cdf(self, points=None):
""" Computes the cumulative density function associated with the Parameter.
Parameters
----------
points : numpy.ndarray, optional
Values of the parameter at which the CDF must be evaluated.
Returns
-------
numpy.ndarray
If ``points!=None``. ndarray containing the cumulative density function evaluated at the points in ``points``.
tuple
If ``points=None``. A tuple (`x`, `cdf`), where `cdf` is the cumulative density function evaluated at the points in `x`.
"""
if points is None:
x = self.distribution.x_range_for_pdf
return x, self.distribution.get_cdf(x)
else:
return self.distribution.get_cdf(points)
[docs]
def get_icdf(self, cdf_values):
""" Computes the inverse cumulative density function associated with the Parameter.
Parameters
----------
cdf_values : numpy.ndarray
Values of the cumulative density function for which its inverse needs to be computed.
Returns
-------
numpy.ndarray
The inverse cumulative density function.
"""
return self.distribution.get_icdf(cdf_values)
[docs]
def get_samples(self, number_of_samples_required):
""" Generates samples from the distribution associated with the Parameter.
Parameters
----------
number_of_samples_required : int
Number of samples that are required.
Returns
-------
numpy.ndarray
The generated samples.
"""
return self.distribution.get_samples(number_of_samples_required)
[docs]
def get_description(self):
""" Provides a description of the Parameter.
Returns
-------
str
A description of the parameter.
"""
return self.distribution.get_description()
[docs]
def get_recurrence_coefficients(self, order=None):
""" Generates the recurrence coefficients.
Parameters
----------
order : int, optional
Order of the recurrence coefficients.
Returns
-------
numpy.ndarray
Array of recurrence coefficients.
"""
return self.distribution.get_recurrence_coefficients(order)
[docs]
def get_jacobi_eigenvectors(self, order=None):
""" Computes the eigenvectors of the Jacobi matrix.
Parameters
----------
order : int
Order of the recurrence coefficients.
Returns
-------
numpy.ndarray
Array of eigenvectors.
"""
if order is None:
order = self.order + 1
JacobiMat = self.get_jacobi_matrix(order)
if order == 1:
V = [1.0]
else:
#D,V = np.linalg.eig(self.get_jacobi_matrix(order))
D, V = sc.linalg.eigh(self.get_jacobi_matrix(order))
idx = D.argsort()[::-1]
eigs = D[idx]
eigVecs = V[:, idx]
#V = np.asmatrix(V) # convert to matrix
#i = np.argsort(D) # get the sorted indices
#i = np.array(i) # convert to array
#V = V[:,i]
return eigVecs
[docs]
def get_jacobi_matrix(self, order=None, ab=None):
""" Computes the Jacobi matrix---a tridiagonal matrix of the recurrence coefficients.
Parameters
----------
order : int
Order of the recurrence coefficients.
Returns
-------
numpy.ndarray
2D array containing the Jacobi matrix.
"""
if order is None and ab is None:
ab = self.get_recurrence_coefficients()
order = self.order + 1
elif ab is None:
ab = self.get_recurrence_coefficients(order)
else:
ab = ab[0:order, :]
order = int(order)
# The case of order 1~
if int(order) == 1:
JacobiMatrix = ab[0, 0]
# For everything else~
else:
JacobiMatrix = np.zeros((int(order), int(order))) # allocate space
JacobiMatrix[0,0] = ab[0,0]
JacobiMatrix[0,1] = np.sqrt(ab[1,1])
k = order - 1
for u in range(1, int(k)):
JacobiMatrix[u,u] = ab[u,0]
JacobiMatrix[u,u-1] = np.sqrt(ab[u,1])
JacobiMatrix[u,u+1] = np.sqrt(ab[u+1,1])
JacobiMatrix[order-1, order-1] = ab[order-1,0]
JacobiMatrix[order-1, order-2] = np.sqrt(ab[order-1,1])
return JacobiMatrix
def _get_orthogonal_polynomial(self, points, order=None):
"""
Private function that evaluates the univariate orthogonal polynomial at quadrature points.
:param Parameter self:
An instance of the Parameter object.
:param numpy.ndarray points:
Points at which the orthogonal polynomial must be evaluated.
:param int order:
Order up to which the orthogonal polynomial must be obtained.
"""
if order is None:
order = self.order + 1
else:
order = order + 1
gridPoints = np.asarray(points).copy()
ab = self.get_recurrence_coefficients(order)
"""
print('Before:')
print(gridPoints)
for q in range(0, gridPoints.shape[0]):
if (gridPoints[q] < self.bounds[0]) or (gridPoints[q] > self.bounds[1]):
grid_flag = 1
if grid_flag == 1:
for r in range(0, gridPoints.shape[0]):
gridPoints[r] = (self.bounds[1] - self.bounds[0]) * ( (gridPoints[r] - self.lower) / (self.upper - self.lower) ) + self.bounds[0]
#print(gridPoints)
print('After:')
print(gridPoints)
"""
orthopoly = np.zeros((order, len(gridPoints))) # create a matrix full of zeros
derivative_orthopoly = np.zeros((order, len(gridPoints)))
dderivative_orthopoly = np.zeros((order, len(gridPoints)))
# Convert the grid points to a numpy array -- simplfy life!
gridPointsII = gridPoints.reshape((-1, 1))
orthopoly[0, :] = 1.0
# Cases
if order == 1: # CHANGED 2/2/18
return orthopoly, derivative_orthopoly, dderivative_orthopoly
orthopoly[1, :] = ((gridPointsII[:, 0] - ab[0, 0]) * orthopoly[0, :]) * (1.0) / (1.0 * np.sqrt(ab[1, 1]))
derivative_orthopoly[1, :] = orthopoly[0, :] / (np.sqrt(ab[1, 1]))
if order == 2: # CHANGED 2/2/18
return orthopoly, derivative_orthopoly, dderivative_orthopoly
if order >= 3: # CHANGED 2/2/18
for u in range(2, order): # CHANGED 2/2/18
# Three-term recurrence rule in action!
orthopoly[u, :] = (((gridPointsII[:, 0] - ab[u - 1, 0]) * orthopoly[u - 1, :]) - np.sqrt(
ab[u - 1, 1]) * orthopoly[u - 2, :]) / (1.0 * np.sqrt(ab[u, 1]))
for u in range(2, order): # CHANGED 2/2/18
# Four-term recurrence formula for derivatives of orthogonal polynomials!
derivative_orthopoly[u,:] = ( ((gridPointsII[:,0] - ab[u-1,0]) * derivative_orthopoly[u-1,:]) - ( np.sqrt(ab[u-1,1]) * derivative_orthopoly[u-2,:] ) + orthopoly[u-1,:] )/(1.0 * np.sqrt(ab[u,1]))
for u in range(2,order):
# Four-term recurrence formula for second derivatives of orthogonal polynomials!
dderivative_orthopoly[u,:] = ( ((gridPointsII[:,0] - ab[u-1,0]) * dderivative_orthopoly[u-1,:]) - ( np.sqrt(ab[u-1,1]) * dderivative_orthopoly[u-2,:] ) + 2.0 * derivative_orthopoly[u-1,:] )/(1.0 * np.sqrt(ab[u,1]))
return orthopoly, derivative_orthopoly, dderivative_orthopoly
def _get_local_quadrature(self, order=None, ab=None):
"""
Returns the 1D quadrature points and weights for the parameter. WARNING: Should not be called under normal circumstances.
:param Parameter self:
An instance of the Parameter class
:param int N:
Number of quadrature points and weights required. If order is not specified, then by default the method will return the number of points defined in the parameter itself.
:return:
A N-by-1 matrix that contains the quadrature points
:return:
A 1-by-N matrix that contains the quadrature weights
"""
if self.endpoints is None:
return get_local_quadrature(self, order, ab)
elif self.endpoints.lower() == 'lower' or self.endpoints.lower() == 'upper':
return get_local_quadrature_radau(self, order, ab)
elif self.endpoints.lower() == 'both':
return get_local_quadrature_lobatto(self, order, ab)
else:
raise(ValueError, 'Error in endpoints specification.')
def get_local_quadrature(self, order=None, ab=None):
# Check for extra input argument!
if order is None:
order = self.order + 1
else:
order = order + 1
if ab is None:
# Get the recurrence coefficients & the jacobi matrix
JacobiMat = self.get_jacobi_matrix(order)
ab = self.get_recurrence_coefficients(order+1)
else:
ab = ab[0:order+1,:]
JacobiMat = self.get_jacobi_matrix(order, ab)
# If statement to handle the case where order = 1
if order == 1:
# Check to see whether upper and lower bound are defined:
if np.isinf(self.distribution.lower) or np.isinf(self.distribution.upper):
p = np.asarray(self.distribution.mean).reshape((1,1))
else:
print('see below!')
print(self.distribution.lower, self.distribution.upper)
print(self.distribution.lower(), self.distribution.upper())
p = np.asarray((self.distribution.upper - self.distribution.lower)/(2.0) + self.distribution.lower).reshape((1,1))
w = [1.0]
else:
# Compute eigenvalues & eigenvectors of Jacobi matrix
#D,V = np.linalg.eig(JacobiMat)
D, V = sc.linalg.eigh(JacobiMat)
#V = np.asmatrix(V) # convert to matrix
#local_points = np.sort(D) # sort by the eigenvalues
#i = np.argsort(D) # get the sorted indices
#i = np.array(i) # convert to array
idx = D.argsort()[::-1]
eigs = D[idx]
eigVecs = V[:, idx]
w = np.linspace(1,order+1,order) # create space for weights
p = np.ones((int(order),1))
for u in range(0, len(idx) ):
w[u] = float(ab[0,1]) * (eigVecs[0,idx[u]]**2) # replace weights with right value
p[u,0] = eigs[u]
#if (p[u,0] < 1e-16) and (-1e-16 < p[u,0]):
# p[u,0] = np.abs(p[u,0])
p = p[::-1]
return p, w
def get_local_quadrature_radau(self, order=None, ab=None):
if self.endpoints.lower() == 'lower':
end0 = self.distribution.lower
elif self.endpoints.lower() == 'upper':
end0 = self.distribution.upper
if order is None:
order = self.order - 1
else:
order = order - 1
N = order
if ab is None:
ab = self.get_recurrence_coefficients(order+1)
else:
ab = ab[0:order+1, :]
p0 = 0.
p1 = 1.
for i in range(0, N+1):
pm1 = p0
p0 = p1
p1 = (end0 - ab[i, 0]) * p0 - ab[i, 1]*pm1
ab[N+1, 0] = end0 - ab[N+1, 1] * p0/p1
return get_local_quadrature(self, order=order+1, ab=ab)
def get_local_quadrature_lobatto(self, order=None, ab=None):
if order is None:
order = self.order - 2
else:
order = order - 2
N = order
endl = self.distribution.lower
endr = self.distribution.upper
if ab is None:
ab = self.get_recurrence_coefficients(order+2)
else:
ab = ab[0:order+2, :]
p0l = 0.
p0r = 0.
p1l = 1.
p1r = 1.
for i in range(0, N+2):
pm1l = p0l
p0l = p1l
pm1r = p0r
p0r = p1r
p1l = (endl - ab[i, 0]) * p0l - ab[i, 1] * pm1l
p1r = (endr - ab[i, 0]) * p0r - ab[i, 1] * pm1r
det = p1l * p0r - p1r * p0l
ab[N+2, 0] = (endl*p1l*p0r-endr*p1r*p0l)/det
ab[N+2, 1] = (endr - endl) * p1l * p1r/det
return get_local_quadrature(self, order=order+2, ab=ab)
