"""
Functional Tolerance Bounds using SRSF
moduleauthor:: J. Derek Tucker <jdtuck@sandia.gov>
"""
import numpy as np
import fdasrsf as fs
import fdasrsf.utility_functions as uf
import fdasrsf.fPCA as fpca
from scipy.stats import chi2
from numpy.linalg import eig
from fdasrsf.boxplots import ampbox, phbox
[docs]def bootTB(f, time, a=0.5, p=.99, B=500, no=5, parallel=True):
"""
This function computes tolerance bounds for functional data containing
phase and amplitude variation using bootstrap sampling
:param f: numpy ndarray of shape (M,N) of N functions with M samples
:param time: vector of size M describing the sample points
:param a: confidence level of tolerance bound (default = 0.05)
:param p: coverage level of tolerance bound (default = 0.99)
:param B: number of bootstrap samples (default = 500)
:param no: number of principal components (default = 5)
:param parallel: enable parallel processing (default = T)
:type f: np.ndarray
:type time: np.ndarray
:rtype: tuple of boxplot objects
:return amp: amplitude tolerance bounds
:rtype out_med: ampbox object
:return ph: phase tolerance bounds
:rtype out_med: phbox object
:return out_med: alignment results
:rtype out_med: fdawarp object
"""
eps = np.finfo(np.double).eps
(M, N) = f.shape
# Align Data
out_med = fs.fdawarp(f, time)
out_med.srsf_align(method="median", parallel=parallel)
# Calculate CI
# a% tolerance bound with p%
fn = out_med.fn
qn = out_med.qn
gam = out_med.gam
q0 = out_med.q0
print("Bootstrap Sampling")
bootlwr_amp = np.zeros((M,B))
bootupr_amp = np.zeros((M,B))
bootlwr_ph = np.zeros((M,B))
bootupr_ph = np.zeros((M,B))
for k in range(B):
out_med.joint_gauss_model(n=100, no=no)
obja = ampbox(out_med)
obja.construct_boxplot(1-p,.3)
objp = phbox(out_med)
objp.construct_boxplot(1-p,.3)
bootlwr_amp[:,k] = obja.Q1a
bootupr_amp[:,k] = obja.Q3a
bootlwr_ph[:,k] = objp.Q1a
bootupr_ph[:,k] = objp.Q3a
# tolerance bounds
boot_amp = np.hstack((bootlwr_amp, bootupr_amp))
f, g, g2 = uf.gradient_spline(time, boot_amp, False)
boot_amp_q = g / np.sqrt(abs(g) + eps)
boot_ph = np.hstack((bootlwr_ph,bootupr_ph))
boot_out = out_med
boot_out.fn = boot_amp
boot_out.qn = boot_amp_q
boot_out.gam = boot_ph
boot_out.rsamps = False
amp = ampbox(boot_out)
amp.construct_boxplot(a, .3)
ph = phbox(boot_out)
ph.construct_boxplot(a, .3)
return amp, ph, out_med
[docs]def pcaTB(f, time, a=0.5, p=.99, no=5, parallel=True):
"""
This function computes tolerance bounds for functional data containing
phase and amplitude variation using fPCA
:param f: numpy ndarray of shape (M,N) of N functions with M samples
:param time: vector of size M describing the sample points
:param a: confidence level of tolerance bound (default = 0.05)
:param p: coverage level of tolerance bound (default = 0.99)
:param no: number of principal components (default = 5)
:param parallel: enable parallel processing (default = T)
:type f: np.ndarray
:type time: np.ndarray
:rtype: tuple of boxplot objects
:return warp: alignment data from time_warping
:return pca: functional pca from jointFPCA
:return tol: tolerance factor
"""
# Align Data
out_warp = fs.fdawarp(f, time)
out_warp.srsf_align( method="median", parallel=parallel)
# Calculate pca
out_pca = fpca.fdajpca(out_warp)
out_pca.calc_fpca(no)
# Calculate TB
tol = mvtol_region(out_pca.coef, a, p, 100000)
return warp, pca, tol
[docs]def mvtol_region(x, alpha, P, B):
"""
Computes tolerance factor for multivariate normal
Krishnamoorthy, K. and Mondal, S. (2006), Improved Tolerance Factors for Multivariate Normal
Distributions, Communications in Statistics - Simulation and Computation, 35, 461–478.
:param x: (M,N) matrix defining N variables of M samples
:param alpha: confidence level
:param P: coverage level
:param B: number of bootstrap samples
:rtype: double
:return tol: tolerance factor
"""
n,p = x.shape
q_squared = chi2.rvs(1, size=(p,B))/n
L = np.zeros((p,B))
for k in range(B):
L[:,k] = eig(rwishart(n-1,p))[0]
c1 = (1+q_squared)/L
c1 = c1.sum()
c2 = (1+2*q_squared)/(L**2)
c2 = c2.sum()
c3 = (1+3*q_squared)/(L**3)
c3 = c3.sum()
a = (c2**3)/(c3**2)
T = (n-1)*(np.sqrt(c2/a) * (chi2.ppf(P,a)-a) + c1)
tol = np.quantile(T,1-alpha)
return tol
[docs]def rwishart(df,p):
"""
Computes a random wishart matrix
:param df: degree of freedom
:param p: number of dimensions
:rtype: double
:return R: matrix
"""
R = np.zeros((p,p))
R = R.flatten()
R[::p+1] = np.sqrt(chi2.rvs(np.arange(df,df-p,-1),size=p))
if p>1:
pseq = np.arange(0,p)
tmp = [np.arange(0,x+1) for x in np.arange(0,p-1)]
R[np.repeat(p*pseq,pseq)+np.concatenate(tmp).ravel()] = np.random.randn(int(p*(p-1)/2))
R = R.reshape((p,p))
R = np.matmul(R.T,R)
return R