关闭

GSL-蒙特卡洛积分

1192人阅读 评论(0) 收藏 举报
分类:

GSL-蒙特卡洛积分  TR:SAN E:VISUALSAN@YAHOO.CN  2011.3.23

---------------------------------------------------------------------------

Gsl中包含有计算N重积分的蒙特卡洛方法,不过只能计算积分上下限确定的多重积分,对于上下限带有参数的积分无能为力,表达式如下:

三种蒙特卡洛积分实现方法分别定义在头文件:

#include <gsl/gsl_monte_miser.h>

#include <gsl/gsl_monte_plain.h>

#include <gsl/gsl_monte_vegas.h>

每组方法的计算思路是一样的:生成控制变量、初始化控制变量、计算积分、释放控制变量。每种方法都将用到随机数产生器。其中需要用户自定义被积分函数:

double (*f)(double * x_array, size_t dim, void * params);

其中x是积分变量数组,dim为积分变量个数,parms为额外传递的参数

对于二次函数f(x,y)=ax2+bxy+cy2,当a=3,b=2,c=1时,下面将定义目标函数:

复制代码
struct my_par
{
double a,b,c;
};
double my_f(double* x_array, size_t dim, void*params)
{
my_par
*mp=(my_par*)params;
if(dim!=2)
{
fprintf(stderr,
"dim!=2");
abort();
}
return x_array[0]*x_array[0]*mp->a+
x_array[
0]*x_array[1]*mp->b+
x_array[
1]*x_array[1]*mp->c;
}

void test_monte()
{
my_par par
={3,2,1};
gsl_monte_function f;
f.dim
=2;
f.f
=my_f;
f.
params=&par;

}
复制代码

下面介绍用法:

  1. plainmonte carlo

     包含头文件:#include <gsl/gsl_monte_plain.h>

     使用步骤:

(1)       gsl_monte_plain_state* gsl_monte_plain_alloc(size_t dim);

               申请一个计算dim重积分的gsl_monte_plain_state类型控制变量

(2)       int gsl_monte_plain_init(gsl_monte_plain_state* state);

               初始化控制变量

(3)       int gsl_monte_plain_integrate (const gsl_monte_function * f,

                           const double xl[], const double xu[],

                           const size_t dim,

                           const size_t calls,

                           gsl_rng * r,

                           gsl_monte_plain_state * state,

                           double *result, double *abserr);

        计算monte carlo积分,参数意义如下:

         :f->被积分函数

         :xl->积分下限数组,共有dim个

         :xu->积分上限数组,共有dim个

         :dim->积分重数

         :calls->目标函数被调用次数,比如500000

         :r->随机数产生器,比如r= gsl_rng_default

         :state->状态空间,由gsl_monte_plain_alloc申请

         :result->计算结果

         :abserr->绝对误差

例子:计算多重积分

复制代码
#include <iostream>
#include
<time.h>
#include
<gsl/gsl_rng.h>
#include
<gsl/gsl_integration.h>
#include
<gsl/gsl_monte.h>
#include
<gsl/gsl_monte_miser.h>
#include
<gsl/gsl_monte_plain.h>
#include
<gsl/gsl_monte_vegas.h>
usingnamespace std;
//visualsan@yahoo.cn SAN NUAA 2011.3.23
struct my_par
{
double a,b,c;
};
double my_f(double* x_array, size_t dim, void*params)
{
my_par
*mp=(my_par*)params;
if(dim!=2)
{
fprintf(stderr,
"dim!=2");
abort();
}
double v= x_array[0]*x_array[0]*mp->a+
x_array[
0]*x_array[1]*mp->b+
x_array[
1]*x_array[1]*mp->c;

return exp(-1/v)/v;
}

void test_monte()
{
my_par par
={3,2,1};
gsl_monte_function f;
f.dim
=2;
f.f
=my_f;
f.
params=&par;

int calls=500000;
double xl[]={0,0},xu[]={3.14,3.14};
double result,er;
gsl_monte_plain_state
*ps=gsl_monte_plain_alloc(2);
const gsl_rng_type*tp=gsl_rng_minstd;
gsl_rng
*pr=gsl_rng_alloc(tp);
gsl_monte_plain_init(ps);
gsl_monte_plain_integrate(
&f,
xl,xu,
2,calls,
pr,ps,
&result,&er);


cout
<<result<<endl;//->0.913118
cout<<er<<endl; //->0.0011711
/*matlab cmd:
fun = @(x,y) exp(-1./(3.*x.*x+2.*x.*y+y.*y))./( 3.*x.*x+2.*x.*y+y.*y );
quad2d(fun,0,3.14,0,3.14)
ans=0.913889482268682
*/

}
void main()
{
test_monte();
}


复制代码
2.    miser monte carlo

     包含头文件:#include <gsl/gsl_monte_ miser.h>

  1. gsl_monte_miser_state* gsl_monte_miser_alloc(size_t dim);申请计算dim重积分的控制变量
  2. int gsl_monte_miser_init(gsl_monte_miser_state* state);初始化积分控制变量
  3. 积分函数:

        int gsl_monte_miser_integrate(gsl_monte_function * f,

                                         const double xl[], const double xh[],

                                         size_t dim, size_t calls,

                                         gsl_rng *r,

                                         gsl_monte_miser_state* state,

                                         double *result, double *abserr);

    4.void gsl_monte_miser_free(gsl_monte_miser_state* state);释放控制变量

miser carlo算法有几个参数需要设置,这些参数定义在gsl_monte_miser_state。

estimate_frac:这个参数在方差计算的递归调用过程中,指明每次函数调用次数占可用函数次数的百分比,默认值为0.1

min_calls:这个参数指明每次方差评估过程中函数调用的最小次数,如果申请函数调用次数n小于min_calls,则n=min_calls。min_calls默认值为16*dim.

min_calls_per_bisection:对分步骤中,函数调用的最小次数,默认值16* min_calls

alpha:指定两个子区域的方差如何组合,默认值为2

Dither:引入变异率,打破被积函数在积分区间的对称性。默认为0,当需要引入时,建议取值0.1

3.vegasmonte carlo

包含头文件gsl_monte_vegas.h

使用步骤:

1.)gsl_monte_vegas_state* gsl_monte_vegas_alloc(size_t dim);申请一个dim重积分的控制变量

2.)int gsl_monte_vegas_init(gsl_monte_vegas_state* state);初始化控制变量

3.)int gsl_monte_vegas_integrate(gsl_monte_function * f,

                                         double xl[], double xu[],

                                         size_t dim, size_t calls,

                                         gsl_rng * r,

                                         gsl_monte_vegas_state *state,

                                         double* result, double* abserr);

4.)void gsl_monte_vegas_free (gsl_monte_vegas_state* state);释放控制变量

下面这个例子是GSL手册上的例子:

计算积分

复制代码
#include <iostream>
#include
<time.h>
#include
<gsl/gsl_math.h>
#include
<gsl/gsl_monte_miser.h>
#include
<gsl/gsl_monte_plain.h>
#include
<gsl/gsl_monte_vegas.h>


usingnamespace std;
//e: visualsan@yahoo.cn tr:SAN ad:NUAA t:2011
constdouble exact=1.3932039296;
double my_f(double* k, size_t dim, void*params)
{
double A =1.0/ ( M_PI * M_PI * M_PI );
return A / ( 1.0- cos(k[0]) * cos(k[1]) * cos(k[2]) );
}
void disp_r(char*t,double r,double er,double ct)
{
printf(
"%s----------------------t=%.1fms \n",t,ct);
printf(
"result=%.6f\n",r);
printf(
"sigma=%.6f\n",er);
printf(
"exact=%.6f\n",exact);
printf(
"error=%.6f = %.1g sigma\n",r-exact,fabs(exact-r) / er);
}
void main()
{
cout
<<"monte carlo sample (visualsan@yahoo.cn tr:SAN )\n";
clock_t t1;
gsl_monte_function f;
f.dim
=3;
f.f
=my_f;


int calls=500000;
double xl[]={0,0,0},xu[]={M_PI,M_PI,M_PI};
double result,er;

gsl_rng_env_setup();
const gsl_rng_type*tp=gsl_rng_minstd;
gsl_rng
*pr=gsl_rng_alloc(tp);



//plain
t1=clock();
gsl_monte_plain_state
*ps=gsl_monte_plain_alloc(3);
gsl_monte_plain_init(ps);
gsl_monte_plain_integrate(
&f,
xl,xu,
3,calls,
pr,ps,
&result,&er);
disp_r(
"plain",result,er, (double)(clock()-t1) );
gsl_monte_plain_free(ps);

//miser
t1=clock();
gsl_monte_miser_state
*pm=gsl_monte_miser_alloc(3);
gsl_monte_miser_init(pm);
gsl_monte_miser_integrate(
&f,
xl,xu,
3,calls,
pr,pm,
&result,&er);
disp_r(
"miser",result,er,(double)(clock()-t1));
gsl_monte_miser_free(pm);

//vegas
t1=clock();
gsl_monte_vegas_state
*pv=gsl_monte_vegas_alloc(3);
gsl_monte_vegas_init(pv);
gsl_monte_vegas_integrate(
&f,
xl,xu,
3,10000,
pr,pv,
&result,&er);
disp_r(
"vegas warm-up",result,er,(double)(clock()-t1));

printf(
"convering......\n");
t1
=clock();
do
{
gsl_monte_vegas_integrate(
&f,
xl,xu,
3,calls/5,
pr,pv,
&result,&er);
printf(
"result=%.6f,sigma=%.6f,chisq/dof=%.1f\n",result,er,pv->chisq);
}
while (fabs(pv->chisq-1.0)>0.5);
gsl_monte_vegas_free(pv);
disp_r(
"vegas final",result,er,(double)(clock()-t1));

}
复制代码

-------------------------------------------------------------------------------------------------

                                                                      SAN VISUALSAN@YAHOO.CN  2011.3.23  NUAA

0
0

查看评论
* 以上用户言论只代表其个人观点,不代表CSDN网站的观点或立场
    个人资料
    • 访问:547551次
    • 积分:7124
    • 等级:
    • 排名:第3270名
    • 原创:7篇
    • 转载:315篇
    • 译文:0篇
    • 评论:60条
    最新评论