Method PCMfunc

SRC/reliability/analysis/analysis/system/PCM.cpp:190–246

}

double
PCM::PCMfunc(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;
	
	rho = rhoin;
	for (i=0; i < n; i++) {
		beta(i) = allbeta(i)*modifier;
		rho(i,i) = beta(i);
	}
	
	if (n == 1)
		return uRV.getCDFvalue( beta(0) );
	else if (n == 2) {
		// check closed-form solution
		CorrelatedStandardNormal phi2(rhoin(1,0));
		double pcf = phi2.getCDF(beta(0),beta(1));
		//opserr << "pcf = " << pcf << " and PCM = " << exp(pf) << endln;
		return pcf;
	} else {
	
		// ���� FIRST CYCLE ���������������� 
		double A1 = uRV.getPDFvalue(rho(1-1,1-1))/(uRV.getCDFvalue(rho(1-1,1-1)));
		double B1 = A1*(rho(1-1,1-1) + A1);
		for (k = 2; k <= n; k++) 
			rho(k-1,1-1) = (rho(k-1,k-1) + rho(1-1,k-1)*A1)/sqrt(1 - rho(1-1,k-1)*rho(1-1,k-1)*B1);

		for (ir = 2; ir <= n - 1; ir++) {
			for (ic = ir + 1; ic <= n; ic++)
				rho(ir-1,ic-1) = (rho(ir-1,ic-1) - rho(1-1,ir-1)*rho(1-1,ic-1)*B1)/sqrt((1-rho(1-1,ir-1)*rho(1-1,ir-1)*B1)*(1 - rho(1-1,ic-1)*rho(1-1,ic-1)*B1));
		}

		// ����- OTHER CYCLES ���������������� 
		for (j = 2; j <= n - 1; j++) {
			A1 = uRV.getPDFvalue(rho(j-1,j-2))/(uRV.getCDFvalue(rho(j-1,j-2)));
			B1 = A1*(rho(j-1,j-2) + A1);
			for (k = j + 1; k <= n; k++) 
				rho(k-1,j-1) = (rho(k-1,j-2) + rho(j-1,k-1)*A1)/sqrt(1 - rho(j-1,k-1)*rho(j-1,k-1)*B1);	
			
			for (ir = j + 1; ir <= n - 1; ir++) {
				for (ic = ir + 1; ic <= n; ic++) 
					rho(ir-1,ic-1) = (rho(ir-1,ic-1) - rho(j-1,ir-1)*rho(j-1,ic-1)*B1)/sqrt((1 - rho(j-1,ir-1)*rho(j-1,ir-1)*B1)*(1 - rho(j-1,ic-1)*rho(j-1,ic-1)*B1));
			}
		}

		// ���� Calculate the product of conditional marginals 
		double pf = log(uRV.getCDFvalue(rho(1-1,1-1))); 
		for (i = 2; i<=n; i++)
			pf = pf + log(uRV.getCDFvalue(rho(i-1,i-2)));
		
		return exp(pf);
	}
}