ADMB Documentation  -a65f1c97
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
vcumdist.cpp
Go to the documentation of this file.
1 /*
2  * $Id$
3  *
4  * Author: David Fournier
5  * Copyright (c) 2008-2012 Regents of the University of California
6  */
11 #include <fvar.hpp>
12 
13 dvariable inv_cumd_norm(const prevariable& x);
14 prevariable& cumd_norm(const prevariable& x);
15 /*
16 double normal_tail_right(const double& x)
17 {
18  const double a3=5;
19  const double a4=9;
20  const double a5=129;
21  double x2=x*x;
22  double z=0.3989422804*exp(-.5*x2);
23  double b1=x2+2;
24  double b2=b1*(x2+4);
25  double b3=b2*(x2+6);
26  double b4=b3*(x2+8);
27  double b5=b4*(x2+10);
28  double tmp=z/x*(1.0 -1.0/b1 +1.0/b2 - a3/b3 +a4/b4 -a5/b5);
29  return tmp;
30 }
31 */
32 
37 dvariable inv_cumd_norm_inner(const prevariable& x)
38 {
39  gradient_structure* gs = gradient_structure::_instance;
40  if (++gs->RETURN_PTR > gs->MAX_RETURN)
41  gs->RETURN_PTR = gs->MIN_RETURN;
42 
43  gs->RETURN_ARRAYS_INCREMENT();
44  const double c0=2.515517;
45  const double c1=0.802853;
46  const double c2=0.010328;
47  const double d1=1.432788;
48  const double d2=0.189269;
49  const double d3=0.001308;
50  if (x<=0 || x>=1.0)
51  {
52  //cerr << "Illegal argument to inv_cumd_norm = " << x << endl;
53  gs->RETURN_ARRAYS_DECREMENT();
54  return 0.0;
55  }
56 
57  if (x<=0.5)
58  {
59  //dvariable t = sqrt(-2.*log(x));
60  //dvariable p=((c2*t+c1)*t+c0)/((((d3*t+d2)*t+d1)*t)+1)-t;
61  //RETURN_ARRAYS_DECREMENT();
62  //return p;
63 
64  double tt = sqrt(-2.*log(value(x)));
65  double u=(c2*tt+c1)*tt+c0;
66  double v=((d3*tt+d2)*tt+d1)*tt+1;
67  double vinv=1.0/v;
68  double pp=u*vinv-tt;
69 
70  //double pp=u*vinv-tt;
71  double dfu=vinv;
72  double dfvinv=u;
73  double dftt=-1.0;
74  //double vinv=1.0/v;
75  double dfv=-vinv*vinv*dfvinv;
76  //double v=((d3*tt+d2)*tt+d1)*tt+1;
77  dftt+=((3*d3*tt+2.0*d2)*tt+d1)*dfv;
78  //double u=(c2*tt+c1)*tt+c0;
79  dftt+=(2.0*c2*tt+c1)*dfu;
80  //double tt = sqrt(-2.*log(value(x)));
81  double dfx=-1.0/(tt*value(x))*dftt;
82 
83  gs->RETURN_ARRAYS_DECREMENT();
84  gs->RETURN_PTR->v->x=pp;
85  gs->GRAD_STACK1->set_gradient_stack(default_evaluation,
86  &(gs->RETURN_PTR->v->x), &(x.v->x),dfx);
87  }
88  else if (x==0.5)
89  {
90  cout << "can't happen" << endl;
91  ad_exit(1);
92  //return 0.0;
93  }
94  else
95  {
96  //dvariable y=1.-x;
97  //dvariable t = sqrt(-2.*log(y));
98 
99  //dvariable p=t-((c2*t+c1)*t+c0)/((((d3*t+d2)*t+d1)*t)+1);
100  //RETURN_ARRAYS_DECREMENT();
101  //return p;
102 
103  double yy=1.-value(x);
104  double tt = sqrt(-2.*log(yy));
105  double u=((c2*tt+c1)*tt+c0);
106  double v=((d3*tt+d2)*tt+d1)*tt+1;
107  double vinv=1/v;
108  double pp=tt-u*vinv;
109 
110  //double pp=tt-u*vinv;
111  double dfu=-vinv;
112  double dfvinv=-u;
113  double dftt=1.0;
114  //double vinv=1.0/v;
115  double dfv=-vinv*vinv*dfvinv;
116  //double v=((d3*tt+d2)*tt+d1)*tt+1;
117  dftt+=((3*d3*tt+2.0*d2)*tt+d1)*dfv;
118  //double u=(c2*tt+c1)*tt+c0;
119  dftt+=(2.0*c2*tt+c1)*dfu;
120  //double tt = sqrt(-2.*log(yy));
121  double dfy=-1.0/(tt*yy)*dftt;
122  //double yy=1.-value(x);
123  double dfx=-dfy;
124 
125  gs->RETURN_ARRAYS_DECREMENT();
126  gs->RETURN_PTR->v->x=pp;
127  gs->GRAD_STACK1->set_gradient_stack(default_evaluation,
128  &(gs->RETURN_PTR->v->x), &(x.v->x),dfx);
129  }
130  return *(gs->RETURN_PTR);
131 }
132 
137 dvariable inv_cumd_norm(const prevariable& x)
138 {
139  dvariable y=inv_cumd_norm_inner(x);
140  if (x>1.e-30)
141  y+=2.50662827*exp(.5*y*y)*(x-cumd_norm(y));
142  return y;
143 }
144 
149 dvariable old_cumd_norm(const prevariable& x)
150 {
151  gradient_structure* gs = gradient_structure::_instance;
152  gs->RETURN_ARRAYS_INCREMENT();
153 
154  const double b1=0.319381530;
155  const double b2=-0.356563782;
156  const double b3=1.781477937;
157  const double b4=-1.821255978;
158  const double b5=1.330274429;
159  const double p=.2316419;
160  if (x>=0)
161  {
162  dvariable u=1./(1+p*x);
163  dvariable y= ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
164  dvariable z=1.0-0.3989422804*exp(-.5*x*x)*y;
165  gs->RETURN_ARRAYS_DECREMENT();
166  return z;
167  }
168  else
169  {
170  dvariable w=-x;
171  dvariable u=1./(1+p*w);
172  dvariable y= ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
173  dvariable z=0.3989422804*exp(-.5*x*x)*y;
174  gs->RETURN_ARRAYS_DECREMENT();
175  return z;
176  }
177 }
178 
184 prevariable& cumd_norm(const prevariable& _x)
185 {
186  gradient_structure* gs = gradient_structure::_instance;
187  if (++gs->RETURN_PTR > gs->MAX_RETURN)
188  gs->RETURN_PTR = gs->MIN_RETURN;
189 
190  dvariable* RETURN_PTR = gs->RETURN_PTR;
191 
192  double x=value(_x);
193  const double b1=0.319381530;
194  const double b2=-0.356563782;
195  const double b3=1.781477937;
196  const double b4=-1.821255978;
197  const double b5=1.330274429;
198  const double b55=b5*5;
199  const double b44=b4*4;
200  const double b33=b3*3;
201  const double b22=b2*2;
202  const double p=.2316419;
203  if (x>=0)
204  {
205  double u=1./(1+p*x);
206  double y= ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
207  double tmp1=-0.3989422804*exp(-.5*x*x);
208  double z=1.0+tmp1*y;
209  RETURN_PTR->v->x=z;
210 
211 
212  // double z=1.0+tmp1*y;
213  double dftmp1=y;
214  double dfy=tmp1;
215 
216  // double tmp1=-0.3989422804*exp(-.5*x*x);
217  double dfx=-tmp1*x*dftmp1;
218 
219  //double y= ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
220  double dfu= ((((b55*u+b44)*u+b33)*u+b22)*u+b1)*dfy;
221 
222  //double u=1./(1+p*x);
223  dfx-=u*u*p*dfu;
224 
225  gs->GRAD_STACK1->set_gradient_stack(default_evaluation,
226  &(RETURN_PTR->v->x), &(_x.v->x), dfx);
227  }
228  else
229  {
230  double w=-x;
231  double u=1./(1+p*w);
232  double y= ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
233  double tmp1=0.3989422804*exp(-.5*x*x);
234  double z=tmp1*y;
235 
236  //double z=tmp1*y;
237  double dftmp1=y;
238  double dfy=tmp1;
239 
240  //double tmp1=0.3989422804*exp(-.5*x*x);
241  double dfx=-tmp1*x*dftmp1;
242 
243  //double y= ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
244  double dfu= ((((b55*u+b44)*u+b33)*u+b22)*u+b1)*dfy;
245 
246  //double u=1./(1+p*w);
247  double dfw=-u*u*p*dfu;
248 
249  //double w=-value(x);
250  dfx-=dfw;
251 
252  RETURN_PTR->v->x=z;
253  gs->GRAD_STACK1->set_gradient_stack(default_evaluation,
254  &(RETURN_PTR->v->x), &(_x.v->x),dfx);
255  }
256  return *RETURN_PTR;
257 }
258 
263 dvar_vector inv_cumd_norm(const dvar_vector& x)
264 {
265  int mmin=x.indexmin();
266  int mmax=x.indexmax();
267  dvar_vector tmp(mmin,mmax);
268  for (int i=mmin;i<=mmax;i++)
269  {
270  tmp(i)=inv_cumd_norm(x(i));
271  }
272  return tmp;
273 }
274 
279 prevariable& bounded_cumd_norm(const prevariable& _x,double beta)
280 {
281  gradient_structure* gs = gradient_structure::_instance;
282  if (++gs->RETURN_PTR > gs->MAX_RETURN)
283  gs->RETURN_PTR = gs->MIN_RETURN;
284 
285  double x=value(_x);
286  const double b1=0.319381530;
287  const double b2=-0.356563782;
288  const double b3=1.781477937;
289  const double b4=-1.821255978;
290  const double b5=1.330274429;
291  const double b55=b5*5;
292  const double b44=b4*4;
293  const double b33=b3*3;
294  const double b22=b2*2;
295  const double p=.2316419;
296  if (x>=0)
297  {
298  double u=1./(1+p*x);
299  double y= ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
300  double tmp1=-0.3989422804*exp(-.5*x*x);
301  double z1=1.0+tmp1*y;
302  double z=beta*(z1-0.5)+0.5;
303  gs->RETURN_PTR->v->x=z;
304 
305  //double z=beta*(z1-0.5)+0.5
306  double dfz1=beta;
307 
308  // double z1=1.0+tmp1*y;
309  double dftmp1=y*dfz1;
310  double dfy=tmp1*dfz1;
311 
312  // double tmp1=-0.3989422804*exp(-.5*x*x);
313  double dfx=-tmp1*x*dftmp1;
314 
315  //double y= ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
316  double dfu= ((((b55*u+b44)*u+b33)*u+b22)*u+b1)*dfy;
317 
318  //double u=1./(1+p*x);
319  dfx-=u*u*p*dfu;
320 
321  gs->GRAD_STACK1->set_gradient_stack(default_evaluation,
322  &(gs->RETURN_PTR->v->x), &(_x.v->x), dfx);
323  }
324  else
325  {
326  double w=-x;
327  double u=1./(1+p*w);
328  double y= ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
329  double tmp1=0.3989422804*exp(-.5*x*x);
330  double z1=tmp1*y;
331  double z=beta*(z1-0.5)+0.5;
332 
333  //double z=beta*(z1-0.5)+0.5;
334  double dfz1=beta;
335 
336  //double z1=tmp1*y;
337  double dftmp1=y*dfz1;
338  double dfy=tmp1*dfz1;
339 
340  //double tmp1=0.3989422804*exp(-.5*x*x);
341  double dfx=-tmp1*x*dftmp1;
342 
343  //double y= ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
344  double dfu= ((((b55*u+b44)*u+b33)*u+b22)*u+b1)*dfy;
345 
346  //double u=1./(1+p*w);
347  double dfw=-u*u*p*dfu;
348 
349  //double w=-value(x);
350  dfx-=dfw;
351 
352  gs->RETURN_PTR->v->x=z;
353  gs->GRAD_STACK1->set_gradient_stack(default_evaluation,
354  &(gs->RETURN_PTR->v->x), &(_x.v->x),dfx);
355  }
356  return *(gs->RETURN_PTR);
357 }
358 
363 dvariable inv_cumd_norm_logistic(const prevariable& x,double p)
364 {
365 #if !defined(OPT_LIB)
366  if (0.0<p || 1.0>p)
367  {
368  cerr << "Error in dvariable inv_cumd_norm_logistic -- illegal p value = "
369  << p << endl;
370  ad_exit(1);
371  }
372 #endif
373  dvariable y=inv_cumd_norm_inner(x);
374  y+=( (1.0-p)*2.50662827*exp(.5*y*y) +
375  p*exp(-x)/square(1.0+exp(-x)) ) * (x-cumd_norm_logistic (y,p));
376  return y;
377 }
378 
383 prevariable& cumd_norm_logistic(const prevariable& _x,double p)
384 {
385  return (1.0-p)*cumd_norm(_x)+p*cumd_logistic(_x);
386 }
prevariable
Base class for dvariable.
Definition: fvar.hpp:1315
old_cumd_norm
dvariable old_cumd_norm(const prevariable &x)
Description not yet available.
Definition: vcumdist.cpp:149
gradient_structure::RETURN_PTR
dvariable * RETURN_PTR
Definition: gradient_structure.h:123
bounded_cumd_norm
double bounded_cumd_norm(const double x, double beta)
Description not yet available.
Definition: cumdist.cpp:119
x
#define x
inv_cumd_norm
double inv_cumd_norm(const double &x)
Description not yet available.
Definition: cumdist.cpp:78
ad_exit
exitptr ad_exit
Definition: gradstrc.cpp:53
gradient_structure::MIN_RETURN
dvariable * MIN_RETURN
Definition: gradient_structure.h:124
dvar_vector
ADMB variable vector.
Definition: fvar.hpp:2172
grad_stack::set_gradient_stack
void set_gradient_stack(void(*func)(void), double *dep_addr, double *ind_addr1=NULL, double mult1=0, double *ind_addr2=NULL, double mult2=0)
Description not yet available.
Definition: fvar.hpp:1045
endl
prnstream & endl(prnstream &)
sqrt
d3_array sqrt(const d3_array &arr3)
Author: David Fournier Copyright (c) 2008-2012 Regents of the University of California.
Definition: d3arr2c.cpp:11
gradient_structure::RETURN_ARRAYS_INCREMENT
void RETURN_ARRAYS_INCREMENT()
Definition: gradstrc.cpp:478
fvar.hpp
Author: David Fournier Copyright (c) 2008-2012 Regents of the University of California.
exp
d3_array exp(const d3_array &arr3)
Returns d3_array results with computed exp from elements in arr3.
Definition: d3arr2a.cpp:28
default_evaluation
void default_evaluation(void)
Description not yet available.
Definition: def_eval.cpp:61
beta
dvariable beta(const prevariable &a, const prevariable &b)
Beta density function.
Definition: vbeta.cpp:18
gradient_structure::_instance
static _THREAD gradient_structure * _instance
Definition: gradient_structure.h:114
prevariable::v
double_and_int * v
pointer to the data
Definition: fvar.hpp:1333
gradient_structure::MAX_RETURN
dvariable * MAX_RETURN
Definition: gradient_structure.h:125
dvar_vector::indexmin
int indexmin() const
Definition: fvar.hpp:2287
w
#define w
inv_cumd_norm_inner
double inv_cumd_norm_inner(double x)
Author: David Fournier Copyright (c) 2008-2012 Regents of the University of California.
Definition: cumdist.cpp:40
gradient_structure::RETURN_ARRAYS_DECREMENT
void RETURN_ARRAYS_DECREMENT()
Definition: gradstrc.cpp:511
cumd_norm_logistic
double cumd_norm_logistic(double _x, double p)
Description not yet available.
Definition: ccumdlog.cpp:78
cumd_logistic
double cumd_logistic(const double &x)
Description not yet available.
Definition: ccumdlog.cpp:13
value
dvector value(const df1_one_vector &v)
Definition: df11fun.cpp:69
gradient_structure::GRAD_STACK1
static _THREAD grad_stack * GRAD_STACK1
Definition: gradient_structure.h:115
gradient_structure
class for things related to the gradient structures, including dimension of arrays, size of buffers, etc.
Definition: gradient_structure.h:107
cumd_norm
double cumd_norm(const double &x)
Culative normal distribution; constant objects.
Definition: cumdist.cpp:90
dvar_vector::indexmax
int indexmax() const
Definition: fvar.hpp:2292
square
double square(const double value)
Return square of value; constant object.
Definition: d3arr4.cpp:16
dvariable
Fundamental data type for reverse mode automatic differentiation.
Definition: fvar.hpp:1518
log
d3_array log(const d3_array &arr3)
Author: David Fournier Copyright (c) 2008-2012 Regents of the University of California.
Definition: d3arr2a.cpp:13
inv_cumd_norm_logistic
double inv_cumd_norm_logistic(double x, double p)
Description not yet available.
Definition: ccumdlog.cpp:39
double_and_int::x
double x
&lt; value of the variable
Definition: fvar.hpp:195
ADMB Generated on Wed Sep 7 2022 00:01:30 for ADMB Documentation by doxygen 1.8.5