NormalRV.cppGo to the documentation of this file.00001 /* ****************************************************************** ** 00002 ** OpenSees - Open System for Earthquake Engineering Simulation ** 00003 ** Pacific Earthquake Engineering Research Center ** 00004 ** ** 00005 ** ** 00006 ** (C) Copyright 2001, The Regents of the University of California ** 00007 ** All Rights Reserved. ** 00008 ** ** 00009 ** Commercial use of this program without express permission of the ** 00010 ** University of California, Berkeley, is strictly prohibited. See ** 00011 ** file 'COPYRIGHT' in main directory for information on usage and ** 00012 ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** 00013 ** ** 00014 ** Developed by: ** 00015 ** Frank McKenna (fmckenna@ce.berkeley.edu) ** 00016 ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** 00017 ** Filip C. Filippou (filippou@ce.berkeley.edu) ** 00018 ** ** 00019 ** Reliability module developed by: ** 00020 ** Terje Haukaas (haukaas@ce.berkeley.edu) ** 00021 ** Armen Der Kiureghian (adk@ce.berkeley.edu) ** 00022 ** ** 00023 ** ****************************************************************** */ 00024 00025 // $Revision: 1.7 $ 00026 // $Date: 2004/08/27 17:51:50 $ 00027 // $Source: /usr/local/cvs/OpenSees/SRC/reliability/domain/distributions/NormalRV.cpp,v $ 00028 00029 00030 // 00031 // Written by Terje Haukaas (haukaas@ce.berkeley.edu) 00032 // 00033 00034 #include <NormalRV.h> 00035 #include <math.h> 00036 #include <string.h> 00037 #include <classTags.h> 00038 #include <OPS_Globals.h> 00039 00040 NormalRV::NormalRV(int passedTag, 00041 double passedMean, 00042 double passedStdv, 00043 double passedStartValue) 00044 :RandomVariable(passedTag, RANDOM_VARIABLE_normal) 00045 { 00046 tag = passedTag ; 00047 mju = passedMean; 00048 sigma = passedStdv; 00049 startValue = passedStartValue; 00050 } 00051 NormalRV::NormalRV(int passedTag, 00052 double passedParameter1, 00053 double passedParameter2, 00054 double passedParameter3, 00055 double passedParameter4, 00056 double passedStartValue) 00057 :RandomVariable(passedTag, RANDOM_VARIABLE_normal) 00058 { 00059 tag = passedTag ; 00060 mju = passedParameter1; 00061 sigma = passedParameter2; 00062 startValue = passedStartValue; 00063 } 00064 NormalRV::NormalRV(int passedTag, 00065 double passedMean, 00066 double passedStdv) 00067 :RandomVariable(passedTag, RANDOM_VARIABLE_normal) 00068 { 00069 tag = passedTag ; 00070 mju = passedMean; 00071 sigma = passedStdv; 00072 startValue = getMean(); 00073 } 00074 NormalRV::NormalRV(int passedTag, 00075 double passedParameter1, 00076 double passedParameter2, 00077 double passedParameter3, 00078 double passedParameter4) 00079 :RandomVariable(passedTag, RANDOM_VARIABLE_normal) 00080 { 00081 tag = passedTag ; 00082 mju = passedParameter1; 00083 sigma = passedParameter2; 00084 startValue = getMean(); 00085 } 00086 00087 00088 NormalRV::~NormalRV() 00089 { 00090 } 00091 00092 00093 void 00094 NormalRV::Print(OPS_Stream &s, int flag) 00095 { 00096 } 00097 00098 00099 double 00100 NormalRV::getPDFvalue(double rvValue) 00101 { 00102 double pi = 3.14159265358979; 00103 return 1 / sqrt ( 2.0 * pi ) * exp ( - 0.5 * pow ( ( ( rvValue - mju ) / sigma ), 2.0 ) ); 00104 } 00105 00106 00107 double 00108 NormalRV::getCDFvalue(double rvValue) 00109 { 00110 double result = 0.5 + errorFunction( ((rvValue-mju)/sigma)/sqrt(2.0) )/2.0; 00111 return result; 00112 } 00113 00114 00115 double 00116 NormalRV::getInverseCDFvalue(double probValue) 00117 { 00118 if (probValue < 0.0 || probValue > 1.0) { 00119 opserr << "WARNING: Illegal probability value input to NormalRV::getInverseCDFvalue()" << endln; 00120 } 00121 double result = getMean() + getStdv() * sqrt(2.0) * inverseErrorFunction(2*probValue-1.0); 00122 return result; 00123 } 00124 00125 00126 const char * 00127 NormalRV::getType() 00128 { 00129 return "NORMAL"; 00130 } 00131 00132 00133 double 00134 NormalRV::getMean() 00135 { 00136 return mju; 00137 } 00138 00139 00140 00141 double 00142 NormalRV::getStdv() 00143 { 00144 return sigma; 00145 } 00146 00147 00148 double 00149 NormalRV::getStartValue() 00150 { 00151 return startValue; 00152 } 00153 00154 double NormalRV::getParameter1() {return mju;} 00155 double NormalRV::getParameter2() {return sigma;} 00156 double NormalRV::getParameter3() {opserr<<"No such parameter in r.v. #"<<tag<<endln; return 0.0;} 00157 double NormalRV::getParameter4() {opserr<<"No such parameter in r.v. #"<<tag<<endln; return 0.0;} 00158 00159 00160 double 00161 NormalRV::errorFunction(double x) 00162 { 00163 00164 // ErrorFunction(x) = 2/sqrt(pi) * integral from 0 to x of exp(-t^2) dt. 00165 00166 double a1,a2,a3,a4,a5; 00167 double b1,b2,b3,b4; 00168 double c1,c2,c3,c4,c5,c6,c7,c8,c9; 00169 double d1,d2,d3,d4,d5,d6,d7,d8; 00170 double p1,p2,p3,p4,p5,p6; 00171 double q1,q2,q3,q4,q5; 00172 double y,z,xnum,xden,del,result; 00173 double pi = 3.14159265358979; 00174 00175 00176 // Evaluate errorFunction for |x| <= 0.46875 00177 if ( fabs(x) <= 0.46875 ) 00178 { 00179 a1 = 3.16112374387056560e00; 00180 a2 = 1.13864154151050156e02; 00181 a3 = 3.77485237685302021e02; 00182 a4 = 3.20937758913846947e03; 00183 a5 = 1.85777706184603153e-1; 00184 b1 = 2.36012909523441209e01; 00185 b2 = 2.44024637934444173e02; 00186 b3 = 1.28261652607737228e03; 00187 b4 = 2.84423683343917062e03; 00188 y = fabs(x); 00189 z = y * y; 00190 xnum = a5 * z; 00191 xden = z; 00192 xnum = (xnum + a1) * z; 00193 xden = (xden + b1) * z; 00194 xnum = (xnum + a2) * z; 00195 xden = (xden + b2) * z; 00196 xnum = (xnum + a3) * z; 00197 xden = (xden + b3) * z; 00198 result = x * (xnum + a4) / (xden + b4); 00199 } 00200 00201 00202 // Evaluate errorFunction for 0.46875 <= |x| <= 4.0 00203 else if ( fabs(x) > 0.46875 && fabs(x) <= 4.0 ) 00204 { 00205 c1 = 5.64188496988670089e-1; 00206 c2 = 8.88314979438837594e00; 00207 c3 = 6.61191906371416295e01; 00208 c4 = 2.98635138197400131e02; 00209 c5 = 8.81952221241769090e02; 00210 c6 = 1.71204761263407058e03; 00211 c7 = 2.05107837782607147e03; 00212 c8 = 1.23033935479799725e03; 00213 c9 = 2.15311535474403846e-8; 00214 d1 = 1.57449261107098347e01; 00215 d2 = 1.17693950891312499e02; 00216 d3 = 5.37181101862009858e02; 00217 d4 = 1.62138957456669019e03; 00218 d5 = 3.29079923573345963e03; 00219 d6 = 4.36261909014324716e03; 00220 d7 = 3.43936767414372164e03; 00221 d8 = 1.23033935480374942e03; 00222 y = fabs(x); 00223 xnum = c9 * y; 00224 xden = y; 00225 xnum = (xnum + c1) * y; 00226 xden = (xden + d1) * y; 00227 xnum = (xnum + c2) * y; 00228 xden = (xden + d2) * y; 00229 xnum = (xnum + c3) * y; 00230 xden = (xden + d3) * y; 00231 xnum = (xnum + c4) * y; 00232 xden = (xden + d4) * y; 00233 xnum = (xnum + c5) * y; 00234 xden = (xden + d5) * y; 00235 xnum = (xnum + c6) * y; 00236 xden = (xden + d6) * y; 00237 xnum = (xnum + c7) * y; 00238 xden = (xden + d7) * y; 00239 result = (xnum + c8) / (xden + d8); 00240 z = floor(y*16)/16; 00241 del = (y-z) * (y+z); 00242 result = exp(-z * z) * exp(-del) * result; 00243 } 00244 00245 00246 // Evaluate erfc for |x| > 4.0 00247 else if ( fabs(x) > 4.0 ) 00248 { 00249 p1 = 3.05326634961232344e-1; 00250 p2 = 3.60344899949804439e-1; 00251 p3 = 1.25781726111229246e-1; 00252 p4 = 1.60837851487422766e-2; 00253 p5 = 6.58749161529837803e-4; 00254 p6 = 1.63153871373020978e-2; 00255 q1 = 2.56852019228982242e00; 00256 q2 = 1.87295284992346047e00; 00257 q3 = 5.27905102951428412e-1; 00258 q4 = 6.05183413124413191e-2; 00259 q5 = 2.33520497626869185e-3; 00260 y = fabs(x); 00261 z = 1 / (y * y); 00262 xnum = p6 * z; 00263 xden = z; 00264 xnum = (xnum + p1) * z; 00265 xden = (xden + q1) * z; 00266 xnum = (xnum + p2) * z; 00267 xden = (xden + q2) * z; 00268 xnum = (xnum + p3) * z; 00269 xden = (xden + q3) * z; 00270 xnum = (xnum + p4) * z; 00271 xden = (xden + q4) * z; 00272 result = z * (xnum + p5) / (xden + q5); 00273 result = (1/sqrt(pi) - result) / y; 00274 z = floor(y*16)/16; 00275 del = (y-z) * (y+z); 00276 result = exp(-z * z) * exp(-del) * result; 00277 } 00278 00279 00280 // Final computations 00281 if ( x > 0.46875 ) 00282 result = (0.5 - result) + 0.5; 00283 if ( x < -0.46875 ) 00284 result = (-0.5 + result) - 0.5; 00285 00286 return result; 00287 } 00288 00289 double 00290 NormalRV::inverseErrorFunction(double y) 00291 { 00292 double a1,a2,a3,a4; 00293 double b1,b2,b3,b4; 00294 double c1,c2,c3,c4; 00295 double d1,d2; 00296 double x,z; 00297 double pi = 3.14159265358979; 00298 00299 // Coefficients in rational approximations. 00300 a1 = 0.886226899; 00301 a2 = -1.645349621; 00302 a3 = 0.914624893; 00303 a4 = 0.140543331; 00304 b1 = -2.118377725; 00305 b2 = 1.442710462; 00306 b3 = -0.329097515; 00307 b4 = 0.012229801; 00308 c1 = -1.970840454; 00309 c2 = -1.624906493; 00310 c3 = 3.429567803; 00311 c4 = 1.641345311; 00312 d1 = 3.543889200; 00313 d2 = 1.637067800; 00314 00315 // Central range 00316 if ( fabs(y) <= 0.7 ) 00317 { 00318 z = y * y; 00319 x = y * (((a4*z+a3)*z+a2)*z+a1) / ((((b4*z+b3)*z+b2)*z+b1)*z+1); 00320 } 00321 00322 // Near end-points of range 00323 if ( y > 0.7 && y < 1 ) 00324 { 00325 z = sqrt(-log((1-y)/2)); 00326 x = (((c4*z+c3)*z+c2)*z+c1) / ((d2*z+d1)*z+1); 00327 } 00328 00329 if ( y < -0.7 && y > -1 ) 00330 { 00331 z = sqrt(-log((1+y)/2)); 00332 x = -(((c4*z+c3)*z+c2)*z+c1) / ((d2*z+d1)*z+1); 00333 } 00334 00335 00336 // Two steps of Newton-Raphson correction to full accuracy. 00337 // Without these steps, erfinv(y) would be about 3 times 00338 // faster to compute, but accurate to only about 6 digits. 00339 x = x - (errorFunction(x) - y) / (2.0/sqrt(pi) * exp(-pow(x,2.0))); 00340 x = x - (errorFunction(x) - y) / (2.0/sqrt(pi) * exp(-pow(x,2.0))); 00341 00342 // Exceptional cases not treated for now 00343 // k = find(y == -1); 00344 // if ~isempty(k), x(k) = -inf*k; end 00345 // k = find(y == 1); 00346 // if ~isempty(k), x(k) = inf*k; end 00347 // k = find(abs(y) > 1); 00348 // if ~isempty(k), x(k) = NaN; end 00349 return x; 00350 00351 } 00352 00353 00354 00355 00356 00357 00358 00359
Generated on Mon Oct 23 15:05:26 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents