Method IPCMfunc

SRC/reliability/analysis/analysis/system/IPCM.cpp:193–300

}

double
IPCM::IPCMfunc(const Vector &allbeta, const Matrix &rhoin, double modifier)
{
	int n = allbeta.Size();
	static NormalRV uRV(1, 0.0, 1.0);
	Vector beta(n);
	Matrix rho(n,n);
	int i,ic,ir,j,k;

	double c1, c2, r, a, b, c21, jprob;
	double orig, newg, modval;
	CorrelatedStandardNormal phi2(0.0);
	
	rho = rhoin;
	for (i=0; i < n; i++) {
		beta(i) = allbeta(i)*modifier;
		rho(i,i) = beta(i);
	}
	
	if (n == 1)
		return 1.0 - uRV.getCDFvalue( beta(0) );
	else if (n == 2) {
		// check closed-form solution
		phi2.setCorrelation(rhoin(1,0));
		double pcf = phi2.getCDF(beta(0),beta(1));
		//opserr << "pcf = " << pcf << " and PCM = " << exp(pf) << endln;
		return 1.0 - pcf;
	} else {
		
		// ���� FIRST CYCLE ���������������� 
		double pdfc1 = uRV.getPDFvalue(rho(1-1,1-1)); 
		double cdfc1 = uRV.getCDFvalue(rho(1-1,1-1)); 
		double A1 = pdfc1/cdfc1;
		double B1 = A1*(rho(1-1,1-1) + A1);
		for (k = 2; k <= n; k++) {
			c1 = -rho(1-1,1-1); 
			c2 = -rho(k-1,k-1);
			r = rho(1-1,k-1); 
			a = pdfc1/(1.0 - cdfc1); 
			b = a*(c1 + a); 
			c21 = (c2 + r*a)/sqrt(1.0 - r*r*b);
			
			// original
			orig = (1.0 - cdfc1)*uRV.getCDFvalue(c21);
			// bootstrap with binormal
			phi2.setCorrelation(r);
			newg = phi2.getCDF(c1,c2);
			
			jprob = uRV.getCDFvalue(c2) - newg;
			// might have a div by 0 problem here
			modval = 1.0 - jprob/cdfc1;
			if ( isnan(modval) ) {
				opserr << "SystemAnalysis::IPCM WARNING illegal norminv value input (nan): " 
					   << "cdfc1 = " << cdfc1 << " and c2 = " << c2 << endln;
			}
			
			rho(k-1,1-1) = uRV.getInverseCDFvalue(modval); 
		}