不均一確率密度分布からのN回サンプリングの最小値について(改変版)



上記の記事を利用して、累積確率密度分布がQ(v)=x^kなる統計量(最小値が0、最大値が1)の期待値、および、それをN回ランダムに繰り返したときに得られる最小の統計量の期待値を求めてみる。

N回ランダムに繰り返したときに、a以下の値が1回もない確率は(1-Q(a))^Nである。少なくとも1回はa以下の値が出る確率はR(a;k,N)=1-(1-Q(a))^Nである。N個の独立な統計量サンプルの中の最も小さい値の期待値はExp(min)=¥int_{min}^{max} v ¥times ¥frac{d}{da}R(a;k,N) dvであるから、上記の記事の通り、Exp(min)=max-¥int_{min}^{max} R(v;k,N) dvである。ここでR(v;k,N)=1-(1-x^k)^Nであり、min=0,max=1であるからExp(min)=1-¥int_0^1 (1-(1-v^k)^N)dv

今、¥int_0^1 1 dv=1であるから、Exp(min)=¥int_0^1 (1-v^k)^N dvまで簡単にできる。この積分の式は、kとNとの関数であるから、改めてF(k,N)=¥int_0^1 (1-v^k)^N dvと定める。式変形により、F(k,N)を、以下のようにNについての漸化式とすることで、解く。

  • F(k,N)=¥int_0^1 (1-v^k)^N dv = ¥int_0^1 (1-v^k)(1-v^k)^{N-1} dv
  • =¥int_0^1 ((1-v-k)^{N-1} -v ¥times v^{k-1}(1-v^k)^{N-1}) dv
  • =F(k,N-1)-¥int_0^1 v ¥times v^{k-1}(1-v^k)^{N-1}dv
  • ここで、((1-v^k)^N)’=-kNv^{k-1}(1-v^k)^{N-1}であることに注意すれば
  • =F(k,N-1)-¥int_0^1 v ¥times (-¥frac{1}{kN}((1-v^k)^N)’)dv
  • 第2項に部分積分の公式をあてはめることによって
  • =F(k,N-1) - ( ¥[v¥times (-(1-v^k)^N)¥]_0^1 -¥int_0^1 -¥frac{1}{kN}(1-v^k)^N dv)
  • =F(k,N-1) - ( 0-0 + ¥frac{1}{kN} F(k,N))
  • ここに
  • F(k,N)=F(k,N-1)-¥frac{1}{kN}F(k,N)なる式が成り立つことが示され、
  • F(k,N)=¥frac{kN}{kN+1}F(k,N-1)なる漸化式が得られた。
  • F(k,0)=¥int_0^1 (1-v^k)^0 dv = ¥[v¥]_0^1であるから
  • F(k,N)=¥prod_{i=1}^N ¥frac{ki}{ki+1}が得られた。
    • k=1,N=1のときF(1,1)=¥frac{1}{2}であり、これは、均一分布(k=1から1回サンプリングしたときに1番小さい値(サンプルの値そのもの)の期待値を確かに表している。均一分布からのN回サンプリングにおける最小値の期待値は¥prod_{i=1}^N¥frac{i}{i+1}=¥frac{1¥times 2 ¥times ... ¥times N}{2¥times 3 ¥times ... ¥times (N+1)}=¥frac{1}{N+1}となり、これは、12月11日の記事と符号している。

エクセルはこちら。ただし、非常に重いです。

Javaソースはこちら。


package StatUtils;

public class ProbabilityDenstiy {

public static void main(String[] args) {
BufferedWriter[] bw = new BufferedWriter[1];
String file = "out2.txt";
bw[0]=null;
int minN = 1000;
int maxN = 1000;//単一のデータがほしいときはminNに一致させる
int N=maxN-minN+1;
double startK = 0.8;
int numK=1;//単一のデータがほしいときは1
double kizamiK=0.1;
double[] K=new double[numK];
K[0] = startK;

//out += "\n";
try{

bw[0] = new BufferedWriter(new FileWriter(file));
//kinjishikiPartial(psTr0,maxn,part);
for(int i=1;i<K.length;i++){
K[i]=K[i-1]+kizamiK;
}
String out="\tk=\t";
for(int i=0;i<numK;i++){
out += K[i] + "\t";
}
System.out.println(out);
out +="\n";
BioBankPerm.InOut4.out3File(bw[0],out);
for(int i=minN;i<maxN+1;i++){
out ="N=\t" + i + "\t";
for(int j=0;j<numK;j++){
double tmp = expMinV1(K[j],i);
out += tmp + "\t";

}
System.out.println("N=" + i + "\t" + out);
out += "\n";
BioBankPerm.InOut4.out3File(bw[0],out);
}

bw[0].close();
}catch(Exception e){

System.out.println(e);

}

}
/*
* y=x^kなるCumulativ probability density functionの
* N回サンプルの最小値の期待値
*/
public static double expMinV1(double k, int N){
double ret = 0;
for(int i=1;i<=N;i++){
ret += Math.log(k*i)-Math.log(k*i+1);
}
ret = Math.pow(Math.E,ret);
return ret;
}
}