(Это nanswer , т.е. не является ответом как таковым. Первоначально он был написан после того, как @con переписал свой вопрос, чтобы включить весь их код, и был шагом к моему Третьему и последнему ответу . Теперь, надеюсь, это полезный ресурс для тех, кто изучает Perl 6.)
Оригинальное вступление к моей первой редакции этого nanswer
Это много кода! :)
До сих пор я думал: вы абсолютно уверены, что на 100% уверены, что не пропустили некоторые входные данные? Кажется, что вероятность того, что вы пропустили данные, гораздо выше, чем у P6, особенно если учесть, что рассчитанное значение является именно тем, которое вы ожидаете получить от следующего правильного результата.
( Обновление Действительно, выяснилось, что проблема была в неправильных входных данных.)
Перевод с комментариями рассматриваемого кода
Остальная часть этого nanswer - построчная «очистка» (не рефакторинг) кода в вопросе. Я имел в виду две цели, вторая из которых является самой важной, поскольку вы, дорогой читатель, читаете это:
Мой перевод эффективно доказал мне и продемонстрировал @con, что я учел каждый бит их кода . Это было сделано, чтобы уменьшить неопределенность в том, где может быть ошибка. Обратите внимание, что большинство моих изменений не могли иметь непосредственного отношения к их ошибке, но я не чувствовал себя комфортно, предполагая что-либо, пока я не переписал.
@ con и мой перевод, вероятно, могут быть полезными для любого, кто изучает Perl 6 . Код @ con является переводом P6 кода Perl 5 . Код P5, в свою очередь, является переводом C кода . И есть другие переводы этого кода на другие языки. Мой код берет перевод @ con и переводит его в более идиоматическую версию с комментариями, объясняющими, почему я изменил их код.
Мой перевод рассматриваемого кода, запускается через tio
Мой перевод кода @ con, за исключением моего комментария (см. Ниже), в tio .
(Мой первоначальный план состоял в том, что мы будем дополнительно исследовать ошибку @ con, основанную на изменении кода в tio и обмене ссылками на измененные версии (версии) tio. Чтобы сделать это, вы / они просто щелкнете внутри значка tio ссылку в вверху (?), чтобы получить ссылку на код в том виде, в котором она была нажата при нажатии на значок ссылки. Tio отлично подходит для возможности писать, редактировать и запускать код, написанный на P6 и других языках, а затем делиться им.)
Мой перевод рассматриваемого кода с комментариями
Я оставил начальную lgamma рутину как есть. Я написал об этом обширный комментарий в сноске, в основном для других читателей, потому что он содержит массу интересных функций 1 :
sub lgamma ( Num(Real) \n --> Num ){
use NativeCall;
sub lgamma (num64 --> num64) is native {}
lgamma( n )
}
Хочу заметить, что код lgamma предназначен для чисел с плавающей запятой (в P6 это типы Num / num...).
sub pvalue (@a, @b) {
Для тех, кто не знает Perls, ключевое слово sub представляет подпрограмму (она же функция).
Этот аргумент принимает два "листовых" ( Positional ) аргумента.
(Когда вы видите @ (либо как "сигил" , например, как в @foo, либо как оператор, например, как в @($bar)), подумайте "список".)
Чтобы ускорить чтение кода и уменьшить количество повторяющихся кодов, я заменил это:
if @a.elems <= 1 {
return 1.0;
}
if @b.elems <= 1 {
return 1.0;
}
с этим:
return 1 if @a | @b <= 1;
Считайте это по-английски как «return 1, если количество элементов в списке @a или списке @b меньше или равно 1».
Я использовал оператор | для создания соединения any .
(Не тратьте время на попытки обдумать теорию о том, что делают узлы и как они это делают. Просто используйте их на практике простым и / или кратким образом. Если они просты, и читайте в простой способ, они великолепны и просто делают очевидные вещи. Если их использование делает код неочевидным в конкретном случае, то подумайте об использовании какой-либо другой конструкции для этого конкретного случая.)
Массив в числовом контексте оценивается по количеству элементов. Числовые операторы типа <= навязывают числовой контекст. Поэтому я бросил .elems.
Числовой контекст и массивы, оценивающие их длину в числовом контексте, являются основными аспектами P6, так что это идиоматическое кодирование, которое подходит для всех, кроме самых основных вводных примеров новичков.
Я переключился с 1.0 на 1.
Я предполагаю, что @con написал 1.0, потому что именно так значение было буквально записано в коде, который они переводили, и / или с намерением, которое представляло значение с плавающей запятой. Но в P6 простые десятичные литералы (без показателя e), такие как 1.0, вместо этого производят рациональные числа.
Использование 1.0 вместо 1 в этом коде заменяет более простой более быстрый тип (целые числа) на более сложный более медленный тип (рациональные числа). Затем значение этого более медленного типа будет принудительно преобразовано (медленнее, чем целое число) в число с плавающей запятой, если оно используется в формулах, в которых любое значение компонента является плавающей запятой. Что произойдет для большинства или всех формул в этом коде, потому что lgamma возвращает число с плавающей запятой.
В целом, P6 работает и читается лучше, если вы оставляете типы переменных и значений неуказанными , оставляя их на усмотрение компилятора, если нет веских причин для их указания. Оставление неопределенных типов снижает когнитивную нагрузку на читателя и повышает гибкость при повторном использовании кода и оптимизации кода компилятором.
Это приводит к дополнительной паре максим:
По умолчанию оставить информацию о типе неявной . Если вы не знаете, имеет ли значение какой тип значения или переменной задано, или знаете, что это не имеет значения, не указывайте это.
Если вы сделаете тип значения, переменной или параметра явным , тогда P6 будет использовать этот тип для принудительного соответствия этому типу, независимо от того, улучшит ли это ваше значение. код или излишне замедляет код или останавливает его вообще.
Если вы абсолютно знаете , что вам нужно использовать другой тип или добавить ограничение типа, чтобы сделать код корректным, тогда все средства продолжаются. Точно так же, если вы уже убедились, что ваш код правильный без более конкретного типа, но вы хотите сделать его быстрее, безопаснее или яснее, продолжайте. Но если вы не знаете, тогда оставьте это общее. P6 был разработан, чтобы делать то, что вы имеете в виду, и будет иметь успех гораздо чаще, чем вы можете себе представить, как однажды сказал Бен.
Чтобы сформулировать этот пункт еще более убедительно, по аналогии с преждевременной оптимизацией , в P6, преждевременная типизация является антишаблоном как для одноразового кода, так и для долгосрочного производственного кода .
my Rat $mean1 = @a.sum / @a.elems;
my Rat $mean2 = @b.sum / @b.elems;
if $mean1 == $mean2 {
return 1.0;
}
стали:
my (\mean_a, \mean_b) = @a.sum / @a, @b.sum / @b;
return 1 if mean_a == mean_b;
Я обычно "режу сигилы" "переменных", если не знаю, что мне понадобится сигил. Если «переменная» на самом деле не меняется, ей, скорее всего, не нужен сигил.
Я переименовал $mean1 в mean_a, потому что он явно соответствует @a.
my Rat $variance1 = 0.0;
my Rat $variance2 = 0.0;
for @a -> $i {
$variance1 += ($mean1 - $i)**2#";" unnecessary for last statement in block
}
for @b -> $i {
$variance2 += ($mean2 - $i)**2
}
стали:
my ($vari_a, $vari_b);
$vari_a += (mean_a - $_)² for @a;
$vari_b += (mean_b - $_)² for @b;
Считайте переменную $_ как "это".
if ($variance1 == 0 && $variance2 == 0) {
return 1.0;
}
$variance1 /= (@a.elems - 1);
$variance2 /= (@b.elems - 1);
стали:
return 1 unless ($vari_a or $vari_b);
$vari_a /= (@a - 1);
$vari_b /= (@b - 1);
В Perls 0 оценивается в логическом тестовом контексте как False. И массив или список, такой как @a, вычисляется в числовом контексте по количеству элементов.
my $WELCH_T_STATISTIC = ($mean1-$mean2)/sqrt($variance1/@a.elems+$variance2/@b.elems);
my $DEGREES_OF_FREEDOM = (($variance1/@a.elems+$variance2/@b.elems)**2)
/
(
($variance1*$variance1)/(@a.elems*@a.elems*(@a.elems-1))+
($variance2*$variance2)/(@b.elems*@b.elems*(@b.elems-1))
);
my $A = $DEGREES_OF_FREEDOM/2;
my $value = $DEGREES_OF_FREEDOM/($WELCH_T_STATISTIC*$WELCH_T_STATISTIC+$DEGREES_OF_FREEDOM);
my Num $beta = lgamma($A)+0.57236494292470009-lgamma($A+0.5);
my Rat $acu = 10**(-15);
my ($ai,$cx,$indx,$ns,$pp,$psq,$qq,$rx,$temp,$term,$xx);
# Check the input arguments.
return $value if $A <= 0.0;# || $q <= 0.0;
стали:
my \WELCH_T_STATISTIC = (mean_a - mean_b)
/ ( $vari_a / @a + $vari_b / @b ).sqrt;
my \DEGREES_OF_FREEDOM = ($vari_a / @a + $vari_b / @b)²
/ ($vari_a² / (@a² * (@a - 1)) + $vari_b² / (@b² * (@b - 1)));
my \A = DEGREES_OF_FREEDOM / 2;
my $value = DEGREES_OF_FREEDOM
/ (WELCH_T_STATISTIC² + DEGREES_OF_FREEDOM);
my \magic-num = 0.57236494292470009;
my \beta = lgamma(A) + magic-num - lgamma(A + 0.5);
my \acu = 1e-15;
# Check the input arguments.
return $value if A <= 0;# || $q <= 0.0;
(Что это за $q в приведенной выше строке?)
Обратите внимание, что я сбросил типы на beta и acu. Значение , присвоенное переменной beta , в любом случае будет Num, поскольку lgamma возвращает Num. Значение, присвоенное acu, также будет Num, потому что использование показателя e в числовом литерале означает, что значение, которое он создает, является Num.
return $value if $value < 0.0 || 1.0 < $value;
# Special cases
return $value if $value == 0.0 || $value == 1.0;
стали:
return $value unless $value ~~ 0^..^1;
Считайте 0^ как «выше нуля» (не включая ноль) и ^1 как «до одного» (не включая единицу).
~~ является оператором «умного совпадения». Возвращает True, если значение справа принимает значение слева.
Таким образом, этот оператор возврата возвращает $value, если $value меньше или равно 0 или больше или равно 1.
В качестве незначительного приближения я переместил объявление my загрузки переменных, которые я пропустил в вышеупомянутом перезаписи, непосредственно перед тем, как они стали релевантными, и добавил $ns тоже:
my ($ai, $cx, $indx, $pp, $psq, $qq, $rx, $temp, $term, $xx, $ns);
$psq = $A + 0.5;
$cx = 1.0 - $value;
if $A < $psq * $value {
($xx, $cx, $pp, $qq, $indx) = ($cx, $value, 0.5, $A, 1);
} else {
($xx, $pp, $qq, $indx) = ($value, $A, 0.5, 0);
}
стали:
$psq = A + 0.5;
$cx = 1 - $value;
($xx, $cx, $pp, $qq, $indx) =
A < $psq * $value
?? ($cx, $value, 0.5, A, 1)
!! ($value, $cx, A, 0.5, 0);
Я переписал условное присвоение группы переменных как троичный , чтобы было легче увидеть, что назначено для чего.
$term = 1.0;
$ai = 1.0;
$value = 1.0;
$ns = $qq + $cx * $psq;
$ns = $ns.Int;
стали:
$term = 1;
$ai = 1;
$value = 1;
$ns = ($qq + $cx * $psq) .Int;
Замена 1.0 с 1 с и объединение выражения, присвоенного $ns, с принуждением .Int.
(Я удалил типы из кода при переводе, и он продолжал вычислять правильные результаты, за исключением того, что удаление приведенного выше Int принуждения сделало код бесконечным. Это то, что в конечном итоге привело меня к поиску в сети, чтобы посмотреть, смогу ли я найти код, который @con переводил. Именно тогда я нашел его на rosettacode.org. Он был явно напечатан как целое число в коде, который я видел, так что, предположительно, он является центральным для обеспечения работы алгоритма.)
#Soper reduction formula.
$rx = $xx / $cx;
$temp = $qq - $ai;
$rx = $xx if $ns == 0;
(не изменился.)
while (True) {
стали:
loop {
$term = $term * $temp * $rx / ( $pp + $ai );
$value = $value + $term;
$temp = $term.abs;
(не изменился.)
if $temp <= $acu && $temp <= $acu * $value {
$value = $value * ($pp * $xx.log + ($qq - 1.0) * $cx.log - $beta).exp / $pp;
$value = 1.0 - $value if $indx;
last;
}
стали:
if $temp <= acu & acu * $value {
$value = $value * ($pp * $xx.log + ($qq - 1) * $cx.log - beta).exp / $pp;
$value = 1 - $value if $indx;
last;
}
На этот раз условие, содержащее соединение (&), читается на английском языке как «если temp меньше или равно значению времени acu и acu».
$ai++;
$ns--;
if 0 <= $ns {
$temp = $qq - $ai;
$rx = $xx if $ns == 0;
} else {
$temp = $psq;
$psq = $psq + 1;
}
}
return $value;
}
Я только что заменил 1.0 на 1.
Теперь проблемный массив. Как я писал в начале, я почти уверен, что вы (или поставщик ваших данных) просто забыли пару строк:
my @array2d =
[27.5,21.0,19.0,23.6,17.0,17.9,16.9,20.1,21.9,22.6,23.1,19.6,19.0,21.7,21.4],
[27.1,22.0,20.8,23.4,23.4,23.5,25.8,22.0,24.8,20.2,21.9,22.1,22.9,20.5,24.4],
[17.2,20.9,22.6,18.1,21.7,21.4,23.5,24.2,14.7,21.8],
[21.5,22.8,21.0,23.0,21.6,23.6,22.5,20.7,23.4,21.8,20.7,21.7,21.5,22.5,23.6,21.5,22.5,23.5,21.5,21.8],
[19.8,20.4,19.6,17.8,18.5,18.9,18.3,18.9,19.5,22.0],
[28.2,26.6,20.1,23.3,25.2,22.1,17.7,27.6,20.6,13.7,23.2,17.5,20.6,18.0,23.9,21.6,24.3,20.4,24.0,13.2],
[30.02,29.99,30.11,29.97,30.01,29.99],
[29.89,29.93,29.72,29.98,30.02,29.98],
[3.0,4.0,1.0,2.1],
[490.2,340.0,433.9],
[<1.0/15.0>, <10.0/62.0>],
[<1.0/10>, <2/50.0>],
[0.010268,0.000167,0.000167],
[0.159258,0.136278,0.122389],
[9/23.0,21/45.0,0/38.0],
[0/44.0,42/94.0,0/22.0];
say @array2d[11][0];
Кто или что говорит, что это правильные ответы? Вы на 100% уверены, что ответ 0.0033... соответствует данным [<1.0/15.0>, <10.0/62.0>],[<1.0/10>, <2/50.0>]?
my @CORRECT_ANSWERS =
0.021378001462867,
0.148841696605327,
0.0359722710297968,
0.090773324285671,
0.0107515611497845,
0.00339907162713746,
0.52726574965384,
0.545266866977794;
И в последнем бите я просто удалил типы и снова использовал верхний индекс для более привлекательного значения плюс показатель (10⁻⁹):
my $i = 0;
my $error = 0;
for @array2d -> @left, @right {
my $pvalue = pvalue(@left, @right);
$error += ($pvalue - @CORRECT_ANSWERS[$i]).abs;
say "$i [" ~ @left.join(',') ~ '] [' ~ @right ~ "] = $pvalue";
if $error > 10⁻⁹ {
say "\$p = $pvalue, but should be @CORRECT_ANSWERS[$i]";
die;
}
# printf("Test sets %u p-value = %.14g\n",$i+1,$pvalue);
$i++
}
printf("the cumulative error is %g\n", $error);
Сноска
1 Люблю эту милую lgamma рутину! (Оказывается Брэд Гилберт написал это .)
sub lgamma ( Num(Real) \n --> Num ){
use NativeCall;
sub lgamma (num64 --> num64) is native {}
lgamma( n )
}
Showcasing:
Слежение за процедурой низкого уровня (C) с процедурой высокого уровня P6 (Perl 6) с тем же именем, чтобы добавить дополнительную обработку перед вызовом функции C.
Явно преобразует тип ввода из более широких чисел типа P6 Real в более узкие Num. Perls «строго типизированы» в соответствии с первоначальным техническим определением термина, но P6 предоставляет дополнительные опции для нескольких других интерпретаций «строгой типизации» и относительно языки менее способны к этим другим интерпретациям "строгой типизации" (например, Perl 5, Python и C). Явное преобразование типов является частью этого изменения в возможностях, которые вводит P6.
Уничтожить сигилу. Это далее обсуждается, где я делал то же самое в другом месте в этом посте.
Лексически оценивая use библиотеки. Внутренняя подпрограмма lgamma и символы, импортированные use Nativecall;, не видны нигде, кроме внутренней подпрограммы lgamma, содержащей ее.
Использование NativeCall (P6 C FFI ), чтобы позволить высокоуровневому коду P6 безупречно отображаться непосредственно в код C, включая автоматическое преобразование из I6E, совместимого с двойным плавающим числом P6 в штучной упаковке Num в распакованном виде машинный тип данных эквивалент num64.
Все в 5 строках! Очень хорошо. :)