3

Добрый день! Мне нужно реализовать один из тестов Кнута. Выбрала критерий сериальной корреляции. Во втором томе Кнута дается формула

формула

Я пытаюсь по ней вычислить этот коэффициент корреляции и он у меня получается при N=50 равный 1,04. Но он не может быть, по модулю, больше единицы!... Помогите, пожалуйста, разобраться...


Кусочек кода:

const int M=N;
unsigned int null_num, first_num, numerator_n, numerator_2, denominator_n, denominator_2;
long double C;
null_num = (a*x+c)%4294967296;
numerator_n=0;
numerator_2 = null_num;
denominator_n = null_num*null_num;

for (i=0; i<M; i++)
{
    first_num = (a*null_num+c)%4294967296;
    numerator_n = numerator_n + null_num*first_num;
    numerator_2 = numerator_2 + first_num;
    denominator_n = denominator_n + first_num*first_num;
    null_num = first_num;
}

numerator_n = numerator_n + null_num*((a*x+c)%4294967296);
numerator_2 = numerator_2 + null_num;
denominator_2 = numerator_2;
denominator_n = denominator_n + null_num*null_num; 
C = (double)((numerator_n*M)-(numerator_2*numerator_2))/(double)((denominator_n*M)-(denominator_2*denominator_2));

printf("Coefficient: %f \n\n", C);
7
  • язык какой? откуда 4294967296?
    – Gedweb
    11 фев 2013 в 21:33
  • 4294967296 - 2^32 четыре байта, 2^(8*sizeof(int))
    – zb'
    11 фев 2013 в 21:41
  • а в формуле он откуда?
    – VladD
    11 фев 2013 в 22:42
  • это только автору ведомо, в формуле никаких делений с остатком нет, думаю защита от переполнения какая-то но я не сишник
    – zb'
    11 фев 2013 в 22:48
  • @777Julia777 значения и типы a,x,c предоставте
    – zb'
    11 фев 2013 в 22:58

1 ответ 1

5

Ваш тест неверно реализован.

Делайте что-то такое:

double sumUV = 0, sumU = 0, sumV = 0, sumU2 = 0, sum V2 = 0;
for (int i = 0; i < n; i++)
{
    double u = GetNextU(i);
    double v = GetNextV(i);

    sumU += u; sumV += v;
    sumU2 += u*u; sumV2 += v*v;
    sumUV += u*v;
}
double coeff = (n * sumUV - sumU * sumV) /
          sqrt((n * sumU2 - sumU * sumU) * (n * sumV2 - sumV * sumV));

Это, к сожалению, достаточно наивный код: из-за неточности при округлении результат может оказаться катастрофически неточен. О подробностях читайте в Кнуте, 4.2.2. Пример после формулы (15) показывает, что при вычислении корня в знаменателе вы можете столкнуться с отрицательным числом, несмотря на то, что это математически невозможно.


Апдейт. С новой формулой у вас V_i = U_{i+1}, вы считаете корреляцию последовательности на себя. Код можно упростить.

double sumU = 0, sumU2 = 0, sumUShifted = 0;
double firstU = GetNextU(0);
double prevU = firstU;
for (int i = 0; i < n; i++)
{
    double nextU = i == (n - 1) ? firstU : GetNextU(i + 1);

    sumU += prevU;
    sumU2 += prevU*prevU;
    sumUShifted += prevU*nextU;

    prevU = nextU;
}
double coeff = (n * sumUShifted - sumU * sumU) / (n * sumU2 - sumU * sumU);
7
  • ну от и ни каких 4294967296
    – Gedweb
    12 фев 2013 в 6:15
  • @GedWeb: 4294967296 уходит в реализацию GetNextU, но там, четно говоря, можно и без неё.
    – VladD
    12 фев 2013 в 10:11
  • 1
    @VladD, а какой тут смысл в auto? Все равно ведь все дальше в double считается. Можно потребовать, чтобы GetNextU() возвращала double. Но, IMHO это не спасет ТС от потери/переполнения при N=50. Недаром приличная часть второго тома отводится вычислениям в рациональных числах.
    – avp
    12 фев 2013 в 16:31
  • @avp: вы правы, сейчас почищу. У ТС они были целыми, но если завести переменную как double, будет правильнее.
    – VladD
    12 фев 2013 в 16:42
  • 2
    @777Julia777: отлично, удачи вам! Но всё же берегитесь накопления неточностей.
    – VladD
    12 фев 2013 в 19:30

Ваш ответ

By clicking “Отправить ответ”, you agree to our terms of service and acknowledge you have read our privacy policy.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.