(&self, mu_2: f64, p_1: f64)
| 99 | } |
| 100 | } |
| 101 | |
| 102 | pub fn iter_4(&self, mu_2: f64, p_1: f64) -> Vec<f64> { |
| 103 | let m_1 = self.moments[0]; |
| 104 | let m_2 = self.moments[1]; |
| 105 | let m_3 = self.moments[2]; |
| 106 | let m_4 = self.moments[3]; |
| 107 | |
| 108 | let mu_1 = (m_1 - (1.0 - p_1) * mu_2) / p_1; |
| 109 | let den_24 = 3.0 * (1.0 - p_1) * (mu_2 - mu_1); |
| 110 | if den_24 == 0.0 { |
| 111 | return vec![]; |
| 112 | } |
| 113 | let sigma_2_squared = (m_3 + 2.0 * p_1 * mu_1.powi(3) + (p_1 - 1.0) * mu_2.powi(3) |
| 114 | - 3.0 * mu_1 * (m_2 + mu_2.powi(2) * (p_1 - 1.0))) |
| 115 | / den_24; |
| 116 | if sigma_2_squared < 0.0 { |
| 117 | return vec![]; |
| 118 | } |
| 119 | let sigma_2 = sigma_2_squared.sqrt(); |
| 120 | |
| 121 | let sigma_1_squared = |
| 122 | ((m_2 - sigma_2.powi(2) - mu_2.powi(2)) / p_1) + sigma_2.powi(2) + mu_2.powi(2) |
| 123 | - mu_1.powi(2); |
| 124 | if sigma_1_squared < 0.0 { |
| 125 | return vec![]; |
| 126 | } |
| 127 | let sigma_1 = sigma_1_squared.sqrt(); |
| 128 | |
| 129 | let p_1_deno = 3.0 * (sigma_1.powi(4) - sigma_2.powi(4)) |
| 130 | + 6.0 * (sigma_1.powi(2) * mu_1.powi(2) - sigma_2.powi(2) * mu_2.powi(2)) |
| 131 | + mu_1.powi(4) |
| 132 | - mu_2.powi(4); |
| 133 | if p_1_deno == 0.0 { |
| 134 | return vec![]; |
| 135 | } |
| 136 | let p_1_new = |
| 137 | (m_4 - 3.0 * sigma_2.powi(4) - 6.0 * sigma_2.powi(2) * mu_2.powi(2) - mu_2.powi(4)) |
| 138 | / p_1_deno; |
| 139 | if !(0.0..=1.0).contains(&p_1_new) { |
| 140 | return vec![]; |
| 141 | } |
| 142 | |
| 143 | vec![mu_1, mu_2, sigma_1, sigma_2, p_1_new] |
| 144 | } |
| 145 | |
| 146 | pub fn iter_5(&self, mu_2: f64, p_1: f64) -> Vec<f64> { |
no outgoing calls