diff --git a/demos/log_normal_wiener_filter.py b/demos/log_normal_wiener_filter.py index 5ab513539d01bb135d5bd95d8c959245cfe73317..5a4eabc8906ed862686ca628d64854e4543ceea4 100644 --- a/demos/log_normal_wiener_filter.py +++ b/demos/log_normal_wiener_filter.py @@ -7,7 +7,7 @@ if __name__ == "__main__": correlation_length = 1. # Typical distance over which the field is correlated field_variance = 2. # Variance of field in position space response_sigma = 0.02 # Smoothing length of response (in same unit as L) - signal_to_noise = 100 # The signal to noise ratio + signal_to_noise = 100 # The signal to noise ratio np.random.seed(43) # Fixing the random seed def power_spectrum(k): # Defining the power spectrum @@ -20,66 +20,65 @@ if __name__ == "__main__": # signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels) signal_space = ift.HPSpace(16) harmonic_space = signal_space.get_default_codomain() - fft = ift.FFTOperator(harmonic_space, target=signal_space) + ht = ift.HarmonicTransformOperator(harmonic_space, target=signal_space) power_space = ift.PowerSpace(harmonic_space) # Creating the mock signal S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum) mock_power = ift.PS_field(power_space, power_spectrum) - mock_signal = fft(ift.power_synthesize(mock_power, real_signal=True)) + mock_signal = ift.power_synthesize(mock_power, real_signal=True) # Setting up an exemplary response mask = ift.Field.ones(signal_space) N10 = int(N_pixels/10) # mask.val[N10*5:N10*9, N10*5:N10*9] = 0. - R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), - exposure=(mask,)) + R = ift.GeometryRemover(signal_space) + R = R*ift.DiagonalOperator(mask) + R = R*ht + R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma) + #R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), + # sensitivity=(mask,)) data_domain = R.target[0] - R_harmonic = R*fft # Setting up the noise covariance and drawing a random noise realization - ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise) - N = ift.DiagonalOperator(ndiag.weight(1)) - noise = ift.Field.from_random(domain=data_domain, random_type='normal', - std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0) - data = R(ift.exp(mock_signal)) + noise + noiseless_data = R(mock_signal) + noise_amplitude = noiseless_data.std()/signal_to_noise + N = ift.DiagonalOperator( + ift.Field.full(data_domain, noise_amplitude**2)) + noise = ift.Field.from_random( + domain=data_domain, random_type='normal', + std=noise_amplitude, mean=0) + data = noiseless_data + noise # Wiener filter m0 = ift.Field.zeros(harmonic_space) ctrl = ift.GradientNormController(tol_abs_gradnorm=0.0001) ctrl2 = ift.GradientNormController(tol_abs_gradnorm=0.1, name="outer") inverter = ift.ConjugateGradient(controller=ctrl) - energy = ift.library.LogNormalWienerFilterEnergy(m0, data, R_harmonic, + energy = ift.library.LogNormalWienerFilterEnergy(m0, data, R, N, S, inverter=inverter) - # minimizer1 = ift.VL_BFGS(controller=ctrl2, max_history_length=20) - minimizer2 = ift.RelaxedNewton(controller=ctrl2) - # minimizer3 = ift.SteepestDescent(controller=ctrl2) + # minimizer = ift.VL_BFGS(controller=ctrl2, max_history_length=20) + minimizer = ift.RelaxedNewton(controller=ctrl2) + # minimizer = ift.SteepestDescent(controller=ctrl2) - # me1 = minimizer1(energy) - me2 = minimizer2(energy) - # me3 = minimizer3(energy) - - # m1 = fft(me1[0].position) - m2 = fft(me2[0].position) - # m3 = fft(me3[0].position) + me = minimizer(energy) + m = ht(me[0].position) # Plotting plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index", "colormap": "Planck-like"} - ift.plot(mock_signal, name="mock_signal.png", **plotdict) + ift.plot(ht(mock_signal), name="mock_signal.png", **plotdict) logdata = np.log(ift.dobj.to_global_data(data.val)).reshape(signal_space.shape) ift.plot(ift.Field(signal_space, val=ift.dobj.from_global_data(logdata)), name="log_of_data.png", **plotdict) - # ift.plot(m1,name='m_LBFGS.png', **plotdict) - ift.plot(m2, name='m_Newton.png', **plotdict) - # ift.plot(m3, name='m_SteepestDescent.png', **plotdict) + ift.plot(m, name='m.png', **plotdict) # Probing the variance class Proby(ift.DiagonalProberMixin, ift.Prober): pass proby = Proby(signal_space, probe_count=1) - proby(lambda z: fft(me2[0].curvature.inverse_times(fft.adjoint_times(z)))) + proby(lambda z: ht(me2[0].curvature.inverse_times(ht.adjoint_times(z)))) sm = ift.FFTSmoothingOperator(signal_space, sigma=0.02) variance = sm(proby.diagonal.weight(-1)) diff --git a/demos/paper_demos/cartesian_wiener_filter.py b/demos/paper_demos/cartesian_wiener_filter.py index 9b4665424be2334a6fa34b622085da6f4ed34000..849d793aa89e93a77bf0f625445ccb2fc9bf760b 100644 --- a/demos/paper_demos/cartesian_wiener_filter.py +++ b/demos/paper_demos/cartesian_wiener_filter.py @@ -2,7 +2,7 @@ import numpy as np import nifty4 as ift if __name__ == "__main__": - signal_to_noise = 1.5 # The signal to noise ratio + signal_to_noise = 0.5 # The signal to noise ratio # Setting up parameters L_1 = 2. # Total side-length of the domain @@ -30,7 +30,7 @@ if __name__ == "__main__": harmonic_domain = ift.DomainTuple.make((harmonic_space_1, harmonic_space_2)) - fft_1 = ift.FFTOperator(harmonic_domain, space=0) + ht_1 = ift.HarmonicTransformOperator(harmonic_domain, space=0) power_space_1 = ift.PowerSpace(harmonic_space_1) mock_power_1 = ift.PS_field(power_space_1, power_spectrum_1) @@ -39,12 +39,12 @@ if __name__ == "__main__": a = 4 * correlation_length_2 * field_variance_2**2 return a / (1 + k * correlation_length_2) ** 2.5 - fft_2 = ift.FFTOperator(mid_domain, space=1) + ht_2 = ift.HarmonicTransformOperator(mid_domain, space=1) power_space_2 = ift.PowerSpace(harmonic_space_2) mock_power_2 = ift.PS_field(power_space_2, power_spectrum_2) - fft = fft_2*fft_1 + ht = ht_2*ht_1 mock_power = ift.Field( (power_space_1, power_space_2), @@ -57,7 +57,7 @@ if __name__ == "__main__": S = ift.DiagonalOperator(diagonal) np.random.seed(10) - mock_signal = fft(ift.power_synthesize(mock_power, real_signal=True)) + mock_signal = ift.power_synthesize(mock_power, real_signal=True) # Setting up a exemplary response N1_10 = int(N_pixels_1/10) @@ -70,49 +70,50 @@ if __name__ == "__main__": mask_2[N2_10*7:N2_10*9] = 0. mask_2 = ift.Field(signal_space_2, ift.dobj.from_global_data(mask_2)) - R = ift.ResponseOperator(signal_domain, - sigma=(response_sigma_1, response_sigma_2), - exposure=(mask_1, mask_2)) + R = ift.GeometryRemover(signal_domain) + R = R*ift.DiagonalOperator(mask_1, signal_domain,spaces=0) + R = R*ift.DiagonalOperator(mask_2, signal_domain,spaces=1) + R = R*ht + R = R * ift.create_harmonic_smoothing_operator(harmonic_domain, 0, + response_sigma_1) + R = R * ift.create_harmonic_smoothing_operator(harmonic_domain, 1, + response_sigma_2) data_domain = R.target - R_harmonic = R*fft + noiseless_data = R(mock_signal) + noise_amplitude = noiseless_data.std()/signal_to_noise # Setting up the noise covariance and drawing a random noise realization - ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise) - N = ift.DiagonalOperator(ndiag.weight(1)) + ndiag = ift.Field.full(data_domain, noise_amplitude**2) + N = ift.DiagonalOperator(ndiag) noise = ift.Field.from_random( domain=data_domain, random_type='normal', - std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0) - data = R(mock_signal) + noise + std=noise_amplitude, mean=0) + data = noiseless_data + noise # Wiener filter - j = R_harmonic.adjoint_times(N.inverse_times(data)) + j = R.adjoint_times(N.inverse_times(data)) ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1) inverter = ift.ConjugateGradient(controller=ctrl) wiener_curvature = ift.library.WienerFilterCurvature( - S=S, N=N, R=R_harmonic, inverter=inverter) + S=S, N=N, R=R, inverter=inverter) m_k = wiener_curvature.inverse_times(j) - m = fft(m_k) + m = ht(m_k) - # Probing the variance - class Proby(ift.DiagonalProberMixin, ift.Prober): - pass - proby = Proby((signal_space_1, signal_space_2), probe_count=1, ncpu=1) - proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z)))) -# sm = SmoothingOperator(signal_space, sigma=0.02) -# variance = sm(proby.diagonal.weight(-1)) - variance = proby.diagonal.weight(-1) - - plot_space = ift.RGSpace((N_pixels_1, N_pixels_2)) - sm = ift.FFTSmoothingOperator(plot_space, sigma=0.03) plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index", "colormap": "Planck-like"} - ift.plot( - ift.log(ift.sqrt(sm(ift.Field(plot_space, val=variance.val.real)))), - name='uncertainty.png', zmin=0., zmax=3., title="Uncertainty map", - **plotdict) - ift.plot(ift.Field(plot_space, val=mock_signal.val.real), - name='mock_signal.png', **plotdict) - ift.plot(ift.Field(plot_space, val=data.val.real), - name='data.png', **plotdict) - ift.plot(ift.Field(plot_space, val=m.val.real), name='map.png', **plotdict) + plot_space = ift.RGSpace((N_pixels_1, N_pixels_2)) + ift.plot(ift.Field(plot_space,val=ht(mock_signal).val), name='mock_signal.png', + **plotdict) + ift.plot(ift.Field(plot_space,val=data.val), name='data.png', **plotdict) + ift.plot(ift.Field(plot_space,val=m.val), name='map.png', **plotdict) + # sampling the uncertainty map + sample_variance = ift.Field.zeros(signal_domain) + sample_mean = ift.Field.zeros(signal_domain) + n_samples = 10 + for i in range(n_samples): + sample = ht(wiener_curvature.generate_posterior_sample()) + m + sample_variance += sample**2 + sample_mean += sample + variance = sample_variance/n_samples - (sample_mean/n_samples)**2 + ift.plot(ift.Field(plot_space, val=ift.sqrt(variance).val), name="uncertainty.png", **plotdict) diff --git a/demos/paper_demos/wiener_filter.py b/demos/paper_demos/wiener_filter.py index b0d682fcd671bb58b34c5b164eab7625661a5166..78cffc3fb7cd0ffbcaa8b98b1104b14c58c3b577 100644 --- a/demos/paper_demos/wiener_filter.py +++ b/demos/paper_demos/wiener_filter.py @@ -6,7 +6,7 @@ if __name__ == "__main__": # Setting up parameters L = 2. # Total side-length of the domain N_pixels = 512 # Grid resolution (pixels per axis) - correlation_length_scale = 1. # Typical distance over which the field is + correlation_length_scale = .2 # Typical distance over which the field is # correlated fluctuation_scale = 2. # Variance of field in position space response_sigma = 0.05 # Smoothing length of response @@ -19,7 +19,7 @@ if __name__ == "__main__": signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels) harmonic_space = signal_space.get_default_codomain() - fft = ift.HarmonicTransformOperator(harmonic_space, target=signal_space) + ht = ift.HarmonicTransformOperator(harmonic_space, target=signal_space) power_space = ift.PowerSpace( harmonic_space, binbounds=ift.PowerSpace.useful_binbounds( harmonic_space, logarithmic=True)) @@ -37,47 +37,43 @@ if __name__ == "__main__": mask = ift.Field(signal_space, ift.dobj.from_global_data(mask)) R = ift.GeometryRemover(signal_space) R = R*ift.DiagonalOperator(mask) - R = R*fft + R = R*ht R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma) data_domain = R.target[0] + noiseless_data = R(mock_signal) + noise_amplitude = noiseless_data.std()/signal_to_noise # Setting up the noise covariance and drawing a random noise realization - ndiag = 1e-8*ift.Field.full(data_domain, fft(mock_signal).var()/signal_to_noise) + ndiag = ift.Field.full(data_domain, noise_amplitude**2) N = ift.DiagonalOperator(ndiag) noise = ift.Field.from_random( domain=data_domain, random_type='normal', - std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0) - data = R(mock_signal) + noise + std=noise_amplitude, mean=0) + data = noiseless_data + noise # Wiener filter j = R.adjoint_times(N.inverse_times(data)) - ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=1e-6) + ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=1e-2) inverter = ift.ConjugateGradient(controller=ctrl) wiener_curvature = ift.library.WienerFilterCurvature( S=S, N=N, R=R, inverter=inverter) m_k = wiener_curvature.inverse_times(j) - m = fft(m_k) + m = ht(m_k) plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index", "colormap": "Planck-like"} - ift.plot(mock_signal, name="mock_signal.png", **plotdict) + ift.plot(ht(mock_signal), name="mock_signal.png", **plotdict) ift.plot(ift.Field(signal_space, val=data.val), name="data.png", **plotdict) ift.plot(m, name="map.png", **plotdict) - # Probing the uncertainty - class Proby(ift.DiagonalProberMixin, ift.Prober): - pass - proby = Proby(harmonic_space, probe_count=1, ncpu=1) - proby(lambda z: wiener_curvature.inverse_times(z)) - sm = ift.FFTSmoothingOperator(signal_space, sigma=0.03) - variance = ift.sqrt(sm(proby.diagonal.weight(-1))) - - # Plotting - plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index", - "colormap": "Planck-like"} - ift.plot(variance, name="uncertainty.png", **plotdict) - ift.plot(mock_signal, name="mock_signal.png", **plotdict) - ift.plot(ift.Field(signal_space, val=data.val), - name="data.png", **plotdict) - ift.plot(m, name="map.png", **plotdict) + # sampling the uncertainty map + sample_variance = ift.Field.zeros(signal_space) + sample_mean = ift.Field.zeros(signal_space) + n_samples = 10 + for i in range(n_samples): + sample = ht(wiener_curvature.generate_posterior_sample()) + m + sample_variance += sample**2 + sample_mean += sample + variance = sample_variance/n_samples - (sample_mean/n_samples)**2 + ift.plot(ift.sqrt(variance), name="uncertainty.png", **plotdict) diff --git a/demos/wiener_filter_easy.py b/demos/wiener_filter_easy.py index 09f946326056ebf2cb3d3ebee96c52b91f8c83bc..2867161e24a63e0bfe5d11581963d4ef86d01b43 100644 --- a/demos/wiener_filter_easy.py +++ b/demos/wiener_filter_easy.py @@ -27,8 +27,8 @@ if __name__ == "__main__": # Set up the geometry s_space = ift.RGSpace([N_pixels, N_pixels], distances=pixel_width) - fft = ift.FFTOperator(s_space) - h_space = fft.target[0] + h_space = s_space.get_default_codomain() + ht = ift.HarmonicTransformOperator(h_space, s_space) p_space = ift.PowerSpace(h_space) # Create mock data @@ -37,19 +37,19 @@ if __name__ == "__main__": sp = ift.PS_field(p_space, pow_spec) sh = ift.power_synthesize(sp, real_signal=True) - ss = fft.inverse_times(sh) - R = ift.FFTSmoothingOperator(s_space, sigma=response_sigma) + R = ht*ift.create_harmonic_smoothing_operator((h_space,),0,response_sigma) - signal_to_noise = 1 - diag = ift.Field(s_space, ss.var()/signal_to_noise).weight(1) - N = ift.DiagonalOperator(diag) + noiseless_data = R(sh) + signal_to_noise = 1. + noise_amplitude = noiseless_data.std()/signal_to_noise + N = ift.DiagonalOperator(ift.Field.full(s_space, noise_amplitude**2)) n = ift.Field.from_random(domain=s_space, random_type='normal', - std=ss.std()/np.sqrt(signal_to_noise), + std=noise_amplitude, mean=0) - d = R(ss) + n + d = noiseless_data + n # Wiener filter @@ -57,8 +57,14 @@ if __name__ == "__main__": IC = ift.GradientNormController(name="inverter", iteration_limit=500, tol_abs_gradnorm=0.1) inverter = ift.ConjugateGradient(controller=IC) - S_inv = fft.adjoint*Sh.inverse*fft - D = (R.adjoint*N.inverse*R + S_inv).inverse + D = (R.adjoint*N.inverse*R + Sh.inverse).inverse # MR FIXME: we can/should provide a preconditioner here as well! D = ift.InversionEnabler(D, inverter) m = D(j) + + plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index", + "colormap": "Planck-like"} + ift.plot(ht(sh), name="mock_signal.png", **plotdict) + ift.plot(ift.Field(s_space, val=d.val), + name="data.png", **plotdict) + ift.plot(ht(m), name="map.png", **plotdict) diff --git a/demos/wiener_filter_via_curvature.py b/demos/wiener_filter_via_curvature.py index 0e3309b748fe63c700f598af1755cab1d45e4270..ce6f1dde5a5af75b7e793599b3d06e875e641370 100644 --- a/demos/wiener_filter_via_curvature.py +++ b/demos/wiener_filter_via_curvature.py @@ -16,14 +16,14 @@ if __name__ == "__main__": # sigma of field in position space sqrt(<|s_x|^2>) field_sigma = 2. * nu.K # smoothing length of response - response_sigma = 0.01*nu.m - # The signal to noise ratio ***CURRENTLY BROKEN*** - signal_to_noise = 70.7 + response_sigma = 0.03*nu.m + # The signal to noise ratio + signal_to_noise = 1 # note that field_variance**2 = a*k_0/4. for this analytic form of power # spectrum def power_spectrum(k): - #RL FIXME: signal_amplitude is not how much signal varies + #RL FIXME: signal_amplitude is not how much signal varies cldim = correlation_length**(2*dimensionality) a = 4/(2*np.pi) * cldim * field_sigma**2 # to be integrated over spherical shells later on @@ -39,7 +39,7 @@ if __name__ == "__main__": signal_space = ift.RGSpace(shape, distances=L/N_pixels) harmonic_space = signal_space.get_default_codomain() - fft = ift.HarmonicTransformOperator(harmonic_space, target=signal_space) + ht = ift.HarmonicTransformOperator(harmonic_space, target=signal_space) power_space = ift.PowerSpace(harmonic_space) # Creating the mock data @@ -53,18 +53,18 @@ if __name__ == "__main__": sensitivity = (1./nu.m)**dimensionality/nu.K R = ift.GeometryRemover(signal_space) R = R*ift.ScalingOperator(sensitivity, signal_space) - R = R*fft + R = R*ht R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma) data_domain = R.target[0] - noise_amplitude = 1./signal_to_noise*field_sigma*sensitivity*((L/N_pixels)**dimensionality) - print "noise amplitude:", noise_amplitude + noiseless_data = R(mock_signal) + noise_amplitude = noiseless_data.std()/signal_to_noise N = ift.DiagonalOperator( ift.Field.full(data_domain, noise_amplitude**2)) noise = ift.Field.from_random( domain=data_domain, random_type='normal', std=noise_amplitude, mean=0) - data = R(mock_signal) + noise + data = noiseless_data + noise # Wiener filter j = R.adjoint_times(N.inverse_times(data)) @@ -75,14 +75,14 @@ if __name__ == "__main__": S=S, N=N, R=R, inverter=inverter) m = wiener_curvature.inverse_times(j) - m_s = fft(m) + m_s = ht(m) sspace2 = ift.RGSpace(shape, distances=L/N_pixels/nu.m) - ift.plot(ift.Field(sspace2, fft(mock_signal).val)/nu.K, name="mock_signal.png") + ift.plot(ift.Field(sspace2, ht(mock_signal).val)/nu.K, name="mock_signal.png") #data = ift.dobj.to_global_data(data.val).reshape(sspace2.shape) #data = ift.Field(sspace2, val=ift.dobj.from_global_data(data)) - ift.plot(ift.Field(sspace2, val=fft(R.adjoint_times(data)).val), name="data.png") - print "msig",np.min(fft(mock_signal).val)/nu.K, np.max(fft(mock_signal).val)/nu.K + ift.plot(ift.Field(sspace2, val=data.val), name="data.png") + print "msig",np.min(ht(mock_signal).val)/nu.K, np.max(ht(mock_signal).val)/nu.K print "map",np.min(m_s.val)/nu.K, np.max(m_s.val)/nu.K ift.plot(ift.Field(sspace2, m_s.val)/nu.K, name="map.png") diff --git a/demos/wiener_filter_via_hamiltonian.py b/demos/wiener_filter_via_hamiltonian.py index 74f5ae90c082dfd8039077e2c19456e54798806b..e3a22db3afbbffeeb51e211093fb61cfa10effd6 100644 --- a/demos/wiener_filter_via_hamiltonian.py +++ b/demos/wiener_filter_via_hamiltonian.py @@ -11,7 +11,7 @@ if __name__ == "__main__": # Define associated harmonic space and harmonic transformation h_space = s_space.get_default_codomain() - fft = ift.HarmonicTransformOperator(h_space, s_space) + ht = ift.HarmonicTransformOperator(h_space, s_space) # Setting up power space p_space = ift.PowerSpace(h_space) @@ -25,27 +25,26 @@ if __name__ == "__main__": # Drawing a sample sh from the prior distribution in harmonic space sp = ift.PS_field(p_space, p_spec) sh = ift.power_synthesize(sp, real_signal=True) - ss = fft.times(sh) # Choosing the measurement instrument - # Instrument = ift.FFTSmoothingOperator(s_space, sigma=0.05) diag = np.ones(s_space.shape) diag[20:80, 20:80] = 0 diag = ift.Field(s_space, ift.dobj.from_global_data(diag)) Instrument = ift.DiagonalOperator(diag) # Adding a harmonic transformation to the instrument - R = Instrument*fft + R = Instrument*ht + noiseless_data = R(sh) signal_to_noise = 1. - ndiag = ift.Field.full(s_space, ss.var()/signal_to_noise) - N = ift.DiagonalOperator(ndiag.weight(1)) + noise_amplitude = noiseless_data.std()/signal_to_noise + N = ift.DiagonalOperator(ift.Field.full(s_space, noise_amplitude**2)) n = ift.Field.from_random(domain=s_space, random_type='normal', - std=ss.std()/np.sqrt(signal_to_noise), + std=noise_amplitude, mean=0) # Creating the mock data - d = R(sh) + n + d = noiseless_data + n j = R.adjoint_times(N.inverse_times(d)) # Choosing the minimization strategy @@ -62,8 +61,8 @@ if __name__ == "__main__": energy, convergence = minimizer(energy) m = energy.position D = energy.curvature - ift.plot(ss, name="signal.png", colormap="Planck-like") - ift.plot(fft(m), name="m.png", colormap="Planck-like") + ift.plot(ht(sh), name="signal.png", colormap="Planck-like") + ift.plot(ht(m), name="m.png", colormap="Planck-like") # sampling the uncertainty map sample_variance = ift.Field.zeros(s_space) @@ -71,7 +70,7 @@ if __name__ == "__main__": n_samples = 50 for i in range(n_samples): - sample = fft(D.generate_posterior_sample() + m) + sample = ht(D.generate_posterior_sample() + m) sample_variance += sample**2 sample_mean += sample sample_mean /= n_samples