|
134 | 134 | "alpha_m_hat = np.sqrt(sw/Delta_prime)\n", |
135 | 135 | "\n", |
136 | 136 | "# Print results\n", |
137 | | - "print(\"Intercept estimate with given uncertainty: ({0:.0f} ± {1:.0f}) mV\".format(c_hat, alpha_c_hat))\n", |
138 | | - "print(\"Slope estimate with given uncertainty: ({0:.2f} ± {1:.2f}) mV/Hz\".format(m_hat, alpha_m_hat))\n", |
| 137 | + "print(\"Intercept estimate with given uncertainty: \"\n", |
| 138 | + " f\"({c_hat:.0f} ± {alpha_c_hat:.0f}) mV\")\n", |
| 139 | + "print(\"Slope estimate with given uncertainty: \"\n", |
| 140 | + " f\"({m_hat:.2f} ± {alpha_m_hat:.2f}) mV/Hz\")\n", |
139 | 141 | "print()\n", |
140 | 142 | "\n", |
141 | 143 | "# Plot data and fit\n", |
|
179 | 181 | "p, V = np.polyfit(frequency, voltage, 1, w=1/alpha_voltage, cov='unscaled')\n", |
180 | 182 | "\n", |
181 | 183 | "# Print results\n", |
182 | | - "print(\"Intercept estimate: ({0:.0f} ± {1:.0f}) mV\".format(p[1], np.sqrt(V[1][1])))\n", |
183 | | - "print(\"Slope estimate: ({0:.2f} ± {1:.2f}) mV/Hz\".format(p[0], np.sqrt(V[0][0])))" |
| 184 | + "print(f\"Intercept estimate: ({p[1]:.0f} ± {np.sqrt(V[1][1]):.0f}) mV\")\n", |
| 185 | + "print(f\"Slope estimate: ({p[0]:.2f} ± {np.sqrt(V[0][0]):.2f}) mV/Hz\")" |
184 | 186 | ] |
185 | 187 | }, |
186 | 188 | { |
|
250 | 252 | "p, V = np.polyfit(load, disp, 1, cov=True)\n", |
251 | 253 | "\n", |
252 | 254 | "# Show optimized fit parameters and uncertainties\n", |
253 | | - "print(\"p[0] = ({0:.3f} ± {1:.3f}) x 1e-7\".format(1e7*p[0], 1e7*np.sqrt(V[0][0])))\n", |
254 | | - "print(\"p[1] = ({0:.1f} ± {1:.1f}) x 1e-3\".format(1e3*p[1], 1e3*np.sqrt(V[1][1])))\n", |
| 255 | + "print(f\"p[0] = ({1e7*p[0]:.3f} ± {1e7*np.sqrt(V[0][0]):.3f}) x 1e-7\")\n", |
| 256 | + "print(f\"p[1] = ({1e3*p[1]:.1f} ± {1e3*np.sqrt(V[1][1]):.1f}) x 1e-3\")\n", |
255 | 257 | "\n", |
256 | 258 | "# Show fit with residuals\n", |
257 | 259 | "from matplotlib.gridspec import GridSpec\n", |
|
319 | 321 | "p2, V2 = np.polyfit(load, disp, 2, cov=True)\n", |
320 | 322 | "\n", |
321 | 323 | "# Show optimized fit parameters and uncertainties\n", |
322 | | - "print(\"p2[0] = ({0:.2f} ± {1:.2f}) x 1e-15\".format(1e15*p2[0], 1e15*np.sqrt(V2[0][0])))\n", |
323 | | - "print(\"p2[1] = ({0:.3f} ± {1:.3f}) x 1e-7\".format(1e7*p2[1], 1e7*np.sqrt(V2[1][1])))\n", |
324 | | - "print(\"p2[2] = ({0:.1f} ± {1:.1f}) x 1e-3\".format(1e3*p2[2], 1e3*np.sqrt(V2[2][2])))\n", |
| 324 | + "print(f\"p2[0] = ({1e15*p2[0]:.2f} ± {1e15*np.sqrt(V2[0][0]):.2f}) x 1e-15\")\n", |
| 325 | + "print(f\"p2[1] = ({1e7*p2[1]:.3f} ± {1e7*np.sqrt(V2[1][1]):.3f}) x 1e-7\")\n", |
| 326 | + "print(f\"p2[2] = ({1e3*p2[2]:.1f} ± {1e3*np.sqrt(V2[2][2]):.1f}) x 1e-3\")\n", |
325 | 327 | "\n", |
326 | 328 | "# Show fit with residuals\n", |
327 | 329 | "from matplotlib.gridspec import GridSpec\n", |
|
353 | 355 | "\n", |
354 | 356 | "# Make the residual plot with a dotted zero line\n", |
355 | 357 | "# Need to adjust ylim, yticks, ylabel for readability\n", |
356 | | - "ax_res = fig.add_subplot(gs[1])\n", |
| 358 | + "fig.add_subplot(gs[1])\n", |
357 | 359 | "plt.plot(load, residuals, 'k.')\n", |
358 | 360 | "plt.plot(load_range, [0,0], 'k:')\n", |
359 | 361 | "plt.xlim(load_range[0], load_range[1])\n", |
|
390 | 392 | "outputs": [], |
391 | 393 | "source": [ |
392 | 394 | "# Define model\n", |
393 | | - "def LCR_model(t, V_bgd, V0, T, phi, tau):\n", |
| 395 | + "def lcr_model(t, V_bgd, V0, T, phi, tau):\n", |
394 | 396 | " return V_bgd + V0*np.cos(2*np.pi*t/T + phi)*np.exp(-t/tau)\n", |
395 | 397 | "\n", |
396 | 398 | "# Define time points\n", |
397 | | - "t = np.linspace(40,950,2000);\n", |
| 399 | + "t = np.linspace(40,950,2000)\n", |
398 | 400 | "\n", |
399 | 401 | "# Define model parameters\n", |
400 | | - "V_bgd = 0.3;\n", |
401 | | - "V0 = 8;\n", |
402 | | - "T = 39;\n", |
403 | | - "phi = 4.5;\n", |
404 | | - "tau = 200;\n", |
| 402 | + "V_bgd = 0.3\n", |
| 403 | + "V0 = 8\n", |
| 404 | + "T = 39\n", |
| 405 | + "phi = 4.5\n", |
| 406 | + "tau = 200\n", |
405 | 407 | "\n", |
406 | 408 | "# Define noise amplitude\n", |
407 | | - "alpha_V = 0.075;\n", |
| 409 | + "alpha_V = 0.075\n", |
408 | 410 | "\n", |
409 | 411 | "# Simulate data\n", |
410 | 412 | "# Compute ideal values\n", |
411 | | - "V_ideal = LCR_model(t, V_bgd, V0, T, phi, tau)\n", |
| 413 | + "V_ideal = lcr_model(t, V_bgd, V0, T, phi, tau)\n", |
412 | 414 | "\n", |
413 | 415 | "# Generate noise with amplitude alpha_V\n", |
414 | 416 | "rg = default_rng(0)\n", |
|
469 | 471 | "\n", |
470 | 472 | "# Plot simulation and preliminary model\n", |
471 | 473 | "plt.plot(t, V_sim, 'k.')\n", |
472 | | - "plt.plot(t, LCR_model(t, V_bgd_init, V0_init, T_init, phi_init, tau_init), 'r-')\n", |
| 474 | + "plt.plot(t, lcr_model(t, V_bgd_init, V0_init, T_init, phi_init, tau_init), 'r-')\n", |
473 | 475 | "plt.xlabel('Time (µs)')\n", |
474 | 476 | "plt.ylabel('Voltage (mV)')\n", |
475 | 477 | "plt.show()" |
|
481 | 483 | "source": [ |
482 | 484 | "### Fit the model to the data and examine the residuals\n", |
483 | 485 | "\n", |
484 | | - "These parameters do not fit the simulated data very well, but they are close enough for the solver to find a best-fit solution that is close to the simulation. We take `alpha_V` as the measurement uncertainty and set `absolute_sigma=True` to tell the solver that it should assume the uncertainty is known independently when calculating the parameter uncertainties. The value of `V0_opt` is negative because `phi_opt` is different from `phi_init` by π; otherwise, all of the best-fit parameter values are within one standard error of the associated simulation input parameter. We show the fit and the residuals as separate plots here because it is usually better to inspect them both at full scale before combining them in a multipanel figure." |
| 486 | + "These parameters do not fit the simulated data very well, but they are close enough for the solver to find a best-fit solution that is close to the simulation. We take `alpha_V` as the measurement uncertainty and set `absolute_sigma=True` to tell the solver that it should assume the uncertainty is known independently when calculating the parameter uncertainties. The value of `V0_opt` is negative because `phi_opt` is different from `phi_init` by π; otherwise, all the best-fit parameter values are within one standard error of the associated simulation input parameter. We show the fit and the residuals as separate plots here because it is usually better to inspect them both at full scale before combining them in a multipanel figure." |
485 | 487 | ] |
486 | 488 | }, |
487 | 489 | { |
|
493 | 495 | "# Fit the model to the data with curve_fit\n", |
494 | 496 | "from scipy.optimize import curve_fit\n", |
495 | 497 | "\n", |
496 | | - "pOpt, pCov = curve_fit(LCR_model, t, V_sim, \n", |
| 498 | + "pOpt, pCov = curve_fit(lcr_model, t, V_sim, \n", |
497 | 499 | " p0=[V_bgd_init, V0_init, T_init, phi_init, tau_init],\n", |
498 | 500 | " sigma=alpha_V*np.ones(t.size),\n", |
499 | 501 | " absolute_sigma=True)\n", |
|
517 | 519 | "alpha_tau = alpha_vec[4]\n", |
518 | 520 | "\n", |
519 | 521 | "# Show optimized fit parameters and uncertainties\n", |
520 | | - "print(\"V = {0:.4f} ± {1:.4f}\".format(V_bgd_opt, alpha_V_bgd))\n", |
521 | | - "print(\"V0 = {0:.2f} ± {1:.2f}\".format(V0_opt, alpha_V0))\n", |
522 | | - "print(\"T = {0:.3f} ± {1:.3f}\".format(T_opt, alpha_T))\n", |
523 | | - "print(\"phi = {0:.3f} ± {1:.3f}\".format(phi_opt, alpha_phi))\n", |
524 | | - "print(\"tau = {0:.1f} ± {1:.1f}\".format(tau_opt, alpha_tau))\n", |
| 522 | + "print(f\"V = {V_bgd_opt:.4f} ± {alpha_V_bgd:.4f}\")\n", |
| 523 | + "print(f\"V0 = {V0_opt:.2f} ± {alpha_V0:.2f}\")\n", |
| 524 | + "print(f\"T = {T_opt:.3f} ± {alpha_T:.3f}\")\n", |
| 525 | + "print(f\"phi = {phi_opt:.3f} ± {alpha_phi:.3f}\")\n", |
| 526 | + "print(f\"tau = {tau_opt:.1f} ± {alpha_tau:.1f}\")\n", |
525 | 527 | "\n", |
526 | 528 | "# Plot data\n", |
527 | 529 | "plt.plot(t, V_sim, 'k.')\n", |
528 | | - "plt.plot(t, LCR_model(t, V_bgd_opt, V0_opt, T_opt, phi_opt, tau_opt), 'r-')\n", |
| 530 | + "plt.plot(t, lcr_model(t, V_bgd_opt, V0_opt, T_opt, phi_opt, tau_opt), 'r-')\n", |
529 | 531 | "plt.xlabel('Time (µs)')\n", |
530 | 532 | "plt.ylabel('Voltage (mV)')\n", |
531 | 533 | "plt.show()\n", |
532 | 534 | "\n", |
533 | 535 | "# Plot the residuals in a separate figure\n", |
534 | | - "plt.plot(t, (V_sim - LCR_model(t, V_bgd_opt, V0_opt, T_opt, phi_opt, tau_opt))/alpha_V, 'k.')\n", |
| 536 | + "plt.plot(t, (V_sim - lcr_model(t, V_bgd_opt, V0_opt, T_opt, phi_opt, tau_opt))/alpha_V, 'k.')\n", |
535 | 537 | "plt.xlabel('Time (µs)')\n", |
536 | 538 | "plt.ylabel(f'Normalized\\nresidual (mV)')\n", |
537 | 539 | "plt.show()" |
|
0 commit comments