Comment detail

重複無し乱数 (Nested Flatten)
標準のrandpermという関数が問のbingoと全く同じ動作をします。投稿のコードはsortがソート前の順序を返すことを使った別解です。
1
2
function r = bingo(n)
[m r] = sort(rand(1,n));

標準で必要な関数があるのにあえて別解を用いるメリットが最初の例では特に無かったので、少し書き換えてk個の列を一度に返せるようにした。標準のrandpermでk個生成するにはループでk回呼ばないといけない。余談だがMATLAB7.4ではデフォルトの乱数生成アルゴリズムとして#2269さんのリンク先でも紹介されているメルセンヌ・ツイスタが採用されている。

1
2
3
4
5
6
7
8
9
function r = bingo(n,k)
% r = bingo(n) retruns a random permutation of 1:n.
% r = bingo(n,k) returns a k-by-n array, each row of which is a random
% permutation of 1:n.
% (ja.doukaku.org Q46)
if nargin<2
    k = 1;
end
[m r] = sort(rand([k n]),2);

最初の投稿でえらそうに別解です、なんて書いたけどMATLABのrandpermは実は.mファイル(つまりMATLABで書かれたスクリプト)でその中身は投稿したコードと全く同じだということに気付いた。情けないのであらためて別解。前と同様第二引数で何個生成するか指定する。sortを使う方法よりずいぶん遅い。次のようにすると計算時間の測定と順序が均等かどうかのおおまかな確認ができる。

tic;r=bingo(100,10000);bar(sum(r,1));toc
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
function r = bingo(n,k)
% r = bingo(n) retruns a random permutation of 1:n.
% r = bingo(n,k) returns a k-by-n array, each row of which is a random
% permutation of 1:n.
% (ja.doukaku.org Q46)
if nargin<2
    k = 1;
end
[s r] = meshgrid(1:k,1:n);
% q = zeros(k,2);
for p = n:-1:2
    q(:,1) = ceil(rand(k,1)*p);
    q(:,2) = p;
    q = q+s(1:2,:)'*n-n;
    r(q) = r(q(:,[2 1]));
end
r = r';

#2264の指摘をみてなるほどと思ったので簡単なテスト。順位相関を使ってみた。

#2285のコード(MATLAB標準のrandpermも同じ方法)は最初に乱数を要求された数列の個数分だけ作り、その順序を返す。MATLABのsortは値が同じ場合は順序を保存するので、たとえば1番目と2番目の乱数が同じ数字だと必ず1が2の前になってしまう。

投稿のコードは乱数の有効桁数(最初の行の最後の引数)を制限し、生成した数列とその位置の順位相関の分布と平均を表示する。有効桁数を1にした場合はかなり偏りが生じるが、3桁くらいでほとんど偏りが無い。順位相関を見る限りでは、普通の精度の乱数を使えばこのアルゴリズムでの偏りを気にする必要は無さそう。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
x = bingo(25,1000,1);
rho = corr((1:size(x,2))',x(1:size(x,1),:)','type','Kendall');
figure(1); hist(rho,-1:0.1:1); axis([-1 1 0 300])
fprintf('Mean = %f, skewness = %f\n', mean(rho), skewness(rho));

function r = bingo(n,k,p)
% r = bingo(n,k) returns a k-by-n array, each row of which is a random
% permutation of 1:n, using random numbers with p significant digits .
[m r] = sort(randr([k n],p),2);

function r = randr(M,p)
r = floor(rand(M).*10^p)./10^p;

Index

Feed

Other

Link

Pathtraq

loading...