| 80 | |
| 81 | |
| 82 | class Hotelling: |
| 83 | |
| 84 | def __init__(self, x, y, dt, eta): |
| 85 | self.x = x |
| 86 | self.y = y |
| 87 | self.sd_x = np.std(x) |
| 88 | self.sd_y = np.std(y) |
| 89 | self.delta_t = dt |
| 90 | self.eta = eta |
| 91 | self.xmin = self.x.min() |
| 92 | self.xmax = self.x.max() |
| 93 | self.ymin = self.y.min() |
| 94 | self.ymax = self.y.max() |
| 95 | s = 100000 |
| 96 | self.hx = (self.xmax-self.xmin)/s |
| 97 | self.hy = (self.ymax-self.ymin)/s |
| 98 | self.flag = 'OK' |
| 99 | try: |
| 100 | self.density, self.potential, self.X, self.Y, self.Z = self.potential_estimation() |
| 101 | self.x_ss_min, self.x_ss_max = self.find_minima() |
| 102 | hessian_min = self.hessian_estimation(self.x_ss_min) |
| 103 | hessian_max = self.hessian_estimation(self.x_ss_max) |
| 104 | |
| 105 | self.covariance_min = np.linalg.inv(hessian_min) |
| 106 | self.covariance_max = np.linalg.inv(hessian_max) |
| 107 | pooled_inverse_covariance_matrix = np.linalg.inv(0.5*(self.covariance_max+self.covariance_min)) |
| 108 | |
| 109 | if self.confidence_level(self.x_ss_max)<0.05 or self.confidence_level(self.x_ss_min)<0.05: |
| 110 | self.hotelling_t_squared = 0 |
| 111 | self.flag = 'LOW_DENSITY' |
| 112 | else: |
| 113 | self.hotelling_t_squared = (self.x_ss_max - self.x_ss_min).dot(pooled_inverse_covariance_matrix.dot(self.x_ss_max - self.x_ss_min)) |
| 114 | if np.linalg.norm(self.x_ss_min-self.x_ss_max)<0.1: |
| 115 | self.hotelling_t_squared = 0 |
| 116 | self.flag = 'TOO_CLOSE' |
| 117 | except np.linalg.LinAlgError as err: |
| 118 | if ( ('singular matrix' in str(err)) or ('Singular matrix' in str(err)) ): |
| 119 | self.hotelling_t_squared = 0 |
| 120 | self.flag = 'SINGULAR_MATRIX' |
| 121 | else: |
| 122 | raise |
| 123 | self.path = None |
| 124 | |
| 125 | def potential_estimation(self): |
| 126 | |
| 127 | X, Y = np.mgrid[self.xmin:self.xmax:200j, self.ymin:self.ymax:200j] |
| 128 | positions = np.vstack([X.ravel(), Y.ravel()]) |
| 129 | values = np.vstack([self.x, self.y]) |
| 130 | kernel = gaussian_kde(values) |
| 131 | Z = np.reshape(kernel(positions).T, X.shape) |
| 132 | |
| 133 | return kernel, kernel.logpdf, X, Y, Z |
| 134 | |
| 135 | def calculate_potential_force(self,path): |
| 136 | |
| 137 | #Initialization |
| 138 | Fv = [] |
| 139 | hx = self.hx |