Population morphology implies a common developmental blueprint for Drosophila motion detectors.

1. from fitter import Fitter
2. import scipy.stats as stats
3. import numpy as np
4. import pandas as pd
6. from .figure_Tools import Subtypes, Subtype_colours
9. class DistributionFitter:
10. def __init__(self, df, DV_col, isExternal, xmin, xmax, n_bins, progress=False):
11. """
12. Core class to fit distributions to subsets of data and summarize results.
13. """
14. self.df = df
15. self.DV_col = DV_col
16. self.isExternal = isExternal
17. self.xmin = xmin
18. self.xmax = xmax
19. self.n_bins = n_bins
20. self.progress = progress
22. self.nTypes = ["T4", "T5"]
23. self.nSubtypes = ["a", "b", "c", "d"]
24. self.output_dict = {t: {s: [] for s in self.nSubtypes} for t in self.nTypes}
26. def fit_distributions(self, distributions=None):
27. """
28. Fit distributions to all type-subtype combinations.
29. """
30. if distributions is None:
31. distributions = ["lognorm", "gamma", "wald", "expon", "fisk", "weibull_min"]
33. for t in self.nTypes:
34. for s in self.nSubtypes:
35. if self.progress:
36. print(f"Fitting to {t}{s}")
37. data = self.df.loc[
38. (self.df.Subtype == t + s)
39. & (self.df.isExternal == self.isExternal),
40. self.DV_col,
41. ].values
42. f = Fitter(
43. data,
44. xmin=self.xmin,
45. xmax=self.xmax,
46. bins=self.n_bins,
47. distributions=distributions,
48. )
49. f.fit()
50. self.output_dict[t][s] = {"Data": data, "Fit": f}
51. return self.output_dict
53. @staticmethod
54. def log_lik(data, f, dist_name):
55. """
56. Compute log-likelihood for given data and distribution.
57. """
58. dist = getattr(stats, dist_name)
59. params = f.fitted_param[dist_name]
61. if dist.shapes:
62. n_shapes = len(dist.shapes.split(","))
63. shapes = params[:n_shapes]
64. else:
65. shapes = ()
66. loc = params[len(shapes)]
67. scale = params[len(shapes) + 1]
69. logpdf_vals = dist.logpdf(data, *shapes, loc=loc, scale=scale)
70. return np.sum(logpdf_vals)
72. def bestFit_summary(self):
73. """
74. Summarize best-fitting distributions based on BIC.
75. """
76. summary_data = {
77. "Type": [],
78. "Subtype": [],
79. "Best_distribution": [],
80. "BIC": [],
81. "Log_likelihood": [],
82. "Parameters": [],
83. }
85. for t in self.nTypes:
86. for s in self.nSubtypes:
87. fit_obj = self.output_dict[t][s]["Fit"]
88. data = self.output_dict[t][s]["Data"]
90. best_dist = fit_obj.summary(Nbest=1, method="bic", plot=False)[
91. "bic"
92. ].index[0]
93. best_bic = fit_obj.summary(Nbest=1, method="bic", plot=False)[
94. "bic"
95. ].values[0]
96. params = fit_obj.fitted_param[best_dist]
97. ll = self.log_lik(data, fit_obj, best_dist)
99. summary_data["Type"].append(t)
100. summary_data["Subtype"].append(s)
101. summary_data["Best_distribution"].append(best_dist)
102. summary_data["BIC"].append(best_bic)
103. summary_data["Log_likelihood"].append(ll)
104. summary_data["Parameters"].append(params)
106. return pd.DataFrame(summary_data)
108. @staticmethod
109. def _plot_pdf(ax, xmin, xmax, num_points, dist_name, params=(), c="r", label=""):
110. """
111. Plot a single PDF curve on an axis.
112. """
113. x = np.linspace(xmin, xmax, num_points)
115. try:
116. dist = getattr(stats, dist_name)
117. except AttributeError:
118. raise ValueError(f"Distribution '{dist_name}' not found in scipy.stats")
120. y = dist.pdf(x, *params)
121. ax.plot(x, y, c=c, label=label)
123. def plot_density(self, ax, dist_name, scale=1):
124. """
125. Plot fitted PDFs and pooled data histogram.
126. """
127. for i, st in enumerate(Subtypes):
128. t = st[:2]
129. s = st[-1]
130. c = Subtype_colours[i]
132. params = self.output_dict[t][s]["Fit"].fitted_param[dist_name]
133. self._plot_pdf(
134. ax, self.xmin, self.xmax, 1000, dist_name, params, c=c, label=st
135. )
137. # Pooled data histogram
138. data = self.df.loc[self.df.isExternal == self.isExternal, self.DV_col].values
139. counts, bins = np.histogram(
140. data, range=(self.xmin, self.xmax), bins=self.n_bins, density=True
141. )
142. ax.bar(
143. x=bins[1:],
144. height=counts,
145. width=(bins[1] - bins[0]) * scale,
146. color="gray",
147. alpha=0.4,
148. label="Observed Pooled data",
149. )
152. def BIC_comparison_plot(ax, dist_fits):
154. index_order = ["gamma", "lognorm", "weibull_min", "fisk", "wald", "expon"]
155. x_labels = [
156. "Gamma",
157. "Log Normal",
158. "Minimum \nWeibull",
159. "Log Logistic",
160. "Wald",
161. "Exponential",
162. ]
163. metric = "bic"
164. types = ["T4", "T5"]
165. subtypes = ["a", "b", "c", "d"]
167. vals = []
168. for i, n_type in enumerate(types):
169. for j, n_subtype in enumerate(subtypes):
170. tmp_df = dist_fits.output_dict[n_type][n_subtype]["Fit"].df_errors.reindex(
171. index_order
172. )
173. vals.append(tmp_df[metric].values)
175. vals = np.column_stack(vals)
177. var = np.std(vals, axis=1, ddof=1)
178. means = vals.mean(axis=1)
180. x = np.arange(1, 7)
182. ax.errorbar(
183. x=x - 0.05,
184. y=means,
185. yerr=var,
186. fmt="o",
187. color="gray",
188. alpha=1,
189. ms = 4,
190. label="mean $\pm$ 1 std",
191. )
193. ax.axhline(
194. y=means.min(), xmin=0, xmax=8, ls="--", c="k", label="Minimum \nMean BIC", alpha = 0.6
195. )
197. # add individual population points
198. for i, subtype in enumerate(Subtypes):
199. ax.scatter(
200. x=x + 0.05,
201. y=vals[:, i],
202. c=Subtype_colours[i],
203. s=8,
204. label=subtype,
205. alpha=1,
206. zorder=100,
207. )
209. ax.set_xticks(ticks=x, labels=x_labels)
210. ax.tick_params(axis="x", labelrotation=45)
212. ax.set_xlabel("Fitted Distribution")
213. ax.set_ylabel("BIC")
215. ax.spines["right"].set_visible(False)
216. ax.spines["top"].set_visible(False)

dist_fit.py at commit 3a1aa1a, no license · at the source