(&self, mu_2: f64, p_1: f64)
| 143 | } |
| 144 | |
| 145 | pub fn iter_5(&self, mu_2: f64, p_1: f64) -> Vec<f64> { |
| 146 | let m_1 = self.moments[0]; |
| 147 | let m_2 = self.moments[1]; |
| 148 | let m_3 = self.moments[2]; |
| 149 | let m_4 = self.moments[3]; |
| 150 | let m_5 = self.moments[4]; |
| 151 | |
| 152 | let mu_1 = (m_1 - (1.0 - p_1) * mu_2) / p_1; |
| 153 | let den_24 = 3.0 * (1.0 - p_1) * (mu_2 - mu_1); |
| 154 | if den_24 == 0.0 { |
| 155 | return vec![]; |
| 156 | } |
| 157 | let sigma_2_squared = (m_3 + 2.0 * p_1 * mu_1.powi(3) + (p_1 - 1.0) * mu_2.powi(3) |
| 158 | - 3.0 * mu_1 * (m_2 + mu_2.powi(2) * (p_1 - 1.0))) |
| 159 | / den_24; |
| 160 | if sigma_2_squared < 0.0 { |
| 161 | return vec![]; |
| 162 | } |
| 163 | let sigma_2 = sigma_2_squared.sqrt(); |
| 164 | |
| 165 | let sigma_1_squared = |
| 166 | ((m_2 - sigma_2.powi(2) - mu_2.powi(2)) / p_1) + sigma_2.powi(2) + mu_2.powi(2) |
| 167 | - mu_1.powi(2); |
| 168 | if sigma_1_squared < 0.0 { |
| 169 | return vec![]; |
| 170 | } |
| 171 | let sigma_1 = sigma_1_squared.sqrt(); |
| 172 | |
| 173 | if (1.0 - p_1) < 1e-4 { |
| 174 | return vec![]; |
| 175 | } |
| 176 | let a_1_squared = 6.0 * sigma_2.powi(4) |
| 177 | + (m_4 |
| 178 | - p_1 |
| 179 | * (3.0 * sigma_1.powi(4) |
| 180 | + 6.0 * sigma_1.powi(2) * mu_1.powi(2) |
| 181 | + mu_1.powi(4))) |
| 182 | / (1.0 - p_1); |
| 183 | if a_1_squared < 0.0 { |
| 184 | return vec![]; |
| 185 | } |
| 186 | let a_1 = a_1_squared.sqrt(); |
| 187 | let mu_2_squared = a_1 - 3.0 * sigma_2.powi(2); |
| 188 | if !mu_2_squared.is_finite() || mu_2_squared < 0.0 { |
| 189 | return vec![]; |
| 190 | } |
| 191 | let mu_2_new = mu_2_squared.sqrt(); |
| 192 | |
| 193 | let a_2 = |
| 194 | 15.0 * sigma_1.powi(4) * mu_1 + 10.0 * sigma_1.powi(2) * mu_1.powi(3) + mu_1.powi(5); |
| 195 | let b_2 = 15.0 * sigma_2.powi(4) * mu_2_new |
| 196 | + 10.0 * sigma_2.powi(2) * mu_2_new.powi(3) |
| 197 | + mu_2_new.powi(5); |
| 198 | if (a_2 - b_2) == 0.0 { |
| 199 | return vec![]; |
| 200 | } |
| 201 | let p_1_new = (m_5 - b_2) / (a_2 - b_2); |
| 202 | if !(0.0..=1.0).contains(&p_1_new) { |
no outgoing calls