diff --git a/cara/tests/test_sampleable_distribution.py b/cara/tests/test_sampleable_distribution.py index b3887763..e3801c7e 100644 --- a/cara/tests/test_sampleable_distribution.py +++ b/cara/tests/test_sampleable_distribution.py @@ -3,8 +3,8 @@ import numpy.testing as npt import pytest from cara.monte_carlo import sampleable -# TODO: seed deterministically the random number generators, here (to -# avoid random issues with tests) +# TODO: seed better the random number generators +np.random.seed(2000) @pytest.mark.parametrize( "mean, std",[ @@ -17,12 +17,14 @@ def test_normal(mean, std): sample_size = 2000000 samples = sampleable.Normal(mean, std).generate_samples(sample_size) histogram, bins = np.histogram(samples,bins=100, density=True) - x = (bins[1:]+bins[:-1])/2 - exact_dist = 1/(np.sqrt(2*np.pi)*std) * np.exp(-((x-mean)/std)**2/2) + selected_bins,selected_histogram = zip(*[(b,h) for b,h in zip( + (bins[1:]+bins[:-1])/2,histogram) if b>=0.25 and b<=1.75]) + exact_dist = 1/(np.sqrt(2*np.pi)*std) * np.exp( + -((np.array(selected_bins)-mean)/std)**2/2) assert len(samples) == sample_size - npt.assert_allclose([samples.mean(), samples.std()], [mean, std], atol=mean/100) - npt.assert_allclose(histogram, exact_dist, atol=exact_dist.mean()/20) + npt.assert_allclose([samples.mean(), samples.std()], [mean, std], rtol=0.01) + npt.assert_allclose(selected_histogram, exact_dist, rtol=0.02) @pytest.mark.parametrize( @@ -37,17 +39,19 @@ def test_lognormal(mean_gaussian, std_gaussian): samples = sampleable.LogNormal(mean_gaussian, std_gaussian ).generate_samples(sample_size) histogram, bins = np.histogram(samples,bins=50, density=True) - x = (bins[1:]+bins[:-1])/2 - exact_dist = ( 1/(x*np.sqrt(2*np.pi)*std_gaussian) * - np.exp(-((np.log(x)-mean_gaussian)/std_gaussian)**2/2) ) + selected_bins,selected_histogram = zip(*[(b,h) for b,h in zip( + (bins[1:]+bins[:-1])/2,histogram) if b>=0.4 and b<=0.6]) + selected_bins = np.array(selected_bins) + exact_dist = ( 1/(selected_bins*np.sqrt(2*np.pi)*std_gaussian) * + np.exp(-((np.log(selected_bins)-mean_gaussian)/std_gaussian)**2/2) ) exact_mean = np.exp(mean_gaussian + std_gaussian**2/2) exact_std = np.sqrt( (np.exp(std_gaussian**2)-1) * np.exp(2*mean_gaussian + std_gaussian**2) ) assert len(samples) == sample_size npt.assert_allclose([samples.mean(), samples.std()], - [exact_mean, exact_std], atol=exact_mean/100) - npt.assert_allclose(histogram, exact_dist, atol=exact_dist.mean()/20) + [exact_mean, exact_std], rtol=0.01) + npt.assert_allclose(selected_histogram, exact_dist, rtol=0.02) @pytest.mark.parametrize( @@ -71,13 +75,36 @@ def test_custom(use_kernel): samples = sampleable.CustomKernel(variable, frequencies, kernel_bandwidth=0.1 ).generate_samples(sample_size) - tolerance = max_function/10 else: samples = sampleable.Custom((0, 10), function, max_function ).generate_samples(sample_size) - tolerance = max_function/50 histogram, bins = np.histogram(samples, bins=100, density=True) - correct_dist = function((bins[1:]+bins[:-1])/2) + selected_bins,selected_histogram = zip(*[(b,h) for b,h in zip( + (bins[1:]+bins[:-1])/2,histogram) if b>=1 and b<=9]) + correct_dist = function(np.array(selected_bins)) assert len(samples) == sample_size - npt.assert_allclose(histogram, correct_dist, atol=tolerance) + npt.assert_allclose(selected_histogram, correct_dist, rtol=0.05) + + +def test_logcustomkernel(): + # test that the sample has approximately the right distribution + # function, for the LogCustomKernel. + # the distribution function is an inverted parabola vs. the log of + # the variable (normalized) + norm = 500/3. + function = lambda x: (-(5 - x)**2 + 25)/norm + sample_size = 2000000 + + log_variable = np.linspace(0.1,9.9,100) + frequencies = function(log_variable) + samples = sampleable.LogCustomKernel(log_variable, frequencies, + kernel_bandwidth=0.1 + ).generate_samples(sample_size) + + histogram, bins = np.histogram(np.log10(samples), bins=100, density=True) + selected_bins,selected_histogram = zip(*[(b,h) for b,h in zip( + (bins[1:]+bins[:-1])/2,histogram) if b>=1 and b<=9]) + correct_dist = function(np.array(selected_bins)) + assert len(samples) == sample_size + npt.assert_allclose(selected_histogram, correct_dist, rtol=0.05)