.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/check_integration_error.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_check_integration_error.py: Check integral error calculation ================================ Generate a fake dataset of an inversion recovery with multiple repeats (φ × t2 × vd × repeats) w/ normally distributed random noise. Check that the following match: - integral w/ error (the canned routine :func:`~pyspecProcScripts.integral_w_errors`) - propagate error based off the programmed σ of the normal distribution - set the error bars based on the standard deviation (along the repeats dimension) of the *real* part of the integral - propagate error based off the variance of the noise in the inactive coherence channels (do this manually inside this script -- should mimic what :func:`~pyspecProcScripts.integral_w_errors` does) .. GENERATED FROM PYTHON SOURCE LINES 18-131 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/images/sphx_glr_check_integration_error_001.png :alt: integration diagnostic :srcset: /auto_examples/images/sphx_glr_check_integration_error_001.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_check_integration_error_002.png :alt: matched filter diagnostic -- signal Energy :srcset: /auto_examples/images/sphx_glr_check_integration_error_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_check_integration_error_003.png :alt: matched filter diagnostic -- time domain :srcset: /auto_examples/images/sphx_glr_check_integration_error_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_check_integration_error_004.png :alt: different types of error :srcset: /auto_examples/images/sphx_glr_check_integration_error_004.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none ---------- logging output to /home/jmfranck/pyspecdata.0.log ---------- shape of all results [(40, 'vd'), (100, 'repeats')] #0 #1 #2 #3 #4 #5 #6 #7 #8 #9 #10 #11 #12 #13 #14 #15 #16 #17 #18 #19 #20 #21 #22 #23 #24 #25 #26 #27 #28 #29 #30 #31 #32 #33 #34 #35 #36 #37 #38 #39 #40 #41 #42 #43 #44 #45 #46 #47 #48 #49 #50 #51 #52 #53 #54 #55 #56 #57 #58 #59 #60 #61 #62 #63 #64 #65 #66 #67 #68 #69 #70 #71 #72 #73 #74 #75 #76 #77 #78 #79 #80 #81 #82 #83 #84 #85 #86 #87 #88 #89 #90 #91 #92 #93 #94 #95 #96 #97 #98 #99 off-pathway std array([1.97220022, 1.92375672, 1.97689265, 1.98747806, 1.90421805, 2.10511932, 2.04874192, 1.97321004, 1.88565543, 2.04797383, 2.00818674, 2.06696218, 2.02617466, 2.02800179, 1.8381604 , 1.9537043 , 2.09145914, 1.98115531, 1.94782724, 1.96258195, 1.85835747, 1.94830219, 2.03456705, 2.05711269, 2.05529398, 1.94220362, 1.94225146, 2.09010836, 2.00995428, 2.10428952, 2.02011849, 1.98634654, 1.9664804 , 1.92916986, 2.00990651, 2.0238763 , 1.86914168, 1.92370829, 2.02447457, 1.99875205]) dimlabels=['vd'] axes={`vd':array([0. , 0.02564103, 0.05128205, 0.07692308, 0.1025641 , 0.12820513, 0.15384615, 0.17948718, 0.20512821, 0.23076923, 0.25641026, 0.28205128, 0.30769231, 0.33333333, 0.35897436, 0.38461538, 0.41025641, 0.43589744, 0.46153846, 0.48717949, 0.51282051, 0.53846154, 0.56410256, 0.58974359, 0.61538462, 0.64102564, 0.66666667, 0.69230769, 0.71794872, 0.74358974, 0.76923077, 0.79487179, 0.82051282, 0.84615385, 0.87179487, 0.8974359 , 0.92307692, 0.94871795, 0.97435897, 1. ]) +/-None} programmed std 2.0 /home/jmfranck/git_repos/pyspecdata/pyspecdata/figlist.py:782: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all axes decorations. plt.gcf().tight_layout() | .. code-block:: Python from numpy import diff, r_, sqrt, real, exp, pi from pyspecdata import ndshape, nddata, init_logging, figlist_var from pyspecProcScripts import integral_w_errors # sphinx_gallery_thumbnail_number = 1 init_logging(level="debug") fl = figlist_var() t2 = nddata(r_[0:1:1024j], "t2") vd = nddata(r_[0:1:40j], "vd") ph1 = nddata(r_[0, 2] / 4.0, "ph1") ph2 = nddata(r_[0:4] / 4.0, "ph2") signal_pathway = {"ph1": 0, "ph2": 1} excluded_pathways = [(0, 0), (0, 3)] # this generates fake clean_data w/ a T₂ of 0.2s # amplitude of 21, just to pick a random amplitude # offset of 300 Hz, FWHM 10 Hz clean_data = ( 21 * (1 - 2 * exp(-vd / 0.2)) * exp(+1j * 2 * pi * 100 * t2 - t2 * 10 * pi) ) clean_data *= exp(signal_pathway["ph1"] * 1j * 2 * pi * ph1) clean_data *= exp(signal_pathway["ph2"] * 1j * 2 * pi * ph2) clean_data["t2":0] *= 0.5 fake_data_noise_std = 2.0 clean_data.reorder(["ph1", "ph2", "vd"]) bounds = (0, 200) # seem reasonable to me result = 0 n_repeats = 100 all_results = ndshape(clean_data) + (n_repeats, "repeats") all_results.pop("t2").pop("ph1").pop("ph2") all_results = all_results.alloc() all_results.setaxis("vd", clean_data.getaxis("vd")) print("shape of all results", ndshape(all_results)) for j in range(n_repeats): data = clean_data.C data.add_noise(fake_data_noise_std) # at this point, the fake data has been generated data.ft(["ph1", "ph2"]) # {{{ usually, we don't use a unitary FT -- this makes it unitary data /= 0.5 * 0.25 # the dt in the integral for both dims data /= sqrt(ndshape(data)["ph1"] * ndshape(data)["ph2"]) # normalization # }}} dt = diff(data.getaxis("t2")[r_[0, 1]]).item() data.ft("t2", shift=True) # {{{ data /= sqrt(ndshape(data)["t2"]) * dt error_pathway = ( set(( (j, k) for j in range(ndshape(data)["ph1"]) for k in range(ndshape(data)["ph2"]) )) - set(excluded_pathways) - set([(signal_pathway["ph1"], signal_pathway["ph2"])]) ) error_pathway = [{"ph1": j, "ph2": k} for j, k in error_pathway] s_int, frq_slice = integral_w_errors( data, signal_pathway, error_pathway, indirect="vd", fl=fl, return_frq_slice=True, ) # }}} manual_bounds = data["ph1", 0]["ph2", 1]["t2":frq_slice] N = ndshape(manual_bounds)["t2"] df = diff(data.getaxis("t2")[r_[0, 1]]).item() manual_bounds.integrate("t2") # N terms that have variance given by fake_data_noise_std**2 each # multiplied by df all_results["repeats", j] = manual_bounds print("#%d" % j) std_off_pathway = ( data["ph1", 0]["ph2", 0]["t2":bounds] .C.run(lambda x: abs(x) ** 2 / 2) # sqrt2 so variance is variance of real .mean_all_but(["t2", "vd"]) .mean("t2") .run(sqrt) ) print( "off-pathway std", std_off_pathway, "programmed std", fake_data_noise_std ) propagated_variance_from_inactive = N * df**2 * std_off_pathway**2 # removed factor of 2 in following, which shouldn't have been there propagated_variance = N * df**2 * fake_data_noise_std**2 fl.next("different types of error") fl.plot(s_int, ".", capsize=6, label="std from int w err", alpha=0.5) manual_bounds.set_error(sqrt(propagated_variance)) fl.plot( manual_bounds, ".", capsize=6, label=r"propagated from programmed variance", alpha=0.5, ) all_results.run(real).mean("repeats", std=True) # by itself, that would give error bars, but the data would be # averaged -- better to put the data in the same position manual_bounds.set_error(all_results.get_error()) # the fact that this matches the previous shows that my sample size is # large enough to give good statistics fl.plot(manual_bounds, ".", capsize=6, label=r"std from repeats", alpha=0.5) manual_bounds.set_error(sqrt(propagated_variance_from_inactive.data)) fl.plot( manual_bounds, ".", capsize=6, label=r"propagated from inactive std", alpha=0.5, ) fl.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 50.853 seconds) .. _sphx_glr_download_auto_examples_check_integration_error.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: check_integration_error.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: check_integration_error.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_