|
37 | 37 | "source": [ |
38 | 38 | "import matplotlib.pyplot as plt\n", |
39 | 39 | "import numpy as np\n", |
40 | | - "from unyt import Jy, nJy\n", |
| 40 | + "from unyt import nJy\n", |
41 | 41 | "\n", |
42 | 42 | "from synference import DepthUncertaintyModel\n", |
43 | 43 | "\n", |
|
158 | 158 | "source": [ |
159 | 159 | "This noise model is implemented in the ```AsinhEmpiricalUncertaintyModel``` class. It works similarly to the ```GeneralEmpiricalUncertaintyModel```, but applies the asinh transformation to the fluxes before fitting the noise model. It only excepts input fluxes in Jy units, and whilst the ```apply_noise``` method will accept the ```out_flux_unit``` argument for consistency with other noise models, the output fluxes will always be in asinh magnitudes.\n", |
160 | 160 | "\n", |
161 | | - "For this model, the asinh softening parameter is not set directly, but in multiples of the median standard deviation of the flux uncertainties in the training data. This is set using the ```f_b_factor``` argument when initializing the model. For example, setting ```f_b_factor=1.0``` will set the softening parameter to the median flux uncertainty, while ```f_b_factor=2.0``` will set it to twice the median flux uncertainty. By default, ```f_b_factor=5.0```, so the softening parameter is set to five times the median flux uncertainty (aka the $5\\sigma$ detection limit)." |
| 161 | + "For this model, the asinh softening parameter is not set directly, but in multiples of the median standard deviation of the flux uncertainties in the training data. This is set using the ```f_b_factor``` argument when initializing the model. For example, setting ```f_b_factor=1.0``` will set the softening parameter to the median flux uncertainty, while ```f_b_factor=2.0``` will set it to twice the median flux uncertainty. By default, ```f_b_factor=5.0```, so the softening parameter is set to five times the median flux uncertainty (aka the $5\\sigma$ detection limit).\n", |
| 162 | + "\n", |
| 163 | + "For this example we setup a more realstic scenario, with a catalogue of 50,000 fake sources with a fixed background level, a fractioanl error and an overall lognormal flux distribution." |
162 | 164 | ] |
163 | 165 | }, |
164 | 166 | { |
|
169 | 171 | "outputs": [], |
170 | 172 | "source": [ |
171 | 173 | "from synference import AsinhEmpiricalUncertaintyModel\n", |
| 174 | + "from synference.utils import f_jy_to_asinh\n", |
| 175 | + "\n", |
| 176 | + "# --- 1. Define Realistic Survey Parameters ---\n", |
| 177 | + "N_SOURCES: int = 50_000\n", |
| 178 | + "\n", |
| 179 | + "# Background noise limit (e.g., 100 uJy).\n", |
| 180 | + "# This dominates the error for faint sources.\n", |
| 181 | + "BACKGROUND_NOISE_JY: float = 0.0001\n", |
| 182 | + "\n", |
| 183 | + "# Fractional/calibration error (e.g., 2%).\n", |
| 184 | + "# This dominates the error for bright sources.\n", |
| 185 | + "FRACTIONAL_ERROR: float = 0.02\n", |
| 186 | + "\n", |
| 187 | + "# Log-normal distribution parameters for \"true\" fluxes.\n", |
| 188 | + "# We set the median flux to be near the noise floor.\n", |
| 189 | + "FLUX_MEDIAN_JY: float = 0.00015\n", |
| 190 | + "FLUX_LOG_SIGMA: float = 1.5 # Width of the log-normal distribution\n", |
| 191 | + "\n", |
| 192 | + "# --- 2. Generate \"True\" Fluxes ---\n", |
| 193 | + "# Use a log-normal distribution: many faint sources, few bright ones\n", |
| 194 | + "true_flux_jy = np.random.lognormal(\n", |
| 195 | + " mean=np.log(FLUX_MEDIAN_JY), sigma=FLUX_LOG_SIGMA, size=N_SOURCES\n", |
| 196 | + ")\n", |
| 197 | + "\n", |
| 198 | + "# --- 3. Calculate \"Ideal\" Error for Each True Flux ---\n", |
| 199 | + "# This is the characteristic error model: sqrt(bg_noise^2 + (frac_err * flux)^2)\n", |
| 200 | + "ideal_error_jy = np.sqrt(BACKGROUND_NOISE_JY**2 + (FRACTIONAL_ERROR * true_flux_jy) ** 2)\n", |
| 201 | + "\n", |
| 202 | + "# --- 4. Simulate \"Observed\" Fluxes by Scattering by the Error ---\n", |
| 203 | + "# The observed flux is the true flux plus Gaussian noise\n", |
| 204 | + "observed_flux_jy = true_flux_jy + np.random.normal(loc=0.0, scale=ideal_error_jy, size=N_SOURCES)\n", |
172 | 205 | "\n", |
173 | | - "# Generate mock data\n", |
174 | | - "fluxes = np.random.uniform(0.1, 100, size=10_000) * Jy\n", |
175 | | - "errors = np.random.normal(0.0, 5.0, size=10_000) * Jy\n", |
| 206 | + "# The \"observed error\" is the error the pipeline *reports*.\n", |
| 207 | + "# We assume the pipeline correctly estimates the ideal error.\n", |
| 208 | + "observed_error_jy = ideal_error_jy\n", |
176 | 209 | "\n", |
| 210 | + "# Plot a flux histogram\n", |
| 211 | + "plt.hist(true_flux_jy, label=\"Fluxes\", bins=100)\n", |
| 212 | + "plt.yscale(\"log\")\n", |
| 213 | + "plt.xlabel(\"Flux [Jy]\")" |
| 214 | + ] |
| 215 | + }, |
| 216 | + { |
| 217 | + "cell_type": "markdown", |
| 218 | + "id": "5dc934bb", |
| 219 | + "metadata": {}, |
| 220 | + "source": [ |
| 221 | + "From this realistic dataset we can generate our `AsinhEmpiricalUncertaintyModel`. Note that some datapoints fall below the softening parameter - these are originally negative fluxes which can be represented by asinh magnitudes." |
| 222 | + ] |
| 223 | + }, |
| 224 | + { |
| 225 | + "cell_type": "code", |
| 226 | + "execution_count": null, |
| 227 | + "id": "2db14a1e", |
| 228 | + "metadata": {}, |
| 229 | + "outputs": [], |
| 230 | + "source": [ |
177 | 231 | "noise_model = AsinhEmpiricalUncertaintyModel(\n", |
178 | | - " observed_phot_jy=fluxes,\n", |
179 | | - " observed_phot_errors_jy=errors,\n", |
| 232 | + " observed_phot_jy=observed_flux_jy,\n", |
| 233 | + " observed_phot_errors_jy=observed_error_jy,\n", |
180 | 234 | ")\n", |
181 | 235 | "\n", |
182 | 236 | "print(noise_model.b)\n", |
183 | 237 | "\n", |
184 | | - "noise_model.plot()" |
| 238 | + "# plot noise_model.b\n", |
| 239 | + "fig, ax = plt.subplots()\n", |
| 240 | + "noise_model.plot(ax=ax)\n", |
| 241 | + "\n", |
| 242 | + "\n", |
| 243 | + "plt.axvline(f_jy_to_asinh(noise_model.b))" |
| 244 | + ] |
| 245 | + }, |
| 246 | + { |
| 247 | + "cell_type": "code", |
| 248 | + "execution_count": null, |
| 249 | + "id": "9b18c966", |
| 250 | + "metadata": {}, |
| 251 | + "outputs": [], |
| 252 | + "source": [ |
| 253 | + "?truncnorm" |
185 | 254 | ] |
186 | 255 | }, |
187 | 256 | { |
|
0 commit comments