| 316 | } |
| 317 | |
| 318 | pub fn most_likely_parameters( |
| 319 | data: &[FitResultRow], |
| 320 | ignore_columns: Option<&[&str]>, |
| 321 | res: usize, |
| 322 | ) -> BTreeMap<String, f64> { |
| 323 | let mut out = BTreeMap::new(); |
| 324 | if data.is_empty() { |
| 325 | return out; |
| 326 | } |
| 327 | let ignored: HashSet<&str> = ignore_columns.unwrap_or(&["error"]).iter().copied().collect(); |
| 328 | |
| 329 | let columns: [(&str, Vec<f64>); 6] = [ |
| 330 | ("mu_1", data.iter().map(|r| r.mu_1).collect()), |
| 331 | ("mu_2", data.iter().map(|r| r.mu_2).collect()), |
| 332 | ("sigma_1", data.iter().map(|r| r.sigma_1).collect()), |
| 333 | ("sigma_2", data.iter().map(|r| r.sigma_2).collect()), |
| 334 | ("p_1", data.iter().map(|r| r.p_1).collect()), |
| 335 | ("error", data.iter().map(|r| r.error).collect()), |
| 336 | ]; |
| 337 | |
| 338 | for (name, vals) in columns { |
| 339 | if ignored.contains(name) { |
| 340 | continue; |
| 341 | } |
| 342 | let min = vals.iter().copied().fold(f64::INFINITY, f64::min); |
| 343 | let max = vals.iter().copied().fold(f64::NEG_INFINITY, f64::max); |
| 344 | if (max - min).abs() < 1e-15 { |
| 345 | out.insert(name.to_string(), round_to_5(min)); |
| 346 | continue; |
| 347 | } |
| 348 | |
| 349 | let n = vals.len() as f64; |
| 350 | let mean = vals.iter().sum::<f64>() / n; |
| 351 | let var = vals.iter().map(|v| (v - mean).powi(2)).sum::<f64>() / n.max(1.0); |
| 352 | let std = var.sqrt().max(1e-12); |
| 353 | let h = (std * n.powf(-1.0 / 5.0)).max(1e-6); |
| 354 | |
| 355 | let steps = res.max(10); |
| 356 | let dx = (max - min) / (steps as f64 - 1.0); |
| 357 | let mut best_x = min; |
| 358 | let mut best_y = f64::NEG_INFINITY; |
| 359 | for i in 0..steps { |
| 360 | let x = min + dx * i as f64; |
| 361 | let y = vals |
| 362 | .iter() |
| 363 | .map(|v| { |
| 364 | let u = (x - v) / h; |
| 365 | (-0.5 * u * u).exp() |
| 366 | }) |
| 367 | .sum::<f64>() |
| 368 | / (n * h * (2.0 * std::f64::consts::PI).sqrt()); |
| 369 | if y > best_y { |
| 370 | best_y = y; |
| 371 | best_x = x; |
| 372 | } |
| 373 | } |
| 374 | out.insert(name.to_string(), round_to_5(best_x)); |
| 375 | } |