カオス法によるq-Gauss生成アルゴリズム

http://ssuopt.amp.i.kyoto-u.ac.jp/akihiro/wiki/index.php?%A5%AB%A5%AA%A5%B9%CB%A1%A4%CB%A4%E8%A4%EBq-Gauss%C0%B8%C0%AE%A5%A2%A5%EB%A5%B4%A5%EA%A5%BA%A5%E0
Last-modified: 2013-10-26 (Sat) 13:43:03 (1972d)

研究内容

カオス法によるq-Gauss生成アルゴリズム

カオス法によるq-Gauss生成アルゴリズムとは次数dの単位円周上二次元写像Q_d(x)とP_d(x,y)と、動径方向の区分線形写像のq-指数関数、q-対数関数を用いた微分同相変換からなる疑似乱数生成アルゴリズムである。この方法を用いることにより独立なq-Gauss分布に従う疑似乱数を生成することが可能である。数学的な詳細は IEEE Transactions on Information Theory, Vol. 59 (2013) pp. 3199--3209 を参照のこと。

サンプルコード

単位円周上二次元写像の次数d=8,動径方向はテント写像(l=2,c=1)を用いている。

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <strings.h>
double qnormal_x,qnormal_y,qnormal_z;
double expq(double q, double w){
  if(q==1.0){
    return(exp(w));
  }
  else{
    return (expl(log(1.0+(1.0-q)*w)/(1.0-q)));
  }
}
double lnq(double q, double w){
  if(q==1.0){
    return(log(w));
  }
  else{
    return ((exp(log(w)*(1.0-q))-1.0)/(1.0-q));
  }
}
void setseed_qnormal(double v0, double z0){
  qnormal_x = sqrt(1-v0*v0); 
  qnormal_y = v0;
  qnormal_z = z0;
}
double Q8(double w, double v){
  return(8*w*v*(((16.0*w*w-24.0)*w*w+10.0)*w*w-1.0));
}
double P8(double w){
  return((((128.0*w*w-256.0)*w*w
         +160.0)*w*w-32.0)*w*w+1.0);
}
double f(double z){
  return(1.0-fabs(1.0-1.99999*z));
}
void qnormal(double q){
  double qq;
  qnormal_y = Q8(qnormal_x,qnormal_y);
  qnormal_x = P8(qnormal_x);
  qq = (q+1.0)/(3.0-q); 
  qnormal_z = f(expq(qq,-qnormal_z*qnormal_z*0.5));
  qnormal_z = sqrt(-2.0*lnq(qq,qnormal_z));
}
int main(int argc, char *argv[]){
  double q,v0,z0,eta,xi;
  int i;
  if(argc != 4){
     printf("%s q v0 z0\n",argv[0]);
     exit(0);
  }
  q = (double)atof(argv[1]);
  v0 = (double)atof(argv[2]);
  z0 = (double)atof(argv[3]);
  setseed_qnormal(v0,z0);
  for(i=0;i<10000;i++){
    qnormal(q);
    xi = qnormal_x*qnormal_z;
    eta = qnormal_y*qnormal_z;
    printf("%lf %lf\n",xi,eta);
  }
}

Front page   Edit Freeze Diff Backup Upload Copy Rename Reload   New List of pages Search Recent changes   Help   RSS of recent changes