Permutation test: Difference between revisions
(Updated second D entry) |
m (→{{header|Wren}}: Minor tidy) |
||
(41 intermediate revisions by 22 users not shown) | |||
Line 9: | Line 9: | ||
was measured and reported in this table. |
was measured and reported in this table. |
||
{| style="text-align: left; width: |
:::::::{| style="text-align: left; width: 25%;" border="4" cellpadding="2" cellspacing="2" |
||
|+ Table of experimental results |
|+ Table of experimental results |
||
|- style="background-color: rgb(255, 204, 255);" |
|- style="background-color: rgb(255, 204, 255);" |
||
Line 44: | Line 44: | ||
* Report the percentage of alternative groupings for which the difference in means is less or equal to the actual experimentally observed difference in means, and the percentage for which it is greater. |
* Report the percentage of alternative groupings for which the difference in means is less or equal to the actual experimentally observed difference in means, and the percentage for which it is greater. |
||
* Note that they should sum to 100%. |
* Note that they should sum to 100%. |
||
Extremely dissimilar values are evidence of an effect not entirely due |
Extremely dissimilar values are evidence of an effect not entirely due |
||
Line 51: | Line 52: | ||
that's easier than loading them at run time. Test your solution on the |
that's easier than loading them at run time. Test your solution on the |
||
data given above. |
data given above. |
||
<br><br> |
|||
=={{header|11l}}== |
|||
{{trans|Kotlin}} |
|||
<syntaxhighlight lang="11l">V data = [85, 88, 75, 66, 25, 29, 83, 39, 97, |
|||
68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
|||
F pick(at, remain, accu, treat) |
|||
I remain == 0 |
|||
R I accu > treat {1} E 0 |
|||
R pick(at - 1, remain - 1, accu + :data[at - 1], treat) + |
|||
(I at > remain {pick(at - 1, remain, accu, treat)} E 0) |
|||
V treat = 0 |
|||
V total = 1.0 |
|||
L(i) 0..8 |
|||
treat += data[i] |
|||
L(i) (19..11).step(-1) |
|||
total *= i |
|||
L(i) (9..1).step(-1) |
|||
total /= i |
|||
V gt = pick(19, 9, 0, treat) |
|||
V le = Int(total - gt) |
|||
print(‘<= : #.6% #.’.format(100 * le / total, le)) |
|||
print(‘ > : #.6% #.’.format(100 * gt / total, gt))</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
<= : 87.197168% 80551 |
|||
> : 12.802832% 11827 |
|||
</pre> |
|||
=={{header|Ada}}== |
=={{header|Ada}}== |
||
< |
<syntaxhighlight lang="ada">with Ada.Text_IO; with Iterate_Subsets; |
||
procedure Permutation_Test is |
procedure Permutation_Test is |
||
Line 115: | Line 150: | ||
New_Line; |
New_Line; |
||
end; |
end; |
||
end Permutation_Test;</ |
end Permutation_Test;</syntaxhighlight> |
||
This solution uses an auxiliary package Iterate_Subsets. Here is the Spec: |
This solution uses an auxiliary package Iterate_Subsets. Here is the Spec: |
||
< |
<syntaxhighlight lang="ada">generic |
||
Subset_Size, More_Elements: Positive; |
Subset_Size, More_Elements: Positive; |
||
package Iterate_Subsets is |
package Iterate_Subsets is |
||
Line 134: | Line 169: | ||
end Iterate_Subsets; |
end Iterate_Subsets; |
||
</syntaxhighlight> |
|||
</lang> |
|||
And here is the implementation: |
And here is the implementation: |
||
< |
<syntaxhighlight lang="ada">package body Iterate_Subsets is |
||
function First return Subset is |
function First return Subset is |
||
Line 171: | Line 206: | ||
end Last; |
end Last; |
||
end Iterate_Subsets;</ |
end Iterate_Subsets;</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 178: | Line 213: | ||
</pre> |
</pre> |
||
=={{header|AWK}}== |
|||
<syntaxhighlight lang="awk"> |
|||
# syntax: GAWK -f PERMUTATION_TEST.AWK |
|||
# converted from C |
|||
BEGIN { |
|||
# "treatment..................control......................" |
|||
n = split("85,88,75,66,25,29,83,39,97,68,41,10,49,16,65,32,92,28,98",data,",") |
|||
for (i=1; i<=n; i++) { # make AWK array look like a C array |
|||
data[i-1] = data[i] |
|||
} |
|||
delete data[n] |
|||
total = 1 |
|||
for (i=0; i<9; i++) { treat += data[i] } |
|||
for (i=19; i>10; i--) { total *= i } |
|||
for (i=9; i>0; i--) { total /= i } |
|||
gt = pick(19,9,0,treat) |
|||
le = total - gt |
|||
printf("<= : %9.6f%% %6d\n",100*le/total,le) |
|||
printf(" > : %9.6f%% %6d\n",100*gt/total,gt) |
|||
exit(0) |
|||
} |
|||
function pick(at,remain,accu,treat) { |
|||
if (!remain) { |
|||
return (accu > treat) ? 1 : 0 |
|||
} |
|||
return pick(at-1,remain-1,accu+data[at-1],treat) + ( (at > remain) ? pick(at-1,remain,accu,treat) : 0 ) |
|||
} |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
<= : 87.197168% 80551 |
|||
> : 12.802832% 11827 |
|||
</pre> |
|||
=={{header|BBC BASIC}}== |
=={{header|BBC BASIC}}== |
||
{{works with|BBC BASIC for Windows}} |
{{works with|BBC BASIC for Windows}} |
||
< |
<syntaxhighlight lang="bbcbasic"> ntreated% = 9 |
||
nplacebo% = 10 |
nplacebo% = 10 |
||
DIM results%(ntreated% + nplacebo% - 1) |
DIM results%(ntreated% + nplacebo% - 1) |
||
Line 218: | Line 286: | ||
N% = (N% + (N% >>> 4)) AND &0F0F0F0F |
N% = (N% + (N% >>> 4)) AND &0F0F0F0F |
||
N% += N% >>> 8 : N% += N% >>> 16 |
N% += N% >>> 8 : N% += N% >>> 16 |
||
= N% AND &7F</ |
= N% AND &7F</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 226: | Line 294: | ||
=={{header|C}}== |
=={{header|C}}== |
||
< |
<syntaxhighlight lang="c">#include <stdio.h> |
||
int data[] = { 85, 88, 75, 66, 25, 29, 83, 39, 97, |
int data[] = { 85, 88, 75, 66, 25, 29, 83, 39, 97, |
||
Line 254: | Line 322: | ||
100 * le / total, le, 100 * gt / total, gt); |
100 * le / total, le, 100 * gt / total, gt); |
||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
Output:<lang><= : 87.197168% 80551 |
Output:<syntaxhighlight lang="text"><= : 87.197168% 80551 |
||
> : 12.802832% 11827</ |
> : 12.802832% 11827</syntaxhighlight> |
||
=={{header|C sharp|C#}}== |
|||
{{trans|Java}} |
|||
<syntaxhighlight lang="cs">using System; |
|||
using System.Collections.Generic; |
|||
namespace PermutationTest { |
|||
class Program { |
|||
static readonly List<int> DATA = new List<int>{ |
|||
85, 88, 75, 66, 25, 29, 83, 39, 97, |
|||
68, 41, 10, 49, 16, 65, 32, 92, 28, 98 |
|||
}; |
|||
static int Pick(int at, int remain, int accu, int treat) { |
|||
if (remain == 0) { |
|||
return (accu > treat) ? 1 : 0; |
|||
} |
|||
return Pick(at - 1, remain - 1, accu + DATA[at - 1], treat) |
|||
+ ((at > remain) ? Pick(at - 1, remain, accu, treat) : 0); |
|||
} |
|||
static void Main() { |
|||
int treat = 0; |
|||
double total = 1.0; |
|||
for (int i = 0; i <= 8; i++) { |
|||
treat += DATA[i]; |
|||
} |
|||
for (int i = 19; i >= 11; i--) { |
|||
total *= i; |
|||
} |
|||
for (int i = 9; i >= 1; --i) { |
|||
total /= i; |
|||
} |
|||
int gt = Pick(19, 9, 0, treat); |
|||
int le = (int) (total - gt); |
|||
Console.WriteLine("<= {0}% {1}", 100.0 * le / total, le); |
|||
Console.WriteLine(" > {0}% {1}", 100.0 * gt / total, gt); |
|||
} |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre><= 87.1971681569205% 80551 |
|||
> 12.8028318430795% 11827</pre> |
|||
=={{header|C++}}== |
|||
This is a translaion of C |
|||
<syntaxhighlight lang="cpp">#include<iostream> |
|||
#include<vector> |
|||
#include<numeric> |
|||
#include<functional> |
|||
class |
|||
{ |
|||
public: |
|||
int64_t operator()(int n, int k){ return partial_factorial(n, k) / factorial(n - k);} |
|||
private: |
|||
int64_t partial_factorial(int from, int to) { return from == to ? 1 : from * partial_factorial(from - 1, to); } |
|||
int64_t factorial(int n) { return n == 0 ? 1 : n * factorial(n - 1);} |
|||
}combinations; |
|||
int main() |
|||
{ |
|||
static constexpr int treatment = 9; |
|||
const std::vector<int> data{ 85, 88, 75, 66, 25, 29, 83, 39, 97, |
|||
68, 41, 10, 49, 16, 65, 32, 92, 28, 98 }; |
|||
int treated = std::accumulate(data.begin(), data.begin() + treatment, 0); |
|||
std::function<int (int, int, int)> pick; |
|||
pick = [&](int n, int from, int accumulated) |
|||
{ |
|||
if(n == 0) |
|||
return accumulated > treated ? 1 : 0; |
|||
else |
|||
return pick(n - 1, from - 1, accumulated + data[from - 1]) + |
|||
(from > n ? pick(n, from - 1, accumulated) : 0); |
|||
}; |
|||
int total = combinations(data.size(), treatment); |
|||
int greater = pick(treatment, data.size(), 0); |
|||
int lesser = total - greater; |
|||
std::cout << "<= : " << 100.0 * lesser / total << "% " << lesser << std::endl |
|||
<< " > : " << 100.0 * greater / total << "% " << greater << std::endl; |
|||
}</syntaxhighlight> |
|||
Output:<syntaxhighlight lang="text"><= : 87.197168% 80551 |
|||
> : 12.802832% 11827</syntaxhighlight> |
|||
=={{header|Common Lisp}}== |
=={{header|Common Lisp}}== |
||
< |
<syntaxhighlight lang="lisp">(defun perm-test (s1 s2) |
||
(let ((more 0) (leq 0) |
(let ((more 0) (leq 0) |
||
(all-data (append s1 s2)) |
(all-data (append s1 s2)) |
||
Line 282: | Line 437: | ||
(format t "<=: ~a ~6f%~% >: ~a ~6f%~%" |
(format t "<=: ~a ~6f%~% >: ~a ~6f%~%" |
||
x (* 100e0 (/ x s)) |
x (* 100e0 (/ x s)) |
||
y (* 100e0 (/ y s))))</ |
y (* 100e0 (/ y s))))</syntaxhighlight>output<syntaxhighlight lang="text"><=: 80551 87.197% |
||
>: 11827 12.803%</ |
>: 11827 12.803%</syntaxhighlight> |
||
=={{header|D}}== |
=={{header|D}}== |
||
< |
<syntaxhighlight lang="d">import std.stdio, std.algorithm, std.array, combinations3; |
||
auto permutationTest(T)(in T[] a, in T[] b) pure nothrow @safe { |
auto permutationTest(T)(in T[] a, in T[] b) pure nothrow @safe { |
||
Line 300: | Line 455: | ||
immutable under = permutationTest(treatmentGroup, controlGroup); |
immutable under = permutationTest(treatmentGroup, controlGroup); |
||
writefln("Under =%6.2f%%\nOver =%6.2f%%", under, 100.0 - under); |
writefln("Under =%6.2f%%\nOver =%6.2f%%", under, 100.0 - under); |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Under = 87.20% |
<pre>Under = 87.20% |
||
Line 307: | Line 462: | ||
Alternative version: |
Alternative version: |
||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight lang="d">void main() @safe { |
||
import std.stdio, std.algorithm, std.range; |
import std.stdio, std.algorithm, std.range; |
||
Line 330: | Line 485: | ||
writefln(" > : %2.2f%% %d", 100.0 * gt / t, gt); |
writefln(" > : %2.2f%% %d", 100.0 * gt / t, gt); |
||
writefln("<= : %2.2f%% %d", 100.0 * le / t, le); |
writefln("<= : %2.2f%% %d", 100.0 * le / t, le); |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> > : 12.80% 11827 |
<pre> > : 12.80% 11827 |
||
<= : 87.20% 80551</pre> |
<= : 87.20% 80551</pre> |
||
=={{header|Delphi}}== |
|||
{{libheader| System.SysUtils}} |
|||
{{Trans|Go}} |
|||
<syntaxhighlight lang="delphi"> |
|||
program Permutation_test; |
|||
{$APPTYPE CONSOLE} |
|||
uses |
|||
System.SysUtils; |
|||
procedure Comb(n, m: Integer; emit: TProc<TArray<Integer>>); |
|||
var |
|||
s: TArray<Integer>; |
|||
last: Integer; |
|||
procedure rc(i, next: Integer); |
|||
begin |
|||
for var j := next to n - 1 do |
|||
begin |
|||
s[i] := j; |
|||
if i = last then |
|||
emit(s) |
|||
else |
|||
rc(i + 1, j + 1); |
|||
end; |
|||
end; |
|||
begin |
|||
SetLength(s, m); |
|||
last := m - 1; |
|||
rc(0, 0); |
|||
end; |
|||
begin |
|||
var tr: TArray<Integer> := [85, 88, 75, 66, 25, 29, 83, 39, 97]; |
|||
var ct: TArray<Integer> := [68, 41, 10, 49, 16, 65, 32, 92, 28, 98]; |
|||
// collect all results in a single list |
|||
var all: TArray<Integer> := concat(tr, ct); |
|||
// compute sum of all data, useful as intermediate result |
|||
var sumAll := 0; |
|||
for var r in all do |
|||
inc(sumAll, r); |
|||
// closure for computing scaled difference. |
|||
// compute results scaled by len(tr)*len(ct). |
|||
// this allows all math to be done in integers. |
|||
var sd := |
|||
function(trc: TArray<Integer>): Integer |
|||
begin |
|||
var sumTr := 0; |
|||
for var x in trc do |
|||
inc(sumTr, all[x]); |
|||
result := sumTr * length(ct) - (sumAll - sumTr) * length(tr); |
|||
end; |
|||
// compute observed difference, as an intermediate result |
|||
var a: TArray<Integer>; |
|||
SetLength(a, length(tr)); |
|||
for var i := 0 to High(a) do |
|||
a[i] := i; |
|||
var sdObs := sd(a); |
|||
// iterate over all combinations. for each, compute (scaled) |
|||
// difference and tally whether leq or gt observed difference. |
|||
var nLe, nGt: Integer; |
|||
comb(length(all), length(tr), |
|||
procedure(c: TArray<Integer>) |
|||
begin |
|||
if sd(c) > sdObs then |
|||
inc(nGt) |
|||
else |
|||
inc(nle); |
|||
end); |
|||
// print results as percentage |
|||
var pc := 100 / (nLe + nGt); |
|||
writeln(format('differences <= observed: %f%%', [nle * pc])); |
|||
writeln(format('differences > observed: %f%%', [ngt * pc])); |
|||
{$IFNDEF UNIX} readln; {$ENDIF} |
|||
end.</syntaxhighlight> |
|||
=={{header|Elixir}}== |
|||
{{trans|Ruby}} |
|||
<syntaxhighlight lang="elixir">defmodule Permutation do |
|||
def statistic(ab, a) do |
|||
sumab = Enum.sum(ab) |
|||
suma = Enum.sum(a) |
|||
suma / length(a) - (sumab - suma) / (length(ab) - length(a)) |
|||
end |
|||
def test(a, b) do |
|||
ab = a ++ b |
|||
tobs = statistic(ab, a) |
|||
{under, count} = Enum.reduce(comb(ab, length(a)), {0,0}, fn perm, {under, count} -> |
|||
if statistic(ab, perm) <= tobs, do: {under+1, count+1}, |
|||
else: {under , count+1} |
|||
end) |
|||
under * 100.0 / count |
|||
end |
|||
defp comb(_, 0), do: [[]] |
|||
defp comb([], _), do: [] |
|||
defp comb([h|t], m) do |
|||
(for l <- comb(t, m-1), do: [h|l]) ++ comb(t, m) |
|||
end |
|||
end |
|||
treatmentGroup = [85, 88, 75, 66, 25, 29, 83, 39, 97] |
|||
controlGroup = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
|||
under = Permutation.test(treatmentGroup, controlGroup) |
|||
:io.fwrite "under = ~.2f%, over = ~.2f%~n", [under, 100-under]</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
under = 87.20%, over = 12.80% |
|||
</pre> |
|||
=={{header|GAP}}== |
=={{header|GAP}}== |
||
< |
<syntaxhighlight lang="gap">a := [85, 88, 75, 66, 25, 29, 83, 39, 97]; |
||
b := [68, 41, 10, 49, 16, 65, 32, 92, 28, 98]; |
b := [68, 41, 10, 49, 16, 65, 32, 92, 28, 98]; |
||
Line 386: | Line 663: | ||
# in order, % less or greater than original diff |
# in order, % less or greater than original diff |
||
PermTest(a, b); |
PermTest(a, b); |
||
[ "87.197", "12.802" ]</ |
[ "87.197", "12.802" ]</syntaxhighlight> |
||
=={{header|FreeBASIC}}== |
|||
{{trans|Phix}} |
|||
<syntaxhighlight lang="freebasic"> |
|||
Dim Shared datos(18) As Integer => {85, 88, 75, 66, 25, 29, 83, 39, 97,_ |
|||
68, 41, 10, 49, 16, 65, 32, 92, 28, 98} |
|||
Function pick(at As Integer, remain As Integer, accu As Integer, treat As Integer) As Integer |
|||
If remain = 0 Then |
|||
If accu > treat Then Return 1 Else Return 0 |
|||
End If |
|||
Dim a As Integer |
|||
If at > remain Then a = pick(at-1, remain, accu, treat) Else a = 0 |
|||
Return pick(at-1, remain-1, accu+datos(at), treat) + a |
|||
End Function |
|||
Dim As Integer treat = 0, le, gt, total = 1, i |
|||
For i = 1 To 9 |
|||
treat += datos(i) |
|||
Next i |
|||
For i = 19 To 11 Step -1 |
|||
total *= i |
|||
Next i |
|||
For i = 9 To 1 Step -1 |
|||
total /= i |
|||
Next i |
|||
gt = pick(19, 9, 0, treat) |
|||
le = total - gt |
|||
Print Using "<= : ##.######% #####"; 100*le/total; le |
|||
Print Using " > : ##.######% #####"; 100*gt/total; gt |
|||
Sleep |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
<= : 92.076035% 85058 |
|||
> : 7.923965% 7320 |
|||
</pre> |
|||
=={{header|Go}}== |
=={{header|Go}}== |
||
A version doing all math in integers until computing final percentages. |
A version doing all math in integers until computing final percentages. |
||
< |
<syntaxhighlight lang="go">package main |
||
import "fmt" |
import "fmt" |
||
Line 461: | Line 781: | ||
} |
} |
||
rc(0, 0) |
rc(0, 0) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 469: | Line 789: | ||
=={{header|Haskell}}== |
=={{header|Haskell}}== |
||
< |
<syntaxhighlight lang="haskell">binomial n m = (f !! n) `div` (f !! m) `div` (f !! (n - m)) |
||
where f = scanl (*) 1 [1..] |
where f = scanl (*) 1 [1..] |
||
Line 489: | Line 809: | ||
[68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
[68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
||
in do putStr "> : "; print r |
in do putStr "> : "; print r |
||
putStr "<=: "; print $ 100 - r</ |
putStr "<=: "; print $ 100 - r</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 497: | Line 817: | ||
Somewhat faster, this goes from top down: |
Somewhat faster, this goes from top down: |
||
< |
<syntaxhighlight lang="haskell">binomial n m = (f !! n) `div` (f !! m) `div` (f !! (n - m)) |
||
where f = scanl (*) 1 [1..] |
where f = scanl (*) 1 [1..] |
||
Line 515: | Line 835: | ||
(a, b) = perms [85, 88, 75, 66, 25, 29, 83, 39, 97] |
(a, b) = perms [85, 88, 75, 66, 25, 29, 83, 39, 97] |
||
[68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
[68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
||
lt = (fromIntegral a) / (fromIntegral b) * 100</ |
lt = (fromIntegral a) / (fromIntegral b) * 100</syntaxhighlight> |
||
In cases where the sample data are a large number of relatively small positive integers, counting number of partial sums is a lot faster: |
In cases where the sample data are a large number of relatively small positive integers, counting number of partial sums is a lot faster: |
||
< |
<syntaxhighlight lang="haskell">combs maxsum len x = foldl f [(0,0,1)] x where |
||
f a n = merge a (map (addNum n) $ filter (\(l,_,_) -> l < len) a) |
f a n = merge a (map (addNum n) $ filter (\(l,_,_) -> l < len) a) |
||
addNum n (a,s,c) |
addNum n (a,s,c) |
||
Line 543: | Line 863: | ||
-- how many combinations are less than current sum |
-- how many combinations are less than current sum |
||
main = print$ permtest [85, 88, 75, 66, 25, 29, 83, 39, 97] |
main = print$ permtest [85, 88, 75, 66, 25, 29, 83, 39, 97] |
||
[68, 41, 10, 49, 16, 65, 32, 92, 28, 98]</ |
[68, 41, 10, 49, 16, 65, 32, 92, 28, 98]</syntaxhighlight> |
||
=={{header|J}}== |
=={{header|J}}== |
||
< |
<syntaxhighlight lang="j">require'stats' |
||
trmt=: 0.85 0.88 0.75 0.66 0.25 0.29 0.83 0.39 0.97 |
trmt=: 0.85 0.88 0.75 0.66 0.25 0.29 0.83 0.39 0.97 |
||
ctrl=: 0.68 0.41 0.1 0.49 0.16 0.65 0.32 0.92 0.28 0.98 |
ctrl=: 0.68 0.41 0.1 0.49 0.16 0.65 0.32 0.92 0.28 0.98 |
||
Line 553: | Line 873: | ||
all=: trmt(#@[ ({. difm }.) |:@([ (comb ~.@,"1 i.@])&# ,) { ,) ctrl |
all=: trmt(#@[ ({. difm }.) |:@([ (comb ~.@,"1 i.@])&# ,) { ,) ctrl |
||
smoutput 'under: ','%',~":100*mean all <: result |
smoutput 'under: ','%',~":100*mean all <: result |
||
smoutput 'over: ','%',~":100*mean all > result</ |
smoutput 'over: ','%',~":100*mean all > result</syntaxhighlight> |
||
Result: |
Result: |
||
<lang>under: 87.1972% |
<syntaxhighlight lang="text">under: 87.1972% |
||
over: 12.8028%</ |
over: 12.8028%</syntaxhighlight> |
||
=={{header|Java}}== |
|||
{{trans|Kotlin}} |
|||
<syntaxhighlight lang="java">public class PermutationTest { |
|||
private static final int[] data = new int[]{ |
|||
85, 88, 75, 66, 25, 29, 83, 39, 97, |
|||
68, 41, 10, 49, 16, 65, 32, 92, 28, 98 |
|||
}; |
|||
private static int pick(int at, int remain, int accu, int treat) { |
|||
if (remain == 0) return (accu > treat) ? 1 : 0; |
|||
return pick(at - 1, remain - 1, accu + data[at - 1], treat) |
|||
+ ((at > remain) ? pick(at - 1, remain, accu, treat) : 0); |
|||
} |
|||
public static void main(String[] args) { |
|||
int treat = 0; |
|||
double total = 1.0; |
|||
for (int i = 0; i <= 8; ++i) { |
|||
treat += data[i]; |
|||
} |
|||
for (int i = 19; i >= 11; --i) { |
|||
total *= i; |
|||
} |
|||
for (int i = 9; i >= 1; --i) { |
|||
total /= i; |
|||
} |
|||
int gt = pick(19, 9, 0, treat); |
|||
int le = (int) (total - gt); |
|||
System.out.printf("<= : %f%% %d\n", 100.0 * le / total, le); |
|||
System.out.printf(" > : %f%% %d\n", 100.0 * gt / total, gt); |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre><= : 87.197168% 80551 |
|||
> : 12.802832% 11827</pre> |
|||
=={{header|jq}}== |
=={{header|jq}}== |
||
{{works with|jq|1.4}} |
{{works with|jq|1.4}} |
||
'''Part 1: Combinations''' |
'''Part 1: Combinations''' |
||
< |
<syntaxhighlight lang="jq"># combination(r) generates a stream of combinations of r items from the input array. |
||
def combination(r): |
def combination(r): |
||
if r > length or r < 0 then empty |
if r > length or r < 0 then empty |
||
Line 568: | Line 925: | ||
( .[1:]|combination(r)) |
( .[1:]|combination(r)) |
||
end; |
end; |
||
</syntaxhighlight> |
|||
</lang> |
|||
'''Part 2: Permutation Test''' |
'''Part 2: Permutation Test''' |
||
< |
<syntaxhighlight lang="jq"># a and b should be arrays: |
||
def permutationTest(a; b): |
def permutationTest(a; b): |
||
Line 595: | Line 952: | ||
| .[1] += 1 ) |
| .[1] += 1 ) |
||
| .[0] * 100.0 / .[1] # under/count |
| .[0] * 100.0 / .[1] # under/count |
||
;</ |
;</syntaxhighlight> |
||
'''Example:''' |
'''Example:''' |
||
< |
<syntaxhighlight lang="jq">def treatmentGroup: [85, 88, 75, 66, 25, 29, 83, 39, 97]; |
||
def controlGroup: [68, 41, 10, 49, 16, 65, 32, 92, 28, 98]; |
def controlGroup: [68, 41, 10, 49, 16, 65, 32, 92, 28, 98]; |
||
permutationTest(treatmentGroup; controlGroup) as $under |
permutationTest(treatmentGroup; controlGroup) as $under |
||
| "% under=\($under)", "% over=\(100 - $under)"</ |
| "% under=\($under)", "% over=\(100 - $under)"</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
$ jq -n -r -f permutation_test.jq |
$ jq -n -r -f permutation_test.jq |
||
Line 607: | Line 964: | ||
% over=12.856957284201869 |
% over=12.856957284201869 |
||
=={{header| |
=={{header|Julia}}== |
||
{{works with|Julia|0.6}} |
|||
<lang mathematica>"<=: " <> ToString[#1] <> " " <> ToString[100. #1/#2] <> "%\n >: " <> |
|||
The primary function for this solution is <tt>permutation_test</tt>, which relies on Julia's <tt>combinations</tt> (from <tt>Combinatorics</tt> module) function to provide all of the possible study arm assignments. <tt>bifurcate</tt> splits the pooled results into "treatment" and "control" groups according to the indices provided by <tt>combinations</tt>. |
|||
'''Functions''' |
|||
<syntaxhighlight lang="julia">using Combinatorics |
|||
meandiff(a::Vector{T}, b::Vector{T}) where T <: Real = mean(a) - mean(b) |
|||
function bifurcate(a::AbstractVector, sel::Vector{T}) where T <: Integer |
|||
x = a[sel] |
|||
asel = trues(length(a)) |
|||
asel[sel] = false |
|||
y = a[asel] |
|||
return x, y |
|||
end |
|||
function permutation_test(treated::Vector{T}, control::Vector{T}) where T <: Real |
|||
effect0 = meandiff(treated, control) |
|||
pool = vcat(treated, control) |
|||
tlen = length(treated) |
|||
plen = length(pool) |
|||
better = worse = 0 |
|||
for subset in combinations(1:plen, tlen) |
|||
t, c = bifurcate(pool, subset) |
|||
if effect0 < meandiff(t, c) |
|||
better += 1 |
|||
else |
|||
worse += 1 |
|||
end |
|||
end |
|||
return better, worse |
|||
end</syntaxhighlight> |
|||
'''Main''' |
|||
<syntaxhighlight lang="julia">const treated = [85, 88, 75, 66, 25, 29, 83, 39, 97] |
|||
const control = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
|||
(better, worse) = permutation_test(treated, control) |
|||
tot = better + worse |
|||
println("Permutation test using the following data:") |
|||
println("Treated: ", treated) |
|||
println("Control: ", control) |
|||
println("\nThere are $tot different permuted groups of these data.") |
|||
@printf("%8d, %5.2f%% showed better than actual results.\n", better, 100 * better / tot) |
|||
print(@sprintf("%8d, %5.2f%% showed equalivalent or worse results.", worse, 100 * worse / tot))</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Permutation test using the following data: |
|||
Treated: [85, 88, 75, 66, 25, 29, 83, 39, 97] |
|||
Control: [68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
|||
There are 92378 different permuted groups of these data. |
|||
11827, 12.80% showed better than actual results. |
|||
80551, 87.20% showed equalivalent or worse results.</pre> |
|||
=={{header|Kotlin}}== |
|||
{{trans|C}} |
|||
<syntaxhighlight lang="scala">// version 1.1.2 |
|||
val data = intArrayOf( |
|||
85, 88, 75, 66, 25, 29, 83, 39, 97, |
|||
68, 41, 10, 49, 16, 65, 32, 92, 28, 98 |
|||
) |
|||
fun pick(at: Int, remain: Int, accu: Int, treat: Int): Int { |
|||
if (remain == 0) return if (accu > treat) 1 else 0 |
|||
return pick(at - 1, remain - 1, accu + data[at - 1], treat) + |
|||
if (at > remain) pick(at - 1, remain, accu, treat) else 0 |
|||
} |
|||
fun main(args: Array<String>) { |
|||
var treat = 0 |
|||
var total = 1.0 |
|||
for (i in 0..8) treat += data[i] |
|||
for (i in 19 downTo 11) total *= i |
|||
for (i in 9 downTo 1) total /= i |
|||
val gt = pick(19, 9, 0, treat) |
|||
val le = (total - gt).toInt() |
|||
System.out.printf("<= : %f%% %d\n", 100.0 * le / total, le) |
|||
System.out.printf(" > : %f%% %d\n", 100.0 * gt / total, gt) |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
<= : 87.197168% 80551 |
|||
> : 12.802832% 11827 |
|||
</pre> |
|||
=={{header|M2000 Interpreter}}== |
|||
{{trans|C}} |
|||
<syntaxhighlight lang="m2000 interpreter"> |
|||
Module Checkit { |
|||
Global data(), treat=0 |
|||
data()=(85, 88, 75, 66, 25, 29, 83, 39, 97,68, 41, 10, 49, 16, 65, 32, 92, 28, 98) |
|||
Function pick(at, remain, accu) { |
|||
If remain Else =If(accu>treat->1,0):Exit |
|||
=pick(at-1,remain-1,accu+data(at-1))+If(at>remain->pick(at-1, remain, accu),0) |
|||
} |
|||
total=1 |
|||
For i=0 to 8 {treat+=data(i)} |
|||
For i=19 to 11 {total*=i} |
|||
For i=9 to 1 {total/=i} |
|||
gt=pick(19,9,0) |
|||
le=total-gt |
|||
Print Format$("<= : {0:1}% {1}", 100*le/total, le) |
|||
Print Format$(" > : {0:1}% {1}", 100*gt/total, gt) |
|||
} |
|||
Checkit |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre><= : 87.2% 80551 |
|||
> : 12.8% 11827 |
|||
</pre> |
|||
Slower version, using a lambda function with a series of inner lambda functions to return each combination at a time. |
|||
<syntaxhighlight lang="m2000 interpreter"> |
|||
Module CheckThis { |
|||
Function CombinationsStep (a, nn) { |
|||
c1=lambda (&f, &a) ->{=car(a) : a=cdr(a) : f=len(a)=0} |
|||
m=len(a) |
|||
c=c1 |
|||
n=m-nn+1 |
|||
p=2 |
|||
while m>n { |
|||
c1=lambda c2=c,n=p, z=(,) (&f, &m) ->{if len(z)=0 then z=cdr(m) |
|||
=cons(car(m),c2(&f, &z)):if f then z=(,) : m=cdr(m) : f=len(m)+len(z)<n |
|||
} |
|||
c=c1 |
|||
p++ |
|||
m-- |
|||
} |
|||
=lambda c, a (&f) ->c(&f, &a) |
|||
} |
|||
treated=(85, 88, 75, 66, 25, 29, 83, 39, 97) |
|||
placebo=(68, 41, 10, 49, 16, 65, 32, 92, 28, 98) |
|||
treat=0 |
|||
m=each(treated): while m {treat+=array(m)} |
|||
total=1 |
|||
for i=len(placebo)+1 to len(placebo) +len(treated):total*=i:next i |
|||
for i=len(placebo)-1 to 1: total/=i:next i |
|||
d=total div 10**int(log(total)) |
|||
k=false |
|||
StepA=CombinationsStep(cons(treated, placebo),len(treated)) |
|||
counter=0 |
|||
gt=0 |
|||
While not k { |
|||
comb=StepA(&k) |
|||
accu=0 |
|||
m=each(comb) |
|||
while m {accu+=array(m)} |
|||
gt+=if(accu>treat->1,0) |
|||
counter++ |
|||
if counter mod d=0 then Print over str$(counter/total," #0.0%"): Refresh 1000 |
|||
} |
|||
print over str$(counter/total," #0.0%") |
|||
print |
|||
lt=total-gt |
|||
print Format$("less or equal={0:1}%, greater={1:1}%, total={2}",lt/total*100, gt/total*100, total) |
|||
} |
|||
CheckThis |
|||
</syntaxhighlight> |
|||
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
|||
<syntaxhighlight lang="mathematica">"<=: " <> ToString[#1] <> " " <> ToString[100. #1/#2] <> "%\n >: " <> |
|||
ToString[#2 - #1] <> " " <> ToString[100. (1 - #1/#2)] <> "%" &[ |
ToString[#2 - #1] <> " " <> ToString[100. (1 - #1/#2)] <> "%" &[ |
||
Count[Total /@ Subsets[Join[#1, #2], {Length@#1}], |
Count[Total /@ Subsets[Join[#1, #2], {Length@#1}], |
||
n_ /; n <= Total@#1], |
n_ /; n <= Total@#1], |
||
Binomial[Length@#1 + Length@#2, Length@#1]] &[{85, 88, 75, 66, 25, |
Binomial[Length@#1 + Length@#2, Length@#1]] &[{85, 88, 75, 66, 25, |
||
29, 83, 39, 97}, {68, 41, 10, 49, 16, 65, 32, 92, 28, 98}]</ |
29, 83, 39, 97}, {68, 41, 10, 49, 16, 65, 32, 92, 28, 98}]</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre><=: 80551 87.1972% |
<pre><=: 80551 87.1972% |
||
>: 11827 12.8028%</pre> |
>: 11827 12.8028%</pre> |
||
=={{header| |
=={{header|Nim}}== |
||
{{trans|C}} |
|||
<lang PicoLisp>(load "@lib/simul.l") # For 'subsets' |
|||
<syntaxhighlight lang="nim">import strformat |
|||
const data = [85, 88, 75, 66, 25, |
|||
(scl 2) |
|||
29, 83, 39, 97, 68, |
|||
41, 10, 49, 16, 65, |
|||
32, 92, 28, 98] |
|||
func pick(at, remain, accu, treat: int): int = |
|||
(de _stat (A) |
|||
if remain == 0: |
|||
(let (LenA (length A) SumA (apply + A)) |
|||
return if accu > treat: 1 else: 0 |
|||
(- |
|||
return pick(at - 1, remain - 1, accu + data[at - 1], treat) + |
|||
(*/ SumA LenA) |
|||
(if at > remain: pick(at - 1, remain, accu, treat) else: 0) |
|||
(de permutationTest (A B) |
|||
(let |
|||
(AB (append A B) |
|||
SumAB (apply + AB) |
|||
LenAB (length AB) |
|||
Tobs (_stat A) |
|||
Count 0 ) |
|||
(*/ |
|||
(sum |
|||
'((Perm) |
|||
(inc 'Count) |
|||
(and (>= Tobs (_stat Perm)) 1) ) |
|||
(subsets (length A) AB) ) |
|||
100.0 |
|||
Count ) ) ) |
|||
var treat = 0 |
|||
(setq |
|||
var le, gt = 0 |
|||
*TreatmentGroup (0.85 0.88 0.75 0.66 0.25 0.29 0.83 0.39 0.97) |
|||
var total = 1.0 |
|||
*ControlGroup (0.68 0.41 0.10 0.49 0.16 0.65 0.32 0.92 0.28 0.98) ) |
|||
for i in countup(0, 8): |
|||
treat += data[i] |
|||
for i in countdown(19, 11): |
|||
total *= float(i) |
|||
for i in countdown(9, 1): |
|||
total /= float(i) |
|||
gt = pick(19, 9, 0, treat) |
|||
(let N (permutationTest *TreatmentGroup *ControlGroup) |
|||
le = int(total - float(gt)) |
|||
(prinl "under = " (round N) "%, over = " (round (- 100.0 N)) "%") )</lang> |
|||
echo fmt"<= : {100.0 * float(le) / total:.6f}% {le}" |
|||
echo fmt" > : {100.0 * float(gt) / total:.6f}% {gt}"</syntaxhighlight> |
|||
{{out}} |
{{out}} |
||
<pre> |
|||
<pre>under = 87.85%, over = 12.15%</pre> |
|||
<= : 87.197168% 80551 |
|||
> : 12.802832% 11827 |
|||
</pre> |
|||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
< |
<syntaxhighlight lang="perl">#!/usr/bin/perl |
||
use warnings; |
use warnings; |
||
use strict; |
use strict; |
||
Line 738: | Line 1,261: | ||
92 |
92 |
||
28 |
28 |
||
98</ |
98</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
equal 313 0.339% |
equal 313 0.339% |
||
Line 744: | Line 1,267: | ||
greater 11827 12.803% |
greater 11827 12.803% |
||
=={{header| |
=={{header|Phix}}== |
||
{{Trans|C}} |
|||
{{works with|Rakudo|2013.05}} |
|||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">data</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">85</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">88</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">75</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">66</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">25</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">83</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">39</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">97</span><span style="color: #0000FF;">,</span> |
|||
<span style="color: #000000;">68</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">41</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">10</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">49</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">16</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">65</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">32</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">92</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">28</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">98</span> <span style="color: #0000FF;">}</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">pick</span><span style="color: #0000FF;">(</span><span style="color: #004080;">int</span> <span style="color: #000000;">at</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">int</span> <span style="color: #000000;">remain</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">int</span> <span style="color: #000000;">accu</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">int</span> <span style="color: #000000;">treat</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">remain</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">accu</span><span style="color: #0000FF;">></span><span style="color: #000000;">treat</span><span style="color: #0000FF;">?</span><span style="color: #000000;">1</span><span style="color: #0000FF;">:</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">pick</span><span style="color: #0000FF;">(</span><span style="color: #000000;">at</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">remain</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">accu</span><span style="color: #0000FF;">+</span><span style="color: #000000;">data</span><span style="color: #0000FF;">[</span><span style="color: #000000;">at</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">treat</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">+</span> |
|||
<span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">at</span><span style="color: #0000FF;">></span><span style="color: #000000;">remain</span><span style="color: #0000FF;">?</span><span style="color: #000000;">pick</span><span style="color: #0000FF;">(</span><span style="color: #000000;">at</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">remain</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">accu</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">treat</span><span style="color: #0000FF;">):</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #004080;">int</span> <span style="color: #000000;">treat</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">le</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">gt</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">total</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">;</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">9</span> <span style="color: #008080;">do</span> <span style="color: #000000;">treat</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">data</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">19</span> <span style="color: #008080;">to</span> <span style="color: #000000;">11</span> <span style="color: #008080;">by</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> <span style="color: #000000;">total</span> <span style="color: #0000FF;">*=</span> <span style="color: #000000;">i</span> <span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">9</span> <span style="color: #008080;">to</span> <span style="color: #000000;">1</span> <span style="color: #008080;">by</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> <span style="color: #000000;">total</span> <span style="color: #0000FF;">/=</span> <span style="color: #000000;">i</span> <span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #000000;">gt</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">pick</span><span style="color: #0000FF;">(</span><span style="color: #000000;">19</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">treat</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">le</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">total</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">gt</span><span style="color: #0000FF;">;</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"<= : %f%% %d\n > : %f%% %d\n"</span><span style="color: #0000FF;">,</span> |
|||
<span style="color: #0000FF;">{</span><span style="color: #000000;">100</span><span style="color: #0000FF;">*</span><span style="color: #000000;">le</span><span style="color: #0000FF;">/</span><span style="color: #000000;">total</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">le</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">100</span><span style="color: #0000FF;">*</span><span style="color: #000000;">gt</span><span style="color: #0000FF;">/</span><span style="color: #000000;">total</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">gt</span><span style="color: #0000FF;">})</span> |
|||
<!--</syntaxhighlight>--> |
|||
{{out}} |
|||
<pre> |
|||
<= : 87.197168% 80551 |
|||
> : 12.802832% 11827 |
|||
</pre> |
|||
=={{header|PicoLisp}}== |
|||
<lang perl6>proto combine (Int, @) {*} |
|||
<syntaxhighlight lang="picolisp">(load "@lib/simul.l") # For 'subsets' |
|||
multi combine (0, @) { [] } |
|||
multi combine ($, []) { () } |
|||
multi combine ($n, [$head, *@tail]) { |
|||
gather { |
|||
take [$head, @$_] for combine($n-1, @tail); |
|||
take [ @$_ ] for combine($n , @tail); |
|||
} |
|||
} |
|||
(scl 2) |
|||
sub stats ( @test, @all ) { |
|||
(([+] @test) / +@test ) - ([+] @all, (@test X* -1)) / (@all - @test) |
|||
} |
|||
(de _stat (A) |
|||
(let (LenA (length A) SumA (apply + A)) |
|||
(- |
|||
(*/ SumA LenA) |
|||
(*/ (- SumAB SumA) (- LenAB LenA)) ) ) ) |
|||
(de permutationTest (A B) |
|||
my @treated = <85 88 75 66 25 29 83 39 97>; |
|||
(let |
|||
my @control = <68 41 10 49 16 65 32 92 28 98>; |
|||
(AB (append A B) |
|||
my @all = @treated, @control; |
|||
SumAB (apply + AB) |
|||
LenAB (length AB) |
|||
Tobs (_stat A) |
|||
Count 0 ) |
|||
(*/ |
|||
(sum |
|||
'((Perm) |
|||
(inc 'Count) |
|||
(and (>= Tobs (_stat Perm)) 1) ) |
|||
(subsets (length A) AB) ) |
|||
100.0 |
|||
Count ) ) ) |
|||
(setq |
|||
my $base = stats( @treated, @all ); |
|||
*TreatmentGroup (0.85 0.88 0.75 0.66 0.25 0.29 0.83 0.39 0.97) |
|||
*ControlGroup (0.68 0.41 0.10 0.49 0.16 0.65 0.32 0.92 0.28 0.98) ) |
|||
my @trials = 0, 0, 0; |
|||
map { @trials[ 1 + ( stats( $_, @all ) <=> $base ) ]++ }, combine( +@treated, @all ); |
|||
say 'Counts: <, =, > ', @trials; |
|||
say 'Less than : %', 100 * @trials[0] / [+] @trials; |
|||
say 'Equal to : %', 100 * @trials[1] / [+] @trials; |
|||
say 'Greater than : %', 100 * @trials[2] / [+] @trials; |
|||
say 'Less or Equal: %', 100 * ( [+] @trials[0,1] ) / [+] @trials;</lang> |
|||
(let N (permutationTest *TreatmentGroup *ControlGroup) |
|||
(prinl "under = " (round N) "%, over = " (round (- 100.0 N)) "%") )</syntaxhighlight> |
|||
{{out}} |
{{out}} |
||
<pre>under = 87.85%, over = 12.15%</pre> |
|||
<pre> |
|||
Counts: <, =, > 80238 313 11827 |
|||
Less than : %86.858343 |
|||
Equal to : %0.338825 |
|||
Greater than : %12.802832 |
|||
Less or Equal: %87.197168 |
|||
</pre> |
|||
=={{header|PureBasic}}== |
=={{header|PureBasic}}== |
||
Line 800: | Line 1,344: | ||
the treatment group or the control group |
the treatment group or the control group |
||
<syntaxhighlight lang="purebasic"> |
|||
<lang PureBasic> |
|||
Define.f meanTreated,meanControl,diffInMeans |
Define.f meanTreated,meanControl,diffInMeans |
||
Line 892: | Line 1,436: | ||
Debug StrF(100*diffLessOrEqual/TotalComBinations,2)+" "+Str(diffLessOrEqual) |
Debug StrF(100*diffLessOrEqual/TotalComBinations,2)+" "+Str(diffLessOrEqual) |
||
Debug StrF(100*diffGreater /TotalComBinations,2)+" "+Str(diffGreater) |
Debug StrF(100*diffGreater /TotalComBinations,2)+" "+Str(diffGreater) |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
Line 899: | Line 1,443: | ||
12.80 11827 |
12.80 11827 |
||
</pre> |
</pre> |
||
=={{header|Python}}== |
=={{header|Python}}== |
||
{{trans|Tcl}} |
{{trans|Tcl}} |
||
< |
<syntaxhighlight lang="python">from itertools import combinations as comb |
||
def statistic(ab, a): |
def statistic(ab, a): |
||
Line 922: | Line 1,465: | ||
controlGroup = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
controlGroup = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
||
under = permutationTest(treatmentGroup, controlGroup) |
under = permutationTest(treatmentGroup, controlGroup) |
||
print("under=%.2f%%, over=%.2f%%" % (under, 100. - under))</ |
print("under=%.2f%%, over=%.2f%%" % (under, 100. - under))</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>under=89.11%, over=10.89%</pre> |
<pre>under=89.11%, over=10.89%</pre> |
||
The above solution does a different thing than the other solutions. I'm not really sure why. If you want to do the same thing as the other solutions, then this is the solution: |
The above solution does a different thing than the other solutions. I'm not really sure why. If you want to do the same thing as the other solutions, then this is the solution: |
||
< |
<syntaxhighlight lang="python">from itertools import combinations as comb |
||
def permutationTest(a, b): |
def permutationTest(a, b): |
||
Line 941: | Line 1,484: | ||
controlGroup = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
controlGroup = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
||
under = permutationTest(treatmentGroup, controlGroup) |
under = permutationTest(treatmentGroup, controlGroup) |
||
print("under=%.2f%%, over=%.2f%%" % (under, 100. - under))</ |
print("under=%.2f%%, over=%.2f%%" % (under, 100. - under))</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>under=87.20%, over=12.80%</pre> |
<pre>under=87.20%, over=12.80%</pre> |
||
=={{header|REXX}}== |
|||
This REXX program is modeled after the C version, with some generalizations and optimization. |
|||
<lang rexx>/*REXX program does a permutation test on N + M subjects (volunteers):*/ |
|||
/* ↑ ↑ */ |
|||
/* │ │ */ |
|||
/* │ └─────control population*/ |
|||
/* └────────treatment population*/ |
|||
n=9 |
|||
data=85 88 75 66 25 29 83 39 97 68 41 10 49 16 65 32 92 28 98 |
|||
w=words(data); m=w-n |
|||
say 'volunteer population given treatment:' right(n,length(w)) |
|||
say ' control population given a placebo:' right(m,length(w)) |
|||
say |
|||
say 'treatment population efficacy % (percentages):' subword(data,1,n) |
|||
say ' control population placebo % (percentages):' subword(data,n+1) |
|||
say |
|||
do v= 0 for w ; #.v=word(data,v+1) ; end |
|||
treat=0; do i= 0 to n-1 ; treat=treat+#.i ; end |
|||
total=1; do j=19 to m+1 by -1 ; total=total*j ; end |
|||
do k= 9 to 1 by -1 ; total=total/k ; end |
|||
gt=picker(n+m, n, 0) |
|||
le=total-gt |
|||
say "<= " format(100*le/total,,3)'%' le /*show 3 decimal places.*/ |
|||
say " > " format(100*gt/total,,3)'%' gt |
|||
exit /*stick a fork in it, we're done.*/ |
|||
/*──────────────────────────────────PICKER subroutine───────────────────*/ |
|||
picker: procedure expose #. treat; parse arg it,rest,eff |
|||
if rest==0 then return eff>treat |
|||
if it>rest then q=picker(it-1, rest, eff) |
|||
else q=0 |
|||
itP=it-1 |
|||
return picker(itP, rest-1, eff+#.itP) + q</lang> |
|||
{{out}} using the default input: |
|||
<pre> |
|||
volunteer population given treatment: 9 |
|||
control population given a placebo: 10 |
|||
treatment population efficacy % (percentages): 85 88 75 66 25 29 83 39 97 |
|||
control population placebo % (percentages): 68 41 10 49 16 65 32 92 28 98 |
|||
<= 87.197% 80551 |
|||
> 12.803% 11827 |
|||
</pre> |
|||
=={{header|R}}== |
=={{header|R}}== |
||
< |
<syntaxhighlight lang="r">permutation.test <- function(treatment, control) { |
||
perms <- combinations(length(treatment)+length(control), |
perms <- combinations(length(treatment)+length(control), |
||
length(treatment), |
length(treatment), |
||
Line 998: | Line 1,497: | ||
p <- mean(rowMeans(perms) <= mean(treatment)) |
p <- mean(rowMeans(perms) <= mean(treatment)) |
||
c(under=p, over=(1-p)) |
c(under=p, over=(1-p)) |
||
}</ |
}</syntaxhighlight> |
||
< |
<syntaxhighlight lang="r">> permutation.test(c(85, 88, 75, 66, 25, 29, 83, 39, 97), |
||
+ c(68, 41, 10, 49, 16, 65, 32, 92, 28, 98)) |
+ c(68, 41, 10, 49, 16, 65, 32, 92, 28, 98)) |
||
under over |
under over |
||
0.8719717 0.1280283 </ |
0.8719717 0.1280283 </syntaxhighlight> |
||
=={{header|Racket}}== |
=={{header|Racket}}== |
||
{{trans|Common Lisp}} |
{{trans|Common Lisp}} |
||
< |
<syntaxhighlight lang="racket">#lang racket/base |
||
(define-syntax-rule (inc! x) |
(define-syntax-rule (inc! x) |
||
Line 1,034: | Line 1,532: | ||
more (real->decimal-string (* 100. (/ more sum)) 2) |
more (real->decimal-string (* 100. (/ more sum)) 2) |
||
leq (real->decimal-string (* 100. (/ leq sum)) 2)))) |
leq (real->decimal-string (* 100. (/ leq sum)) 2)))) |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
Line 1,040: | Line 1,538: | ||
<=: 80551 87.20% |
<=: 80551 87.20% |
||
>: 11827 12.80% |
>: 11827 12.80% |
||
</pre> |
|||
=={{header|Raku}}== |
|||
(formerly Perl 6) |
|||
The use of <code>.race</code> to allow concurrent calculations means that multiple 'workers' will be updating <code>@trials</code> simultaneously. To avoid race conditions, the <code>⚛++</code> operator is used, which guarantees safe updates without the use of locks. That is turn requires declaring that array as being composed of <code>atomicint</code>. |
|||
{{works with|Rakudo|2018.09}} |
|||
<syntaxhighlight lang="raku" line>sub stats ( @test, @all ) { |
|||
([+] @test / +@test) - ([+] flat @all, (@test X* -1)) / @all - @test |
|||
} |
|||
my int @treated = <85 88 75 66 25 29 83 39 97>; |
|||
my int @control = <68 41 10 49 16 65 32 92 28 98>; |
|||
my int @all = flat @treated, @control; |
|||
my $base = stats( @treated, @all ); |
|||
my atomicint @trials[3] = 0, 0, 0; |
|||
@all.combinations(+@treated).race.map: { @trials[ 1 + ( stats( $_, @all ) <=> $base ) ]⚛++ } |
|||
say 'Counts: <, =, > ', @trials; |
|||
say 'Less than : %', 100 * @trials[0] / [+] @trials; |
|||
say 'Equal to : %', 100 * @trials[1] / [+] @trials; |
|||
say 'Greater than : %', 100 * @trials[2] / [+] @trials; |
|||
say 'Less or Equal: %', 100 * ( [+] @trials[0,1] ) / [+] @trials;</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Counts: <, =, > 80238 313 11827 |
|||
Less than : %86.858343 |
|||
Equal to : %0.338825 |
|||
Greater than : %12.802832 |
|||
Less or Equal: %87.197168 |
|||
</pre> |
|||
=={{header|REXX}}== |
|||
This REXX program is modeled after the '''C''' version, with some generalizations and optimization added. |
|||
<syntaxhighlight lang="rexx">/*REXX program performs a permutation test on N + M subjects (volunteers): */ |
|||
/* ↑ ↑ */ |
|||
/* │ │ */ |
|||
/* │ └─────control population. */ |
|||
/* └────────treatment population. */ |
|||
n= 9 /*define the number of the control pop.*/ |
|||
data= 85 88 75 66 25 29 83 39 97 68 41 10 49 16 65 32 92 28 98 |
|||
w= words(data); m= w - n /*w: volunteers + control population*/ |
|||
L= length(w) /*L: used to align some output numbers*/ |
|||
say '# of volunteers & control population: ' w |
|||
say 'volunteer population given treatment: ' right(n, L) |
|||
say ' control population given a placebo: ' right(m, L) |
|||
say |
|||
say 'treatment population efficacy % (percentages): ' subword(data, 1, n) |
|||
say ' control population placebo % (percentages): ' subword(data, n+1 ) |
|||
say |
|||
do v= 0 for w ; #.v= word(data, v+1) ; end |
|||
treat= 0; do i= 0 to n-1 ; treat= treat + #.i ; end |
|||
tot= 1; do j= w to m+1 by -1 ; tot= tot * j ; end |
|||
do k=w%2 to 1 by -1 ; tot= tot / k ; end |
|||
GT= picker(n+m, n, 0) /*compute the GT value from PICKER func*/ |
|||
LE= tot - GT /* " " LE " via subtraction.*/ |
|||
say "<= " format(100 * LE / tot, ,3)'%' LE /*display number with 3 decimal places.*/ |
|||
say " > " format(100 * GT / tot, ,3)'%' GT /* " " " " " " */ |
|||
exit /*stick a fork in it, we're all done. */ |
|||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
|||
picker: procedure expose #. treat; parse arg it,rest,eff /*get the arguments.*/ |
|||
if rest==0 then return eff > treat /*is REST = to zero?*/ |
|||
if it>rest then q= picker(it-1, rest, eff) /*maybe recurse. */ |
|||
else q= 0 |
|||
itP= it - 1 /*set temporary var.*/ |
|||
return picker(itP, rest - 1, eff + #.itP) + q /*recurse. */</syntaxhighlight> |
|||
{{out|output|text= when using the default input:}} |
|||
<pre> |
|||
# of volunteers & control population: 19 |
|||
volunteer population given treatment: 9 |
|||
control population given a placebo: 10 |
|||
treatment population efficacy % (percentages): 85 88 75 66 25 29 83 39 97 |
|||
control population placebo % (percentages): 68 41 10 49 16 65 32 92 28 98 |
|||
<= 87.197% 80551 |
|||
> 12.803% 11827 |
|||
</pre> |
</pre> |
||
=={{header|Ruby}}== |
=={{header|Ruby}}== |
||
{{trans|Python}} |
{{trans|Python}} |
||
< |
<syntaxhighlight lang="ruby">def statistic(ab, a) |
||
sumab, suma = ab.inject(:+).to_f, a.inject(:+).to_f |
sumab, suma = ab.inject(:+).to_f, a.inject(:+).to_f |
||
suma / a.size - (sumab - suma) / (ab.size - a.size) |
suma / a.size - (sumab - suma) / (ab.size - a.size) |
||
Line 1,063: | Line 1,643: | ||
controlGroup = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
controlGroup = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
||
under = permutationTest(treatmentGroup, controlGroup) |
under = permutationTest(treatmentGroup, controlGroup) |
||
puts "under=%.2f%%, over=%.2f%%" % [under, 100 - under]</ |
puts "under=%.2f%%, over=%.2f%%" % [under, 100 - under]</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,069: | Line 1,649: | ||
under=87.20%, over=12.80% |
under=87.20%, over=12.80% |
||
</pre> |
</pre> |
||
=={{header|Rust}}== |
|||
<syntaxhighlight lang="rust"> |
|||
fn main() { |
|||
let treatment = vec![85, 88, 75, 66, 25, 29, 83, 39, 97]; |
|||
let control = vec![68, 41, 10, 49, 16, 65, 32, 92, 28, 98]; |
|||
let mut data_set = control.clone(); |
|||
data_set.extend_from_slice(&treatment); |
|||
let greater = combinations(treatment.iter().sum(), treatment.len() as i64, &data_set) as f64; |
|||
let lesser = combinations(control.iter().sum(), control.len() as i64, &data_set) as f64; |
|||
let total = binomial(data_set.len() as i64, treatment.len() as i64) as f64; |
|||
println!("<= : {}%", (lesser / total * 100.0)); |
|||
println!(" > : {}%", (greater / total * 100.0)); |
|||
} |
|||
fn factorial(x: i64) -> i64 { |
|||
let mut product = 1; |
|||
for a in 1..(x + 1) { |
|||
product *= a; |
|||
} |
|||
product |
|||
} |
|||
fn binomial(n: i64, k: i64) -> i64 { |
|||
let numerator = factorial(n); |
|||
let denominator = factorial(k) * factorial(n - k); |
|||
numerator / denominator |
|||
} |
|||
fn combinations(total: i64, number: i64, data: &[i64]) -> i64 { |
|||
if total < 0 { |
|||
return binomial(data.len() as i64, number); |
|||
} |
|||
if number == 0 { |
|||
return 0; |
|||
} |
|||
if number > data.len() as i64 { |
|||
return 0; |
|||
} |
|||
if number == data.len() as i64 { |
|||
if total < data.iter().sum() { |
|||
return 1; |
|||
} else { |
|||
return 0; |
|||
} |
|||
} |
|||
let tail = &data[1..]; |
|||
combinations(total - data[0], number - 1, &tail) + combinations(total, number, &tail) |
|||
} |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
<= : 86.8583428954946% |
|||
> : 12.80283184307952% |
|||
</pre> |
|||
=={{header|Scala}}== |
|||
===Imperative version (Ugly, side effects)=== |
|||
{{Out}}Best seen running in your browser either by [https://scalafiddle.io/sf/69o8tJX/0 ScalaFiddle (ES aka JavaScript, non JVM)] or [https://scastie.scala-lang.org/jpUyNPetRXabZkEKc68FAA Scastie (remote JVM)]. |
|||
<syntaxhighlight lang="scala">object PermutationTest extends App { |
|||
private val data = |
|||
Array(85, 88, 75, 66, 25, 29, 83, 39, 97, 68, 41, 10, 49, 16, 65, 32, 92, 28, 98) |
|||
private var (total, treat) = (1.0, 0) |
|||
private def pick(at: Int, remain: Int, accu: Int, treat: Int): Int = { |
|||
if (remain == 0) return if (accu > treat) 1 else 0 |
|||
pick(at - 1, remain - 1, accu + data(at - 1), treat) + |
|||
(if (at > remain) pick(at - 1, remain, accu, treat) else 0) |
|||
} |
|||
for (i <- 0 to 8) treat += data(i) |
|||
for (j <- 19 to 11 by -1) total *= j |
|||
for (g <- 9 to 1 by -1) total /= g |
|||
private val gt = pick(19, 9, 0, treat) |
|||
private val le = (total - gt).toInt |
|||
println(f"<= : ${100.0 * le / total}%f%% ${le}%d") |
|||
println(f" > : ${100.0 * gt / total}%f%% ${gt}%d") |
|||
}</syntaxhighlight> |
|||
=={{header|Seed7}}== |
=={{header|Seed7}}== |
||
< |
<syntaxhighlight lang="seed7">$ include "seed7_05.s7i"; |
||
include "float.s7i"; |
include "float.s7i"; |
||
Line 1,109: | Line 1,778: | ||
writeln("<= : " <& 100.0 * flt(le) / flt(total) digits 6 <& "% " <& le); |
writeln("<= : " <& 100.0 * flt(le) / flt(total) digits 6 <& "% " <& le); |
||
writeln(" > : " <& 100.0 * flt(gt) / flt(total) digits 6 <& "% " <& gt); |
writeln(" > : " <& 100.0 * flt(gt) / flt(total) digits 6 <& "% " <& gt); |
||
end func;</ |
end func;</syntaxhighlight> |
||
{{out}} |
|||
<= : 87.197168% 80551 |
|||
> : 12.802832% 11827 |
|||
=={{header|Sidef}}== |
|||
{{trans|Ruby}} |
|||
<syntaxhighlight lang="ruby">func statistic(ab, a) { |
|||
var(sumab, suma) = (ab.sum, a.sum) |
|||
suma/a.size - ((sumab-suma) / (ab.size-a.size)) |
|||
} |
|||
func permutationTest(a, b) { |
|||
var ab = (a + b) |
|||
var tobs = statistic(ab, a) |
|||
var under = (var count = 0) |
|||
ab.combinations(a.len, {|*perm| |
|||
statistic(ab, perm) <= tobs && (under += 1) |
|||
count += 1 |
|||
}) |
|||
under * 100 / count |
|||
} |
|||
var treatmentGroup = [85, 88, 75, 66, 25, 29, 83, 39, 97] |
|||
var controlGroup = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
|||
var under = permutationTest(treatmentGroup, controlGroup) |
|||
say ("under=%.2f%%, over=%.2f%%" % (under, 100 - under))</syntaxhighlight> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
under=87.20%, over=12.80% |
|||
<= : 87.197168% 80551 |
|||
> : 12.802832% 11827 |
|||
</pre> |
</pre> |
||
=={{header|Tcl}}== |
=={{header|Tcl}}== |
||
< |
<syntaxhighlight lang="tcl">package require Tcl 8.5 |
||
# Difference of means; note that the first list must be the concatenation of |
# Difference of means; note that the first list must be the concatenation of |
||
Line 1,164: | Line 1,857: | ||
set count [tcl::mathfunc::double $count] |
set count [tcl::mathfunc::double $count] |
||
list [expr {$overcount / $count}] [expr {$undercount / $count}] |
list [expr {$overcount / $count}] [expr {$undercount / $count}] |
||
}</ |
}</syntaxhighlight> |
||
Demonstration code: |
Demonstration code: |
||
< |
<syntaxhighlight lang="tcl">set treatmentGroup {0.85 0.88 0.75 0.66 0.25 0.29 0.83 0.39 0.97} |
||
set controlGroup {0.68 0.41 0.10 0.49 0.16 0.65 0.32 0.92 0.28 0.98} |
set controlGroup {0.68 0.41 0.10 0.49 0.16 0.65 0.32 0.92 0.28 0.98} |
||
lassign [permutationTest $treatmentGroup $controlGroup] over under |
lassign [permutationTest $treatmentGroup $controlGroup] over under |
||
puts [format "under=%.2f%%, over=%.2f%%" [expr {$under*100}] [expr {$over*100}]]</ |
puts [format "under=%.2f%%, over=%.2f%%" [expr {$under*100}] [expr {$over*100}]]</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>under=86.90%, over=13.10%</pre> |
<pre>under=86.90%, over=13.10%</pre> |
||
=={{header|Ursala}}== |
=={{header|Ursala}}== |
||
< |
<syntaxhighlight lang="ursala">#import std |
||
#import nat |
#import nat |
||
#import flo |
#import flo |
||
Line 1,189: | Line 1,882: | ||
#show+ |
#show+ |
||
t = --* *-'%'@lrNCC printf/$'%0.2f' times/$100. f(treatment_group,control_group)</ |
t = --* *-'%'@lrNCC printf/$'%0.2f' times/$100. f(treatment_group,control_group)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
12.80% |
12.80% |
||
87.20%</pre> |
87.20%</pre> |
||
=={{header|Wren}}== |
|||
{{trans|C}} |
|||
{{libheader|Wren-fmt}} |
|||
<syntaxhighlight lang="wren">import "./fmt" for Fmt |
|||
var data = [85, 88, 75, 66, 25, 29, 83, 39, 97, 68, 41, 10, 49, 16, 65, 32, 92, 28, 98] |
|||
var pick // recursive |
|||
pick = Fn.new { |at, remain, accu, treat| |
|||
if (remain == 0) return (accu > treat) ? 1 : 0 |
|||
return pick.call(at-1, remain-1, accu + data[at-1], treat) + |
|||
((at > remain) ? pick.call(at-1, remain, accu, treat) : 0) |
|||
} |
|||
var treat = 0 |
|||
var total = 1 |
|||
for (i in 0..8) treat = treat + data[i] |
|||
for (i in 19..11) total = total * i |
|||
for (i in 9..1) total = total / i |
|||
var gt = pick.call(19, 9, 0, treat) |
|||
var le = (total - gt).truncate |
|||
Fmt.print("<= : $f\% $d", 100 * le / total, le) |
|||
Fmt.print(" > : $f\% $d", 100 * gt / total, gt)</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
<= : 87.197168% 80551 |
|||
> : 12.802832% 11827 |
|||
</pre> |
|||
=={{header|XPL0}}== |
|||
{{trans|C}} |
|||
<syntaxhighlight lang "XPL0">include xpllib; \For Print |
|||
int Data; |
|||
func Pick(At, Remain, Accu, Treat); |
|||
int At, Remain, Accu, Treat; |
|||
[if Remain = 0 then return if Accu > Treat then 1 else 0; |
|||
return Pick(At-1, Remain-1, Accu + Data(At-1), Treat) + |
|||
(if At > Remain then Pick(At-1, Remain, Accu, Treat) else 0); |
|||
]; |
|||
int Treat, GT, LE, I; |
|||
real Total; |
|||
[Treat:= 0; Total:= 1.; |
|||
Data:= [85, 88, 75, 66, 25, 29, 83, 39, 97, |
|||
68, 41, 10, 49, 16, 65, 32, 92, 28, 98]; |
|||
for I:= 0 to 8 do Treat:= Treat + Data(I); |
|||
for I:= 19 downto 11 do Total:= Total * float(I); |
|||
for I:= 9 downto 1 do Total:= Total / float(I); |
|||
GT:= Pick(19, 9, 0, Treat); |
|||
LE:= (fix(Total) - GT); |
|||
Print("<= : %f% %d\n", 100. * float(LE) / Total, LE); |
|||
Print(" > : %f% %d\n", 100. * float(GT) / Total, GT); |
|||
]</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
<= : 87.19717% 80551 |
|||
> : 12.80283% 11827 |
|||
</pre> |
|||
=={{header|zkl}}== |
=={{header|zkl}}== |
||
A solution that is not going to scale gracefully at all. |
A solution that is not going to scale gracefully at all. |
||
{{trans|D}} |
{{trans|D}} |
||
< |
<syntaxhighlight lang="zkl">fcn permutationTest(a,b){ |
||
ab := a.extend(b); |
ab := a.extend(b); |
||
tObs := a.sum(0); |
tObs := a.sum(0); |
||
Line 1,209: | Line 1,964: | ||
controlGroup := T(68, 41, 10, 49, 16, 65, 32, 92, 28, 98); |
controlGroup := T(68, 41, 10, 49, 16, 65, 32, 92, 28, 98); |
||
under := permutationTest(treatmentGroup, controlGroup); |
under := permutationTest(treatmentGroup, controlGroup); |
||
println("Under =%6.2f%%\nOver =%6.2f%%".fmt(under, 100.0 - under));</ |
println("Under =%6.2f%%\nOver =%6.2f%%".fmt(under, 100.0 - under));</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
Latest revision as of 12:30, 24 January 2024
You are encouraged to solve this task according to the task description, using any language you may know.
A new medical treatment was tested on a population of volunteers, with each volunteer randomly assigned either to a group of treatment subjects, or to a group of control subjects.
Members of the treatment group were given the treatment, and members of the control group were given a placebo. The effect of the treatment or placebo on each volunteer was measured and reported in this table.
Table of experimental results Treatment group Control group 85 68 88 41 75 10 66 49 25 16 29 65 83 32 39 92 97 28 98
Write a program that performs a permutation test to judge whether the treatment had a significantly stronger effect than the placebo.
- Do this by considering every possible alternative assignment from the same pool of volunteers to a treatment group of size and a control group of size (i.e., the same group sizes used in the actual experiment but with the group members chosen differently), while assuming that each volunteer's effect remains constant regardless.
- Note that the number of alternatives will be the binomial coefficient .
- Compute the mean effect for each group and the difference in means between the groups in every case by subtracting the mean of the control group from the mean of the treatment group.
- Report the percentage of alternative groupings for which the difference in means is less or equal to the actual experimentally observed difference in means, and the percentage for which it is greater.
- Note that they should sum to 100%.
Extremely dissimilar values are evidence of an effect not entirely due
to chance, but your program need not draw any conclusions.
You may assume the experimental data are known at compile time if
that's easier than loading them at run time. Test your solution on the
data given above.
11l
V data = [85, 88, 75, 66, 25, 29, 83, 39, 97,
68, 41, 10, 49, 16, 65, 32, 92, 28, 98]
F pick(at, remain, accu, treat)
I remain == 0
R I accu > treat {1} E 0
R pick(at - 1, remain - 1, accu + :data[at - 1], treat) +
(I at > remain {pick(at - 1, remain, accu, treat)} E 0)
V treat = 0
V total = 1.0
L(i) 0..8
treat += data[i]
L(i) (19..11).step(-1)
total *= i
L(i) (9..1).step(-1)
total /= i
V gt = pick(19, 9, 0, treat)
V le = Int(total - gt)
print(‘<= : #.6% #.’.format(100 * le / total, le))
print(‘ > : #.6% #.’.format(100 * gt / total, gt))
- Output:
<= : 87.197168% 80551 > : 12.802832% 11827
Ada
with Ada.Text_IO; with Iterate_Subsets;
procedure Permutation_Test is
type Group_Type is array(Positive range <>) of Positive;
Treat_Group: constant Group_Type := (85, 88, 75, 66, 25, 29, 83, 39, 97);
Ctrl_Group: constant Group_Type := (68, 41, 10, 49, 16, 65, 32, 92, 28, 98);
package Iter is new Iterate_Subsets(Treat_Group'Length, Ctrl_Group'Length);
Full_Group: constant Group_Type(1 .. Iter.All_Elements)
:= Treat_Group & Ctrl_Group;
function Mean(S: Iter.Subset) return Float is
Sum: Natural := 0;
begin
for I in S'Range loop
Sum := Sum + Full_Group(S(I));
end loop;
return Float(Sum)/Float(S'Length);
end Mean;
package FIO is new Ada.Text_IO.Float_IO(Float);
T_Avg: Float := Mean(Iter.First);
S_Avg: Float;
S: Iter.Subset := Iter.First;
Equal: Positive := 1; -- Mean(Iter'First) = Mean(Iter'First)
Higher: Natural := 0;
Lower: Natural := 0;
begin -- Permutation_Test;
-- first, count the subsets with a higher, an equal or a lower mean
loop
Iter.Next(S);
S_Avg := Mean(S);
if S_Avg = T_Avg then
Equal := Equal + 1;
elsif S_Avg >= T_Avg then
Higher := Higher + 1;
else
Lower := Lower + 1;
end if;
exit when Iter.Last(S);
end loop;
-- second, output the results
declare
use Ada.Text_IO;
Sum: Float := Float(Higher + Equal + Lower);
begin
Put("Less or Equal: ");
FIO.Put(100.0*Float(Lower+Equal) / Sum, Fore=>3, Aft=>1, Exp=>0);
Put(Integer'Image(Lower+Equal));
New_Line;
Put("More: ");
FIO.Put(100.0*Float(Higher) / Sum, Fore=>3, Aft=>1, Exp=>0);
Put(Integer'Image(Higher));
New_Line;
end;
end Permutation_Test;
This solution uses an auxiliary package Iterate_Subsets. Here is the Spec:
generic
Subset_Size, More_Elements: Positive;
package Iterate_Subsets is
All_Elements: Positive := Subset_Size + More_Elements;
subtype Index is Integer range 1 .. All_Elements;
type Subset is array (1..Subset_Size) of Index;
-- iterate over all subsets of size Subset_Size
-- from the set {1, 2, ..., All_Element}
function First return Subset;
procedure Next(S: in out Subset);
function Last(S: Subset) return Boolean;
end Iterate_Subsets;
And here is the implementation:
package body Iterate_Subsets is
function First return Subset is
S: Subset;
begin
for I in S'Range loop
S(I) := I;
end loop;
return S;
end First;
procedure Next(S: in out Subset) is
I: Natural := S'Last;
begin
if S(I) < Index'Last then
S(I) := S(I) + 1;
else
while S(I-1)+1 = S(I) loop
I := I - 1;
end loop;
S(I-1) := S(I-1) + 1;
for J in I .. S'Last loop
S(J) := S(J-1) + 1;
end loop;
end if;
return;
end Next;
function Last(S: Subset) return Boolean is
begin
return S(S'First) = Index'Last-S'Length+1;
end Last;
end Iterate_Subsets;
- Output:
Less or Equal: 87.2 80551 More: 12.8 11827
AWK
# syntax: GAWK -f PERMUTATION_TEST.AWK
# converted from C
BEGIN {
# "treatment..................control......................"
n = split("85,88,75,66,25,29,83,39,97,68,41,10,49,16,65,32,92,28,98",data,",")
for (i=1; i<=n; i++) { # make AWK array look like a C array
data[i-1] = data[i]
}
delete data[n]
total = 1
for (i=0; i<9; i++) { treat += data[i] }
for (i=19; i>10; i--) { total *= i }
for (i=9; i>0; i--) { total /= i }
gt = pick(19,9,0,treat)
le = total - gt
printf("<= : %9.6f%% %6d\n",100*le/total,le)
printf(" > : %9.6f%% %6d\n",100*gt/total,gt)
exit(0)
}
function pick(at,remain,accu,treat) {
if (!remain) {
return (accu > treat) ? 1 : 0
}
return pick(at-1,remain-1,accu+data[at-1],treat) + ( (at > remain) ? pick(at-1,remain,accu,treat) : 0 )
}
- Output:
<= : 87.197168% 80551 > : 12.802832% 11827
BBC BASIC
ntreated% = 9
nplacebo% = 10
DIM results%(ntreated% + nplacebo% - 1)
results%() = 85, 88, 75, 66, 25, 29, 83, 39, 97, \ REM treated group
\ 68, 41, 10, 49, 16, 65, 32, 92, 28, 98 : REM placebo group
greater% = 0
FOR comb% = 0 TO 2^(ntreated%+nplacebo%)-1
IF FNnbits(comb%) = ntreated% THEN
tsum% = 0 : psum% = 0
FOR b% = 0 TO ntreated%+nplacebo%-1
IF comb% AND 2^b% THEN
tsum% += results%(b%)
ELSE
psum% += results%(b%)
ENDIF
NEXT
meandiff = tsum%/ntreated% - psum%/nplacebo%
IF comb% = 2^ntreated% - 1 THEN
actual = meandiff
ELSE
greater% -= meandiff > actual
groups% += 1
ENDIF
ENDIF
NEXT
percent = 100 * greater%/groups%
PRINT "Percentage groupings <= actual experiment: "; 100 - percent
PRINT "Percentage groupings > actual experiment: "; percent
END
DEF FNnbits(N%)
N% -= N% >>> 1 AND &55555555
N% = (N% AND &33333333) + (N% >>> 2 AND &33333333)
N% = (N% + (N% >>> 4)) AND &0F0F0F0F
N% += N% >>> 8 : N% += N% >>> 16
= N% AND &7F
- Output:
Percentage groupings <= actual experiment: 87.1970296 Percentage groupings > actual experiment: 12.8029704
C
#include <stdio.h>
int data[] = { 85, 88, 75, 66, 25, 29, 83, 39, 97,
68, 41, 10, 49, 16, 65, 32, 92, 28, 98 };
int pick(int at, int remain, int accu, int treat)
{
if (!remain) return (accu > treat) ? 1 : 0;
return pick(at - 1, remain - 1, accu + data[at - 1], treat) +
( at > remain ? pick(at - 1, remain, accu, treat) : 0 );
}
int main()
{
int treat = 0, i;
int le, gt;
double total = 1;
for (i = 0; i < 9; i++) treat += data[i];
for (i = 19; i > 10; i--) total *= i;
for (i = 9; i > 0; i--) total /= i;
gt = pick(19, 9, 0, treat);
le = total - gt;
printf("<= : %f%% %d\n > : %f%% %d\n",
100 * le / total, le, 100 * gt / total, gt);
return 0;
}
Output:
<= : 87.197168% 80551
> : 12.802832% 11827
C#
using System;
using System.Collections.Generic;
namespace PermutationTest {
class Program {
static readonly List<int> DATA = new List<int>{
85, 88, 75, 66, 25, 29, 83, 39, 97,
68, 41, 10, 49, 16, 65, 32, 92, 28, 98
};
static int Pick(int at, int remain, int accu, int treat) {
if (remain == 0) {
return (accu > treat) ? 1 : 0;
}
return Pick(at - 1, remain - 1, accu + DATA[at - 1], treat)
+ ((at > remain) ? Pick(at - 1, remain, accu, treat) : 0);
}
static void Main() {
int treat = 0;
double total = 1.0;
for (int i = 0; i <= 8; i++) {
treat += DATA[i];
}
for (int i = 19; i >= 11; i--) {
total *= i;
}
for (int i = 9; i >= 1; --i) {
total /= i;
}
int gt = Pick(19, 9, 0, treat);
int le = (int) (total - gt);
Console.WriteLine("<= {0}% {1}", 100.0 * le / total, le);
Console.WriteLine(" > {0}% {1}", 100.0 * gt / total, gt);
}
}
}
- Output:
<= 87.1971681569205% 80551 > 12.8028318430795% 11827
C++
This is a translaion of C
#include<iostream>
#include<vector>
#include<numeric>
#include<functional>
class
{
public:
int64_t operator()(int n, int k){ return partial_factorial(n, k) / factorial(n - k);}
private:
int64_t partial_factorial(int from, int to) { return from == to ? 1 : from * partial_factorial(from - 1, to); }
int64_t factorial(int n) { return n == 0 ? 1 : n * factorial(n - 1);}
}combinations;
int main()
{
static constexpr int treatment = 9;
const std::vector<int> data{ 85, 88, 75, 66, 25, 29, 83, 39, 97,
68, 41, 10, 49, 16, 65, 32, 92, 28, 98 };
int treated = std::accumulate(data.begin(), data.begin() + treatment, 0);
std::function<int (int, int, int)> pick;
pick = [&](int n, int from, int accumulated)
{
if(n == 0)
return accumulated > treated ? 1 : 0;
else
return pick(n - 1, from - 1, accumulated + data[from - 1]) +
(from > n ? pick(n, from - 1, accumulated) : 0);
};
int total = combinations(data.size(), treatment);
int greater = pick(treatment, data.size(), 0);
int lesser = total - greater;
std::cout << "<= : " << 100.0 * lesser / total << "% " << lesser << std::endl
<< " > : " << 100.0 * greater / total << "% " << greater << std::endl;
}
Output:
<= : 87.197168% 80551
> : 12.802832% 11827
Common Lisp
(defun perm-test (s1 s2)
(let ((more 0) (leq 0)
(all-data (append s1 s2))
(thresh (apply #'+ s1)))
(labels
((recur (data sum need avail)
(cond ((zerop need) (if (>= sum thresh)
(incf more)
(incf leq)))
((>= avail need)
(recur (cdr data) sum need (1- avail))
(recur (cdr data) (+ sum (car data)) (1- need) (1- avail))))))
(recur all-data 0 (length s1) (length all-data))
(cons more leq))))
(let* ((a (perm-test '(68 41 10 49 16 65 32 92 28 98)
'(85 88 75 66 25 29 83 39 97)))
(x (car a))
(y (cdr a))
(s (+ x y)))
(format t "<=: ~a ~6f%~% >: ~a ~6f%~%"
x (* 100e0 (/ x s))
y (* 100e0 (/ y s))))
output
<=: 80551 87.197%
>: 11827 12.803%
D
import std.stdio, std.algorithm, std.array, combinations3;
auto permutationTest(T)(in T[] a, in T[] b) pure nothrow @safe {
immutable tObs = a.sum;
auto combs = combinations!false(a ~ b, a.length);
immutable under = combs.count!(perm => perm.sum <= tObs);
return under * 100.0 / combs.length;
}
void main() {
immutable treatmentGroup = [85, 88, 75, 66, 25, 29, 83, 39, 97];
immutable controlGroup = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98];
immutable under = permutationTest(treatmentGroup, controlGroup);
writefln("Under =%6.2f%%\nOver =%6.2f%%", under, 100.0 - under);
}
- Output:
Under = 87.20% Over = 12.80%
Alternative version:
void main() @safe {
import std.stdio, std.algorithm, std.range;
immutable treatment = [85, 88, 75, 66, 25, 29, 83, 39, 97];
immutable control = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98];
immutable both = treatment ~ control;
immutable sTreat = treatment.sum;
T pick(T)(in size_t at, in size_t remain, in T accu) pure nothrow @safe @nogc {
if (remain == 0)
return accu > sTreat;
return pick(at - 1, remain - 1, accu + both[at - 1]) +
(at > remain ? pick(at - 1, remain, accu) : 0);
}
alias mul = reduce!q{a * b};
immutable t = mul(1.0, iota(both.length, treatment.length + 1, -1))
.reduce!q{a / b}(iota(treatment.length, 0, -1));
immutable gt = pick(both.length, treatment.length, 0);
immutable le = cast(int)(t - gt);
writefln(" > : %2.2f%% %d", 100.0 * gt / t, gt);
writefln("<= : %2.2f%% %d", 100.0 * le / t, le);
}
- Output:
> : 12.80% 11827 <= : 87.20% 80551
Delphi
program Permutation_test;
{$APPTYPE CONSOLE}
uses
System.SysUtils;
procedure Comb(n, m: Integer; emit: TProc<TArray<Integer>>);
var
s: TArray<Integer>;
last: Integer;
procedure rc(i, next: Integer);
begin
for var j := next to n - 1 do
begin
s[i] := j;
if i = last then
emit(s)
else
rc(i + 1, j + 1);
end;
end;
begin
SetLength(s, m);
last := m - 1;
rc(0, 0);
end;
begin
var tr: TArray<Integer> := [85, 88, 75, 66, 25, 29, 83, 39, 97];
var ct: TArray<Integer> := [68, 41, 10, 49, 16, 65, 32, 92, 28, 98];
// collect all results in a single list
var all: TArray<Integer> := concat(tr, ct);
// compute sum of all data, useful as intermediate result
var sumAll := 0;
for var r in all do
inc(sumAll, r);
// closure for computing scaled difference.
// compute results scaled by len(tr)*len(ct).
// this allows all math to be done in integers.
var sd :=
function(trc: TArray<Integer>): Integer
begin
var sumTr := 0;
for var x in trc do
inc(sumTr, all[x]);
result := sumTr * length(ct) - (sumAll - sumTr) * length(tr);
end;
// compute observed difference, as an intermediate result
var a: TArray<Integer>;
SetLength(a, length(tr));
for var i := 0 to High(a) do
a[i] := i;
var sdObs := sd(a);
// iterate over all combinations. for each, compute (scaled)
// difference and tally whether leq or gt observed difference.
var nLe, nGt: Integer;
comb(length(all), length(tr),
procedure(c: TArray<Integer>)
begin
if sd(c) > sdObs then
inc(nGt)
else
inc(nle);
end);
// print results as percentage
var pc := 100 / (nLe + nGt);
writeln(format('differences <= observed: %f%%', [nle * pc]));
writeln(format('differences > observed: %f%%', [ngt * pc]));
{$IFNDEF UNIX} readln; {$ENDIF}
end.
Elixir
defmodule Permutation do
def statistic(ab, a) do
sumab = Enum.sum(ab)
suma = Enum.sum(a)
suma / length(a) - (sumab - suma) / (length(ab) - length(a))
end
def test(a, b) do
ab = a ++ b
tobs = statistic(ab, a)
{under, count} = Enum.reduce(comb(ab, length(a)), {0,0}, fn perm, {under, count} ->
if statistic(ab, perm) <= tobs, do: {under+1, count+1},
else: {under , count+1}
end)
under * 100.0 / count
end
defp comb(_, 0), do: [[]]
defp comb([], _), do: []
defp comb([h|t], m) do
(for l <- comb(t, m-1), do: [h|l]) ++ comb(t, m)
end
end
treatmentGroup = [85, 88, 75, 66, 25, 29, 83, 39, 97]
controlGroup = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98]
under = Permutation.test(treatmentGroup, controlGroup)
:io.fwrite "under = ~.2f%, over = ~.2f%~n", [under, 100-under]
- Output:
under = 87.20%, over = 12.80%
GAP
a := [85, 88, 75, 66, 25, 29, 83, 39, 97];
b := [68, 41, 10, 49, 16, 65, 32, 92, 28, 98];
# Compute a decimal approximation of a rational
Approx := function(x, d)
local neg, a, b, n, m, s;
if x < 0 then
x := -x;
neg := true;
else
neg := false;
fi;
a := NumeratorRat(x);
b := DenominatorRat(x);
n := QuoInt(a, b);
a := RemInt(a, b);
m := 10^d;
s := "";
if neg then
Append(s, "-");
fi;
Append(s, String(n));
n := Size(s) + 1;
Append(s, String(m + QuoInt(a*m, b)));
s[n] := '.';
return s;
end;
PermTest := function(a, b)
local c, d, p, q, u, v, m, n, k, diff, all;
p := Size(a);
q := Size(b);
v := Concatenation(a, b);
n := p + q;
m := Binomial(n, p);
diff := Sum(a)/p - Sum(b)/q;
all := [1 .. n];
k := 0;
for u in Combinations(all, p) do
c := List(u, i -> v[i]);
d := List(Difference(all, u), i -> v[i]);
if Sum(c)/p - Sum(d)/q > diff then
k := k + 1;
fi;
od;
return [Approx((1 - k/m)*100, 3), Approx(k/m*100, 3)];
end;
# in order, % less or greater than original diff
PermTest(a, b);
[ "87.197", "12.802" ]
FreeBASIC
Dim Shared datos(18) As Integer => {85, 88, 75, 66, 25, 29, 83, 39, 97,_
68, 41, 10, 49, 16, 65, 32, 92, 28, 98}
Function pick(at As Integer, remain As Integer, accu As Integer, treat As Integer) As Integer
If remain = 0 Then
If accu > treat Then Return 1 Else Return 0
End If
Dim a As Integer
If at > remain Then a = pick(at-1, remain, accu, treat) Else a = 0
Return pick(at-1, remain-1, accu+datos(at), treat) + a
End Function
Dim As Integer treat = 0, le, gt, total = 1, i
For i = 1 To 9
treat += datos(i)
Next i
For i = 19 To 11 Step -1
total *= i
Next i
For i = 9 To 1 Step -1
total /= i
Next i
gt = pick(19, 9, 0, treat)
le = total - gt
Print Using "<= : ##.######% #####"; 100*le/total; le
Print Using " > : ##.######% #####"; 100*gt/total; gt
Sleep
- Output:
<= : 92.076035% 85058 > : 7.923965% 7320
Go
A version doing all math in integers until computing final percentages.
package main
import "fmt"
var tr = []int{85, 88, 75, 66, 25, 29, 83, 39, 97}
var ct = []int{68, 41, 10, 49, 16, 65, 32, 92, 28, 98}
func main() {
// collect all results in a single list
all := make([]int, len(tr)+len(ct))
copy(all, tr)
copy(all[len(tr):], ct)
// compute sum of all data, useful as intermediate result
var sumAll int
for _, r := range all {
sumAll += r
}
// closure for computing scaled difference.
// compute results scaled by len(tr)*len(ct).
// this allows all math to be done in integers.
sd := func(trc []int) int {
var sumTr int
for _, x := range trc {
sumTr += all[x]
}
return sumTr*len(ct) - (sumAll-sumTr)*len(tr)
}
// compute observed difference, as an intermediate result
a := make([]int, len(tr))
for i, _ := range a {
a[i] = i
}
sdObs := sd(a)
// iterate over all combinations. for each, compute (scaled)
// difference and tally whether leq or gt observed difference.
var nLe, nGt int
comb(len(all), len(tr), func(c []int) {
if sd(c) > sdObs {
nGt++
} else {
nLe++
}
})
// print results as percentage
pc := 100 / float64(nLe+nGt)
fmt.Printf("differences <= observed: %f%%\n", float64(nLe)*pc)
fmt.Printf("differences > observed: %f%%\n", float64(nGt)*pc)
}
// combination generator, copied from combination task
func comb(n, m int, emit func([]int)) {
s := make([]int, m)
last := m - 1
var rc func(int, int)
rc = func(i, next int) {
for j := next; j < n; j++ {
s[i] = j
if i == last {
emit(s)
} else {
rc(i+1, j+1)
}
}
return
}
rc(0, 0)
}
- Output:
differences <= observed: 87.197168% differences > observed: 12.802832%
Haskell
binomial n m = (f !! n) `div` (f !! m) `div` (f !! (n - m))
where f = scanl (*) 1 [1..]
permtest treat ctrl = (fromIntegral less) / (fromIntegral total) * 100
where
total = binomial (length avail) (length treat)
less = combos (sum treat) (length treat) avail
avail = ctrl ++ treat
combos total n a@(x:xs)
| total < 0 = binomial (length a) n
| n == 0 = 0
| n > length a = 0
| n == length a = fromEnum (total < sum a)
| otherwise = combos (total - x) (n - 1) xs
+ combos total n xs
main = let r = permtest
[85, 88, 75, 66, 25, 29, 83, 39, 97]
[68, 41, 10, 49, 16, 65, 32, 92, 28, 98]
in do putStr "> : "; print r
putStr "<=: "; print $ 100 - r
- Output:
> : 12.80283184307952 <=: 87.19716815692048
Somewhat faster, this goes from top down:
binomial n m = (f !! n) `div` (f !! m) `div` (f !! (n - m))
where f = scanl (*) 1 [1..]
perms treat ctrl = (less,total) where
total = binomial (length ctrl + length treat) (length treat)
less = length $ filter (<= sum treat)
$ sums (treat ++ ctrl) (length treat)
sums x n
| l < n || n < 0 = []
| n == 0 = [0]
| l == n = [sum x]
| otherwise = [a + b | i <- [0..n], a <- sums left i, b <- sums right (n - i)]
where (l, l1) = (length x, l `div` 2)
(left, right) = splitAt l1 x
main = print $ (lt, 100 - lt) where
(a, b) = perms [85, 88, 75, 66, 25, 29, 83, 39, 97]
[68, 41, 10, 49, 16, 65, 32, 92, 28, 98]
lt = (fromIntegral a) / (fromIntegral b) * 100
In cases where the sample data are a large number of relatively small positive integers, counting number of partial sums is a lot faster:
combs maxsum len x = foldl f [(0,0,1)] x where
f a n = merge a (map (addNum n) $ filter (\(l,_,_) -> l < len) a)
addNum n (a,s,c)
-- anything larger than maxsum is as good as infinity
| s + n > maxsum = (a+1, maxsum + 1, c)
| otherwise = (a+1, s+n, c)
merge a [] = a
merge [] a = a
merge a@((a1,a2,a3):as) b@((b1,b2,b3):bs)
| a1 == b1 && a2 == b2 = (a1,a2,a3+b3):merge as bs
| a1 < b1 || (a1 == b1 && a2 < b2) = (a1,a2,a3):merge as b
| otherwise = (b1,b2,b3):merge a bs
permtest a b = (lt, ge) where
lt = sum $ map (\(a,b,c) -> if a == la && b < sa then c else 0)
$ combs sa la (a++b)
ge = (binomial (la + lb) la) - lt
(sa, la, lb) = (sum a, length a, length b)
binomial n m = (f !! n) `div` (f !! m) `div` (f !! (n - m))
where f = scanl (*) 1 [1..]
-- how many combinations are less than current sum
main = print$ permtest [85, 88, 75, 66, 25, 29, 83, 39, 97]
[68, 41, 10, 49, 16, 65, 32, 92, 28, 98]
J
require'stats'
trmt=: 0.85 0.88 0.75 0.66 0.25 0.29 0.83 0.39 0.97
ctrl=: 0.68 0.41 0.1 0.49 0.16 0.65 0.32 0.92 0.28 0.98
difm=: -&mean
result=: trmt difm ctrl
all=: trmt(#@[ ({. difm }.) |:@([ (comb ~.@,"1 i.@])&# ,) { ,) ctrl
smoutput 'under: ','%',~":100*mean all <: result
smoutput 'over: ','%',~":100*mean all > result
Result:
under: 87.1972%
over: 12.8028%
Java
public class PermutationTest {
private static final int[] data = new int[]{
85, 88, 75, 66, 25, 29, 83, 39, 97,
68, 41, 10, 49, 16, 65, 32, 92, 28, 98
};
private static int pick(int at, int remain, int accu, int treat) {
if (remain == 0) return (accu > treat) ? 1 : 0;
return pick(at - 1, remain - 1, accu + data[at - 1], treat)
+ ((at > remain) ? pick(at - 1, remain, accu, treat) : 0);
}
public static void main(String[] args) {
int treat = 0;
double total = 1.0;
for (int i = 0; i <= 8; ++i) {
treat += data[i];
}
for (int i = 19; i >= 11; --i) {
total *= i;
}
for (int i = 9; i >= 1; --i) {
total /= i;
}
int gt = pick(19, 9, 0, treat);
int le = (int) (total - gt);
System.out.printf("<= : %f%% %d\n", 100.0 * le / total, le);
System.out.printf(" > : %f%% %d\n", 100.0 * gt / total, gt);
}
}
- Output:
<= : 87.197168% 80551 > : 12.802832% 11827
jq
Part 1: Combinations
# combination(r) generates a stream of combinations of r items from the input array.
def combination(r):
if r > length or r < 0 then empty
elif r == length then .
else ( [.[0]] + (.[1:]|combination(r-1))),
( .[1:]|combination(r))
end;
Part 2: Permutation Test
# a and b should be arrays:
def permutationTest(a; b):
def normalize(a;b): # mainly to avoid having to compute $sumab
(a|add) as $sa
| (b|add) as $sb
| (($sa + $sb)/((a|length) + (b|length))) as $avg
| [(a | map(.-$avg)), (b | map(.-$avg))];
# avg(a) - avg(b) (assuming ab==a+b and avg(ab) is 0)
def statistic(ab; a):
(a | add) as $suma
# (ab|add) should be 0, by normalization
| ($suma / (a|length)) +
($suma / ((ab|length) - (a|length)));
normalize(a;b)
| (a + b) as $ab # pooled observations
| .[0] as $a | .[1] as $b
| statistic($ab; $a) as $t_observed # observed difference in means
| reduce ($ab|combination($a|length)) as $perm # for each combination...
([0,0]; # state: [under,count]
if statistic($ab; $perm) <= $t_observed then .[0] += 1 else . end
| .[1] += 1 )
| .[0] * 100.0 / .[1] # under/count
;
Example:
def treatmentGroup: [85, 88, 75, 66, 25, 29, 83, 39, 97];
def controlGroup: [68, 41, 10, 49, 16, 65, 32, 92, 28, 98];
permutationTest(treatmentGroup; controlGroup) as $under
| "% under=\($under)", "% over=\(100 - $under)"
- Output:
$ jq -n -r -f permutation_test.jq % under=87.14304271579813 % over=12.856957284201869
Julia
The primary function for this solution is permutation_test, which relies on Julia's combinations (from Combinatorics module) function to provide all of the possible study arm assignments. bifurcate splits the pooled results into "treatment" and "control" groups according to the indices provided by combinations.
Functions
using Combinatorics
meandiff(a::Vector{T}, b::Vector{T}) where T <: Real = mean(a) - mean(b)
function bifurcate(a::AbstractVector, sel::Vector{T}) where T <: Integer
x = a[sel]
asel = trues(length(a))
asel[sel] = false
y = a[asel]
return x, y
end
function permutation_test(treated::Vector{T}, control::Vector{T}) where T <: Real
effect0 = meandiff(treated, control)
pool = vcat(treated, control)
tlen = length(treated)
plen = length(pool)
better = worse = 0
for subset in combinations(1:plen, tlen)
t, c = bifurcate(pool, subset)
if effect0 < meandiff(t, c)
better += 1
else
worse += 1
end
end
return better, worse
end
Main
const treated = [85, 88, 75, 66, 25, 29, 83, 39, 97]
const control = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98]
(better, worse) = permutation_test(treated, control)
tot = better + worse
println("Permutation test using the following data:")
println("Treated: ", treated)
println("Control: ", control)
println("\nThere are $tot different permuted groups of these data.")
@printf("%8d, %5.2f%% showed better than actual results.\n", better, 100 * better / tot)
print(@sprintf("%8d, %5.2f%% showed equalivalent or worse results.", worse, 100 * worse / tot))
- Output:
Permutation test using the following data: Treated: [85, 88, 75, 66, 25, 29, 83, 39, 97] Control: [68, 41, 10, 49, 16, 65, 32, 92, 28, 98] There are 92378 different permuted groups of these data. 11827, 12.80% showed better than actual results. 80551, 87.20% showed equalivalent or worse results.
Kotlin
// version 1.1.2
val data = intArrayOf(
85, 88, 75, 66, 25, 29, 83, 39, 97,
68, 41, 10, 49, 16, 65, 32, 92, 28, 98
)
fun pick(at: Int, remain: Int, accu: Int, treat: Int): Int {
if (remain == 0) return if (accu > treat) 1 else 0
return pick(at - 1, remain - 1, accu + data[at - 1], treat) +
if (at > remain) pick(at - 1, remain, accu, treat) else 0
}
fun main(args: Array<String>) {
var treat = 0
var total = 1.0
for (i in 0..8) treat += data[i]
for (i in 19 downTo 11) total *= i
for (i in 9 downTo 1) total /= i
val gt = pick(19, 9, 0, treat)
val le = (total - gt).toInt()
System.out.printf("<= : %f%% %d\n", 100.0 * le / total, le)
System.out.printf(" > : %f%% %d\n", 100.0 * gt / total, gt)
}
- Output:
<= : 87.197168% 80551 > : 12.802832% 11827
M2000 Interpreter
Module Checkit {
Global data(), treat=0
data()=(85, 88, 75, 66, 25, 29, 83, 39, 97,68, 41, 10, 49, 16, 65, 32, 92, 28, 98)
Function pick(at, remain, accu) {
If remain Else =If(accu>treat->1,0):Exit
=pick(at-1,remain-1,accu+data(at-1))+If(at>remain->pick(at-1, remain, accu),0)
}
total=1
For i=0 to 8 {treat+=data(i)}
For i=19 to 11 {total*=i}
For i=9 to 1 {total/=i}
gt=pick(19,9,0)
le=total-gt
Print Format$("<= : {0:1}% {1}", 100*le/total, le)
Print Format$(" > : {0:1}% {1}", 100*gt/total, gt)
}
Checkit
- Output:
<= : 87.2% 80551 > : 12.8% 11827
Slower version, using a lambda function with a series of inner lambda functions to return each combination at a time.
Module CheckThis {
Function CombinationsStep (a, nn) {
c1=lambda (&f, &a) ->{=car(a) : a=cdr(a) : f=len(a)=0}
m=len(a)
c=c1
n=m-nn+1
p=2
while m>n {
c1=lambda c2=c,n=p, z=(,) (&f, &m) ->{if len(z)=0 then z=cdr(m)
=cons(car(m),c2(&f, &z)):if f then z=(,) : m=cdr(m) : f=len(m)+len(z)<n
}
c=c1
p++
m--
}
=lambda c, a (&f) ->c(&f, &a)
}
treated=(85, 88, 75, 66, 25, 29, 83, 39, 97)
placebo=(68, 41, 10, 49, 16, 65, 32, 92, 28, 98)
treat=0
m=each(treated): while m {treat+=array(m)}
total=1
for i=len(placebo)+1 to len(placebo) +len(treated):total*=i:next i
for i=len(placebo)-1 to 1: total/=i:next i
d=total div 10**int(log(total))
k=false
StepA=CombinationsStep(cons(treated, placebo),len(treated))
counter=0
gt=0
While not k {
comb=StepA(&k)
accu=0
m=each(comb)
while m {accu+=array(m)}
gt+=if(accu>treat->1,0)
counter++
if counter mod d=0 then Print over str$(counter/total," #0.0%"): Refresh 1000
}
print over str$(counter/total," #0.0%")
print
lt=total-gt
print Format$("less or equal={0:1}%, greater={1:1}%, total={2}",lt/total*100, gt/total*100, total)
}
CheckThis
Mathematica/Wolfram Language
"<=: " <> ToString[#1] <> " " <> ToString[100. #1/#2] <> "%\n >: " <>
ToString[#2 - #1] <> " " <> ToString[100. (1 - #1/#2)] <> "%" &[
Count[Total /@ Subsets[Join[#1, #2], {Length@#1}],
n_ /; n <= Total@#1],
Binomial[Length@#1 + Length@#2, Length@#1]] &[{85, 88, 75, 66, 25,
29, 83, 39, 97}, {68, 41, 10, 49, 16, 65, 32, 92, 28, 98}]
- Output:
<=: 80551 87.1972% >: 11827 12.8028%
Nim
import strformat
const data = [85, 88, 75, 66, 25,
29, 83, 39, 97, 68,
41, 10, 49, 16, 65,
32, 92, 28, 98]
func pick(at, remain, accu, treat: int): int =
if remain == 0:
return if accu > treat: 1 else: 0
return pick(at - 1, remain - 1, accu + data[at - 1], treat) +
(if at > remain: pick(at - 1, remain, accu, treat) else: 0)
var treat = 0
var le, gt = 0
var total = 1.0
for i in countup(0, 8):
treat += data[i]
for i in countdown(19, 11):
total *= float(i)
for i in countdown(9, 1):
total /= float(i)
gt = pick(19, 9, 0, treat)
le = int(total - float(gt))
echo fmt"<= : {100.0 * float(le) / total:.6f}% {le}"
echo fmt" > : {100.0 * float(gt) / total:.6f}% {gt}"
- Output:
<= : 87.197168% 80551 > : 12.802832% 11827
Perl
#!/usr/bin/perl
use warnings;
use strict;
use List::Util qw{ sum };
sub means {
my @groups = @_;
return map sum(@$_) / @$_, @groups;
}
sub following {
my $pattern = shift;
my $orig_count = grep $_, @$pattern;
my $count;
do {
my $i = $#{$pattern};
until (0 > $i) {
$pattern->[$i] = $pattern->[$i] ? 0 : 1;
last if $pattern->[$i];
--$i;
}
$count = grep $_, @$pattern;
} until $count == $orig_count or not $count;
undef @$pattern unless $count;
}
my @groups;
my $i = 0;
while (<DATA>) {
chomp;
$i++, next if /^$/;
push @{ $groups[$i] }, $_;
}
my @orig_means = means(@groups);
my $orig_cmp = $orig_means[0] - $orig_means[1];
my $pattern = [ (0) x @{ $groups[0] },
(1) x @{ $groups[1] }
];
my @cmp = (0) x 3;
while (@$pattern) {
my @perms = map { my $g = $_;
[ (@{ $groups[0] }, @{ $groups[1] } ) [ grep $pattern->[$_] == $g, 0 .. $#{$pattern} ] ];
} 0, 1;
my @means = means(@perms);
$cmp[ ($means[0] - $means[1]) <=> $orig_cmp ]++;
} continue {
following($pattern);
}
my $all = sum(@cmp);
my $length = length $all;
for (0, -1, 1) {
printf "%-7s %${length}d %6.3f%%\n",
(qw(equal greater less))[$_], $cmp[$_], 100 * $cmp[$_] / $all;
}
__DATA__
85
88
75
66
25
29
83
39
97
68
41
10
49
16
65
32
92
28
98
- Output:
equal 313 0.339% less 80238 86.858% greater 11827 12.803%
Phix
constant data = {85, 88, 75, 66, 25, 29, 83, 39, 97, 68, 41, 10, 49, 16, 65, 32, 92, 28, 98 } function pick(int at, int remain, int accu, int treat) if remain=0 then return iff(accu>treat?1:0) end if return pick(at-1, remain-1, accu+data[at], treat) + iff(at>remain?pick(at-1, remain, accu, treat):0) end function int treat = 0, le, gt atom total = 1; for i=1 to 9 do treat += data[i] end for for i=19 to 11 by -1 do total *= i end for for i=9 to 1 by -1 do total /= i end for gt = pick(19, 9, 0, treat) le = total - gt; printf(1,"<= : %f%% %d\n > : %f%% %d\n", {100*le/total, le, 100*gt/total, gt})
- Output:
<= : 87.197168% 80551 > : 12.802832% 11827
PicoLisp
(load "@lib/simul.l") # For 'subsets'
(scl 2)
(de _stat (A)
(let (LenA (length A) SumA (apply + A))
(-
(*/ SumA LenA)
(*/ (- SumAB SumA) (- LenAB LenA)) ) ) )
(de permutationTest (A B)
(let
(AB (append A B)
SumAB (apply + AB)
LenAB (length AB)
Tobs (_stat A)
Count 0 )
(*/
(sum
'((Perm)
(inc 'Count)
(and (>= Tobs (_stat Perm)) 1) )
(subsets (length A) AB) )
100.0
Count ) ) )
(setq
*TreatmentGroup (0.85 0.88 0.75 0.66 0.25 0.29 0.83 0.39 0.97)
*ControlGroup (0.68 0.41 0.10 0.49 0.16 0.65 0.32 0.92 0.28 0.98) )
(let N (permutationTest *TreatmentGroup *ControlGroup)
(prinl "under = " (round N) "%, over = " (round (- 100.0 N)) "%") )
- Output:
under = 87.85%, over = 12.15%
PureBasic
Given a treatment group with [n=9] and a control group with [m=10]. The numbers [x] from [1] to [1<<(n+m)] exhaust the possible states.
Any bit-String of Length [n+m] containing [n=9] "1's" is a Valid bit String, as tested by: IsValidBitString(x,n+m,n).
Then we can use these bits to Select whether a particular index For our array should be assigned to: the treatment group or the control group
Define.f meanTreated,meanControl,diffInMeans
Define.f actualmeanTreated,actualmeanControl,actualdiffInMeans
Dim poolA(19)
poolA(1) =85 ; first 9 the treated
poolA(2) =88
poolA(3) =75
poolA(4) =66
poolA(5) =25
poolA(6) =29
poolA(7) =83
poolA(8) =39
poolA(9) =97
poolA(10) =68 ; last 10 the control
poolA(11) =41
poolA(12) =10
poolA(13) =49
poolA(14) =16
poolA(15) =65
poolA(16) =32
poolA(17) =92
poolA(18) =28
poolA(19) =98
Procedure.i IsValidBitString(x,pool,treated)
Protected c,i
For i=1 to pool
mask=1<<(i-1)
If mask&x:c+1:EndIf
Next
If c=treated :ProcedureReturn x
Else :ProcedureReturn 0
EndIf
EndProcedure
treated=9
control=10
pool =treated+control
; actual Experimentally observed difference in means
For i=1 to Treated
sumTreated+poolA(i)
Next
For i=Treated+1 to Treated+Control
sumControl+poolA(i)
Next
actualmeanTreated=sumTreated /Treated
actualmeanControl=sumControl /Control
actualdiffInMeans=actualmeanTreated-actualmeanControl
; exhaust the possibilites
For x=1 to 1<<pool
; Valid? i.e. are there 9 "1's" ?
If IsValidBitString(x,pool,treated)
TotalComBinations+1:sumTreated=0:sumControl=0
; separate the groups
For i=pool to 1 Step -1
mask=1<<(i-1):idx=pool-i+1
If mask&x
sumTreated+poolA(idx)
Else
sumControl+poolA(idx)
EndIf
Next
meanTreated=sumTreated /Treated
meanControl=sumControl /Control
diffInMeans=meanTreated-meanControl
; gather the statistics
If (diffInMeans)<=(actualdiffInMeans)
diffLessOrEqual+1
Else
diffGreater+1
EndIf
EndIf
Next
; show our results
; cw(StrF(100*diffLessOrEqual/TotalComBinations,2)+" "+Str(diffLessOrEqual))
; cw(StrF(100*diffGreater /TotalComBinations,2)+" "+Str(diffGreater))
Debug StrF(100*diffLessOrEqual/TotalComBinations,2)+" "+Str(diffLessOrEqual)
Debug StrF(100*diffGreater /TotalComBinations,2)+" "+Str(diffGreater)
- Output:
87.20 80551 12.80 11827
Python
from itertools import combinations as comb
def statistic(ab, a):
sumab, suma = sum(ab), sum(a)
return ( suma / len(a) -
(sumab -suma) / (len(ab) - len(a)) )
def permutationTest(a, b):
ab = a + b
Tobs = statistic(ab, a)
under = 0
for count, perm in enumerate(comb(ab, len(a)), 1):
if statistic(ab, perm) <= Tobs:
under += 1
return under * 100. / count
treatmentGroup = [85, 88, 75, 66, 25, 29, 83, 39, 97]
controlGroup = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98]
under = permutationTest(treatmentGroup, controlGroup)
print("under=%.2f%%, over=%.2f%%" % (under, 100. - under))
- Output:
under=89.11%, over=10.89%
The above solution does a different thing than the other solutions. I'm not really sure why. If you want to do the same thing as the other solutions, then this is the solution:
from itertools import combinations as comb
def permutationTest(a, b):
ab = a + b
Tobs = sum(a)
under = 0
for count, perm in enumerate(comb(ab, len(a)), 1):
if sum(perm) <= Tobs:
under += 1
return under * 100. / count
treatmentGroup = [85, 88, 75, 66, 25, 29, 83, 39, 97]
controlGroup = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98]
under = permutationTest(treatmentGroup, controlGroup)
print("under=%.2f%%, over=%.2f%%" % (under, 100. - under))
- Output:
under=87.20%, over=12.80%
R
permutation.test <- function(treatment, control) {
perms <- combinations(length(treatment)+length(control),
length(treatment),
c(treatment, control),
set=FALSE)
p <- mean(rowMeans(perms) <= mean(treatment))
c(under=p, over=(1-p))
}
> permutation.test(c(85, 88, 75, 66, 25, 29, 83, 39, 97),
+ c(68, 41, 10, 49, 16, 65, 32, 92, 28, 98))
under over
0.8719717 0.1280283
Racket
#lang racket/base
(define-syntax-rule (inc! x)
(set! x (add1 x)))
(define (permutation-test control-gr treatment-gr)
(let ([both-gr (append control-gr treatment-gr)]
[threshold (apply + control-gr)]
[more 0]
[leq 0])
(let loop ([data both-gr] [sum 0] [needed (length control-gr)] [available (length both-gr)])
(cond [(zero? needed) (if (>= sum threshold)
(inc! more)
(inc! leq))]
[(>= available needed) (loop (cdr data) sum needed (sub1 available))
(loop (cdr data) (+ sum (car data)) (sub1 needed) (sub1 available))]
[else (void)]))
(values more leq)))
(let-values ([(more leq) (permutation-test '(68 41 10 49 16 65 32 92 28 98)
'(85 88 75 66 25 29 83 39 97))])
(let ([sum (+ more leq)])
(printf "<=: ~a ~a%~n>: ~a ~a%~n"
more (real->decimal-string (* 100. (/ more sum)) 2)
leq (real->decimal-string (* 100. (/ leq sum)) 2))))
- Output:
<=: 80551 87.20% >: 11827 12.80%
Raku
(formerly Perl 6)
The use of .race
to allow concurrent calculations means that multiple 'workers' will be updating @trials
simultaneously. To avoid race conditions, the ⚛++
operator is used, which guarantees safe updates without the use of locks. That is turn requires declaring that array as being composed of atomicint
.
sub stats ( @test, @all ) {
([+] @test / +@test) - ([+] flat @all, (@test X* -1)) / @all - @test
}
my int @treated = <85 88 75 66 25 29 83 39 97>;
my int @control = <68 41 10 49 16 65 32 92 28 98>;
my int @all = flat @treated, @control;
my $base = stats( @treated, @all );
my atomicint @trials[3] = 0, 0, 0;
@all.combinations(+@treated).race.map: { @trials[ 1 + ( stats( $_, @all ) <=> $base ) ]⚛++ }
say 'Counts: <, =, > ', @trials;
say 'Less than : %', 100 * @trials[0] / [+] @trials;
say 'Equal to : %', 100 * @trials[1] / [+] @trials;
say 'Greater than : %', 100 * @trials[2] / [+] @trials;
say 'Less or Equal: %', 100 * ( [+] @trials[0,1] ) / [+] @trials;
- Output:
Counts: <, =, > 80238 313 11827 Less than : %86.858343 Equal to : %0.338825 Greater than : %12.802832 Less or Equal: %87.197168
REXX
This REXX program is modeled after the C version, with some generalizations and optimization added.
/*REXX program performs a permutation test on N + M subjects (volunteers): */
/* ↑ ↑ */
/* │ │ */
/* │ └─────control population. */
/* └────────treatment population. */
n= 9 /*define the number of the control pop.*/
data= 85 88 75 66 25 29 83 39 97 68 41 10 49 16 65 32 92 28 98
w= words(data); m= w - n /*w: volunteers + control population*/
L= length(w) /*L: used to align some output numbers*/
say '# of volunteers & control population: ' w
say 'volunteer population given treatment: ' right(n, L)
say ' control population given a placebo: ' right(m, L)
say
say 'treatment population efficacy % (percentages): ' subword(data, 1, n)
say ' control population placebo % (percentages): ' subword(data, n+1 )
say
do v= 0 for w ; #.v= word(data, v+1) ; end
treat= 0; do i= 0 to n-1 ; treat= treat + #.i ; end
tot= 1; do j= w to m+1 by -1 ; tot= tot * j ; end
do k=w%2 to 1 by -1 ; tot= tot / k ; end
GT= picker(n+m, n, 0) /*compute the GT value from PICKER func*/
LE= tot - GT /* " " LE " via subtraction.*/
say "<= " format(100 * LE / tot, ,3)'%' LE /*display number with 3 decimal places.*/
say " > " format(100 * GT / tot, ,3)'%' GT /* " " " " " " */
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
picker: procedure expose #. treat; parse arg it,rest,eff /*get the arguments.*/
if rest==0 then return eff > treat /*is REST = to zero?*/
if it>rest then q= picker(it-1, rest, eff) /*maybe recurse. */
else q= 0
itP= it - 1 /*set temporary var.*/
return picker(itP, rest - 1, eff + #.itP) + q /*recurse. */
- output when using the default input:
# of volunteers & control population: 19 volunteer population given treatment: 9 control population given a placebo: 10 treatment population efficacy % (percentages): 85 88 75 66 25 29 83 39 97 control population placebo % (percentages): 68 41 10 49 16 65 32 92 28 98 <= 87.197% 80551 > 12.803% 11827
Ruby
def statistic(ab, a)
sumab, suma = ab.inject(:+).to_f, a.inject(:+).to_f
suma / a.size - (sumab - suma) / (ab.size - a.size)
end
def permutationTest(a, b)
ab = a + b
tobs = statistic(ab, a)
under = count = 0
ab.combination(a.size) do |perm|
under += 1 if statistic(ab, perm) <= tobs
count += 1
end
under * 100.0 / count
end
treatmentGroup = [85, 88, 75, 66, 25, 29, 83, 39, 97]
controlGroup = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98]
under = permutationTest(treatmentGroup, controlGroup)
puts "under=%.2f%%, over=%.2f%%" % [under, 100 - under]
- Output:
under=87.20%, over=12.80%
Rust
fn main() {
let treatment = vec![85, 88, 75, 66, 25, 29, 83, 39, 97];
let control = vec![68, 41, 10, 49, 16, 65, 32, 92, 28, 98];
let mut data_set = control.clone();
data_set.extend_from_slice(&treatment);
let greater = combinations(treatment.iter().sum(), treatment.len() as i64, &data_set) as f64;
let lesser = combinations(control.iter().sum(), control.len() as i64, &data_set) as f64;
let total = binomial(data_set.len() as i64, treatment.len() as i64) as f64;
println!("<= : {}%", (lesser / total * 100.0));
println!(" > : {}%", (greater / total * 100.0));
}
fn factorial(x: i64) -> i64 {
let mut product = 1;
for a in 1..(x + 1) {
product *= a;
}
product
}
fn binomial(n: i64, k: i64) -> i64 {
let numerator = factorial(n);
let denominator = factorial(k) * factorial(n - k);
numerator / denominator
}
fn combinations(total: i64, number: i64, data: &[i64]) -> i64 {
if total < 0 {
return binomial(data.len() as i64, number);
}
if number == 0 {
return 0;
}
if number > data.len() as i64 {
return 0;
}
if number == data.len() as i64 {
if total < data.iter().sum() {
return 1;
} else {
return 0;
}
}
let tail = &data[1..];
combinations(total - data[0], number - 1, &tail) + combinations(total, number, &tail)
}
- Output:
<= : 86.8583428954946% > : 12.80283184307952%
Scala
Imperative version (Ugly, side effects)
- Output:
Best seen running in your browser either by ScalaFiddle (ES aka JavaScript, non JVM) or Scastie (remote JVM).
object PermutationTest extends App {
private val data =
Array(85, 88, 75, 66, 25, 29, 83, 39, 97, 68, 41, 10, 49, 16, 65, 32, 92, 28, 98)
private var (total, treat) = (1.0, 0)
private def pick(at: Int, remain: Int, accu: Int, treat: Int): Int = {
if (remain == 0) return if (accu > treat) 1 else 0
pick(at - 1, remain - 1, accu + data(at - 1), treat) +
(if (at > remain) pick(at - 1, remain, accu, treat) else 0)
}
for (i <- 0 to 8) treat += data(i)
for (j <- 19 to 11 by -1) total *= j
for (g <- 9 to 1 by -1) total /= g
private val gt = pick(19, 9, 0, treat)
private val le = (total - gt).toInt
println(f"<= : ${100.0 * le / total}%f%% ${le}%d")
println(f" > : ${100.0 * gt / total}%f%% ${gt}%d")
}
Seed7
$ include "seed7_05.s7i";
include "float.s7i";
const array integer: treatmentGroup is [] (85, 88, 75, 66, 25, 29, 83, 39, 97);
const array integer: controlGroup is [] (68, 41, 10, 49, 16, 65, 32, 92, 28, 98);
const array integer: both is treatmentGroup & controlGroup;
const func integer: pick (in integer: at, in integer: remain, in integer: accu, in integer: treat) is func
result
var integer: picked is 0;
begin
if remain = 0 then
picked := ord(accu > treat);
else
picked := pick(at - 1, remain - 1, accu + both[at], treat);
if at > remain then
picked +:= pick(at - 1, remain, accu, treat);
end if;
end if;
end func;
const proc: main is func
local
var integer: experimentalResult is 0;
var integer: treat is 0;
var integer: total is 1;
var integer: le is 0;
var integer: gt is 0;
var integer: i is 0;
begin
for experimentalResult range treatmentGroup do
treat +:= experimentalResult;
end for;
total := 19 ! 10; # Binomial coefficient
gt := pick(19, 9, 0, treat);
le := total - gt;
writeln("<= : " <& 100.0 * flt(le) / flt(total) digits 6 <& "% " <& le);
writeln(" > : " <& 100.0 * flt(gt) / flt(total) digits 6 <& "% " <& gt);
end func;
- Output:
<= : 87.197168% 80551 > : 12.802832% 11827
Sidef
func statistic(ab, a) {
var(sumab, suma) = (ab.sum, a.sum)
suma/a.size - ((sumab-suma) / (ab.size-a.size))
}
func permutationTest(a, b) {
var ab = (a + b)
var tobs = statistic(ab, a)
var under = (var count = 0)
ab.combinations(a.len, {|*perm|
statistic(ab, perm) <= tobs && (under += 1)
count += 1
})
under * 100 / count
}
var treatmentGroup = [85, 88, 75, 66, 25, 29, 83, 39, 97]
var controlGroup = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98]
var under = permutationTest(treatmentGroup, controlGroup)
say ("under=%.2f%%, over=%.2f%%" % (under, 100 - under))
- Output:
under=87.20%, over=12.80%
Tcl
package require Tcl 8.5
# Difference of means; note that the first list must be the concatenation of
# the two lists (because this is cheaper to work with).
proc statistic {AB A} {
set sumAB [tcl::mathop::+ {*}$AB]
set sumA [tcl::mathop::+ {*}$A]
expr {
$sumA / double([llength $A]) -
($sumAB - $sumA) / double([llength $AB] - [llength $A])
}
}
# Selects all k-sized combinations from a list.
proc selectCombinationsFrom {k l} {
if {$k == 0} {return {}} elseif {$k == [llength $l]} {return [list $l]}
set all {}
set n [expr {[llength $l] - [incr k -1]}]
for {set i 0} {$i < $n} {} {
set first [lindex $l $i]
incr i
if {$k == 0} {
lappend all $first
} else {
foreach s [selectCombinationsFrom $k [lrange $l $i end]] {
lappend all [list $first {*}$s]
}
}
}
return $all
}
# Compute the permutation test value and its complement.
proc permutationTest {A B} {
set whole [concat $A $B]
set Tobs [statistic $whole $A]
set undercount 0
set overcount 0
set count 0
foreach perm [selectCombinationsFrom [llength $A] $whole] {
set t [statistic $whole $perm]
incr count
if {$t <= $Tobs} {incr undercount} else {incr overcount}
}
set count [tcl::mathfunc::double $count]
list [expr {$overcount / $count}] [expr {$undercount / $count}]
}
Demonstration code:
set treatmentGroup {0.85 0.88 0.75 0.66 0.25 0.29 0.83 0.39 0.97}
set controlGroup {0.68 0.41 0.10 0.49 0.16 0.65 0.32 0.92 0.28 0.98}
lassign [permutationTest $treatmentGroup $controlGroup] over under
puts [format "under=%.2f%%, over=%.2f%%" [expr {$under*100}] [expr {$over*100}]]
- Output:
under=86.90%, over=13.10%
Ursala
#import std
#import nat
#import flo
treatment_group = <85,88,75,66,25,29,83,39,97>
control_group = <68,41,10,49,16,65,32,92,28,98>
f = # returns the fractions of alternative mean differences above and below the actual
float~*; -+
vid^~G(plus,~&)+ (not fleq@rlX)*|@htX; ~~ float+ length,
minus*+ mean^~*C/~& ^DrlrjXS(~&l,choices)^/-- length@l+-
#show+
t = --* *-'%'@lrNCC printf/$'%0.2f' times/$100. f(treatment_group,control_group)
- Output:
12.80% 87.20%
Wren
import "./fmt" for Fmt
var data = [85, 88, 75, 66, 25, 29, 83, 39, 97, 68, 41, 10, 49, 16, 65, 32, 92, 28, 98]
var pick // recursive
pick = Fn.new { |at, remain, accu, treat|
if (remain == 0) return (accu > treat) ? 1 : 0
return pick.call(at-1, remain-1, accu + data[at-1], treat) +
((at > remain) ? pick.call(at-1, remain, accu, treat) : 0)
}
var treat = 0
var total = 1
for (i in 0..8) treat = treat + data[i]
for (i in 19..11) total = total * i
for (i in 9..1) total = total / i
var gt = pick.call(19, 9, 0, treat)
var le = (total - gt).truncate
Fmt.print("<= : $f\% $d", 100 * le / total, le)
Fmt.print(" > : $f\% $d", 100 * gt / total, gt)
- Output:
<= : 87.197168% 80551 > : 12.802832% 11827
XPL0
include xpllib; \For Print
int Data;
func Pick(At, Remain, Accu, Treat);
int At, Remain, Accu, Treat;
[if Remain = 0 then return if Accu > Treat then 1 else 0;
return Pick(At-1, Remain-1, Accu + Data(At-1), Treat) +
(if At > Remain then Pick(At-1, Remain, Accu, Treat) else 0);
];
int Treat, GT, LE, I;
real Total;
[Treat:= 0; Total:= 1.;
Data:= [85, 88, 75, 66, 25, 29, 83, 39, 97,
68, 41, 10, 49, 16, 65, 32, 92, 28, 98];
for I:= 0 to 8 do Treat:= Treat + Data(I);
for I:= 19 downto 11 do Total:= Total * float(I);
for I:= 9 downto 1 do Total:= Total / float(I);
GT:= Pick(19, 9, 0, Treat);
LE:= (fix(Total) - GT);
Print("<= : %f% %d\n", 100. * float(LE) / Total, LE);
Print(" > : %f% %d\n", 100. * float(GT) / Total, GT);
]
- Output:
<= : 87.19717% 80551 > : 12.80283% 11827
zkl
A solution that is not going to scale gracefully at all.
fcn permutationTest(a,b){
ab := a.extend(b);
tObs := a.sum(0);
combs := Utils.Helpers.pickNFrom(a.len(),ab); // 92,378
under := combs.reduce('wrap(sum,perm){ sum+(perm.sum(0) <= tObs) },0);
100.0 * under / combs.len();
}
treatmentGroup := T(85, 88, 75, 66, 25, 29, 83, 39, 97);
controlGroup := T(68, 41, 10, 49, 16, 65, 32, 92, 28, 98);
under := permutationTest(treatmentGroup, controlGroup);
println("Under =%6.2f%%\nOver =%6.2f%%".fmt(under, 100.0 - under));
- Output:
Under = 87.20% Over = 12.80%
- Clarified and Needing Review
- Programming Tasks
- Probability and statistics
- 11l
- Ada
- AWK
- BBC BASIC
- C
- C sharp
- C++
- Common Lisp
- D
- Delphi
- System.SysUtils
- Elixir
- GAP
- FreeBASIC
- Go
- Haskell
- J
- Java
- Jq
- Julia
- Kotlin
- M2000 Interpreter
- Mathematica
- Wolfram Language
- Nim
- Perl
- Phix
- PicoLisp
- PureBasic
- Python
- R
- Racket
- Raku
- REXX
- Ruby
- Rust
- Scala
- Seed7
- Sidef
- Tcl
- Ursala
- Wren
- Wren-fmt
- XPL0
- Zkl