ADMB Documentation  11.1.2192
 All Classes Files Functions Variables Typedefs Friends Defines
nmonte.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: nmonte.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <admodel.h>
00008 
00009 void generate_actual_multivariate_mixture(int nvar, const dvector& a,
00010   const dvector& b, dmatrix& ch,long int& iseed, const double& lprob,
00011   const dvector& w);
00012 
00013 double cumd_mixture_02(const double& x);
00014 
00015 void generate_actual_multivariate_cauchy(int nvar, const dvector& a,
00016   const dvector& b, dmatrix& ch,long int& iseed, const double& lprob,
00017   const dvector& w);
00018 double log_likelihood_mixture_02(const double& x);
00019 double log_likelihood_mixture(const double& x);
00020 double inv_cumd_norm(const double& x);
00021 double inv_cumd_norm_ln(const double& x);
00022 double inv_cumd_norms(const double& x);
00023 double cumd_norm(const double& x);
00024 double ln_normal_tail_right(const double& x);
00025 double ln_normal_tail_left(const double& x);
00026 double ln_normal_tail(const double& x);
00027 double better_rand(long int& iseed);
00028 void get_bounded_mixture(double x1,double x2, const double& y,
00029   const double& density, long int& iseed);
00030 void get_bounded_cauchy(double x1,double x2, const double& y,
00031   const double& density, long int& iseed);
00032 void get_bounded_normal(double x1,double x2, const double& y,
00033   const double& density, long int& iseed);
00034 
00035 void get_bounded_normal_virtual(double x1,double x2, const double& _y,
00036   double& _log_density);
00037 
00038 dvector bounded_multivariate_normal(int nvar, const dvector& a1,
00039   const dvector& b1, dmatrix& ch,long int& iseed);
00040 
00041 void generate_virtual_multivariate(int nvar, const dvector& a,
00042   const dvector& b, const dmatrix& ch, const double& lprob, const dvector& eps);
00043 
00044 void generate_actual_multivariate(int nvar, const dvector& a, const dvector& b,
00045   const dmatrix& ch, long int& iseed, const double& lprob, const dvector& w);
00046 
00047 dvector bounded_robust_multivariate_normal(int nvar, const dvector& a1,
00048   dvector& b1, const dmatrix& ch, const dmatrix& ch3, const dmatrix& chinv,
00049   const dmatrix& ch3inv, double contaminant,long int& iseed,
00050   const double& lprob, const double& _lprob3, double& log_tprob,
00051   const int& _outflag)
00052 {
00053   double& lprob3=(double&) _lprob3;
00054   int& outflag=(int&) _outflag;
00055   dvector w(1,nvar);
00056   dvector a(1,nvar);
00057   dvector b(1,nvar);
00058   w.initialize();
00059   double r = better_rand(iseed);
00060   if (r>contaminant)
00061   {
00062     outflag=0;
00063     a=a1;
00064     b=b1;
00065     generate_actual_multivariate(nvar,a,b,ch,iseed,lprob,w);
00066     a=a1;
00067     b=b1;
00068     generate_virtual_multivariate(nvar,a,b,ch3,lprob3,ch3inv*w);
00069   }
00070   else
00071   {
00072     outflag=1;
00073     a=a1;
00074     b=b1;
00075     generate_actual_multivariate(nvar,a,b,ch3,iseed,lprob3,w);
00076     a=a1;
00077     b=b1;
00078     generate_virtual_multivariate(nvar,a,b,ch,lprob,chinv*w);
00079   }
00080   double tmpa=(lprob+.5*nvar);
00081   double tmpb=(lprob3-nvar*log(3.0)+.5*nvar);
00082   lprob3=lprob3-nvar*log(3.0);
00083   //tprob=(1.-contaminant)*exp(tmpa)+contaminant*exp(tmpb);
00084   if (tmpa>tmpb)
00085   {
00086     //tprob=exp(tmpa)*
00087       //( (1.-contaminant)+contaminant*exp(tmpb-tmpa) );
00088     log_tprob=tmpa+log((1.-contaminant)+contaminant*exp(tmpb-tmpa));
00089   }
00090   else
00091   {
00092     log_tprob=tmpb+log((1.-contaminant)*exp(tmpa-tmpb)+contaminant);
00093   }
00094   return w;
00095 }
00096 
00097 dvector bounded_multivariate_cauchy(int nvar, const dvector& a1,
00098   dvector& b1, const dmatrix& _ch,long int& iseed, const double& lprob,
00099   double& log_tprob, const int& _outflag)
00100 {
00101   dmatrix& ch=(dmatrix&)_ch;
00102   int& outflag=(int&) _outflag;
00103   dvector w(1,nvar);
00104   dvector a(1,nvar);
00105   dvector b(1,nvar);
00106   w.initialize();
00107   {
00108     outflag=0;
00109     a=a1;
00110     b=b1;
00111     generate_actual_multivariate_cauchy(nvar,a,b,ch,iseed,lprob,w);
00112   }
00113   double tmpa=(lprob+.5*nvar);
00114   log_tprob=tmpa;
00115   return w;
00116 }
00117 
00118 void generate_actual_multivariate_cauchy(int nvar, const dvector& _a,
00119   const dvector& _b, dmatrix& ch,long int& iseed, const double& _lprob,
00120   const dvector& _w)
00121 {
00122   double& lprob=(double&) _lprob;
00123   dvector& a=(dvector&) _a;
00124   dvector& b=(dvector&) _b;
00125   dvector& w=(dvector&) _w;
00126   double ah;
00127   double bl;
00128   dvector y(1,nvar);
00129   lprob=0;
00130   double log_density=0.0;;
00131   for (int i=1;i<=nvar;i++)
00132   {
00133     ah=a(i)/ch(i,i);
00134     bl=b(i)/ch(i,i);
00135 
00136     get_bounded_cauchy(ah,bl,y(i),log_density,iseed);
00137     lprob += log_density;
00138     for (int j=i;j<=nvar;j++)
00139     {
00140       double tmp=y(i)*ch(j,i);
00141       w(j)+=tmp;
00142       a(j)-=tmp;
00143       b(j)-=tmp;
00144     }
00145   }
00146 }
00147 
00148 dvector bounded_multivariate_mixture(int nvar, const dvector& a1,
00149   dvector& b1, const dmatrix& _ch,long int& iseed, const double& lprob,
00150   const int& _outflag)
00151 {
00152   dmatrix& ch=(dmatrix&)_ch;
00153   int& outflag=(int&) _outflag;
00154   dvector w(1,nvar);
00155   dvector a(1,nvar);
00156   dvector b(1,nvar);
00157   w.initialize();
00158   {
00159     outflag=0;
00160     a=a1;
00161     b=b1;
00162     generate_actual_multivariate_mixture(nvar,a,b,ch,iseed,lprob,w);
00163   }
00164   return w;
00165 }
00166 
00167 void generate_actual_multivariate_mixture(int nvar, const dvector& _a,
00168   const dvector& _b, dmatrix& ch,long int& iseed, const double& _lprob,
00169   const dvector& _w)
00170 {
00171   double& lprob=(double&) _lprob;
00172   dvector& a=(dvector&) _a;
00173   dvector& b=(dvector&) _b;
00174   dvector& w=(dvector&) _w;
00175   double ah;
00176   double bl;
00177   dvector y(1,nvar);
00178   lprob=0;
00179   double log_density=0.0;;
00180   for (int i=1;i<=nvar;i++)
00181   {
00182     ah=a(i)/ch(i,i);
00183     bl=b(i)/ch(i,i);
00184 
00185     get_bounded_mixture(ah,bl,y(i),log_density,iseed);
00186     lprob += log_density;
00187     for (int j=i;j<=nvar;j++)
00188     {
00189       double tmp=y(i)*ch(j,i);
00190       w(j)+=tmp;
00191       a(j)-=tmp;
00192       b(j)-=tmp;
00193     }
00194   }
00195 }
00196 
00197 void generate_actual_multivariate(int nvar, const dvector& _a,
00198   const dvector& _b, const dmatrix& ch, long int& iseed, const double& _lprob,
00199   const dvector& _w)
00200 {
00201   double& lprob=(double&) _lprob;
00202   dvector& a=(dvector&) _a;
00203   dvector& b=(dvector&) _b;
00204   dvector& w=(dvector&) _w;
00205   double ah;
00206   double bl;
00207   dvector y(1,nvar);
00208   lprob=0;
00209   double log_density=0.0;
00210   for (int i=1;i<=nvar;i++)
00211   {
00212     ah=a(i)/ch(i,i);
00213     bl=b(i)/ch(i,i);
00214 
00215     get_bounded_normal(ah,bl,y(i),log_density,iseed);
00216     lprob += log_density;
00217     for (int j=i;j<=nvar;j++)
00218     {
00219       double tmp=y(i)*ch(j,i);
00220       w(j)+=tmp;
00221       a(j)-=tmp;
00222       b(j)-=tmp;
00223     }
00224   }
00225 }
00226 
00227 void generate_virtual_multivariate(int nvar, const dvector& _a,
00228   const dvector& _b, const dmatrix& ch, const double& _lprob,
00229   const dvector& eps)
00230 {
00231   double& lprob=(double&) _lprob;
00232   dvector& a=(dvector&) _a;
00233   dvector& b=(dvector&) _b;
00234   double ah;
00235   double bl;
00236   dvector& y=(dvector&)(eps);
00237   lprob=0;
00238   double log_density;
00239   for (int i=1;i<=nvar;i++)
00240   {
00241     ah=a(i)/ch(i,i);
00242     bl=b(i)/ch(i,i);
00243     get_bounded_normal_virtual(ah,bl,y(i),log_density);
00244     lprob -= log_density;
00245     for (int j=i;j<=nvar;j++)
00246     {
00247       double tmp=y(i)*ch(j,i);
00248       a(j)-=tmp;
00249       b(j)-=tmp;
00250     }
00251   }
00252 }
00253 
00254 void get_bounded_cauchy(double x1,double x2, const double& _y,
00255   const double& _log_density, long int& iseed)
00256 {
00257   double& log_density=(double&) _log_density;
00258   double& y=(double&) _y;
00259   //double lp1;
00260   //double lp2;
00261   //double v;
00262   double w;
00263   double u = better_rand(iseed);
00264   {
00265     double lower=cumd_cauchy(x1);
00266     double upper=cumd_cauchy(x2);
00267     w=-log(upper-lower);
00268     u=u*.9998+.0001;
00269     y = inv_cumd_cauchy(u*(upper-lower)+lower);
00270   }
00271   log_density=w-log(1+0.5*y*y);
00272 }
00273 
00274 void get_bounded_mixture(double x1,double x2, const double& _y,
00275   const double& _log_density, long int& iseed)
00276 {
00277   double& y=(double&) _y;
00278   double& log_density=(double&) _log_density;
00279   //double lp1;
00280   //double lp2;
00281   //double v;
00282   double w;
00283   double u = better_rand(iseed);
00284   {
00285     //double lower=cumd_mixture_02(x1);
00286     //double upper=cumd_mixture_02(x2);
00287     double lower=cumd_mixture(x1);
00288     double upper=cumd_mixture(x2);
00289     w=-log(upper-lower);
00290     u=u*.9998+.0001;
00291     //y = inv_cumd_mixture_02(u*(upper-lower)+lower);
00292     y = inv_cumd_mixture(u*(upper-lower)+lower);
00293   }
00294   //log_density=w+log_likelihood_mixture_02(y);
00295   log_density=w+log_likelihood_mixture(y);
00296 }
00297 
00298 void get_bounded_normal(double x1,double x2, const double& _y,
00299   const double& _log_density, long int& iseed)
00300 {
00301   double& y=(double&) _y;
00302   double& log_density=(double&) _log_density;
00303   double lp1;
00304   double lp2;
00305   double v;
00306   double w;
00307   double u = better_rand(iseed);
00308   if (x2<-5.0)
00309   {
00310     lp1=ln_normal_tail(x1);
00311     lp2=ln_normal_tail(x2);
00312     v=exp(lp1-lp2);
00313     w=-lp2-log(1.0-v);
00314     y = inv_cumd_norm_ln( lp2 + log(u*(1.0-v)+v) );
00315   }
00316   else if (x1>5.0)
00317   {
00318     lp1=ln_normal_tail(x2);
00319     lp2=ln_normal_tail(x1);
00320     v=exp(lp1-lp2);
00321     w=-lp2-log(1.0-v);
00322     y = inv_cumd_norm_ln( lp2 + log(u*(1.0-v)+v) );
00323     y = -y;
00324   }
00325   else
00326   {
00327     double lower=cumd_norm(x1);
00328     double upper=cumd_norm(x2);
00329     w=-log(upper-lower);
00330     u=u*.9998+.0001;
00331     y = inv_cumd_norm(u*(upper-lower)+lower);
00332   }
00333   log_density=w-0.5*y*y;
00334 }
00335 
00336 void get_bounded_normal_virtual(double x1,double x2, const double& _y,
00337   double& _log_density)
00338 {
00339   double& y=(double&) _y;
00340   double& log_density=(double&) _log_density;
00341   double lp1;
00342   double lp2;
00343   double v;
00344   double w;
00345   if (x2<-5.0)
00346   {
00347     lp1=ln_normal_tail(x1);
00348     lp2=ln_normal_tail(x2);
00349     v=exp(lp1-lp2);
00350     w=-lp2-log(1.0-v);
00351   }
00352   else if (x1>5.0)
00353   {
00354     lp1=ln_normal_tail(x2);
00355     lp2=ln_normal_tail(x1);
00356     v=exp(lp1-lp2);
00357     w=-lp2-log(1.0-v);
00358   }
00359   else
00360   {
00361     double lower=cumd_norm(x1);
00362     double upper=cumd_norm(x2);
00363     w=-log(upper-lower);
00364   }
00365   log_density=w-0.5*y*y;
00366 }
00367 
00368 double ln_normal_tail(const double& x)
00369 {
00370   if (x<0) return ln_normal_tail_left(x);
00371   return ln_normal_tail_right(x);
00372 }
00373 
00374 double ln_normal_tail_left(const double& x)
00375 {
00376   if (x>-2.0)
00377   {
00378     cerr << "arugument of ln_normal_tail_left must be < -2.0 you have "
00379          << x << endl;
00380     exit(1);
00381   }
00382   return ln_normal_tail_right(-x);
00383 }
00384 
00385 double ln_normal_tail_right(const double& x)
00386 {
00387   if (x<2.0)
00388   {
00389     cerr << "arugument of ln_normal_tail_right must be > 2.0 you have "
00390          << x << endl;
00391     exit(1);
00392   }
00393   const double a3=5;
00394   const double a4=9;
00395   const double a5=129;
00396   double x2=x*x;
00397   double lz=log(0.3989422804)-.5*x2;
00398   double b1=x2+2;
00399   double b2=b1*(x2+4);
00400   double b3=b2*(x2+6);
00401   double b4=b3*(x2+8);
00402   double b5=b4*(x2+10);
00403   double tmp=lz-log(x) +
00404     log(1.0 -1.0/b1 +1.0/b2 - a3/b3 +a4/b4 -a5/b5);
00405   return tmp;
00406 }
00407 
00408 /*
00409 Same as linad99/cumdist.cpp -> double inv_cumd_norm_inner(const double& x)
00410 double inv_cumd_norm(const double& x)
00411 {
00412   const double c0=2.515517;
00413   const double c1=0.802853;
00414   const double c2=0.010328;
00415   const double d1=1.432788;
00416   const double d2=0.189269;
00417   const double d3=0.001308;
00418   if (x<=0 || x>=1.0)
00419   {
00420     cerr << "Illegal argument to inv_cumd_norm = " << x << endl;
00421     return 0;
00422   }
00423 
00424   if (x<0.5)
00425   {
00426     double t = sqrt(-2.*log(x));
00427     double p=((c2*t+c1)*t+c0)/((((d3*t+d2)*t+d1)*t)+1)-t;
00428     return p;
00429   }
00430   else if (x==0.5)
00431   {
00432     return 0.0;
00433   }
00434   else
00435   {
00436     double y=1.-x;
00437     double t = sqrt(-2.*log(y));
00438     double p=t-((c2*t+c1)*t+c0)/((((d3*t+d2)*t+d1)*t)+1);
00439     return p;
00440   }
00441 }
00442 */
00443 
00444 double inv_cumd_norm_ln(const double& x)
00445 {
00446   const double c0=2.515517;
00447   const double c1=0.802853;
00448   const double c2=0.010328;
00449   const double d1=1.432788;
00450   const double d2=0.189269;
00451   const double d3=0.001308;
00452   if (x>=0)
00453   {
00454     cerr << "Illegal argument to inv_cumd_norm = " << x << endl;
00455     return 0;
00456   }
00457 
00458   if (x<log(0.5))
00459   {
00460     double t = sqrt(-2.*x);
00461     double p=((c2*t+c1)*t+c0)/((((d3*t+d2)*t+d1)*t)+1)-t;
00462     return p;
00463   }
00464   else
00465   {
00466     double y;
00467     if (x<-1.e-10)
00468     {
00469       y=log(1.-exp(x));
00470     }
00471     else
00472     {
00473       y=log(-x)*log(1.0+x*(0.5+0.1666666666667*x));
00474     }
00475     double t = sqrt(-2.*y);
00476     double p=t-((c2*t+c1)*t+c0)/((((d3*t+d2)*t+d1)*t)+1);
00477     return p;
00478   }
00479 }
00480 
00481 /*
00482 Same as in linad99/cumdist.cpp
00483 double cumd_norm(const double& x)
00484 {
00485   const double b1=0.319381530;
00486   const double b2=-0.356563782;
00487   const double b3=1.781477937;
00488   const double b4=-1.821255978;
00489   const double b5=1.330274429;
00490   const double p=.2316419;
00491   if (x>30.0)
00492   {
00493     return 1.0;
00494   }
00495   else if (x<-30.0)
00496   {
00497     return 0.0;
00498   }
00499   else if (x>=0)
00500   {
00501     double u=1./(1+p*x);
00502     double y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00503     double z=1.0-0.3989422804*exp(-.5*x*x)*y;
00504     return z;
00505   }
00506   else
00507   {
00508     double w=-x;
00509     double u=1./(1+p*w);
00510     double y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00511     double z=0.3989422804*exp(-.5*x*x)*y;
00512     return z;
00513   }
00514 }
00515 */
00516 
00517 double cumd_mixture_02(const double& x)
00518 {
00519   double w=0.0;
00520   double w1=0.0;
00521   double w2=0.0;
00522   w1=cumd_norm(x);
00523   w2=cumd_cauchy(x);
00524   w=0.98*w1+.02*w2;
00525   return w;
00526 }
00527 
00528 
00529 double cumd_mixture(const double& x)
00530 {
00531   double w=0.0;
00532   double w1=0.0;
00533   double w2=0.0;
00534   w1=cumd_norm(x);
00535   w2=cumd_cauchy(x);
00536   w=0.95*w1+.05*w2;
00537   return w;
00538 }
00539 
00540 double log_likelihood_mixture(const double& x)
00541 {
00542   double y=0.5*x*x;
00543   double w=log(0.379*exp(-y)+0.01254/(1.+y));
00544   return w;
00545 }
00546 
00547 
00548 double log_likelihood_mixture_02(const double& x)
00549 {
00550   double y=0.5*x*x;
00551   double w=log(0.391*exp(-y)+0.004502/(1.+y));
00552   return w;
00553 }
00554 
00555 
00556 double inv_cumd_mixture(const double& zz)
00557 {
00558   if (zz<=0.5)
00559   {
00560     if (zz> 0.030328178)  // 2.02
00561     {
00562       const double beta=0.2000006361;
00563       const double a1= -1.20100758;
00564       const double a2=  0.705759703;
00565       const double a3= -0.3969207118;
00566       const double a4=  0.1013877547;
00567       const double b1=  0.4064582431;
00568       const double b2= -1.313226944;
00569       const double b3= -0.4745760236;
00570       const double b4=  0.8704844718;
00571       double t=2*zz;
00572       double x=-log(t);
00573       if (x>0) x=pow(x,beta);
00574       double pred=inv_cumd_cauchy(zz);
00575       double u=(((a4*x+a3)*x+a2)*x+a1)*x+1.0;
00576       double v=(((b4*x+b3)*x+b2)*x+b1)*x+1.0;
00577       return u/v*pred;
00578     }
00579     else if (zz>0.002235963385) // -10.0
00580     {
00581       const double beta=1.391943399;
00582       const double a1=-0.3960836641;
00583       const double a2=0.06458378977;
00584       const double a3=-0.005347047973;
00585       const double a4=0.0001938683488;
00586       const double b1=0.1252881802;
00587       const double b2=0.01855936157;
00588       const double b3=-0.01768441436;
00589       const double b4=0.001537604957;
00590       double t=2*zz;
00591       double x=-log(t);
00592       if (x>0) x=pow(x,beta);
00593       double pred=inv_cumd_cauchy(zz);
00594       double u=(((a4*x+a3)*x+a2)*x+a1)*x+1.0;
00595       double v=(((b4*x+b3)*x+b2)*x+b1)*x+1.0;
00596       return u/v*pred;
00597     }
00598     else
00599     {
00600       return inv_cumd_cauchy(20.*zz);
00601     }
00602   }
00603   else
00604   {
00605     return -inv_cumd_mixture(1-zz);
00606   }
00607 }
00608 
00609 double inv_cumd_mixture_02(const double& zz)
00610 {
00611   if (zz<=0.5)
00612   {
00613     if (zz> 0.030328178)  // 2.02
00614     {
00615       const double beta=0.20000048;
00616       const double a1=-2.1053015;
00617       const double a2=2.4746767;
00618       const double a3=-1.7449831;
00619       const double a4=0.49304616;
00620       const double b1=-0.27804868;
00621       const double b2=-1.391687;
00622       const double b3=1.0477872;
00623       const double b4=-0.09850461;
00624       double t=2*zz;
00625       double x=-log(t);
00626       if (x>0) x=pow(x,beta);
00627       double pred=inv_cumd_cauchy(zz);
00628       double u=(((a4*x+a3)*x+a2)*x+a1)*x+1.0;
00629       double v=(((b4*x+b3)*x+b2)*x+b1)*x+1.0;
00630       return u/v*pred;
00631     }
00632     else if (zz>0.002235963385) // -10.0
00633     {
00634       const double beta=1.2392131;
00635       const double a1=-0.42498166;
00636       const double a2=0.07016042;
00637       const double a3=-0.0053910956;
00638       const double a4=0.0001646762;
00639       const double b1=0.18071484;
00640       const double b2=0.00029307283;
00641       const double b3=-0.014820049;
00642       const double b4=0.0013084549;
00643       double t=2*zz;
00644       double x=-log(t);
00645       if (x>0) x=pow(x,beta);
00646       double pred=inv_cumd_cauchy(zz);
00647       double u=(((a4*x+a3)*x+a2)*x+a1)*x+1.0;
00648       double v=(((b4*x+b3)*x+b2)*x+b1)*x+1.0;
00649       return u/v*pred;
00650     }
00651     else
00652     {
00653       return inv_cumd_cauchy(50.*zz);
00654     }
00655   }
00656   else
00657   {
00658     return -inv_cumd_mixture_02(1-zz);
00659   }
00660 }
00661 
00662 /*
00663 dvariable inv_cumd_norm(const prevariable& x)
00664 {
00665   const double c0=2.515517;
00666   const double c1=0.802853;
00667   const double c2=0.010328;
00668   const double d1=1.432788;
00669   const double d2=0.189269;
00670   const double d3=0.001308;
00671   if (x<=0 || x>=1.0)
00672   {
00673     cerr << "Illegal argument to inv_cumd_norm = " << x << endl;
00674     return 0;
00675   }
00676 
00677   if (x<0.5)
00678   {
00679     dvariable t = sqrt(-2.*log(x));
00680     dvariable p=((c2*t+c1)*t+c0)/((((d3*t+d2)*t+d1)*t)+1)-t;
00681     return p;
00682   }
00683   else if (x==0.5)
00684   {
00685     return 0.0;
00686   }
00687   else
00688   {
00689     dvariable y=1.-x;
00690     dvariable t = sqrt(-2.*log(y));
00691     dvariable p=t-((c2*t+c1)*t+c0)/((((d3*t+d2)*t+d1)*t)+1);
00692     return p;
00693   }
00694 }
00695 */