人工的なデータを作成するアルゴリズム

質問

人工的なデータを作成するアルゴリズムを探しています。
条件は，
1. データの個数はN=500個を3セット
2. それぞれのセットは1から500までの整数を並べ替えたもの
3. セット間の相関係数は，セット1とセット2が0.7，セット1とセット3が0.6，セット2とセット3が0.5とします。
このうち条件1と条件2は絶対です。条件3はできるだけ近ければよいとします。
どのようにして作成すれば良いでしょうか？

回答No.1

ask-it-aurora

探索すべき対象が非常に多く、期待している相関係数は比較的大きいです。そこでリスト [1, 2, ..., 500] をふたつに分けひっくり返した形のものだけを探してみました。
セット1 = x = [1, 2, ..., 500],
セット2 = y = [22, 23, ..., 500, 1, 2, ..., 21],
セット3 = z = [470, 471, ..., 500, 1, 2, ..., 469]
を考えると相関係数はおよそ
r(x, y) ~ 0.76,
r(x, z) ~ 0.65,
r(y, z) ~ 0.44
となります。これくらいではいかが？

以下、計算に使ったpythonのコード（どうもインデントが正しく表示できないようですが、理解する上で特に問題はないでしょう）。計算はラップトップで一分ほどでした。
from scipy.stats import pearsonr

N = 500
a = range(1, N+1)

def r(x, y): return pearsonr(x, y)[0]

def rev(n): return a[n:] + a[:n]

## find date set
score = float("inf")
x = rev(0)
for j in range(N):
y = rev(j)
score1 = (r(x, y) - 0.7)**2
for k in range(N):
z = rev(k)
score2 = (r(x, z) - 0.6)**2 + (r(y, z) - 0.5)**2
if score1 + score2 < score:
score = score1 + score2
j0 = j; k0 = k

## result
x = rev(0); y = rev(j0); z = rev(k0)
print j0, k0
print r(x, y), r(x, z), r(y, z)

お礼

ありがとうございます。
結局のところ，以下のようにしました。9セット分のデータを作成しても0.3秒ですので，十分に実用に耐えます。
簡単にやり方を説明すると，
1.一様乱数を必要なセット数だけ発生させる。
2.各データセットを与えられた相関があるように変換する。
3.変換後のデータセットをそれぞれ順位づけして1から500までの整数に変換する。
相関係数は
1.0
0.6，1.0
0.2，0.1，1.0
0.1，0.2，0.7，1.0
0.1，0.2，0.6，0.6，1.0
0.1，0.2，0.1，0.1，0.1，1.0
0.1，0.2，0.1，0.2，0.2，0.6，1.0
0.1，0.1，0.1，0.1，0.1，0.1，0.1，1.0
0.1，0.1，0.1，0.1，0.1，0.1，0.1，0.6，1.0
を目標にしていますが，実際に得られた相関係数は以下のとおりです。
1.000
0.585，1.000
0.231，0.097，1.000
0.191，0.219，0.726，1.000
0.159，0.217，0.570，0.591，1.000
0.107，0.170，0.075，0.031，0.082，1.000
0.059，0.178，0.040，0.112，0.152，0.565，1.000
0.174，0.178，0.168，0.166，0.176，0.127，0.120，1.000
0.153，0.203，0.152，0.093，0.133，0.124，0.100，0.591，1.000

#include<iostream>
#include<iomanip>
#include<algorithm>
#include<random>
using std::begin;
using std::end;
using std::cout;
using std::endl;
inline double sqr(double v) {return v*v;}
double pearsonr(double *x,double *y,int n)
{
double sx=0.0,sy=0.0,sxy=0.0,sxx=0.0,syy=0.0;
int i;
for (i=0;i<n;++i) {
sx+=x[i];
sy+=y[i];
}
double mx=sx/n,my=sy/n;
for (i=0;i<n;++i) {
sxy += (x[i]-mx)*(y[i]-my);
sxx += sqr(x[i]-mx);
syy += sqr(y[i]-my);
}
return sxy/sqrt(sxx*syy);
}
int main() {
int i,j,k;
const int NSET=9;
const int N=500;
double x[NSET][N],y[NSET][N];
double r[NSET][NSET],r1[NSET][NSET],s;
double r0[NSET][NSET]={{1.0},
{0.6,1.0},
{0.2,0.1,1.0},
{0.1,0.2,0.7,1.0},
{0.1,0.2,0.6,0.6,1.0},
{0.1,0.2,0.1,0.1,0.1,1.0},
{0.1,0.2,0.1,0.2,0.2,0.6,1.0},
{0.1,0.1,0.1,0.1,0.1,0.1,0.1,1.0},
{0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.6,1.0},
};
for (i=0;i<NSET;++i) {
s=0.0;
for (j=0;j<=i;++j) {
if (j<i) {
r1[i][j]=r0[i][j];
for (k=0;k<j;++k) r1[i][j]-=r1[j][k]*r1[i][k];
if (0<j) r1[i][j]/=r1[j][j];
s+=sqr(r1[i][j]);
}
else r1[i][j]=sqrt(1-s);
}
}
// 一様分布で作成
std::random_device rnd;
std::mt19937 mt(rnd());
std::uniform_real_distribution<> rand01(0.0, 1.0);
for (i=0;i<NSET;++i) {
for (j=0;j<N;++j) {
x[i][j]=rand01(mt);
}
}
// 相関付ける
for (i=0;i<NSET;++i) {
for (j=0;j<N;++j) {
y[i][j]=0.0;
for (k=0;k<=i;++k) {
y[i][j] += x[k][j]*r1[i][k];
}
}
}
struct R {
double val;
int rank;
int id;
} data[N];
for (i=0;i<NSET;++i) {
for (j=0;j<N;++j) {
data[j].id = j;
data[j].val = y[i][j];
}
// valでソートする
auto val_more = [](const R& x, const R& y) { return x.val < y.val; };
std::sort(begin(data), end(data), val_more);
// 順位づけ
for (j=0;j<N;++j) {
data[j].rank = static_cast<int>(std::distance(begin(data), std::lower_bound(begin(data), end(data), data[j], val_more)))+1;
}
// idでソートする
auto id_more = [](const R& x, const R& y) { return x.id < y.id; };
std::sort(begin(data), end(data), id_more);
//
for (j=0;j<N;++j) x[i][j]= data[j].rank;
}
// 順位付けに従った値になったので，相関係数を検証する
for (i=0;i<NSET;++i) {
for (j=0;j<NSET;++j) {
r[i][j]=pearsonr(x[i],x[j],N);
}
}
cout << "r=";
for (i=0;i<NSET;++i) {
for (j=0;j<=i;++j) cout << "\t" << std::fixed << std::setprecision(3) << r[i][j];
cout << endl;
}
cout << endl;
cout << resetiosflags(std::ios_base::floatfield);
cout << "x=";
for (i=0;i<N;++i) {
for (j=0;j<NSET;++j) cout << "\t" << x[j][i];
cout << endl;
}
cout << endl;
cout << endl;
} 復活

補足

以下のように書いてみると約0.3秒で答えがでます。しかし，実際にやりたかったことは3セットではなく9セットなのです。同じように書いてみても実際的な計算時間では収まりません。4セットですら約38秒かかります。
その上，9セットは大まかには2-3-2-2のサブセットに分かれていて，与えられる相関係数はサブセット内は大きく0.63くらいですが，サブセット間は小さくて約0.07です。
なにか画期的な方法はないものでしょうか？

#include<iostream>
#include<cmath>
inline double sqr(double v) {return v*v;}
double pearsonr(double *x,double *y,int n)
{
double sx=0.0,sy=0.0,sxy=0.0,sxx=0.0,syy=0.0;
int i;
for (i=0;i<n;++i) {
sx+=x[i];
sy+=y[i];
}
double mx=sx/n,my=sy/n;
for (i=0;i<n;++i) {
sxy += (x[i]-mx)*(y[i]-my);
sxx += sqr(x[i]-mx);
syy += sqr(y[i]-my);
}
return sxy/sqrt(sxx*syy);
}
int main() {
double x[1000];
double r[3][3],rr[3][3],s[3];
double r0[3][3]={{1.0},
{0.7,1.0},
{0.6,0.5,1.0},
};
int n=500;
int i,j,k[3],kk[3];
for (i=0;i<n;++i) {
x[i]=x[i+n]=i+1;
}
double ss=999;
kk[0]=k[0]=0;
for (k[1]=1;k[1]<n;++k[1]) {
r[1][0] = pearsonr(x+k[1],x+k[0],n);
s[1]=sqr(r[1][0]-r0[1][0]);
if (ss<s[1]) continue;
for (k[2]=1;k[2]<n;++k[2]) {
for (i=0;i<2;++i) r[2][i] = pearsonr(x+k[2],x+k[i],n);
s[2]=s[1];
for (i=0;i<2;++i) s[2]+=sqr(r[2][i]-r0[2][i]);
if (ss>s[2]) {
ss=s[2];
for (i=1;i<3;++i) {
kk[i]=k[i];
for (j=0;j<i;++j) rr[i][j]=r[i][j];
}
}
}
}
std::cout << "kk=";
for (i=1;i<3;++i) std::cout << " " << kk[i];
std::cout << std::endl;
std::cout << "r=";
for (i=1;i<3;++i) {
for (j=0;j<i;++j) std::cout << " " << rr[i][j];
std::cout << std::endl;
}
std::cout << std::endl;
}