Conditions | 6 |
Total Lines | 76 |
Code Lines | 42 |
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 |
||
67 | def OU_fit(self, bins=None, rang=None, a=1e-5, b=100): |
||
68 | """Description goes here.""" |
||
69 | |||
70 | # TODO: make a generic fit method |
||
71 | res = dict() |
||
72 | |||
73 | # estimate sigma |
||
74 | try: |
||
75 | # dy = np.diff(self.values) |
||
76 | # dt = np.diff(self.times) |
||
77 | sigma_est = (np.nanmean(np.diff(self.values) ** 2 / self.values[:-1] ** 2 / np.diff(self.times))) ** 0.5 |
||
78 | res["sigma_est"] = sigma_est |
||
79 | except Exception as e: |
||
80 | raise Exception(f"Could not estimate sigma: {e}") |
||
81 | |||
82 | # plot histogram |
||
83 | if bins is None: |
||
84 | bins = np.int(self.values.size ** 0.5 / 1.5) # bins='auto' |
||
85 | if rang is None: |
||
86 | rang = (np.percentile(self.values, 0), np.percentile(self.values, 99)) |
||
87 | p, x = np.histogram(self.values, density=True, bins=bins, range=rang) # bins='sqrt') |
||
88 | x = (x + np.roll(x, -1))[:-1] / 2.0 |
||
89 | |||
90 | plt.subplots() |
||
91 | plt.hist(self.values, density=True, alpha=0.75, bins=bins, range=rang) |
||
92 | plt.plot(x, p, "r-", alpha=0.5) |
||
93 | anchored_text = AnchoredText( |
||
94 | f"mean {np.mean(self.values):.2g} \n " |
||
95 | "median {np.median(self.values):.2g} \n " |
||
96 | "mode {scipy_stats.mode(self.values)[0][0]:.2g} \n " |
||
97 | "std {np.std(self.values):.2g} \n " |
||
98 | "var {np.var(self.values):.2g}", |
||
99 | loc="upper right", |
||
100 | ) |
||
101 | plt.gca().add_artist(anchored_text) |
||
102 | |||
103 | # curve_fit |
||
104 | try: |
||
105 | popt, pcov = scipy_optimize.curve_fit(f=self.pdf, xdata=x, ydata=p) |
||
106 | # print('curve_fit: (l, mu)') |
||
107 | # print('popt: ') |
||
108 | # print(popt) |
||
109 | # print('pcov: ') |
||
110 | # print(np.sqrt(np.diag(pcov))) |
||
111 | x_c = np.linspace(1e-5, 1.1 * np.max(x), 1000) |
||
112 | plt.plot(x_c, self.pdf(x_c, *popt), "k-", label="curve_fit", alpha=0.8) |
||
113 | res["curve_fit"] = (popt, np.sqrt(np.diag(pcov))) |
||
114 | except Exception as e: |
||
115 | raise Exception(f"Some error fitting with curve_fit {e}") |
||
116 | |||
117 | # TODO: place outside as a helper class |
||
118 | # and understand how it is used |
||
119 | # is it mu really needed ? |
||
120 | # mu here is different than self.mu |
||
121 | # fit pdf with MLE |
||
122 | class OU(scipy_stats.rv_continuous): |
||
123 | def _pdf(self, x, l, mu): |
||
124 | return (l * mu) ** (1 + l) / scipy_special.gamma(1 + l) * np.exp(-l * mu / x) / x ** (l + 2) |
||
125 | |||
126 | try: |
||
127 | fit = OU(a=a, b=np.percentile(self.values, b)).fit(self.values, 1, 1, floc=0, fscale=1) |
||
128 | # print('MLE fit: (l, mu)') |
||
129 | # print(fit) |
||
130 | x_c = np.linspace(0, 1.1 * np.max(x), 1000) |
||
131 | plt.plot(x_c, self.pdf(x_c, fit[0], fit[1]), "k-.", label="MLE", alpha=0.8) |
||
132 | res["MLE_fit"] = fit[:-2] |
||
133 | except Exception as e: |
||
134 | raise Exception(f"Some error fitting with MLE {e}") |
||
135 | plt.legend(loc="lower right") |
||
136 | plt.show() |
||
137 | |||
138 | # estimate theta |
||
139 | res["th_est1"] = fit[0] * sigma_est ** 2 / 2 |
||
140 | res["th_est2"] = popt[0] * sigma_est ** 2 / 2 |
||
141 | |||
142 | return res |
||
143 | |||
204 |