/** zipfUTest.cxxtest ---
*
* Copyright (C) 2019, 2020 Linas Vepstas
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License v3 as
* published by the Free Software Foundation and including the exceptions
* at http://opencog.org/wiki/Licenses
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program; if not, write to:
* Free Software Foundation, Inc.,
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*/
#include <climits>
#include <opencog/util/zipf.h>
using namespace opencog;
class zipfUTest : public CxxTest::TestSuite
{
std::random_device rd;
std::mt19937 gen;
public:
// zipfUTest() : gen(rd()) {}
zipfUTest() : gen(2) {}
// Check that the generated distribution is normally distributed.
void verify(const std::vector<size_t>& pdf, size_t ndraw,
int n, double s, double q, double sigma)
{
int len = n + 1;
// Create the expected distribution
double norm = 0.0;
std::vector<double> expect;
expect.resize(len);
for (int i = len-1; 0<i; i--)
{
double term = std::pow (i+q, -s);
expect[i] = term;
norm += term;
}
// Normalize
for (int i = 1; i <len; i++)
expect[i] *= 1.0 / norm;
// Also normalize the pdf.
size_t cnt = 0;
for (int i = 1; i <len; i++)
cnt += pdf[i];
TS_ASSERT_EQUALS (cnt, ndraw);
double fcnt = (double) cnt;
// Expected standard deviation. 1/(2pi sqrt(cnt))
double erms = 1.0 / sqrt(2.0 * 3.14159 * fcnt);
// Take the difference. We expect the difference
// to be zero. Also adjust the variance; we expect
// a normal distribution.
double mean = 0.0;
double meansq = 0.0;
for (int i = 1; i <len; i++)
{
// #define GRAPH
#ifdef GRAPH
printf("%d %g %g %g\n", i, expect[i],
pdf[i]/fcnt, (expect[i] - pdf[i]/fcnt)/erms);
#endif
expect[i] -= pdf[i] / fcnt;
// The tail of the distribution is very accurate.
// The head, not so much. Multiplying by sqrt seems
// to turn the distribution into a normal distribution.
// I don't know the theory; this seems to work, though.
if (abs(q) < 0.1)
expect[i] *= pow(i, 0.5*s);
else
expect[i] *= pow((i+q)/q, 0.5*s);
expect[i] /= erms;
mean += expect[i];
meansq += expect[i] * expect[i];
}
mean *= 1.0 / n;
meansq *= 1.0 / n;
double rms = sqrt(meansq);
printf("Dist: n=%d s=%4.3f q=%4.2f cnt=%zu "
"mean=%6.4g rms=%5.3g\n",
n, s, q, cnt, mean, rms);
TS_ASSERT_LESS_THAN(mean, 2.0);
// Hmm. There is some non-trivial expression for nrms
// that I do not know. Its close to 1.0 but seems to
// depend on s and depend on .. what ? log(n) ???
// TS_ASSERT_DELTA(nrms, 1.0, allow);
// Actually test. We're going to use sigma as the max allowed
// deviation. So six-sigma allows the unit tests to pass most
// of the time.
for (int i = 1; i < len; i++)
{
// Should have approx == 1.0, approximately.
TS_ASSERT_DELTA(expect[i], 0.0, sigma);
}
// #define PRINT_ERROR_DISTRIBUTION
#ifdef PRINT_ERROR_DISTRIBUTION
// Dump the error distribution to a file, for graphing.
std::vector<int> edf;
edf.resize(200);
for (int i = 0; i <200; i++)
edf[i] = 0;
for (int i = 1; i < len; i++)
{
int bin = (100.0 * expect[i] / 6.0) + 100;
if (0>bin) bin=0;
if (199 < bin) bin=199;
edf[bin] += 1;
}
printf ("\n#The error distribution:\n#\n");
for (int i = 0; i <200; i++)
printf("%d %f %d\n", i, 6.0*(i-100.0)/100.0, edf[i]);
#endif // PRINT_ERROR_DISTRIBUTION
}
// Draw from the Zipf distribution per given parameters.
void check_zipf(size_t ndraw, int n, double s, double q, double sigma)
{
zipf_distribution<> zipf(n, s, q);
printf("Init rejection-inversion: n=%d s=%4.3f q=%4.2f draw=%zu\n",
n, s, q, ndraw);
std::vector<size_t> pdf;
pdf.resize(n+1);
for (int i = 0; i <= n; i++)
pdf[i] = 0;
for (size_t s = 0; s < ndraw; s++)
{
unsigned int draw = zipf(gen);
TS_ASSERT_LESS_THAN_EQUALS(1, draw);
TS_ASSERT_LESS_THAN_EQUALS(draw, n);
pdf[draw] = pdf[draw] + 1;
}
verify(pdf, ndraw, n, s, q, sigma);
printf("\n");
}
// Draw from the Zipf distribution per given parameters.
void check_table(size_t ndraw, int n, double s, double q)
{
zipf_table_distribution<> zipf(n, s);
printf("Init table: n=%d s=%4.3f q=%4.2f draw=%zu\n",
n, s, q, ndraw);
std::vector<size_t> pdf;
pdf.resize(n+1);
for (int i = 0; i <= n; i++)
pdf[i] = 0;
for (size_t s = 0; s < ndraw; s++)
{
unsigned int draw = zipf(gen);
TS_ASSERT_LESS_THAN_EQUALS(1, draw);
TS_ASSERT_LESS_THAN_EQUALS(draw, n);
pdf[draw] ++;
}
verify(pdf, ndraw, n, s, q, 6.0);
printf("\n");
}
// Test the Zipf distribution, 300 bins, for exponent s=1.
// This is a pretty small test, as such things go.
void test_zipf()
{
printf("\n");
// Six sigma should do the trick.
double sigma = 7.0;
check_zipf(1623, 30, 1.0, 0.0, sigma);
check_zipf(1623000, 300, 1.0, 0.0, sigma);
check_zipf(1623000, 3000, 1.0, 0.0, sigma);
check_zipf(5623000, 30000, 1.0, 0.0, sigma);
check_zipf(1623000, 30, 0.2, 0.0, sigma);
check_zipf(1623000, 30, 0.6, 0.0, sigma);
check_zipf(1623000, 30, 0.8, 0.0, sigma);
check_zipf(1623000, 30, 1.2, 0.0, sigma);
check_zipf(1623000, 30, 1.5, 0.0, sigma);
check_zipf(1623000, 30, 2.0, 0.0, sigma);
check_zipf(1623000, 30, 3.1, 0.0, sigma);
check_zipf(1623000, 30, 4.5, 0.0, sigma);
check_zipf(1623000, 30, 8.2, 0.0, sigma);
check_zipf(1623000, 30, 1.0 - 1e-3, 0.0, sigma);
check_zipf(1623000, 30, 1.0 - 1e-6, 0.0, sigma);
check_zipf(1623000, 30, 1.0 - 1e-9, 0.0, sigma);
check_zipf(1623000, 30, 1.0 - 1e-12, 0.0, sigma);
check_zipf(1623000, 30, 1.0 + 1e-3, 0.0, sigma);
check_zipf(1623000, 30, 1.0 + 1e-6, 0.0, sigma);
check_zipf(1623000, 30, 1.0 + 1e-9, 0.0, sigma);
check_zipf(1623000, 30, 1.0 + 1e-12, 0.0, sigma);
check_zipf(1623000, 3021, 0.01, 0.0, sigma);
check_zipf(1623000, 3070, 0.53, 0.0, sigma);
check_zipf(1623000, 2999, 1.03, 0.0, sigma);
check_zipf(2623010, 213, 2.31, 0.0, sigma);
check_zipf(1623000, 30, 1, -0.49, sigma);
check_zipf(1623000, 30, 1, -0.4, sigma);
check_zipf(1623000, 30, 1, -0.1, sigma);
check_zipf(1623000, 30, 1, 0.1, sigma);
check_zipf(1623000, 30, 1.1, -0.49, sigma);
check_zipf(1623000, 30, 1.1, -0.4, sigma);
check_zipf(1623000, 30, 1.1, -0.1, sigma);
check_zipf(1623000, 30, 1.1, 0.1, sigma);
// I don't understand the theory of the error
// distribution for q != 0 but this seems to work.
check_zipf(1623000, 30, 1.2, 0.4, 1.2*sigma);
check_zipf(1623000, 30, 1.5, 0.4, 1.5*sigma);
check_zipf(1623000, 30, 2.0, 0.4, 1.5*sigma);
check_zipf(1623000, 30, 3.1, 0.4, 3.1*sigma);
check_zipf(1623000, 30, 4.5, 0.4, 4.5*sigma);
check_zipf(1623000, 30, 8.2, 0.4, 8.2*sigma);
check_zipf(1623000, 30, 1.2, 4, sigma);
check_zipf(1623000, 30, 1.5, 4, sigma);
check_zipf(1623000, 30, 2.0, 4, sigma);
check_zipf(1623000, 30, 3.1, 4, sigma);
check_zipf(1623000, 30, 4.5, 4, 9.5);
check_zipf(1623000, 30, 8.2, 4, 24.0);
check_zipf(1623000, 30, 1.2, 44, sigma);
check_zipf(1623000, 30, 1.5, 44, sigma);
check_zipf(1623000, 30, 2.0, 44, sigma);
check_zipf(1623000, 30, 3.1, 44, sigma);
check_zipf(1623000, 30, 4.5, 44, sigma);
check_zipf(1623000, 30, 8.2, 44, sigma);
check_zipf(1623000, 30, 1.2, 44e4, sigma);
check_zipf(1623000, 30, 1.5, 44e4, sigma);
check_zipf(1623000, 30, 2.0, 44e4, sigma);
check_zipf(1623000, 30, 3.1, 44e4, sigma);
check_zipf(1623000, 30, 4.5, 44e4, sigma);
check_zipf(1623000, 30, 8.2, 44e4, sigma);
check_table(1623000, 310, 1.0, 0.0);
check_table(1623000, 312, 1.2, 0.0);
check_table(1623000, 217, 0.8, 0.0);
}
};