mirror of
https://github.com/TeachBooks/template.git
synced 2026-11-12 22:25:52 +00:00
247 lines
No EOL
8.6 KiB
Python
247 lines
No EOL
8.6 KiB
Python
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
from scipy.special import erfc
|
|
from scipy import stats
|
|
import seaborn as sns
|
|
|
|
def process_negative_values(q, n, a):
|
|
"""Prevents problems with negative values in sqrt of C_C0 function."""
|
|
N_q = sum(q<0)
|
|
N_n = sum(n<0)
|
|
N_a = sum(a<0)
|
|
set_value = 0.001
|
|
if N_q > 0:
|
|
for i in range(len(q)):
|
|
if q[i]<0:
|
|
q[i] = set_value
|
|
print(f'q has {N_q} negative values that have been set to {set_value:.1e}.')
|
|
if N_n > 0:
|
|
for i in range(len(n)):
|
|
if n[i]<0:
|
|
n[i] = set_value
|
|
print(f'n has {N_n} negative values that have been set to {set_value:.1e}.')
|
|
if N_a > 0:
|
|
for i in range(len(a)):
|
|
if a[i]<0:
|
|
a[i] = set_value
|
|
print(f'a has {N_a} negative values that have been set to {set_value:.1e}.')
|
|
return q, n, a
|
|
|
|
class RandomVariableDistribution():
|
|
"""Define a N or LN probability to use in MUDE week 2.7 exercise.
|
|
|
|
Usage:
|
|
-----
|
|
Create an instance of the class
|
|
by default, the random variable is ~N(0,1)
|
|
Specify distribution using 1st and 2nd moments for N or LN:
|
|
--> use method define_distribution to change it!
|
|
(the method prints a summary of the change)
|
|
for example, copy/paste the following:
|
|
test = RandomVariableDistribution()
|
|
test.print_status()
|
|
test.define_distribution(mean=5, stddev=1, dist_type='LN')
|
|
An instance will include the following attributes:
|
|
- distribution moment and (scipy) parameters shape, location, scale
|
|
- dist_object: distribution object (scipy.stats)
|
|
- get_sample: method for getting random samples from the distribution
|
|
- plot_pdf_or_cdf:
|
|
"""
|
|
|
|
def __init__(self):
|
|
self.mean = 0.
|
|
self.stddev = 1.
|
|
self.dist_type = 'N'
|
|
|
|
def print_status(self):
|
|
"""Returns current parameter values and distribution type."""
|
|
print(f'The mean is {self.mean:.3e}')
|
|
print(f'The standard deviation is {self.stddev:.3e}')
|
|
if self.dist_type == 'N':
|
|
distribution_name = 'Normal'
|
|
elif self.dist_type == 'LN':
|
|
distribution_name = 'Lognormal'
|
|
print('The distribution is ' + distribution_name)
|
|
|
|
def define_distribution(self,
|
|
mean=None,
|
|
stddev=None,
|
|
dist_type=None,
|
|
print_new_dist=False):
|
|
"""Compute distribution parameters and RV object from moments.
|
|
|
|
Usage:
|
|
-----
|
|
Optional arguments will override current dist:
|
|
--> default distribution is standard normal, ~N(0,1)
|
|
mean: float
|
|
stddev: float
|
|
dist_type: 'N' or 'LN'
|
|
print_new_dist=True to print distribution changes.
|
|
"""
|
|
if mean is None:
|
|
mean = self.mean
|
|
if stddev is None:
|
|
stddev = self.stddev
|
|
if dist_type is None:
|
|
dist_type = self.dist_type
|
|
if dist_type == 'N':
|
|
shape = None
|
|
location = mean
|
|
scale = stddev
|
|
dist_object = stats.norm(loc=location, scale=scale)
|
|
elif dist_type == 'LN':
|
|
"""
|
|
Finds distribution parameters of LN-distributed X,
|
|
which are equivalent to the mean and standard deviation
|
|
of the log of X: Y = ln(X), then converts to the
|
|
parameterization used by scipy.stats.lognorm.
|
|
"""
|
|
stddev_y = np.sqrt(np.log(1 + (stddev/mean)**2))
|
|
mean_y = np.log(mean) - 0.5*stddev_y**2
|
|
|
|
shape = stddev_y
|
|
location = 0
|
|
scale = np.exp(mean_y)
|
|
dist_object = stats.lognorm(s=shape,
|
|
loc=location,
|
|
scale=scale)
|
|
else:
|
|
raise ValueError("You can only use N or LN for self.dist_type.")
|
|
self.mean = mean
|
|
self.stddev = stddev
|
|
self.dist_type = dist_type
|
|
self.shape = shape
|
|
self.location = location
|
|
self.scale = scale
|
|
self.dist_object = dist_object
|
|
self.get_sample = lambda size : self.dist_object.rvs(size=size)
|
|
if print_new_dist:
|
|
self.print_status()
|
|
|
|
def plot_pdf_or_cdf(self, plot_type='pdf', plot_limits=None):
|
|
"""Plots the pdf or cdf (specify choice with plot_type).
|
|
|
|
Usage:
|
|
-----
|
|
No inputs required, but there are optional inputs:
|
|
plot_type='cdf' to plot cdf (pdf is default value)
|
|
plot_limits: list with 2 values for x-axis
|
|
"""
|
|
if plot_limits is None:
|
|
if self.dist_type == 'N':
|
|
plot_limits = [self.mean - 3*self.stddev,
|
|
self.mean + 3*self.stddev]
|
|
elif self.dist_type == 'LN':
|
|
'Avoid negative or zero value for lognormal dist.'
|
|
lower_value = self.mean - 3*self.stddev
|
|
if lower_value <= 0:
|
|
lower_value = 0.01
|
|
plot_limits = [lower_value,
|
|
self.mean + 4*self.stddev]
|
|
else:
|
|
raise ValueError("You can only use N or LN for self.dist_type.")
|
|
|
|
plot_x = np.linspace(plot_limits[0], plot_limits[1], 1000)
|
|
if plot_type == 'pdf':
|
|
plot_y = self.dist_object.pdf(plot_x)
|
|
y_label = 'PDF [-]'
|
|
elif plot_type == 'cdf':
|
|
plot_y = self.dist_object.cdf(plot_x)
|
|
y_label = 'CDF [-]'
|
|
else:
|
|
raise ValueError("You must specify pdf or cdf for plot_type.")
|
|
fontsize = 15
|
|
plt.figure()
|
|
plt.plot(plot_x, plot_y, 'k-', linewidth=3)
|
|
plt.xlabel('Random Variable, X', fontsize=fontsize)
|
|
plt.ylabel(y_label, fontsize=fontsize)
|
|
plt.show()
|
|
|
|
def pdf_of_function_of_RV(cases, case_C=True, case_name=None, plot_limits=None):
|
|
"""Plot the PDF of a sampled function of random variables (C_C0 or Z).
|
|
|
|
Usage:
|
|
-----
|
|
cases: a list of arrays
|
|
case_C=True (default) assumes C_C0 is given. False is for Z.
|
|
case_name: list of strings with the labels for each array.
|
|
e.g., ['Original', 'High sigma_a', 'Low sigma_a']
|
|
(default will return Case i for all arrays)
|
|
plot_limits: list of two values that override default limits
|
|
for x-axis.
|
|
"""
|
|
bw_adjust = 3 # adjusts smoothness of Seaborn kernel density contours
|
|
labels = []
|
|
|
|
plt.figure()
|
|
|
|
for i in range(len(cases)):
|
|
if case_name == None:
|
|
labels.append(f'Case {i+1}')
|
|
else:
|
|
labels.append(case_name[i])
|
|
sns.kdeplot(cases[i], label=labels[i], ax=plt.gca(), bw_adjust=bw_adjust)
|
|
plt.legend()
|
|
if case_C:
|
|
plt.xlabel('Concentration ratio, $C/C_0$ [$-$]')
|
|
if plot_limits == None:
|
|
plt.xlim(0,1)
|
|
else:
|
|
plt.xlim(plot_limits)
|
|
else:
|
|
plt.xlabel('Limit-state function, $Z=\mathrm{threshold}-C/C_0$ [$-$]')
|
|
if plot_limits == None:
|
|
plt.xlim(-2,1)
|
|
else:
|
|
plt.xlim(plot_limits)
|
|
plt.ylabel('PDF [$]$')
|
|
plt.grid()
|
|
plt.show()
|
|
|
|
def monte_carlo(Z):
|
|
"""Monte Carlo analysis for sampled function of random variables.
|
|
|
|
Returns a summary of the analysis and a plot of pf and c.o.v..
|
|
The only input is an array of the sampled limit state-function.
|
|
Assumes that failure is when the limit-state function Z<0.
|
|
"""
|
|
N = len(Z)
|
|
Nf = np.sum(Z<0)
|
|
|
|
if Nf == 0:
|
|
raise ValueError('No failure realizations in the sample!')
|
|
elif Nf == N:
|
|
raise ValueError('All realizations in the sample failed!')
|
|
|
|
pf = Nf/N
|
|
|
|
N_f_c = np.cumsum(Z<0)
|
|
N_i = np.arange(1,len(N_f_c)+1,1)
|
|
p_f_c = N_f_c / N_i
|
|
with np.errstate(divide='ignore'):
|
|
V = 1/(np.sqrt(N_i*p_f_c))
|
|
|
|
print(f'MONTE CARLO RESULTS\n'
|
|
f'--------------------------------------------------\n'
|
|
f' N = {N}, number of samples\n'
|
|
f' Nf = {Nf}, number of failed realizations\n'
|
|
f' pf = {pf:.3e}, estimated failure probability\n'
|
|
f' V = {V[-1]:.3f}, c.o.v. of estimated pf\n')
|
|
|
|
plt.loglog(N_i,p_f_c)
|
|
plt.xlim(10,N)
|
|
plt.ylim(10**(np.floor(np.log10(pf)) - 1),1)
|
|
plt.grid()
|
|
plt.xlabel('Realization No.')
|
|
plt.ylabel('Estimated failure probability, $\hat{p}_f$')
|
|
plt.show()
|
|
|
|
plt.plot(N_i,V)
|
|
plt.xscale('log')
|
|
plt.grid()
|
|
plt.xlabel('time step')
|
|
plt.ylabel('c.o.v.')
|
|
plt.show()
|
|
|
|
return Nf,pf |