| Conditions | 8 |
| Total Lines | 101 |
| Code Lines | 55 |
| Lines | 0 |
| Ratio | 0 % |
| Changes | 0 | ||
Small methods make your code easier to understand, in particular if combined with a good name. Besides, if your method is small, finding a good name is usually much easier.
For example, if you find yourself adding comments to a method's body, this is usually a good sign to extract the commented part to a new method, and use the comment as a starting point when coming up with a good name for this new method.
Commonly applied refactorings include:
If many parameters/temporary variables are present:
| 1 | # Licensed under a 3-clause BSD style license - see LICENSE |
||
| 63 | def OU_fit(self, bins=None, rang=None, a=1e-5, b=100): |
||
| 64 | """Fit the signal to an OU stochastic process, using several statistical approaches. |
||
| 65 | |||
| 66 | This function tries to fit the signal to an OU stochastic |
||
| 67 | process using both basic curve fitting and the Maximum |
||
| 68 | Likelihood Estimation method, and returns some plots of |
||
| 69 | the signals and its properties, and the estimated parameters. |
||
| 70 | """ |
||
| 71 | |||
| 72 | # TODO: make a generic fit method |
||
| 73 | res = dict() |
||
| 74 | |||
| 75 | y = self.values |
||
| 76 | t = self.times |
||
| 77 | |||
| 78 | dy = np.diff(y) |
||
| 79 | dt = np.diff(t) |
||
| 80 | |||
| 81 | # estimate sigma |
||
| 82 | try: |
||
| 83 | sigma_est = (np.nanmean(dy ** 2 / y[:-1] ** 2 / dt)) ** 0.5 |
||
| 84 | res["sigma_est"] = sigma_est |
||
| 85 | except Exception as e: |
||
| 86 | log.error(f"Could not estimate sigma: {e}") |
||
| 87 | |||
| 88 | # plot histogram |
||
| 89 | if bins is None: |
||
| 90 | bins = np.int(y.size ** 0.5 / 1.5) # bins='auto' |
||
| 91 | if rang is None: |
||
| 92 | rang = (np.percentile(y, 0), np.percentile(y, 99)) |
||
| 93 | |||
| 94 | p, x = np.histogram(y, density=True, bins=bins, range=rang) # bins='sqrt') |
||
| 95 | x = (x + np.roll(x, -1))[:-1] / 2.0 |
||
| 96 | |||
| 97 | plt.subplots() |
||
| 98 | |||
| 99 | plt.hist(y, density=True, alpha=0.75, bins=bins, range=rang) |
||
| 100 | plt.plot(x, p, "r-", alpha=0.5) |
||
| 101 | |||
| 102 | anchored_text = AnchoredText( |
||
| 103 | f"mean {np.mean(y):.2g} \n " |
||
| 104 | f"median {np.median(y):.2g} \n " |
||
| 105 | f"mode {scipy_stats.mode(y)[0][0]:.2g} \n " |
||
| 106 | f"std {np.std(y):.2g} \n " |
||
| 107 | f"var {np.var(y):.2g}", |
||
| 108 | loc="upper right", |
||
| 109 | ) |
||
| 110 | plt.gca().add_artist(anchored_text) |
||
| 111 | |||
| 112 | # curve_fit |
||
| 113 | try: |
||
| 114 | popt, pcov = scipy_optimize.curve_fit(f=self.pdf, xdata=x, ydata=p) |
||
| 115 | # print('curve_fit: (l, mu)') |
||
| 116 | # print('popt: ') |
||
| 117 | # print(popt) |
||
| 118 | # print('pcov: ') |
||
| 119 | # print(np.sqrt(np.diag(pcov))) |
||
| 120 | x_c = np.linspace(1e-5, 1.1 * np.max(x), 1000) |
||
| 121 | plt.plot(x_c, self.pdf(x_c, *popt), "k-", label="curve_fit", alpha=0.8) |
||
| 122 | res["curve_fit"] = (popt, np.sqrt(np.diag(pcov))) |
||
| 123 | except Exception as e: |
||
| 124 | log.error(f"Some error fitting with curve_fit {e}") |
||
| 125 | |||
| 126 | # fit pdf with MLE |
||
| 127 | # TODO: place outside as a helper class |
||
| 128 | # mu here is different than self.mu |
||
| 129 | class OU(scipy_stats.rv_continuous): |
||
| 130 | def _pdf(self, x, l, mu): |
||
| 131 | return ( |
||
| 132 | (l * mu) ** (1 + l) |
||
| 133 | / scipy_special.gamma(1 + l) |
||
| 134 | * np.exp(-l * mu / x) |
||
| 135 | / x ** (l + 2) |
||
| 136 | ) |
||
| 137 | |||
| 138 | try: |
||
| 139 | fit = OU(a=a, b=np.percentile(y, b)).fit(y, 1, 1, floc=0, fscale=1) |
||
| 140 | # print('MLE fit: (l, mu)') |
||
| 141 | # print(fit) |
||
| 142 | x_c = np.linspace(0, 1.1 * np.max(x), 1000) |
||
| 143 | plt.plot(x_c, self.pdf(x_c, fit[0], fit[1]), "k-.", label="MLE", alpha=0.8) |
||
| 144 | res["MLE_fit"] = fit[:-2] |
||
| 145 | except Exception as e: |
||
| 146 | log.error(f"Some error fitting with MLE {e}") |
||
| 147 | |||
| 148 | plt.legend(loc="lower right") |
||
| 149 | # plt.show() |
||
| 150 | |||
| 151 | # estimate theta (from curve_fit) |
||
| 152 | try: |
||
| 153 | res["th_est1"] = fit[0] * sigma_est ** 2 / 2 |
||
| 154 | except NameError as e: |
||
| 155 | res["th_est1"] = None |
||
| 156 | |||
| 157 | # estimate theta (from MLE) |
||
| 158 | try: |
||
| 159 | res["th_est2"] = popt[0] * sigma_est ** 2 / 2 |
||
| 160 | except NameError as e: |
||
| 161 | res["th_est2"] = None |
||
| 162 | |||
| 163 | return res |
||
| 164 | |||
| 254 |