我正在从Perl5学习Perl6.
为了进行编译,我将发布整个程序:
sub lgamma ( Num(Real) \n --> Num ){
use NativeCall;
sub lgamma (num64 --> num64) is native {}
lgamma( n )
}
sub pvalue (@a, @b) {
if @a.elems <= 1 {
return 1.0;
}
if @b.elems <= 1 {
return 1.0;
}
my Rat $mean1 = @a.sum / @a.elems;
my Rat $mean2 = @b.sum / @b.elems;
if $mean1 == $mean2 {
return 1.0;
}
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
}
if ($variance1 == 0 && $variance2 == 0) {
return 1.0;
}
$variance1 /= (@a.elems - 1);
$variance2 /= (@b.elems - 1);
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;
return $value if $value < 0.0 || 1.0 < $value;
# Special cases
return $value if $value == 0.0 || $value == 1.0;
$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);
}
$term = 1.0;
$ai = 1.0;
$value = 1.0;
$ns = $qq + $cx * $psq;
$ns = $ns.Int;
#Soper reduction formula.
$rx = $xx / $cx;
$temp = $qq - $ai;
$rx = $xx if $ns == 0;
while (True) {
$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;
}
$ai++;
$ns--;
if 0 <= $ns {
$temp = $qq - $ai;
$rx = $xx if $ns == 0;
} else {
$temp = $psq;
$psq = $psq + 1.0;
}
}
return $value;
}
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],#0.
[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];
my @CORRECT_ANSWERS = (0.021378001462867,
0.148841696605327,
0.0359722710297968,
0.090773324285671,
0.0107515611497845,
0.00339907162713746,
0.52726574965384,
0.545266866977794);
my UInt $i = 0;
my Real $error = 0.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**-9 {
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);
这个子阵列的不同之处在于它具有"/"用于除法.如何让Perl6 for-loop评估这个子数组？
编辑:我正在努力构建一个最小的工作示例.我发布了整个代码,以便编译.
(事实证明,这个答案完全错过了病弱的@con.但是我不会删除它,因为它收集了一些希望有用的关于理性数字的链接.)
为什么数组在声明中跳过计算值？
事实并非如此.
我正在学习Perl6 ......
1.3数学和数学中的小数数字和Perl 6中的数字.一些解释:
P6 doc的Numerics页面.(由Zoffix撰写.感谢您为P6 Zoffix所做的大量优秀工作.)
我回答了问题"Perl 6的性能是否因使用十进制数的有理数而受到影响".
......来自Perl5
关于它的一些东西属于Perl 5到Perl 6指南文档之一.你愿意开一个新的doc问题吗？
[1.0/15.0, 10.0/62.0],#this isn't evaluated
[1.0/10, 2/50.0],#neither is this
他们都被评估了.
在数学和P6中,字面值1.0是小数是理性的.foo / bar也是理性的.
(当然,它是如果foo和bar均为整数或有理数和所述结果中的任分母仍然是64位或更少,或之一foo或者bar是一个任意精度FatRat 理性.)
但是,Perl6似乎不像我在这里指出的那样.
你还没有解释你所看到的让你觉得P6不喜欢它们的情况.
很有可能的是,根据布拉德吉尔伯特的评论,你看到的是价值,比方说,<1/15>而不是0.066667.前者是P6字面值,100%准确表示1/15.为了100%准确地表示在十进制中,显示器必须在末尾0.06666666666666...具有一个...或一些这样的表示无限重复的最终结果6.但<1/15>表示相同数量和较短,可以说是如此简单dd和.perl使用的<1/15>形式来代替.
当我声明2D数组时,如何让Perl6计算这样的表达式？
你不需要做任何事情.它正在评估它们.:)
(这是一个nanswer,即不是答案本身.它最初是在@con重写他们的问题以包括他们的所有代码之后编写的,并且朝着我的第三个也是最后一个答案迈出了一步.现在希望这对于那些学习Perl 6的人来说是一个有用的资源. .)
这是一大堆代码!:)
到目前为止,我的想法是:你是否绝对,积极地,100%肯定你还没有错过一些输入数据？你跳过数据似乎比P6更有可能,特别是考虑到计算出的值正是你期望的下一个正确结果.
(更新确实证明问题是输入数据不正确.)
这个nanswer的其余部分是问题中代码的逐行"清理"(不是重构).我有两个目标,第二个是你最重要的目标,亲爱的读者,请阅读:
我的翻译有效地向我证明并向@con证明我已经考虑了他们的所有代码.这旨在减少关于bug可能位置的不确定性.请注意,我的大部分更改都可能与他们的错误没有直接关系,但在我完成重写之前,我感觉不舒服.
@con的代码,以及我的翻译,对于任何学习Perl 6的人来说都可能是有用的.@con的代码是Perl 5代码的P6翻译.P5代码又是C代码的翻译.并且该代码还有其他语言的其他翻译.我的代码采用@con的翻译并将其翻译成更具惯用性的版本,其中包含解释为什么我更改了代码.
我对@ con的代码的翻译,减去我的评论(见下文),在tio中.
(我最初的计划是,我们会根据修改tio中的代码并共享修改后的tio版本的链接来进一步探索@con的错误.为此,你/他们只需点击顶部的链接图标( )单击链接图标时获取代码链接.Tio非常适合编写,编辑和运行用P6和其他语言编写的代码然后共享它.)
我lgamma按原样离开了开门程序.我在脚注中写了很多关于它的评论,主要是针对其他读者,因为它包含了大量有趣的功能1:
sub lgamma ( Num(Real) \n --> Num ){
use NativeCall;
sub lgamma (num64 --> num64) is native {}
lgamma( n )
}
我注意到lgamma代码是关于浮点数的(在P6中是Num/ num...types).
sub pvalue (@a, @b) {
对于那些不了解Perls的人,sub关键字引入了一个子程序(又称函数).
这个需要两个"listy"(位置)参数.
(当你看到@(或者作为"sigil",例如在@foo操作员中或作为操作员,例如在中@($bar))时,请考虑"list".)
为了加快代码读取速度并减少重复代码,我将其替换为:
if @a.elems <= 1 {
return 1.0;
}
if @b.elems <= 1 {
return 1.0;
}
有了这个:
return 1 if @a | @b <= 1;
用英语阅读" 1如果列表@a或列表中的元素数@b小于或等于1则返回".
我用|运算符构造了一个any连接点.
(Don't waste time trying to wrap your head around the theory of what junctions do, and how they do it. Just use them in practice in simple and/or succinct ways. If they're kept simple, and read in a simple fashion, they're great and just do the obvious thing. If using them makes code not obvious in a particular case, then consider using some other construct for that particular case.)
An array in numeric context evaluates to its number of elements. Numeric operators like <= impose numeric context. So I dropped the .elems.
Numeric context, and arrays evaluating to their length in numeric context, are basic aspects of P6 so this is idiomatic coding that's appropriate for all but the most basic introductory newbie examples.
I switched from 1.0 to 1.
My guess is that @con wrote 1.0 because that was how the value was literally written in code they were translating and/or with the intent it represented a floating point value. But in P6 plain decimal literals (without an e exponent) like 1.0 produce rational numbers instead.
Using 1.0 instead of 1 in this code replaces a simpler faster type (integers) with a more complicated slower type (rationals). A value of this slower type will then be forced to convert (more slowly than an integer) into a floating point number when it's used in formulae in which any component value is floating point. Which will happen for most or all the formulae in this code because lgamma returns a floating point number.
More generally, P6 works and reads best if you leave the types of variables and values unspecified, leaving them for the compiler to figure out, unless there's a compelling reason to specify them. Leaving types unspecified reduces cognitive load for the reader and increases flexibility for reuse of the code and for optimization of the code by the compiler.
This leads to a complementary pair of maxims:
By default, leave type information implicit. If you don't know if it matters what type a value or variable is given, or know that it doesn't matter, then don't specify it.
If you make the type of a value, variable or parameter explicit, then P6 will use that type forcing compliance with that type, regardless of whether that improves your code or unnecessarily slows code down or stops it altogether.
If you absolutely know that you need to use a different type or add a type constraint to make code correct, then by all means go ahead. Likewise, if you've already verified that your code is correct without the more specific type but you want to make it faster, safer, or clearer, go ahead. But if you don't know, then leave it generic. P6 has been designed to do what you mean and will succeed far more often than you can possibly imagine, as Ben once said.
To state this point even more forcefully, in analogy with premature optimization, in P6, premature typing is an anti-pattern in both one-off code and long term production code.
my Rat $mean1 = @a.sum / @a.elems;
my Rat $mean2 = @b.sum / @b.elems;
if $mean1 == $mean2 {
return 1.0;
}
became:
my (\mean_a, \mean_b) = @a.sum / @a, @b.sum / @b;
return 1 if mean_a == mean_b;
I typically "slash sigils" of "variables" unless I know I'll need a sigil. If a "variable" doesn't actually vary, it likely doesn't need a sigil.
I've renamed $mean1 to mean_a because it clearly corresponds to @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
}
became:
my ($vari_a, $vari_b);
$vari_a += (mean_a - $_)² for @a;
$vari_b += (mean_b - $_)² for @b;
Read the variable $_ as "it".
if ($variance1 == 0 && $variance2 == 0) {
return 1.0;
}
$variance1 /= (@a.elems - 1);
$variance2 /= (@b.elems - 1);
became:
return 1 unless ($vari_a or $vari_b);
$vari_a /= (@a - 1);
$vari_b /= (@b - 1);
In Perls, 0 evaluates in a boolean test context as False. And an array or list like @a evaluates in a numeric context to its element count.
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;
became:
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;
(What's that $q about in the above line?)
Note that I dropped the types on beta and acu. The value assigned to the beta variable will be a Num anyway because lgamma returns a Num. The value assigned to acu will also be a Num because use of the e exponent in a number literal means the value it constructs is a Num.
return $value if $value < 0.0 || 1.0 < $value;
# Special cases
return $value if $value == 0.0 || $value == 1.0;
became:
return $value unless $value ~~ 0^..^1;
Read the 0^ as "above zero" (not including zero) and the ^1 as "up to one" (not including one).
The ~~ is the "smart match" operator. It returns True if the value on its right accepts the value on its left.
So this return statement returns $value if $value is less than or equal to 0 or greater than or equal to 1.
As a minor tidy up I moved a my declaration of a load of variables that I dropped out in the above rewrite to this point, immediately before they became relevant, and added $ns too:
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);
}
became:
$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);
I rewrote the conditional assignment of a bunch of variables as a ternary so it was easier to see what got assigned to what.
$term = 1.0;
$ai = 1.0;
$value = 1.0;
$ns = $qq + $cx * $psq;
$ns = $ns.Int;
became:
$term = 1;
$ai = 1;
$value = 1;
$ns = ($qq + $cx * $psq) .Int;
Replacing 1.0s with 1s and combining the expression assigned to $ns with the .Int coercion.
(I stripped types from the code as I translated and it continued to calculate the correct results except that removing the above Int coercion made the code infiniloop. This is what finally led me to search the net to see if I could find the code @con was translating. That's when I found it on rosettacode.org. It was explicitly typed as an integer in the code I saw so presumably it's central to ensuring the algorithm works.)
#Soper reduction formula.
$rx = $xx / $cx;
$temp = $qq - $ai;
$rx = $xx if $ns == 0;
(Didn't change.)
while (True) {
became:
loop {
$term = $term * $temp * $rx / ( $pp + $ai );
$value = $value + $term;
$temp = $term.abs;
(Didn't change.)
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;
}
became:
if $temp <= acu & acu * $value {
$value = $value * ($pp * $xx.log + ($qq - 1) * $cx.log - beta).exp / $pp;
$value = 1 - $value if $indx;
last;
}
This time the condition containing the junction (&) reads in English as "if temp is less than or equal to both acu and acu times value".
$ai++;
$ns--;
if 0 <= $ns {
$temp = $qq - $ai;
$rx = $xx if $ns == 0;
} else {
$temp = $psq;
$psq = $psq + 1;
}
}
return $value;
}
I just replaced 1.0 with 1.
Now the problematic array. As I wrote at the start, I'm pretty sure you have (or the supplier of your data has) just forgotten a couple lines:
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];
Who or what says these are the correct answers? Are you 100% sure the 0.0033... answer goes with the [<1.0/15.0>, <10.0/62.0>],[<1.0/10>, <2/50.0>] data?
my @CORRECT_ANSWERS =
0.021378001462867,
0.148841696605327,
0.0359722710297968,
0.090773324285671,
0.0107515611497845,
0.00339907162713746,
0.52726574965384,
0.545266866977794;
And in the last bit I just removed types and once again used a superscript for a prettier value-plus-exponent (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 Love that nice lgamma routine! (Turns out Brad Gilbert wrote it.)
sub lgamma ( Num(Real) \n --> Num ){
use NativeCall;
sub lgamma (num64 --> num64) is native {}
lgamma( n )
}
Showcasing:
Shadowing a low level (C) routine with a high level P6 (Perl 6) routine with the same name to add further processing before calling the C function.
Explicitly converting the input type from the broader P6 type Real numbers to the narrower Num. Perls are "strongly typed" per the original technical definition of the term, but P6 provides additional options per several other interpretations of "strong typing" and relative to languages less capable of these other interpretations of "strong typing" (like Perl 5, Python, and C). Explicit type conversion is part of this shift in capability that P6 introduces.
Slashing the sigil. This is further discussed where I've done the same elsewhere in this post.
Lexically scoping use of a library. The inner lgamma routine, and the symbols imported by use Nativecall;, are not visible anywhere but inside the outer lgamma routine containing it.
Using NativeCall (the P6 C FFI) to allow high level P6 code to sweetly map directly to C code including automatic conversion from P6's IEEE double float compatible boxed type Num to the unboxed machine data type equivalent num64.
All in 5 lines! Very nice. :)
.oO（他说，希望如此。有没有人写过一个 SO 问题的四个答案？！？）
如果您交换数据中的这些行：
[<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],
相反：
[0.010268,0.000167,0.000167],
[0.159258,0.136278,0.122389],
[<1.0/15.0>, <10.0/62.0>],
[<1.0/10>, <2/50.0>],
那么你的程序的输出（至少是我重写你的代码 nanswer 的版本）是：
0.159258
0 [27.5,21,19,23.6,17,17.9,16.9,20.1,21.9,22.6,23.1,19.6,19,21.7,21.4] [27.1 22 20.8 23.4 23.4 23.5 25.8 22 24.8 20.2 21.9 22.1 22.9 20.5 24.4] = 0.02137800146286709
1 [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 23 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] = 0.14884169660532756
2 [19.8,20.4,19.6,17.8,18.5,18.9,18.3,18.9,19.5,22] [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 23.9 21.6 24.3 20.4 24 13.2] = 0.035972271029797116
3 [30.02,29.99,30.11,29.97,30.01,29.99] [29.89 29.93 29.72 29.98 30.02 29.98] = 0.09077332428566681
4 [3,4,1,2.1] [490.2 340 433.9] = 0.010751561149784494
5 [0.010268,0.000167,0.000167] [0.159258 0.136278 0.122389] = 0.003399071627137453
6 [1.0/15.0,10.0/62.0] [1.0/10 2/50.0] = 0.5272657496538401
7 [0.391304,0.466667,0] [0 0.446809 0] = 0.5452668669777938
the cumulative error is 5.50254e-15
回想起来，这显然是错误的。2020 年事后诸葛亮等等。:)
| 归档时间： |
|
| 查看次数： |
172 次 |
| 最近记录： |