Conditions | 8 |
Total Lines | 97 |
Code Lines | 52 |
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 |
||
75 | def OU_fit(self, bins=None, rang=None, a=1e-5, b=100): |
||
76 | """Fit the signal to an OU stochastic process, using several statistical approaches. |
||
77 | |||
78 | This function tries to fit the signal to an OU stochastic |
||
79 | process using both basic curve fitting and the Maximum |
||
80 | Likelihood Estimation method, and returns some plots of |
||
81 | the signals and its properties, and the estimated parameters. |
||
82 | """ |
||
83 | |||
84 | # TODO: make a generic fit method |
||
85 | res = dict() |
||
86 | |||
87 | y = self.values |
||
88 | t = self.times |
||
89 | |||
90 | dy = np.diff(y) |
||
91 | dt = np.diff(t) |
||
92 | |||
93 | # estimate sigma |
||
94 | try: |
||
95 | sigma_est = (np.nanmean(dy ** 2 / y[:-1] ** 2 / dt)) ** 0.5 |
||
96 | res["sigma_est"] = sigma_est |
||
97 | except Exception as e: |
||
98 | log.error(f"Could not estimate sigma: {e}") |
||
99 | |||
100 | # plot histogram |
||
101 | if bins is None: |
||
102 | bins = np.int(y.size ** 0.5 / 1.5) # bins='auto' |
||
103 | if rang is None: |
||
104 | rang = (np.percentile(y, 0), np.percentile(y, 99)) |
||
105 | |||
106 | p, x = np.histogram(y, density=True, bins=bins, range=rang) # bins='sqrt') |
||
107 | x = (x + np.roll(x, -1))[:-1] / 2.0 |
||
108 | |||
109 | plt.subplots() |
||
110 | |||
111 | plt.hist(y, density=True, alpha=0.75, bins=bins, range=rang) |
||
112 | plt.plot(x, p, "r-", alpha=0.5) |
||
113 | |||
114 | anchored_text = AnchoredText( |
||
115 | f"mean {np.mean(y):.2g} \n " |
||
116 | f"median {np.median(y):.2g} \n " |
||
117 | f"mode {scipy_stats.mode(y)[0][0]:.2g} \n " |
||
118 | f"std {np.std(y):.2g} \n " |
||
119 | f"var {np.var(y):.2g}", |
||
120 | loc="upper right", |
||
121 | ) |
||
122 | plt.gca().add_artist(anchored_text) |
||
123 | |||
124 | # curve_fit |
||
125 | try: |
||
126 | popt, pcov = scipy_optimize.curve_fit(f=self.pdf, xdata=x, ydata=p) |
||
127 | # print('curve_fit: (l, mu)') |
||
128 | # print('popt: ') |
||
129 | # print(popt) |
||
130 | # print('pcov: ') |
||
131 | # print(np.sqrt(np.diag(pcov))) |
||
132 | x_c = np.linspace(1e-5, 1.1 * np.max(x), 1000) |
||
133 | plt.plot(x_c, self.pdf(x_c, *popt), "k-", label="curve_fit", alpha=0.8) |
||
134 | res["curve_fit"] = (popt, np.sqrt(np.diag(pcov))) |
||
135 | except Exception as e: |
||
136 | log.error(f"Some error fitting with curve_fit {e}") |
||
137 | |||
138 | # fit pdf with MLE |
||
139 | # TODO: place outside as a helper class |
||
140 | # mu here is different than self.mu |
||
141 | class OU(scipy_stats.rv_continuous): |
||
142 | def _pdf(self, x, l, mu): |
||
143 | return (l * mu) ** (1 + l) / scipy_special.gamma(1 + l) * np.exp(-l * mu / x) / x ** (l + 2) |
||
144 | |||
145 | try: |
||
146 | fit = OU(a=a, b=np.percentile(y, b)).fit(y, 1, 1, floc=0, fscale=1) |
||
147 | # print('MLE fit: (l, mu)') |
||
148 | # print(fit) |
||
149 | x_c = np.linspace(0, 1.1 * np.max(x), 1000) |
||
150 | plt.plot(x_c, self.pdf(x_c, fit[0], fit[1]), "k-.", label="MLE", alpha=0.8) |
||
151 | res["MLE_fit"] = fit[:-2] |
||
152 | except Exception as e: |
||
153 | log.error(f"Some error fitting with MLE {e}") |
||
154 | |||
155 | plt.legend(loc="lower right") |
||
156 | |||
157 | plt.show() |
||
158 | |||
159 | # estimate theta (from curve_fit) |
||
160 | try: |
||
161 | res["th_est1"] = fit[0] * sigma_est ** 2 / 2 |
||
162 | except NameError as e: |
||
163 | res["th_est1"] = None |
||
164 | |||
165 | # estimate theta (from MLE) |
||
166 | try: |
||
167 | res["th_est2"] = popt[0] * sigma_est ** 2 / 2 |
||
168 | except NameError as e: |
||
169 | res["th_est2"] = None |
||
170 | |||
171 | return res |
||
172 | |||
267 |