using System; using System.IO; // ======================================================================================================================= // Class name: TSTForm4 // // Description : General procedures to calculate the index ... // // Comments : Most of the variables are given as a parameter // in stead of properties. this is only done for // compability with the old Form Call. // New values necessary should be given as properties. // The zfunction has a predefined call and has the name // GG. ( kept from the form routine ) // The use of this function is the same as a property // thus : theForm.GG := Zfunction (with zfunction being // the zFunction in the program // // ======================================================================================================================= namespace Deltares.Stability.Calculation.Inner { // type // TGiveProgress = procedure(AProgress: Double) of object; // Onderstaande is nu input bij FORM4 // Const NM = 100; // NC = 10; // NE = 25; public class TDForm4 { public TZfunctiePointerType GG = null; public int IVC = 0; public double Variantie = 0; private double[] FABLZ; private double[] FEIVAL; private double[] FEZA; private int[] FIH; private int[] FISOR; private int[] FIUG; private int FIerror = 0; private int FMaxIterMonteCarlo = 0; private double FMinBeta = 0; private double[] FMinEX; private double[] FMinU; private bool FMonteCarloExecuted = false; private double[] FRZR; private double[] FSZA; private double[] FZ; private double FZInit = 0; private double FZMinimumFound = 0; public int MaxiterMonteCarlo { get { return FMaxIterMonteCarlo; } set { FMaxIterMonteCarlo = value; } } public bool MonteCarloExecuted { get { return FMonteCarloExecuted; } } public double[] MinEX { get { return FMinEX; } set { FMinEX = value; } } public double[] MinU { get { return FMinU; } set { FMinU = value; } } public double ZMinimumFound { get { return FZMinimumFound; } set { FZMinimumFound = value; } } public double MinBeta { get { return FMinBeta; } set { FMinBeta = value; } } public int Ierror { get { return FIerror; } } public void FORM4(int NM, int[] IV, double[] Ex, double[] Sx, ref double[][] Vp, ref double[][] Cov, int NC, double[][] Eivec, int NE, ref double stuurpar, ref double Beta, ref double[] X, ref double[] U, ref double[] zes, ref int ierror, ref int Izs) { // // double FX = 0; double DFX = 0; double[] EXA; double[] SXA; double[] ABL; double[] EZ; double[] SZ; int iz = 0; int n1 = 0; double RMAX = 0; double r1 = 0; double r2 = 0; double theta = 0; double dtheta = 0; double fm = 0; int i = 0; int j = 0; int ij = 0; int k = 0; int l = 0; int ls = 0; int il = 0; int ks = 0; double EZS = 0; double SZS = 0; double SZJ = 0; double FD = 0; double FV = 0; double FDI = 0; double FVI = 0; double XV = 0; double XVI = 0; double XA = 0; double ZST = 0; double BET = 0; double FXOld = 0; double LocStuurpar = 0; double FPower = 0; double SABL = 0; int KOMP = 0; int KOMPI = 0; int IRO = 0; int IROT = 0; int IRHO = 0; // ZES: Tarr3; bool ready = false; bool finish = false; FPower = 2; DFX = 1; FXOld = 1; FMonteCarloExecuted = false; if ((NC == 1)) { IRHO = 0; } else { IRHO = 1; } InitLocalArrays(NM); EXA = new double[NM]; SXA = new double[NM]; ABL = new double[NM]; EZ = new double[NM]; SZ = new double[NM]; finish = false; iz = 0; // IZS = IMPORTANT VARIABLE DETERMINING THE NO. OF ITERAIONS if (Izs == 0) { Izs = 20; } n1 = NM + 1; ierror = 0; // PRUEFUNG DES STOCHASTISCHEN MODELLS YPRUEF(IRHO, NM, NC, IV, Vp, Cov, Sx, Ex, ref ierror); if ((ierror <= 0)) { KOMP = 0; FZ[n1 - 1] = 0.0; FEZA[n1 - 1] = 0.0; FSZA[n1 - 1] = 0.0; if ((Beta <= FormCalculation.CEMIN2)) { for (i = 1; i <= NM; i ++) { X[i - 1] = Ex[i - 1]; } } for (i = 1; i <= NM; i ++) { ABL[i - 1] = 1.0; } IRO = IRHO; if ((NM == 1)) { IRO = 0; } // VERTEILUNGSPARAMETER YTXZUY(NM, IV, Ex, Sx, ref EXA, ref SXA, ref Vp, X, ref Cov, NC, IRO); // (* // { Make covariance positive definite, if not so } // if (NC > 1) then // begin // If (Not IsMatrixPositiveDefinite(COV) ) then // MakeMatrixPositiveDefinite(COV); // end; // *) for (i = 1; i <= NM; i ++) { FEIVAL[i - 1] = SXA[i - 1]; } if ((IRO != 0)) { // ERMITTLUNG VON EIGENWERTEN UND EIGENVEKTOREN BEI KORRELIERTEN VARIABLEN if ((NE < NM)) { ierror = 143000; Foutafhandeling(IRO, NM, IV, Sx, EXA, SXA, ref Cov, NC, iz, Izs, ref ierror); finish = true; } if (!finish) { Yacobi(1, NM, ref Cov, NC, ref FEIVAL, NM, ref Eivec, NE, ref IROT, ref ierror, ref EZ, ref SZ); RMAX = FEIVAL[0]; for (i = 2; i <= NM; i ++) { if ((RMAX < FEIVAL[i - 1])) { RMAX = FEIVAL[i - 1]; } for (j = 1; j <= i; j ++) { Cov[j - 1][i - 1] = Cov[i - 1][j - 1]; } } for (i = 1; i <= NM; i ++) { if ((FEIVAL[i - 1] < 0.0)) { if ((FEIVAL[i - 1]/RMAX >= -FormCalculation.COGE14)) { FEIVAL[i - 1] = 0.0; } else { ierror = 145000; Foutafhandeling(IRO, NM, IV, Sx, EXA, SXA, ref Cov, NC, iz, Izs, ref ierror); finish = true; } if (!finish) { FEIVAL[i - 1] = Math.Sqrt(Convert.ToDouble(FEIVAL[i - 1])); } } } } } if (!finish) { // VARIABLENREIHENFOLGE AUFSTEIGend NACH WIEDERHOLUNGSZAHLEN YSORT(NM, ref Vp, ref FISOR); // VERHAELTNISSE UND GRUPPEN DER WIEDERHOLUNGSZAHLEN r2 = Vp[4][FISOR[0] - 1]; FRZR[0] = r2; FIUG[0] = 1; FIUG[n1 - 1] = n1; if ((NM > 1)) { for (i = 2; i <= NM; i ++) { // R1 := VP[5,FISOR[i]]; r1 = Vp[4][FISOR[i - 1] - 1]; FRZR[i - 1] = r1/r2; r2 = r1; FIUG[i - 1] = 0; if ((Math.Abs(FRZR[i - 1] - 1.0) > FormCalculation.CFIX4)) { FIUG[i - 1] = i; } } } // Call the external procedure, provided by the user GG(NM, Ex, theta, ref ierror, ref fm); FZInit = fm; if ((ierror > 0)) { Foutafhandeling(IRO, NM, IV, Sx, Ex, SXA, ref Cov, NC, iz, Izs, ref ierror); finish = true; } if (!finish) { fm = Math.Abs(fm); if ((fm < 1.0)) { fm = 1.0; } XV = 0.0; // CopyExToMinEx // Begin OF ITERATION LocStuurpar = stuurpar; do { iz = iz + 1; // STCalcProgressDlg.UpdateProgressLiftProb(true,Iz); // VI := V1**(1./FLOAT(IZ)) // LocStuurpar := LocStuurpar * FPower; // (* Smallstep := false; *) // eventuele verklein loop // repeat // BERECHNUNG DER PARTIELLEN ABLEITUNGEN ierror = 1; GG(NM, X, theta, ref ierror, ref FX); // WriteLn(LFErr, ' x0 x1 en fx ', x[0]: 12: 5, X[1]: 12: 5, Fx: 12: 5); if ((ierror > 0)) { Foutafhandeling(IRO, NM, IV, Sx, Ex, SXA, ref Cov, NC, iz, Izs, ref ierror); finish = true; } if (!finish) { DFX = Math.Abs(FX)/fm; YGGS(NM, X, IV, ref theta, 1, FX, SZ, ref dtheta, ref ABL, ref Cov, NC, ref ierror); if ((ierror > 0)) { Foutafhandeling(IRO, NM, IV, Sx, EXA, SXA, ref Cov, NC, iz, Izs, ref ierror); finish = true; } if (!finish) { SABL = 0.0; for (i = 1; i <= NM; i ++) { SABL = SABL + ABL[i - 1]*ABL[i - 1]; // TRANSFORMATION DER LOGNORMALEN BV UND IHRER ABLEITUNG // IN EINE NORMALVERTEILTE BV MIT ENTSPRECHendER ABLEITUNG if (((IV[i - 1] == 3) || (IV[i - 1] == 10))) { ABL[i - 1] = ABL[i - 1]*(X[i - 1] - Vp[2][i - 1]); X[i - 1] = Math.Log(Convert.ToDouble(X[i - 1] - Vp[2][i - 1])); } FABLZ[i - 1] = ABL[i - 1]; } if ((SABL < FormCalculation.CEMIN3)) { ierror = 1; finish = true; } if (!finish) { if ((IRO > 0)) { Yatmuv(ref Eivec, ABL, FABLZ, NE, 100, NM); } i = n1; j = i; do { // BEGIN OF COMPOSITION i = i - 1; ready = (i == 0); if (!ready) { // ERMITTLUNG DES ANFANGS- UND DES endINDEXES // EINER GRUPPE VON VARIABLEN GLEICHER WIEDERHOLUNGSZAHL j = j - 1; while ((FIUG[j - 1] != j) && (j > 0)) { j -= 1; } ij = i + j; EZS = FEZA[i]; SZS = FSZA[i]; for (k = j; k <= i; k ++) { l = ij - k; ls = FISOR[l - 1]; il = IV[ls - 1]; if (((il != 0) && (il != 2) && (il != 3) && (il != 10))) { YPARTR(NM, ls, IV[ls - 1], X[ls - 1], Vp, ref EXA[ls - 1], ref SXA[ls - 1], ref FD, ref FV); FEIVAL[ls - 1] = SXA[ls - 1]; if ((IRO > 0)) { Cov[ls - 1][ls - 1] = SXA[ls - 1]*SXA[ls - 1]; } } FZ[l - 1] = FZ[l + 1 - 1] + ABL[ls - 1]*X[ls - 1]; EZS = EZS + ABL[ls - 1]*EXA[ls - 1]; // UNABHAENGIGE VARIABLEN if ((IRO == 0)) { SZS = SZS + Math.Pow(ABL[ls - 1]*SXA[ls - 1], 2.0); } } if ((IRO > 0)) { // ABHAENGIGE VARIABLEN // De compiler schreeuwt "variable 'k' may be undefined after loop" // Ik weet niet waarom, maar kan geen kwaad?. 19990928,HenkJan. for (k = j; k <= i; k ++) { ks = FISOR[k - 1]; SZ[k - 1] = 0.0; for (l = j; l <= i; l ++) { ls = FISOR[l - 1]; SZ[k - 1] = SZ[k - 1] + ABL[ls - 1]*Cov[ls - 1][ks - 1]; } SZS = SZS + SZ[k - 1]*ABL[ks - 1]; } } SZJ = Math.Sqrt(SZS); // BERECHNUNG VON ALPHA U BEZOGEN AUF DIE BISHER ERFASSTE(N) // GRUPPE(N) for (k = j; k <= i; k ++) { ks = FISOR[k - 1]; U[ks - 1] = 0.0; if ((SZJ > FormCalculation.COGE14)) { U[ks - 1] = FEIVAL[ks - 1]*FABLZ[ks - 1]/SZJ; } } i = j; EZ[j - 1] = EZS; SZ[j - 1] = SZS; FEZA[j - 1] = EZS; FSZA[j - 1] = SZS; if (((Math.Abs(1.0 - FRZR[j - 1]) >= FormCalculation.COGE8) && (SZJ >= FormCalculation.COGE8))) { // BERUECKSICHTIGUNG DER ZUFAELLIGEN FOLGE YNORM(FZ[j - 1], EZ[j - 1], SZJ, ref FDI, ref FVI, ref KOMPI); YBINEX(2, FZ[j - 1], FRZR[j - 1], 1.0, FDI, FVI, KOMPI, ref FD, ref FV, ref KOMP); YPARTR(NM, j, 0, FZ[j - 1], Vp, ref FEZA[j - 1], ref FSZA[j - 1], ref FD, ref FV); if ((KOMP == 1)) { FEZA[j - 1] = 2.0*FZ[j - 1] - FEZA[j - 1]; } FSZA[j - 1] = FSZA[j - 1]*FSZA[j - 1]; } } // end OF COMPOSITION and loop i } while (!(ready)); FZ[0] = FZ[0] - FX; Beta = (FEZA[0] - FZ[0])/Math.Sqrt(Convert.ToDouble(FSZA[0])); XV = 0.0; ZST = Beta; do { // BEGIN OF DECOMPOSITION i = i + 1; ready = (i > NM); if (!ready) { j = i; // ERMITTLUNG DES ANFANGS- UND endINDEXES // EINER GRUPPE VON VARIABLEN GLEICHER WIEDERHOLUNGSZAHL j = j + 1; while ((FIUG[j - 1] != j)) { j ++; } j = j - 1; BET = 0.0; if ((SZ[i - 1] > FormCalculation.COGE14)) { BET = (EZ[i - 1] - FZ[i - 1])/SZ[i - 1]; } FZ[j + 1 - 1] = FEZA[j + 1 - 1] - BET*FSZA[j + 1 - 1]; for (k = i; k <= j; k ++) { ks = FISOR[k - 1]; XA = X[ks - 1]; if ((IRO == 0)) { // UNABHAENGIGE VARIABLEN X[ks - 1] = ABL[ks - 1]*SXA[ks - 1]*SXA[ks - 1]; } else { // ABHAENGIGE VARIABLEN X[ks - 1] = 0.0; for (l = i; l <= j; l ++) { ls = FISOR[l - 1]; X[ks - 1] = X[ks - 1] + Cov[ks - 1][ls - 1]*ABL[ls - 1]; } } // BESTIMMUNG DER NEUEN KOORDINATEN IM U-RAUM UND IM RAUM // DER BASISVARIABLEN X U[ks - 1] = -U[ks - 1]*ZST; X[ks - 1] = EXA[ks - 1] - BET*X[ks - 1]; XVI = X[ks - 1] - XA; X[ks - 1] = (1.0 - LocStuurpar)*XA + LocStuurpar*X[ks - 1]; if ((IV[ks - 1] == 3) || (IV[ks - 1] == 10)) { X[ks - 1] = Math.Exp(Convert.ToDouble(X[ks - 1])) + Vp[2][ks - 1]; } if ((SXA[ks - 1] > 1.0)) { XVI = XVI/SXA[ks - 1]; } XV = XV + XVI*XVI; } // REDUKTION VON ZST GEMAESS DEM GEWICHT DER RESTGRUPPE(N) XA = ZST; ZST = 0.0; if ((SZ[i - 1] > FormCalculation.COGE14)) { ZST = XA*Math.Sqrt(Convert.ToDouble(FSZA[j]/SZ[i - 1])); } i = j; } } while (!(ready)); // end OF DECOMPOSITION XV = Math.Sqrt(XV)/NM; // Print the last three iterations // if (IZ > IZS - 3) then // begin // WriteLn(lferr, ' Iteration ', IZ: 5, ' Array X ='); // for i := 1 to Nm do // begin // WriteLn(lferr, ' ', x[i - 1]: 15: 3); // if i mod 5 = 0 then WriteLn(lferr); // end; // end; } } } // Keep track of minimal dfx for use if iteration doe not succeed if ((iz == 1)) { FZMinimumFound = DFX; CopyMinimalResults(Beta, U, X); } else { if ((Math.Abs(FZMinimumFound) > Math.Abs(DFX))) { FZMinimumFound = DFX; CopyMinimalResults(Beta, U, X); } } if ((iz > 1) && (Math.Abs(DFX) > Math.Abs(FXOld))) { LocStuurpar = LocStuurpar/FPower; } else { LocStuurpar = LocStuurpar*0.7*FPower; // max in step removed LocStuurpar = Math.Min(LocStuurpar, 2.0); } // until (iz = 1) or ( abs(DFX) < abs(FXOld) * 1.1 ) or (LocStuurpar < 0.00001); // IMPORTANT CONVERGENCE CRITERIA FOR 'FORM4' ! FXOld = DFX; // GiveProgress(IZ / IZS); } while (!((iz >= Izs) || ((DFX <= FormCalculation.CFIX4) && (XV <= FormCalculation.CFIX2)))); // Test for minimum z value // or finish or smallstep or STCalcProgressDlg.Aborted; // End of iteration // if STCalcProgressDlg.Aborted then finish := true; // (* If (not finish ) then // {indien niet aan het afbreek criterium wordt voldaan // en er geen fout is dan een monte carlo simulatie } // if (( Iz >= Izs ) or Smallstep ) then // begin // idum := -13; // FMonteCarloExecuted := true; // MonteCarlo(Nm,cov,Ex,Sx,X,U,Zes,IV,Beta,LFErr); // end; *) // (* procedure FORM4(NM: Integer; IV: TDynamic1DIntArr; // Ex,Sx : TDynamic1DDoubleArr; var Vp,Cov : TDynamic2DDoubleArr; // NC : Integer; Eivec : TDynamic2DDoubleArr; // NE : Integer; Var v1,Beta : Double; // var X,U : TDynamic1DDoubleArr; var zes : Tarr3; var ierror,Izs : Integer);*) // if( not finish) and ( not FMonteCarloExecuted) then // begin // STandARDAUSDRUCK GG(NM, X, theta, ref ierror, ref FX); if ((ierror > 0) || (iz >= Izs)) { Foutafhandeling(IRO, NM, IV, Sx, Ex, SXA, ref Cov, NC, iz, Izs, ref ierror); } // WriteLn(lferr, ' betrouwbaarheidsindex =', BETA: 8: 3, ' met bijbehorende faalkans', // YPHIN(-BETA): 12); zes[0] = FZ[0]; zes[1] = FEZA[0]; zes[2] = Math.Sqrt(Convert.ToDouble(FSZA[0])); } // end; } } // geef geheugen van dynamische arrays weer vrij FreeLocalArrays(); EXA = null; SXA = null; ABL = null; EZ = null; SZ = null; } public double NormalDistribution(double Argm) { double result; double Value; const double Cnst = 2.506628274631; // sqrt(2*Pi) double Arg2; double ArgA; int Term; int Sign; Value = 0; Arg2 = Argm*Argm; ArgA = Math.Abs(Argm); Term = 15; // supplies accuracy of 4.7(-10) in X=2.6 Sign = 2*(Term%2) - 1; if (ArgA > 2.6) { do { Value = Term/(ArgA + Value); Term -= 1; } while (!(Term == 0)); Value = Math.Exp(-Arg2/2)/Cnst/(ArgA + Value); if (Argm > 0) { Value = 1 - Value; } } else { do { Value = Arg2/(2 + (1 + Sign*Value)/Term); Term -= 1; Sign = -Sign; } while (!(Term == 0)); Value = 0.5 + Math.Exp(-Arg2/2)/Cnst*Argm/(1 - Value); } result = Value; return result; } public void KNOT(ref double Average, ref double Stdv, double AMU, double ZIGMA, double A) { int Komp = 0; double Fd = 0; double Fv = 0; double Phi = 0; double GedExp = 0; double Ex = 0; double Ex2 = 0; YNORM(A, AMU, ZIGMA, ref Fd, ref Fv, ref Komp); if ((Komp == 1)) { Phi = Fv; } else { Phi = 1.0 - Fv; } GedExp = Math.Exp(-0.5*(A - AMU)*(A - AMU)/(ZIGMA*ZIGMA)); if ((Math.Abs(Phi) < 1.0E-10)) { Phi = 1.0E-10; } Ex = (ZIGMA*GedExp/Math.Sqrt(2.0*Math.PI) + AMU*Phi)/Phi; Ex2 = ((ZIGMA*(A + AMU)/Math.Sqrt(2.0*Math.PI))*GedExp + (ZIGMA*ZIGMA + AMU*AMU)*Phi)/Phi; Average = Ex; Stdv = Math.Sqrt(Ex2 - Ex*Ex); } public void MonteCarlo(int Nm, double[] Cov, double[] Ex, double[] Sx, ref double[] X, ref double[] U, ref double[] zes, int[] IV, ref double Beta, ref Stream lferr) { // (* // var // i, j: Integer; // Theta: Double; // Z: Double; // NResult: Integer; // PFaal: Double; // Loop: Integer; // Kdfz: Double; // kdfmax: Double; // Muz: Double; // SigZ: Double; // ZSom: Double; // ZSom2: Double; // Ontwerp: TDynamic1DDoubleArr; // Valdone: TDynamic1DIntArr; // Correl: Boolean; // *) // (* // { initialisatie } // ZSom := 0; // ZSom2 := 0; // PFaal := 1; // NResult := 0; // KdfMax := 0; // Loop := 0; // Theta := 0; // SetLength(Ontwerp, Nm); // SetLength(Valdone, nm); // repeat // Inc(loop); // for i := Low(Valdone) to High(Valdone) do // valdone[i] := 0; // // for i := Low(x) to High(x) do // if (valdone[i] = 0) then // begin // correl := False; // j := i; // repeat // if (i < High(x)) then // begin // Inc(j); // if (High(cov) >= j) then // if (High(cov[j]) >= i) then // correl := (Abs(cov[j, i]) > 1.E-8); // if correl then // begin // GetCorrelRealisation(Iv[j], iv[i], Ex[j], Ex[i], Sx[j], Sx[i], // X[j], X[i], Cov[j, i]); // Valdone[i] := 1; // Valdone[j] := 1; // end; // end; // until (j = High(x)) or (Correl); // // if (not Correl) then // begin // GetSingleRealisation(Iv[i], Ex[i], Sx[i], x[i]); // valdone[i] := 1; // end; // end; // // GG(Nm, X, THETA, FIERror, Z); // ZSom := Zsom + Z; { voor berekenen van verwachting} // ZSom2 := ZSom2 + Z * Z; { voor berekenen stdv } // // if (Z < 0) then // begin // Inc(NResult); // { Bepaal de gezamelijke kansdichtheid van de variabelen } // KdfZ := CombinedKdf(Nm, Ex, Sx, X, Iv); // {Bewaar waarden van grootste kdf } // if Kdfz > KdfMax then // begin // { vul de ontwerpwaarden Cases en bewaar de nieuwe kdf} // for i := Low(X) to High(X) do // Ontwerp[i] := X[i]; // kdfMax := kdfz; // end; // end; // if NResult > 1 then // PFaal := NResult / Loop; // // until (PFaal < 1.E-7) or (Loop = FMaxiterMonteCarlo); // // begin // if (NResult = 0) then // begin // { beta afleiden uit z } // Muz := Zsom / Loop; // SigZ := Sqrt((Zsom2 - ZSom * Zsom / loop) / (Loop - 1)); // PFaal := NormalDistribution(-Muz / SigZ); // ZES[1] := PFaal; // ZES[2] := Muz; // ZES[3] := SigZ; // Beta := -NormalDistrInverse(Pfaal); // end // else // begin // ZES[1] := pFaal; // ZES[2] := 0; // ZES[3] := 0; // // { van de faalkans een Beta maken } // Beta := -NormalDistrInverse(Pfaal); // for i := Low(X) to High(X) do // begin // {de ontwerp waarden in X zetten } // X[i] := Ontwerp[i]; // {Nu de alpha's hier de U's invullen } // U[i] := getAlfa(IV[i], Ex[i], Sx[i], X[i], Beta); // end; // end; // end; // { weggooien tijdelijke pointer } // Ontwerp := nil; // ValDone := nil; // *) } // ======================================================================================================================= // Description: InitLocalArrays(NM) // Initialises the lo0cal dynamic arrays; // // Date ID Modification // Created // ======================================================================================================================= protected void InitLocalArrays(int NM) { FIH = new int[NM]; FZ = new double[NM + 1]; FEZA = new double[NM + 1]; FSZA = new double[NM + 1]; FRZR = new double[NM]; FEIVAL = new double[NM]; FABLZ = new double[NM]; FISOR = new int[NM]; FIUG = new int[NM + 1]; } // ======================================================================================================================= // Description: FreeLocalArrays // Realeases memory of the local dynamic arrays; // // Date ID Modification // Created // ======================================================================================================================= protected void FreeLocalArrays() { FIH = null; FZ = null; FEZA = null; FSZA = null; FRZR = null; FEIVAL = null; FABLZ = null; FISOR = null; FIUG = null; } // ======================================================================================================================= // Date ID Modification // 2008-10-13 Best Created // ======================================================================================================================= protected void CopyMinimalResults(double ABeta, double[] AU, double[] AX) { int LNDim; int i; FMinBeta = ABeta; LNDim = AX.Length; FMinEX = new double[LNDim]; FMinU = new double[LNDim]; for (i = 0; i < LNDim; i ++) { MinEX[i] = AX[i]; MinU[i] = AU[i]; } } // procedure FORM4 (V1,BETA : double; error,IZS: Integer); // ======================================================================================================================= // Description : // Algoritm to calculate reliability index Beta, given the mechanical // model (by meand of procedure GG) and the stochastic model // (by means of the variables in the parameter list ). // // BEI BEDARF AN HOEHERER DIMENSION ALS N = 100 MUESSEN DIE STATISCHEN // FELDVEREINBARUNGEN UND DIE VARIABLE NST (ZIEMLICH AM ANFANG DER // WERTZUWEISUNGEN) AUF DEN NEUEN WERT N BZW N+1 GEBRACHT WERDEN. // // DER COMMONBLOCK CFORM KANN DURCH VEREINBARUNG IM AUFRUFendEN // PROGRAMM DAZU BENUTZT WERDEN, INFORMATION UEBER AKTUELLE WERTE // DER DARIN VEREINBARTEN FELDER ZU ERLANGEN, INSBESONDERE DIE // PARTIELLEN ABLEITUNGEN ABL[i] = D(GG)/D(X[i]) // // Date ID Modification // Created // 2008-01-16 Best Maximum in locstuurpar removed, Increase in stuurpar is smaller // 2009-01-06 Best MakeMatrixPositiveDefinite only if there is correlation // ======================================================================================================================= // (**) // Onderstaande routine is gecreeerd ivm de vele jumps naar regel 600 in de FORM4 // routine. protected void Foutafhandeling(int IRO, int nm, int[] IV, double[] SX, double[] EXA, double[] SXA, ref double[][] COV, int NC, int iz, int izs, ref int Ierror) { if ((IRO > 0)) { YTYZUX(nm, IV, SX, EXA, SXA, ref COV, NC); } // FEHLERBEHANDLUNG if (((iz >= izs) || (Ierror != 0))) { if ((Ierror == 0)) { Ierror = 300000; } else { if ((Ierror < 1000)) { Ierror = 200000 + Ierror; } } } } // ======================================================================================================================= // Description: Gamma function to replace the NAG lib routine // // Date ID Modification // 1992-12 best Created // ======================================================================================================================= protected double Gamma(double ZZ) { double result; double P; double Z; double X; double G; P = 1.0; Z = ZZ; while ((Z >= 2)) { Z = Z - 1.0; P = P*Z; } if ((Z < 1.0)) { P = P/Z; Z = Z + 1; } // Nu is 1 < Z < 2 // volgens Abramovitz geldt dan : X = Z - 1; G = (((((((3.5868343E-2*X - 1.93527818E-1)*X + 4.82199394E-1)*X - 7.56704078E-1)*X + 9.18206857E-1)*X - 8.97056937E-1)*X + 9.88205891E-1)*X - 5.77191652E-1)*X + 1.0; result = P*G; return result; } // procedure Kopf(IRHO: Integer; NM: Integer); // ======================================================================================================================= // Description: Yacobi transformatie // JACOBI-VERFAHREN // IVEC = 0: NUR EIGENWERTE (EW) // = 1: AUCH EIGENVEKTOREN (EV) // BENOETIGT WIRD NUR DIE OBERE HAELFTE EINSCHLIESSLICH // DIAGONALE DER SYMM. MATRIX. DIESER TEIL DER MATRIX IST // ANSCHLIESSend ZERSTOERT. // DIAGONALE UND UNTERER TEIL BLEIBEN UNVERAendERT // // Date ID Modification // Created // ======================================================================================================================= protected void Yacobi(int IVEC, int N, ref double[][] A, int NA, ref double[] EW, int __New, ref double[][] EV, int NEV, ref int IROT, ref int IERror, ref double[] B, ref double[] Z) { int i; int im1; int j; int k; int l; int lp1; int n1; int K1; int Km1; int Lm1; double SM; double tresh; double g; double h; double t; double theta; double c; double s; double tau; bool Ready; IERror = 0; if (((N > NA) || (N > __New))) { IERror = 1; // regel 1 } else { // regel 2 if ((IVEC != 0)) { if ((N > NEV)) { IERror = 1; // regel 1 } else { EV[0][0] = 1.0; for (i = 2; i <= N; i ++) { im1 = i - 1; for (j = 1; j <= im1; j ++) { EV[i - 1][j - 1] = 0.0; EV[j - 1][i - 1] = 0.0; } EV[i - 1][i - 1] = 1.0; } } } // regel 10 if ((IERror == 0)) { for (i = 1; i <= N; i ++) { B[i - 1] = A[i - 1][i - 1]; EW[i - 1] = B[i - 1]; Z[i - 1] = 0.0; } IROT = 0; i = 0; do { // regel 14 i = i + 1; SM = 0.0; n1 = N - 1; for (k = 1; k <= n1; k ++) { K1 = k + 1; for (l = K1; l <= N; l ++) { SM = SM + Math.Abs(A[k - 1][l - 1]); } } Ready = (SM < FormCalculation.CEMIN); // SM return. if (!Ready) { tresh = 0.0; if ((i < 4)) { tresh = 0.2*SM/(N*N); } for (k = 1; k <= n1; k ++) { K1 = k + 1; for (l = K1; l <= N; l ++) { g = 100.0*Math.Abs(A[k - 1][l - 1]); if (((i >= 5) && ((Math.Abs(EW[k - 1]) + g) == Math.Abs(EW[k - 1])) && ((Math.Abs(EW[l - 1]) + g) == Math.Abs(EW[l - 1])))) { A[k - 1][l - 1] = 0; // -> 70 } else { if ((Math.Abs(A[k - 1][l - 1]) > tresh)) { h = EW[l - 1] - EW[k - 1]; if (((Math.Abs(h) + g) == Math.Abs(h))) { t = A[k - 1][l - 1]/h; } else { theta = 0.5*h/A[k - 1][l - 1]; t = 1.0/(Math.Abs(theta) + Math.Sqrt(1.0 + theta*theta)); if ((theta < 0)) { t = -t; } } // regel 24 c = 1.0/Math.Sqrt(1.0 + t*t); s = t*c; tau = s/(1.0 + c); h = t*A[k - 1][l - 1]; Z[k - 1] = Z[k - 1] - h; Z[l - 1] = Z[l - 1] + h; EW[k - 1] = EW[k - 1] - h; EW[l - 1] = EW[l - 1] + h; A[k - 1][l - 1] = 0.0; Km1 = k - 1; if ((Km1 >= 1)) { for (j = 1; j <= Km1; j ++) { g = A[j - 1][k - 1]; h = A[j - 1][l - 1]; A[j - 1][k - 1] = g - s*(h + g*tau); A[j - 1][l - 1] = h + s*(g - h*tau); } } Lm1 = l - 1; // regel 31 if ((Lm1 >= K1)) { for (j = K1; j <= Lm1; j ++) { g = A[k - 1][j - 1]; h = A[j - 1][l - 1]; A[k - 1][j - 1] = g - s*(h + g*tau); A[j - 1][l - 1] = h + s*(g - h*tau); } } // regel 41 lp1 = l + 1; if ((N >= lp1)) { for (j = lp1; j <= N; j ++) { g = A[k - 1][j - 1]; h = A[l - 1][j - 1]; A[k - 1][j - 1] = g - s*(h + g*tau); A[l - 1][j - 1] = h + s*(g - h*tau); } } // regel 51 if ((IVEC != 0)) { for (j = 1; j <= N; j ++) { g = EV[j - 1][k - 1]; h = EV[j - 1][l - 1]; EV[j - 1][k - 1] = g - s*(h + g*tau); EV[j - 1][l - 1] = h + s*(g - h*tau); } } // regel 69 IROT = IROT + 1; } // IF ABS } // IF } } for (k = 1; k <= N; k ++) { B[k - 1] = B[k - 1] + Z[k - 1]; EW[k - 1] = B[k - 1]; Z[k - 1] = 0.0; } Ready = (i >= 50); } } while (!(Ready)); } } } // ************************************************************ // ************************************************************ // ************************************************************ // NB: A is een matrix die wordt meegegeven als array. // (Is het mogelijk in Pascal een 2D-array van ongedefinieerde grootte // mee te geven als argument?) // Dit geeft dus de omzetting A[j-1,i-1] -> A[j-1*NA+i-1]; protected void Yatmuv(ref double[][] A, double[] X, double[] Y, int NA, int NX, int NY) { // C MULTIPLIKATION Y = A(TRANSPONIERT) * X // DIMENSION A(NA,NA), X(NX), Y(NY) int i; int j; int N; N = NA; if ((NX < N)) { N = NX; } if ((NY < N)) { N = NY; } for (i = 1; i <= N; i ++) { Y[i - 1] = 0.0; for (j = 1; j <= N; j ++) { Y[i - 1] = Y[i - 1] + A[j - 1][i - 1]*X[j - 1]; } } } // ****************************************************************** // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected void YBINEX(int IT, double X, double RN, double P, double FDI, double FVI, int KOMPI, ref double FD, ref double FV, ref int KOMP) { // (*C // C BINOMIALFOLGE VON RN UNABH�NGIGEN ZUFALLLSVARIABLEN X // C MIT AUFTRETENSWAHRSCHEINLICHKEIT P. // C X NIMMT MIT WAHRSCHEINLICHKEIT (1-P) DEN WERT 0 AN UND MIT // C WAHRSCHEINLICHKEIT P EINEN WERT AUS EINER GEGEBENEN VERTEILUNG F. // C HIERBEI GELTEN FVI = F(X) UND FDI = D(X) (DICHTE). // C // C IT STEUERT, OB DIE VERTEILUNG DES MAXIMUMS ODER DIE DES MINIMUMS // C BERECHNET WIRD. // C // C 0 < RN // C 0 < P <= 1 // C // C // C IT = 1 VERTEILUNG DES MAXIMUMS // C // C (P*FVI) ** RN X < 0 // C FV = // C (1-P+P*FVI) ** RN X >= 0 // C // C RN * (P*FVI) ** (RN-1) * P * FDI X < 0 // C FD = // C RN * (1-P+P*FVI) ** (RN-1) * P * FDI X > 0 // C // C // C IT = 2 VERTEILUNG DES MINIMUMS // C // C 1 - (1-P*FVI) ** RN X < 0 // C FV = // C 1 - (P-P*FVI) ** RN X >= 0 // C // C RN * (1-P*FVI) ** (RN-1) * P * FDI X < 0 // C FD = // C RN * (P-P*FVI) ** (RN-1) * P * FDI X > 0 // C // C EINGANG: IT, X, RN, P, FDI, FVI, KOMPI // C // C AUSGANG: FD, FV, KOMP // C // C BENOETIGTE UP: YREIHE, YNPOT // C // // C // C VERTEILUNG // (**) double f; double g; double h; double ik; double pn; FD = 0.0; FV = 0.0; KOMP = 0; f = P*FVI; if ((IT == 1)) { if ((X < 0)) { if ((KOMPI == 1)) { h = P - f; } else { h = f; } // -> 70 } else { if ((KOMPI == 1)) { h = 1.0 - f; } else { h = 1.0 - P + f; } // -> 70 } } else { // IT = 2 KOMP = 1; if ((X < 0)) { if ((KOMPI == 1)) { h = 1.0 - P + f; } else { h = 1.0 - f; } // -> 70 } else { if ((KOMPI == 1)) { h = f; } else { h = P - f; } // -> 70 } } // 70 if ((h != 0.0)) { h = RN*Math.Log(h); if ((h <= -0.1)) { if ((h > -FormCalculation.CEXP2)) { FV = Math.Exp(h); } // -> 200 } else { if ((IT == 1)) { KOMP = 1; if ((KOMPI == 1)) { if ((X < 0)) { g = YNPOT(RN, -FVI); h = RN*Math.Log(P); pn = 0.0; if ((h > -FormCalculation.CEXP2)) { pn = Math.Exp(h); } if ((pn >= 0.9)) { FV = 1.0 - pn + pn*g; // -> 200 } else { KOMP = 0; FV = pn - pn*g; } // -> 200 } else { // x >= 0 FV = -YNPOT(RN, -f); } // -> 200 } else { FV = -YREIHE(h); } // -> 200 } else { // IT = 2 KOMP = 0; if ((KOMPI == 0)) { if ((X >= 0.0)) { g = -YNPOT(RN, -FVI); h = RN*Math.Log(P); pn = 0.0; if ((h > -FormCalculation.CEXP2)) { pn = Math.Exp(h); } if ((pn >= 0.9)) { FV = 1.0 - pn + pn*g; // -> 200 } else { KOMP = 0; FV = pn - pn*g; } // -> 200 } else { // x < 0 FV = -YNPOT(RN, -f); } // -> 200 } else { FV = -YREIHE(h); } // -> 200 } } } // 200 // C // C DICHTE FD = FDI*P*RN; f = FV; ik = IT + KOMP; if ((ik == 2)) { f = 1.0 - f; } FD = FD*f; if ((FD != 0.0)) { ik = IT + KOMPI; f = P*FVI; if ((ik == 2)) { h = 1.0; if (((IT == 1) && (X < 0.0)) || ((IT == 2) && (X >= 0))) { h = P; } f = h - f; } else { if (!(((ik == 1) && (X < 0.0)) || ((ik == 3) && (X >= 0)))) { f = 1.0 - P + f; } } FD = FD/f; } } // ======================================================================================================================= // Description : // // Date ID Modification // Created // ======================================================================================================================= protected void YCHI2(double X0, double RN, ref double FD, ref double FV, ref int KOMP) { // (*C // C ZENTRALE C H I Q U A D R A T V E R T E I L U N G // C // C EINGANG: // C X0 WERT DER ZUFALLSVARIABLEN X // C RN FREIHEITSGRAD VON X // C // C AUSGANG: // C FD WERT DER DICHTEFUNKTION VON X AN DER STELLE X0 // C FV P (X < X0) BEI KOMP=0 // C P (X > X0) BEI KOMP=1 // C KOMP SIEHE FV // C // C BENOETIGTES UP: ALGAMA AUS BIBLIOTHEK IMSL // C (LOG.NAT. DER VOLLST. GAMMAFUNKTION) // C(**) double sg; double sga; double a; double y; double a_s; double gn2; double b; double se; double dfv; double db; double rn2; double egn2; double en; double hexp; KOMP = 0; FD = 0.0; FV = 0.0; if ((RN > 0.0) && (X0 >= 0.0)) { if ((X0 <= 0.0)) { if ((Math.Abs(RN - 2.0) < FormCalculation.COGE8)) { FD = 0.5; } // --> return } else { // regel 1 y = X0/2.0; rn2 = RN/2.0; egn2 = Math.Log(Gamma(rn2)); gn2 = egn2; a_s = (rn2 - 1.0)*Math.Log(y) - y - gn2; db = X0/RN; b = Math.Log(db); a = a_s + b; FD = 0.0; if ((a_s > -FormCalculation.CEXP2)) { FD = Math.Exp(a_s)/2.0; } se = 1.0; sg = 0.0; en = RN; do { sga = sg; en = en + 2.0; se = X0*se/en; sg = sg + se; } while (!((sga == sg))); hexp = 0.0; if ((a > -FormCalculation.CEXP2)) { hexp = Math.Exp(a); } dfv = (1.0 + sg)*hexp; if ((dfv >= 0.9)) { dfv = 1.0 - dfv; KOMP = 1; } FV = dfv; } } } // ======================================================================================================================= // Description : // // Date ID Modification // Created // ======================================================================================================================= protected void YCONV(int N, int[] IV, double[] EX, double[] SX, double[][] VP, ref double[] X, double C0, double[] C, ref double BETA, ref double[] ZES, ref int IERror) { // (*C // C VEREINFACHTE VERSION DES FALTUNGSALGORITHMUS FORM4 // C // C EINSCHRAENKUNGEN GEGENUEBER FORM4 // C 1. NUR LINEARKOMBINATION FX = C0 + SUMME(C[i-1]*X[i-1]) // C ALS MECHANISCHES MODELL // C 2. ALLE VARIABLEN UNABHAENGIG // C 3. KEINE LASTWIEDERHOLUNGEN // C // C EINGANG: // C N, IV, EX, SX, VP, EVT. X (JE NACH BETA), NAUS, BETA // C WIE BEI FORM4 // C C0, C(N) KOEFFIZIENTEN DES MECHANISCHEN MODELLS // C // C AUSGANG: // C X, BETA, ZES, IER // C WIE BEI FORM4 // C (**) // DIMENSION EX(N), SX(N), X(N), VP(10,N), C(N), ZES(3), IV(N) // COMMON /CCONV/EXAC(100), SXAC(100) // C int IZ = 0; int IZS = 0; int i = 0; double ABT = 0; double SZ = 0; double EZ = 0; var EXAC = new double[100 - 1 + 1]; var SXAC = new double[100 - 1 + 1]; double FD = 0; double FV = 0; double XV = 0; double XA = 0; double BTV = 0; double ABET = 0; // wat doen met NAUS???? IERror = 0; IZ = 0; IZS = 20; BETA = 0.0; for (i = 1; i <= N; i ++) { YPARB(N, i, IV[i - 1], EX[i - 1], SX[i - 1], ref VP); } do { // regel 50... IZ = IZ + 1; ABT = BETA; SZ = 0.0; EZ = 0.0; for (i = 1; i <= N; i ++) { YPARTR(N, i, IV[i - 1], X[i - 1], VP, ref EXAC[i - 1], ref SXAC[i - 1], ref FD, ref FV); EZ = EZ + C[i - 1]*EXAC[i - 1]; //@ Unsupported property or method(A): 'Power' SZ = SZ + Math.Pow(C[i - 1]*SXAC[i - 1], 2); } SZ = Math.Sqrt(SZ); BETA = -(-C0 - EZ)/SZ; XV = 0.0; for (i = 1; i <= N; i ++) { XA = X[i - 1]; X[i - 1] = EXAC[i - 1] - C[i - 1]*Math.Pow(SXAC[i - 1], 2.0)*BETA/SZ; XA = X[i - 1] - XA; if ((SXAC[i - 1] > 1.0)) { XA = XA/SXAC[i - 1]; } XV = XV + XA*XA; } XV = Math.Sqrt(XV)/N; BTV = Math.Abs(BETA - ABT); ABET = Math.Abs(BETA); if ((ABET > 1.0)) { BTV = BTV/ABET; } // if (IZ > IZS-3) and (NAUS < 99) then YFAUS (NAUS, N, X,4H X ); MAG DEZE REGEL WEG?? // van de repeat~until met regel 50... } while (!(((IZ >= IZS) || ((BTV <= FormCalculation.CFIX4) && (XV <= FormCalculation.CFIX2))))); // NAUS-gebeuren negeren... // (*if (NAUS <> 99) then // begin // FX := CO; // for i:= 1 to N do FX := FX + C[i-1]*X[i-1]; // // WriteLn (NAUS,1040) FX, IZ // YFAUS (NAUS, N, X,4H X-* ); // PF := YPHIN(-BETA); // // WriteLn (NAUS,1050) BETA, PF // end;(**) ZES[1] = C0; ZES[2] = EZ; ZES[3] = SZ; if ((IZ >= IZS)) { IERror = 100000; } } // 1030 FORMAT ('1' ,3X,E12.3) // 1040 FORMAT ('0','G (X STERN) =' ,E10.3, // . 'ANZ. IT. = ',I2) // 1050 FORMAT ('0', 'BETA =',F8.3,'PF =',1G12.3) // end // *************************************************************** // (**) // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected void YEXPO(double X, double A, double RL, ref double FD, ref double FV, ref int KOMP) { // (*C // C E X P O N E N T I A L V E R T E I L U N G // C // C FV = 1 - EXP ( - RL * ( X - A ) ) // C // C FD = RL * EXP ( - RL * ( X - A ) ) // C // C MIT A = EX - SX // C RL = 1 / SX // C UND X > A // C // C BENOETIGTES UP: YREIHE // C (**) double Y; double H; KOMP = 0; FD = 0.0; FV = 0.0; Y = X - A; if ((Y >= 0.0)) { H = -Y*RL; FD = RL*Math.Exp(H); if ((Y > 0.0)) { if ((Y < (1.0/RL))) { FV = -YREIHE(H); } else { FV = Math.Exp(H); KOMP = 1; } } } } // ======================================================================================================================= // Description : // // Date ID Modification // Created // ======================================================================================================================= protected void YFDFV(int N, int I, int IV, double X0, double[][] VP, ref double FD, ref double FV, ref int KOMP) { // C // C DICHTE- UND VERTEILUNGSFUNKTION DER VARIABLEN X // C VOM VERTEILUNGSTYP IV AN DER STELLE X0 // C // C EINGANG: // C N SPALTENZAHL DER 10-ZEILIGEN MATRIX VP // C = ZAHL DER VEREINBARTEN VARIABLEN // C I INDEX DER VARIABLEN X[i-1], FUER DIE FD UND FV ZU // C BERECHNEN SIND // C IV VERTEILUNGSTYP DER VARIABLEN X[i-1]; ZAHLENSCHLUESSEL // C SIEHE HEFT 43 / ANHANG F // C X0 WERT DER ZUFALLSVARIABLEN X[i-1] // C VP VERTEILUNGSPARAMETER; BESTIMMUNG VORAB BEISPIELSWEISE // C DURCH UP YPARB // C // C AUSGANG: // C FD WERT DER DICHTEFUNKTION VON X AN DER STELLE X0 // C FV P (XX0) BEI KOMP=1 // C KOMP SIEHE FV // C // C BENOETIGTE UP: YUNI, YNORM, YLNORM, YEXPO, YGAMVD, YGUMB, // C YFRECH, YWEIB, YV96, YV97, YV98, YV99, // C YBINEX, YGAMAX, YGAMIN, YPFEX // C // C ANMERKUNG: BETA-VERTEILUNG IST NOCH NICHT IMPLEMENTIERT. // DIMENSION VP(10,N) int IIV; double FDI; double FVI; int KOMPI; int IT; bool Eruit; // wordt gebruikt ipv de exit-commando's. Eruit = false; KOMP = 0; if (((IV == 302) || (IV == 303) || (IV == 310))) { YGAMAX(X0, VP[0][I - 1], VP[1][I - 1], VP[2][I - 1], VP[6][I - 1], IV, ref FD, ref FV, ref KOMP); Eruit = true; } if (!Eruit) { if (((IV == 402) || (IV == 403) || (IV == 410))) { // YGAMIN (X0, VP[1,I-1], VP[2,I-1], VP[3,I-1], VP[7,I-1], IV, FD, FV, KOMP); YGAMIN(X0, VP[0][I - 1], VP[1][I - 1], VP[2][I - 1], VP[6][I - 1], IV, ref FD, ref FV, ref KOMP); Eruit = true; } } if (!Eruit) { IIV = IV - Convert.ToInt32(IV/100)*100; switch (IIV) { case 96: // afgeknotte verdeling YV96(X0, VP[0][I - 1], VP[1][I - 1], VP[2][I - 1], VP[3][I - 1], ref FD, ref FV, ref KOMP); break; case 97: break; case 98: break; case 99: break; case 1: // YV97 (X0, VP[1,I-1], VP[2,I-1], VP[3,I-1], VP[4,I-1],FD, FV, KOMP); 19990925,HenkJanFaber. // YV98 (X0, VP[1,I-1], VP[2,I-1], VP[3,I-1], VP[4,I-1],FD, FV, KOMP); // YV99 (X0, VP[1,I-1], VP[2,I-1], VP[3,I-1], VP[4,I-1],FD, FV, KOMP); YUNI(X0, VP[0][I - 1], VP[1][I - 1], ref FD, ref FV, ref KOMP); break; case 2: YNORM(X0, VP[0][I - 1], VP[1][I - 1], ref FD, ref FV, ref KOMP); break; case 3: VP[2][I - 1] = 0.0; YLNORM(X0, VP[0][I - 1], VP[1][I - 1], VP[2][I - 1], ref FD, ref FV, ref KOMP); break; case 4: YEXPO(X0, VP[0][I - 1], VP[1][I - 1], ref FD, ref FV, ref KOMP); break; case 5: YGAMVD(X0, VP[0][I - 1], VP[1][I - 1], ref FD, ref FV, ref KOMP); break; case 6: IV = 2; // (* deze optie bestaat nog niet wordt normale verdeling van gemaakt*) YNORM(X0, VP[0][I - 1], VP[1][I - 1], ref FD, ref FV, ref KOMP); break; case 7: YGUMB(X0, VP[0][I - 1], VP[1][I - 1], ref FD, ref FV, ref KOMP); break; case 8: YFRECH(X0, VP[0][I - 1], VP[1][I - 1], VP[2][I - 1], ref FD, ref FV, ref KOMP); break; case 9: YWEIB(X0, VP[0][I - 1], VP[1][I - 1], VP[2][I - 1], ref FD, ref FV, ref KOMP); break; case 10: YLNORM(X0, VP[0][I - 1], VP[1][I - 1], VP[2][I - 1], ref FD, ref FV, ref KOMP); break; } if ((IV >= 100)) { FDI = FD; FVI = FV; KOMPI = KOMP; IT = 1; if ((IV <= 500)) { if ((IV > 200)) { IT = 2; } YBINEX(IT, X0, VP[6][I - 1], VP[7][I - 1], FDI, FVI, KOMPI, ref FD, ref FV, ref KOMP); } else { if ((IV > 600)) { IT = 2; } YPFEX(IT, VP[6][I - 1], FDI, FVI, KOMPI, ref FD, ref FV, ref KOMP); } } } } // ******************************************************************** // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected void YFRECH(double X, double U, double RK, double EPS, ref double FD, ref double FV, ref int KOMP) { // (*C // C TYP II F R E C H E T GROESSTWERTE // C // C FV = EXP ( - ( (U-EPS) / (X-EPS) ** RK ) X > EPS // C // C FD = RK / (U-EPS) * ( (U-EPS) / (X-EPS) ) ** ( RK + 1 ) * FV // C // C BENOETIGTES UP: YREIHE // C(**) double H; double Y; double V; KOMP = 0; FD = 0.0; FV = 0.0; Y = X - EPS; V = U - EPS; if ((Y > 0.0)) { H = Math.Pow((V/Y), RK); if ((H <= FormCalculation.CEXP2)) { FV = Math.Exp(-H); FD = RK*H*FV/Y; if ((X > U)) { KOMP = 1; FV = -YREIHE(-H); } } } } // ======================================================================================================================= // Description : // ASYMPTOTISCHE VERTEILUNG DES MAXIMUMS // EINES NORMALPROZESSES (IV=2) // BZW. EINES LOGNORMALPROZESSES (IV=3) // // FV = FVI * EXP ( -AG * FDI / FVI ) // // FD = ( 1 + AG * ( U + FDI / FVI ) ) * FDI // * EXP ( -AG * FDI / FVI ) / SSTERN // // MIT AG = ALPHA * GAMMA // U = (X-EX)/SX (IV=2) BZW. = DLOG((X-TAU)/EXL)/SXL (IV=3) // SSTERN = SX (IV=2) BZW. = (X-TAU)*SXL (IV=3) // // BENOETIGTE UP: YNORM, YREIHE // // Date ID Modification // Created // ======================================================================================================================= protected void YGAMAX(double X, double EX, double SX, double TAU, double AG, int IV, ref double FD, ref double FV, ref int KOMP) { int J = 0; double U = 0; double XMT = 0; double DIVert = 0; double C = 0; double EGA = 0; double FVI = 0; double H = 0; double DV = 0; double FDI = 0; KOMP = 0; J = IV - Convert.ToInt32(IV/100)*100; if ((J <= 2)) { U = (X - EX)/SX; DIVert = SX; } else { XMT = X - TAU; U = Math.Log(XMT/EX)/SX; DIVert = XMT*SX; } YNORM(U, 0.0, 1.0, ref FDI, ref FVI, ref KOMP); H = FVI; if ((KOMP != 0)) { U = -U; H = 1.0 - H; } DV = FDI/H; C = -AG*DV; EGA = Math.Exp(C); FD = EGA*FDI*(1.0 + AG*(U + DV))/DIVert; FV = EGA*H; if ((KOMP != 0)) { if ((FV <= 0.5)) { KOMP = 0; } else { FVI = EGA*FVI; FV = FVI - YREIHE(C); } } } // ======================================================================================================================= // Description : // ASYMPTOTISCHE VERTEILUNG DES MINIMUMS // EINES NORMALPROZESSES (IV=2) // BZW EINES LOGNORMALPROZESSES (IV=3) // // WEITERE BESCHREIBUNG ANALOG ZU YGAMAX // // BENOETIGTES UP: YGAMAX // // Date ID Modification // Created // ======================================================================================================================= protected void YGAMIN(double X, double EX, double SX, double TAU, double AG, int IV, ref double FD, ref double FV, ref int KOMP) { int J; double U; double XS; XS = 0.0; J = IV - Convert.ToInt32(IV/100)*100; if ((J <= 2)) { U = (EX - X)/SX; } else { XS = X - TAU; U = Math.Log(EX/XS)/SX; } YGAMAX(U, 0.0, 1.0, 0.0, AG, 2, ref FD, ref FV, ref KOMP); FD = FD/SX; // Compiler geeft warning op onderstaande regel -> negeren. if ((J > 2)) { FD = FD/XS; } KOMP = 1 - KOMP; } // ======================================================================================================================= // Description : // G A M M A V E R T E I L U N G // // EINGANG: // X0 WERT DER ZUFALLSVARIABLEN X (IM DEFINITIONSBEREICH // 0 <= X0 < +UNendLICH) // RK = V*V // RL = V/SX // MIT V = EX/SX // EX ERWARTUNGSWERT VON X // UND SX STANDARDABWEICHUNG VON X // // AUSGANG: // FD DICHTE AN DER STELLE X0 // FV P (X < X0) KOMP=0 // P (X > X0) KOMP=1 // KOMP SIEHE FV // // METHODE: BERECHNUNG MIT HILFE DER ZENTRALEN CHIQUADRATVERTEILUNG // // BENOETIGTE UP: YCHI2 SOWIE MITTELBAR UEBER YCHI2 ALGAMA AUS IMSL // // Date ID Modification // Created // ======================================================================================================================= protected void YGAMVD(double X0, double RK, double RL, ref double FD, ref double FV, ref int KOMP) { double Z; double A; KOMP = 0; FD = 0.0; FV = 0.0; if ((X0 >= 0.0)) { Z = 2.0*X0*RL; A = 2.0*RK; YCHI2(Z, A, ref FD, ref FV, ref KOMP); FD = 2.0*FD*RL; } } // ======================================================================================================================= // Description : // // Date ID Modification // Created // ======================================================================================================================= protected void YGGS(int N, double[] X, int[] IV, ref double THETA, int IGGS, double FX, double[] HX, ref double ABLT, ref double[] ABL, ref double[][] ABL2, int NA, ref int IERror) { // GGFunctie: TZfunctiePointerType; // EXTERNAL GG // DIMENSION X(N), ABL(N), ABL2(NA,NA), HX(N), IV(N) // (*C // C NUMERISCHE BERECHNUNG VON PARTIELLEN ABLEITUNGEN // C DER FUNKTION GG := GG (X, THETA) // C // C EINGANG: // C N DIMENSION // C X VEKTOR DER AKTUELLEN WERTE VON X // C IV VEKTOR MIT DIMENSION N // C BEI IV[i-1]=0 WERDEN ABL[i-1]=ABL2[i,j-1]=0 GESETZT, // C WOMIT DIE BERECHNUNG EINGESPART WERDEN KANN. // C THETA EINZELNE VARIABLE IN GG // C IGGS STEUERGROESSE; BESTIMMT, WELCHE ABLEITUNGEN // C BERECHNET WERDEN. // C IGGS = 0: ABLT = D(GG)/D(THETA) // C IGGS = 1: ABL[i-1] = D(GG)/D(X[i-1]) // C IGGS = 2: ABL[i-1] WIE BEI IGGS=1 // C ABL2[i,j-1] = D2(GG)/(D(X[i-1])*D(X[j-1])) // C FX FUNKTIONSWERT VON GG AN DER STELLE (X, THETA) // C NA DIMENSION DER QUADRATISCHEN MATRIX ABL2 // C NA >= 1, WENN IGGS < 2 // C NA >= N, WENN IGGS = 2 // C // C HILFSGROESSE: // C HX VEKTOR MIT DIMENSION N // C // C AUSGANG: // C ABLT SIEHE IGGS // C ABL SIEHE IGGS // C ABL2 SIEHE IGGS // C IER FEHLERPARAMETER AUS GG; IER > 0 BEWIRKT // C ABBRUCH VON YGGS // C // C METHODE: // C SIEHE HEFT 43/ABSCHN. 3.5.3 // C // C(**) double EPS = 0; double XV = 0; double XH = 0; double XJ = 0; double HF = 0; double FO = 0; int i = 0; int i1 = 0; int j = 0; bool Eruit = false; // Wordt gebruikt ipv de exit-commando's.. Eruit = false; IERror = 0; EPS = FormCalculation.COGE27; if ((IGGS <= 0)) { // else goto 10 // C // C BERECHNUNG VON ABLT XV = Math.Abs(THETA); if ((XV < 1.0)) { XV = 1.0; } XH = XV*EPS; XV = THETA; if ((ABLT < 0.0)) { XH = -XH; } THETA = THETA + XH; GG(N, X, THETA, ref IERror, ref FO); // WriteLn(LFErr, ' x0 x1 en fx(fo) ', x[0]: 12: 5, X[1]: 12: 5, FO: 12: 5); if ((IERror == 0)) { XH = THETA - XV; HF = FO - FX; THETA = XV; if ((Math.Abs(HF) <= FormCalculation.CEMIN3)) { ABLT = 0.0; // --> return } else { if ((Math.Abs(XH) <= FormCalculation.CEMIN3)) { ABLT = FormCalculation.CEMAX4*signUnit(XH)*signUnit(HF); // --> return } else { ABLT = HF/XH; } // --> return } } // --> return // C // C BERECHNUNG VON ABL UND ABL2 : FUNKTIONSWERTE // 10 } else { // IGGS > 0 for (i = 1; i <= N; i ++) { if ((IV[i - 1] == 0)) { ABL[i - 1] = 0.0; if ((IGGS != 1)) { for (j = 1; j <= i; j ++) { ABL2[j - 1][i - 1] = 0.0; } } } } for (i = 1; i <= N; i ++) { if ((IV[i - 1] != 0)) { XV = Math.Abs(X[i - 1]); if ((XV < 0.01)) { XV = 0.01; } XH = XV*EPS; XV = X[i - 1]; if ((ABL[i - 1] < 0.0)) { XH = -XH; } X[i - 1] = X[i - 1] + XH; GG(N, X, THETA, ref IERror, ref ABL[i - 1]); // WriteLn(LFErr, ' x0 x1 en fx(ABL[i]) ', x[0]: 12: 5, X[1]: 12: 5, ABL[i - 1]: 12: // 5); if ((IERror > 0)) { Eruit = true; } if (!Eruit) { XH = X[i - 1] - XV; X[i - 1] = XV; HX[i - 1] = XH; } } } if (!Eruit) { if ((IGGS != 1)) { // C // C BERECHNUNG VON ABL2 for (i = 1; i <= N; i ++) { if ((IV[i - 1] != 0)) { XV = X[i - 1]; X[i - 1] = X[i - 1] + HX[i - 1]; i1 = i - 1; if ((i1 != 0)) { for (j = 1; j <= i1; j ++) { if ((IV[j - 1] != 0)) { XJ = X[j - 1]; X[j - 1] = X[j - 1] + HX[j - 1]; GG(N, X, THETA, ref IERror, ref ABL2[j - 1][i - 1]); // WriteLn(LFErr, ' x0 x1 en fx(ABL2[j,i])) ', x[0]: 12: 5, X[1]: 12: 5, // ABL2[j - 1, i - 1]: 12: 5); if ((IERror > 0)) { Eruit = true; } if (!Eruit) { X[j - 1] = XJ; } } } if (!Eruit) { X[i - 1] = XV - HX[i - 1]; GG(N, X, THETA, ref IERror, ref ABL2[i - 1][i - 1]); // WriteLn(LFErr, ' x0 x1 en fx(ABL2[I,I])) ', x[0]: 12: 5, X[1]: 12: 5, ABL2[I // - 1, I - 1]: 12: 5); if ((IERror > 0)) { Eruit = true; } X[i - 1] = XV; } } } } if (!Eruit) { for (i = 1; i <= N; i ++) { if ((IV[i - 1] != 0)) { i1 = i - 1; if ((i1 != 0)) { for (j = 1; j <= i1; j ++) { if ((IV[j - 1] != 0)) { ABL2[j - 1][i - 1] = (ABL2[j - 1][i - 1] - ABL[i - 1] - ABL[j - 1] + FX)/(HX[i - 1]*HX[j - 1]); } } } ABL2[i - 1][i - 1] = (ABL2[i - 1][i - 1] + ABL[i - 1] - 2.0*FX)/(HX[i - 1]*HX[i - 1]); } } } } } // BERECHNUNG VON ABL if (!Eruit) { for (i = 1; i <= N; i ++) { if ((IV[i - 1] != 0)) { HF = ABL[i - 1] - FX; if ((Math.Abs(HF) <= FormCalculation.CEMIN3)) { ABL[i - 1] = 0.0; // --> 90 } else { if ((Math.Abs(HX[i - 1]) <= FormCalculation.CEMIN3)) { ABL[i - 1] = FormCalculation.CEMAX4*signUnit(HX[i - 1])*signUnit(HF); } else { ABL[i - 1] = HF/HX[i - 1]; } } } } } } } // ************************************************************ // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected void YGUMB(double X, double U, double A, ref double FD, ref double FV, ref int KOMP) { // (*C // C TYP I G U M B E L GROESSTWERTE // C // C FV = EXP ( - EXP ( - A * ( X - U ) ) ) // C // C FD = A * EXP ( - A * ( X - U ) - EXP ( - A * ( X - U ) ) ) // C // C BENOETIGTES UP: YREIHE // C(**) double H; double EH; FD = 0.0; FV = 0.0; KOMP = 0; H = -A*(X - U); if ((X <= U)) { if ((H <= FormCalculation.CEXP3)) { EH = Math.Exp(H); FV = Math.Exp(-EH); FD = A*Math.Exp(H - EH); } } else { KOMP = 1; if ((H >= -FormCalculation.CEXP2)) { EH = Math.Exp(H); FV = -YREIHE(-EH); FD = A*Math.Exp(H - EH); } } } // ************************************************************ // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected void YLNORM(double X0, double ELN, double SLN, double TAU, ref double FD, ref double FV, ref int KOMP) { // (*C // C L O G N O R M A L V E R T E I L U N G // C // C EINGANG: // C X0 WERT DER ZUFALLSVARIABLEN X // C ELN = (EX-TAU) / DSQRT(1+V*V) // C SLN = DSQRT(DLOG(1+V*V)) // C TAU UNTERE GRENZE DES DEFINITIONSBEREICHES VON X // C MIT V = SX/(EX-TAU) // C EX ERWARTUNGSWERT VON X // C UND SX STANDARDABWEICHUNG VON X // C // C AUSGANG: // C FD WERT DER DICHTEFUNKTION VON X AN DER STELLE X0 // C FV P (X < X0) KOMP=0 // C P (X > X0) KOMP=1 // C KOMP SIEHE FV // C // C METHODE: X WIRD AUF NORMALVERTEILUNG TRANSFORMIERT // C // C BENOETIGTES UP: YPHIN // C(**) double H; KOMP = 0; FD = 0.0; FV = 0.0; H = X0 - TAU; if ((H > 0.0)) { H = H/ELN; if ((H > FormCalculation.CEMIN)) { H = Math.Log(H)/SLN; if ((H > 0.0)) { KOMP = 1; H = -H; } FV = YPHIN(H); H = H*H/2.0; if ((H <= FormCalculation.CEXP2)) { FD = Math.Exp(-H)/((X0 - TAU)*SLN*2.506628275); } } } } // ************************************************************ protected double YNINVF(double P) { double result; // (*C // C INVERTIERTES STANDARDNORMALVERTEILUNGSINTEGRAL // C X := YNINVF (P) = INVERSE (PHI (X)) // C // C METHODE: // C YNINVG LIEFERT DEI ANFANGSLOESUNG, DIE MIT NEWTONVERFAHREN // C VERBESSERT WIRD. // C // C BENOETIGTE UP: YNINVG, YPHIN // C(**) double H; double EPS; double FXS; double FX; H = 0.0; if ((P != 0.5)) { H = -FormCalculation.CEXP1; if ((P != 0.0)) { EPS = FormCalculation.COGE3; H = YNINVG(P); FXS = YPHIN(H + EPS); FX = YPHIN(H); while (((Math.Abs(FX/P - 1.0) >= FormCalculation.COGE23) && ((FX - FXS) != 0.0))) { EPS = (FX - P)*EPS/(FXS - FX); H = H - EPS; FXS = FX; FX = YPHIN(H); } } } // 10 result = H; return result; } // ************************************************************ // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected double YNINVG(double P) { double result; // (*C // C INVERTIERTES STANDARDNORMALVERTEILUNGSINTEGRAL // C X = YNINVG (P) = INVERSE (PHI (X)) // C // C METHODE: FORMEL 26.2.23 (ABRAMOWITZ/SEGUN) // C // C ANMERKUNG: MAESSIGE GENAUIGKEIT // C(**) double A; double H; A = P; H = 1.0; if ((P > 0.5)) { A = 1.0 - P; H = -1.0; } if ((A <= 0.0)) { A = FormCalculation.CEMIN; } A = Math.Sqrt(-2.0*Math.Log(A)); A = A - ((0.010328*A + 0.802853)*A + 2.515517)/(((0.001308*A + 0.189269)*A + 1.432788)*A + 1.0); if ((H > 0.0)) { A = -A; } result = A; return result; } // ************************************************************ // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected double YNOR1G(double X0) { double result; // (*C ******************************************************************** // C * * // C *** BERECHNUNG DER EINDIMENSIONALEN NORMALVERTEILUNG MIT HOHER * // C *** GENAUIGKEIT * // C * * // C * EINGABE: * // C * X0 WERT FUER DEN DIE VERTEILUNG BERECHNET WERDEN SOLL * // C * * // C * AUSGABE: * // C * YNOR1G WERT DER EINDIMENSIONALEN NORMALVERTEILUNG * // C * * // C ******************************************************************** // C(**) int J; int N; double XA; double X2; double PI; double FAK; double A0; double A1; double A2; double AI; double AN; double B0; double B1; double B2; double BN; double F1; double F2; double Alpha; double YNor12; double Ynor11; // FNOR11, FNOR12: double; bool Eruit; // wordt gebruikt ipv de exit-commando's.. Eruit = false; Ynor11 = 0.0; result = 0.0; if ((X0 >= -FormCalculation.CEXP1)) { result = 0.5; XA = Math.Abs(X0); if ((XA >= FormCalculation.COGEN)) { result = 1.0; if ((X0 <= FormCalculation.CONE1)) { PI = FormCalculation.COMPI; X2 = XA*XA; FAK = Math.Exp(-X2*0.5)/Math.Sqrt(2.0*PI); if ((XA < 4.0)) { // C *** BERECHNUNG NACH FORMEL 26.2.15 (ABRAMOWITZ/SEGUN) A0 = 0.0; A1 = XA; B0 = 1.0; B1 = 1.0; F1 = XA; N = 2; J = -1; AN = J*X2*(N - 1.0); BN = 2.0*N - 1.0; A2 = A1*BN + A0*AN; B2 = B1*BN + B0*AN; F2 = A2/B2; while ((Math.Abs((F2 - F1)/F1) >= FormCalculation.COGE8)) { A0 = A1; B0 = B1; B1 = B2; A1 = A2; F1 = F2; N = N + 1; J = -J; AN = J*X2*(N - 1.0); BN = 2.0*N - 1.0; A2 = A1*BN + A0*AN; B2 = B1*BN + B0*AN; F2 = A2/B2; } // 500 Ynor11 = 0.5 - FAK*F2; if ((XA <= 3.9)) { if ((X0 > 0.0)) { Ynor11 = 1.0 - Ynor11; } result = Ynor11; Eruit = true; } } if (!Eruit) { // C *** BERECHNUNG NACH FORMEL 26.2.14 (ABRAMOWITZ/SEGUN) // 100 B0 = 1.0; B1 = XA; A0 = 0.0; A1 = 1.0; F1 = 1.0/XA; AI = 1.0; // 110 A2 = XA*A1 + AI*A0; B2 = XA*B1 + AI*B0; F2 = A2/B2; while ((Math.Abs((F2 - F1)/F1) >= FormCalculation.COGE8)) { A0 = A1; B0 = B1; A1 = A2; B1 = B2; F1 = F2; AI = AI + 1.0; A2 = XA*A1 + AI*A0; B2 = XA*B1 + AI*B0; F2 = A2/B2; } // 510 YNor12 = FAK*F2; if ((XA >= 4.0)) { if ((X0 > 0.0)) { YNor12 = 1.0 - YNor12; } result = YNor12; } else { Alpha = (4.0 - XA)*10.0; A0 = (1.0 - Alpha)*YNor12 + Alpha*Ynor11; if ((X0 > 0.0)) { A0 = 1.0 - A0; } result = A0; } } } } } return result; } // ======================================================================================================================= // Method : NormalDistribution(Argm: Double) : Double; // // Member of class : // // Acces : // // Description: Returns Gaussian Probability Function; // accuracy 4.7(-10) // ERROR IS PROPORTIONAL WITH Z(y) // // Pre-condition : // // Post-condition : // // Method History : Abramowitz and Stegun, 26.2.14/15 // // Date ID Modification // 1996-12-09 Sel Created // ======================================================================================================================= // ======================================================================================================================= // Method : NormalDistrInverse(Argm: Double) : Double; // // Member of class : // // Acces : // // Description: Returns Inverse Gaussian Probability Function; // accuracy 7.2(-10) // // Pre-condition : Argument Range between 0 and 1 // // Post-condition : For Outcome MaxDouble, check Argument // // Comment : Abramowitz and Stegun, 26.2.23; // Taylor Series correction with two terms; // Outcome is between -MaxDouble and +MaxDouble, // which practically means between -10 and +10; // MaxDouble = 1.7E+308; If this value is obtained, // it probably indicates an invalid argument range. // // Method History : // Date ID Modification // 1996-12-09 Sel Created // 1997-03-06 Sel MaxDouble added to indicate invalid argument range // ======================================================================================================================= protected double NormalDistrInverse(double Argm) { double result; double Temp; double value; value = 0; if ((Argm <= 0) || (Argm >= 1)) { if (Argm <= 0) { value = -10000; } if (Argm >= 1) { value = 10000; } } else { if (Argm < 0.5) { Temp = Math.Sqrt(-2*Math.Log(Argm)); } else { Temp = Math.Sqrt(-2*Math.Log(1 - Argm)); } value = Temp - (2.515517 + Temp*(0.802853 + Temp*0.010328))/(1 + Temp*(1.432788 + Temp*(0.189269 + Temp*0.001308))); if (Argm < 0.5) { value = -value; } Temp = 2*value*(Argm - NormalDistribution(value)); Temp = 1 - Math.Sqrt(1 - Temp*Math.Sqrt(2*Math.PI)*Math.Exp(value*value/2)); value = value + Temp/value; } result = value; return result; } // ************************************************************ // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected void YNORM(double X0, double EX, double SX, ref double FD, ref double FV, ref int KOMP) { // (*// N O R M A L V E R T E I L U N G // // // // EINGANG: // // X0 WERT DER ZUFALLSVARIABLEN X // // EX ERWARTUNGSWERT VON X // // SX STandARDABWEICHUNG VON X // // // // ***** Y N O R M IS Y P H I N AS A procedure. // // AUSGANG: // // FD WERT DER DICHTEFUNKTION VON X AN DER STELLE X0 // // FV P (X < X0) BEI KOMP=0 // // P (X > X0) BEI KOMP=1 // // KOMP SIEHE FV // // // // METHODE: // // BERECHNUNG DES NORMALVERTEILUNGSINTEGRALS NACH FORMEL 26.2.17 // // (ABRAMOWITZ/SEGUN) // (**) double Y; double R; double T; FD = 0.0; FV = 0.0; KOMP = 0; if ((SX <= 0.0)) { if ((X0 == EX)) { FD = 1.0; } if ((X0 >= EX)) { FV = 1.0; } } else { Y = (X0 - EX)/SX; if ((Y > 0.0)) { Y = -Y; KOMP = 1; } Y = -Y; R = Y; T = 0.0; if ((R <= FormCalculation.CEXP1)) { R = Math.Exp(-R*R/2.0)/2.506628275; FD = R/SX; T = 0.5; if ((Y != 0.0)) { T = 1.0/(1.0 + 0.2316419*Y); T = ((((1.330274429*T - 1.821255978)*T + 1.781477937)*T - 0.356563782)*T + 0.31938153)*T; T = R*T; } } FV = T; } } // ************************************************************ // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected double YNORMD(double X0, double EX, double SX) { double result; // (*C // C DICHTEFUNKTION EINER NORMALVERTEILTEN ZUFALLSVARIABLEN X // C MIT MITTELWERT EX UND STandARDABWEICHUNG SX AN DER STELLE X0 // C(**) double H; result = 0.0; if ((SX <= 0.0)) { if ((X0 == EX)) { result = 1.0; } } else { H = Math.Abs((X0 - EX)/SX); if ((H <= FormCalculation.CEXP1)) { result = Math.Exp(-H*H/2.0)/(2.506628275*SX); } } return result; } // ************************************************************ // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected double YNPOT(double RN, double Y) { double result; // (*C // C *****EVALUATION OF ( (1+Y)**RN -1. ) WITH POWER-SERIES : // C **** SU := S0 + S0*F1 + (S0*F1)*F2 + ..... // C ***** FOR : -1. < = Y < = +1. and // C ***** RN > 0. // C // C ***** N = INT(RN) if RN = FLOAT(Integer) ON INPUT, IN THIS CASE // C ***** THERE ARE ONLY N (MAX = 100 !!) ITERATIONS. // C(**) int N1; int N; int I; double SI; double SU; double RI; double FI; bool Eruit; // Er uit wordt gebruikt ipv 'exit' in de repeat-until. Eruit = false; N1 = 100; if ((RN == Convert.ToInt32(RN))) { N1 = Convert.ToInt32(RN); } N = Int32.MinValue; // N := MIN0(N1,100); I = 0; SI = 1.0; // C *****BELOW MACHINE PRECISION THERE IS NO NEED FOR ITERATION : SU = RN*Y; if ((Math.Abs(SU) >= FormCalculation.COGEN)) { // C *****POWER-SERIES : SU = 0.0; do { I = I + 1; RI = I; FI = ((RN - RI + 1.0)/RI)*Y; SI = SI*FI; SU = SU + SI; if ((I >= N)) { Eruit = true; } } while (!((Math.Abs(SI/SU) <= FormCalculation.COGE8) || Eruit)); } // 2 result = SU; return result; } // ************************************************************ // // // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected void YPARB(int N, int I, int IVT, double E, double S, ref double[][] VP) { // (*C // C ERMITTLUNG VON VERTEILUNGSPARAMETERN AUS ERWARTUNGSWERT, STANDARD- // C ABWEICHUNG UND EVT. DEFINITIONSBEREICHSGRENZEN // C DIE VERTEILUNGS-UP WERDEN IM ALLGEMEINEN MIT DIESEN PARAMETERN // C UND NICHT MIT DEN MOMENTEN AUFGERUFEN, DA DIE FORMELN MEIST AUF // C AUF DIE PARAMETER ZUGESCHNITTEN SIND. // C // C EINGANG: // C N SPALTENZAHL DER 10-ZEILIGEN MATRIX VP // C := ZAHL DER VEREINBARTEN ZUFALLSVARIABLEN // C I INDEX DER VARIABLEN X[i-1], FUER DIE DIE PARAMETER ZU // C BERECHNEN SIND // C IVT VERTEILUNGSTYP DER VARIABLEN // C ZAHLENSCHLUESSEL SIEHE HEFT 43 / ANHANG F // C E ERWARTUNGSWERT // C S STANDARDABWEICHUNG // C VP NUR BEI IVT=6,8,9 ODER 10 // C VP(3,I) = UNTERE GRENZE DES DEFINITIONSBEREICHES // C NUR BEI IVT = 6 // C VP(4,I) = OBERE GRENZE DES DEFINITIONSBEREICHES // C // C AUSGANG: // C VP BESETZT WERDEN VP(1,I) UND VP(2,I) // C ZUM FORMELMAESSIGEN ZUSAMMENHANG ZWISCHEN DEN PARAMETERN // C UND E UND S SIEHE HEFT 43 / ANHANG F / TABELLE F1 // C // C BENOETIGTE UP: YREFA // C YVGAM // C GAMMA (VOLLST. GAMMAFUNKTION AUS BIBL. IMSL) // C(**) // De Integer IVC en double Variantie zijn global gemaakt.. int IV; double H; double F; double G; double X0; double EPS; int MAXIT; IV = IVT; if ((IV > 99)) { IV = IV - Convert.ToInt32(Convert.ToInt32(IV/100)*100); } if (((IV < 96) || (IV > 99))) { switch (IV) { case 1: H = S*Math.Sqrt(3.0); VP[0][I - 1] = E - H; VP[1][I - 1] = E + H; break; case 2: VP[0][I - 1] = E; VP[1][I - 1] = S; break; case 3: VP[2][I - 1] = 0.0; // 31 H = S/(E - VP[2][I - 1]); H = 1.0 + H*H; VP[0][I - 1] = (E - VP[2][I - 1])/Math.Sqrt(H); VP[1][I - 1] = Math.Sqrt(Math.Log(H)); break; case 4: VP[0][I - 1] = E - S; VP[1][I - 1] = 1.0/S; break; case 5: H = E/S; VP[0][I - 1] = H*H; VP[1][I - 1] = H/S; break; case 6: H = S/(VP[3][I - 1] - VP[2][I - 1]); G = (E - VP[2][I - 1])/(VP[3][I - 1] - VP[2][I - 1]); VP[0][I - 1] = ((1.0 - G)*G/(H*H) - 1.0)*G; VP[1][I - 1] = (1.0/G - 1.0)*VP[0][I - 1]; break; case 7: VP[1][I - 1] = 1.28254983/S; VP[0][I - 1] = E - 0.4500534689*S; break; case 8: F = -1.0; // 91 Variantie = S/(E - VP[2][I - 1]); IVC = IV; EPS = FormCalculation.COGE27*S; X0 = 1.0/Variantie; if (((IV == 8) && (X0 < 2.1))) { X0 = 2.1; } if (((IV == 9) && (X0 < 0.01))) { X0 = 0.01; } MAXIT = 20; FIerror = 0; YREFA(YVGAM, ref VP[1][I - 1], X0, EPS, MAXIT, ref FIerror); VP[0][I - 1] = (E - VP[2][I - 1])/Gamma(1.0 + F/VP[1][I - 1]) + VP[2][I - 1]; break; case 9: F = 1.0; // 91 Variantie = S/(E - VP[2][I - 1]); IVC = IV; EPS = FormCalculation.COGE27*S; X0 = 1.0/Variantie; if (((IV == 8) && (X0 < 2.1))) { X0 = 2.1; } if (((IV == 9) && (X0 < 0.01))) { X0 = 0.01; } MAXIT = 20; FIerror = 0; // YREFA (YVGAM, VP[2,I-1], X0, EPS, MAXIT, IERror); // VP[1,I-1] := (E-VP[3,I-1]) / GAMMA(1.+F/VP[2,I-1]) + VP[3,I-1]; YREFA(YVGAM, ref VP[1][I - 1], X0, EPS, MAXIT, ref FIerror); VP[0][I - 1] = (E - VP[2][I - 1])/Gamma(1.0 + F/VP[1][I - 1]) + VP[2][I - 1]; break; case 10: // 31 H = S/(E - VP[2][I - 1]); H = 1.0 + H*H; VP[0][I - 1] = (E - VP[3][I - 1])/Math.Sqrt(H); VP[1][I - 1] = Math.Sqrt(Math.Log(H)); break; } } } // *********************************************************** protected void YPARTR(int N, int I, int IV, double X0, double[][] VP, ref double EA, ref double SA, ref double FD, ref double FV) { // (*C // C VERTEILUNGSAPPROXIMATION // C YPARTR BESTIMMT ERWARTUNGSWERT EA UND STANDARDABWEICHUNG SA // C EINER NORMALVERTEILUNG SO, DASS AN DER STELLE X=X0 SICH NACH // C DER NORMALVERTEILUNG DIE GLEICHEN WERTE FUER DICHTE UND VERTEI- // C LUNGSFUNKTION ERGEBEN WIE NACH DER NICHTNORMALEN VERTEILUNG VON // C TYP IV VON X. // C // C EINGANG: // C N SPALTENZAHL DER 10-ZEILIGEN MATRIX VP // C = ZAHL DER VEREINBARTEN VARIABLEN // C I INDEX DER VARIABLEN X[i-1], FUER DIE EA UND SA ZU BESTIMMEN // C SIND // C IV VERTEILUNGSTYP DER VARIABLEN X[i-1] // C ZAHLENSCHLUESSEL SIEHE HEFT 43 / ANHANG F // C BEI IV=0 WERDEN FD UND FV ZU EINGAENGEN, GLEICHZEITIG // C WERDEN ALLE and EREN EINGAENGE BEDEUTUNGSLOS // C X0 WERT DER VARIABLEN X[i-1] // C VP VERTEILUNGSPARAMETER. VORHERIGE BESTIMMUNG BEISPIELSWEISE // C MIT YPARB // C // C AUSGANG: // C EA ERWARTUNGSWERT DER APPROX. NORMALVERTEILUNG // C SA STANDARDABWEICHUNG DER APPROX. NORMALVERTEILUNG // C FD // C UND FV WERTE DER DICHTEFUNKTION UND DER VERTEILUNGSFUNKTION VON X // C AN DER STELLE X0. SIEHE AUCH EINGANG IV. // C // C METHODE: // C SIEHE HEFT 14 ODER 29 // C WENN YPARTR (AUS YFDFV) FUER FD UND/ODER FV WERTE VON NULL // C BZW. FUER FV AUCH VON EINS ERHAELT, KANN KEINE APPROXIMATION // C DURCHGEFUEHRT WERDEN. BEI IV=0 WERDEN DANN WILLKUERLICHE WERTE // C GESETZT, BEI IV>0 WIRD VERSUCHT DURCH VERSCHIEBEN DER STELLE X0 // C ZU EINER BRAUCHBAREN LOESUNG ZU GELANGEN. // C // C BENOETIGTE UP: YFDFV, YNINVF, YNORMD // C(**) // DIMENSION VP(10,N) // C int ISHIFT; int KOMP; double WEG; double FVS; bool klaar; // Om de breaks in de repeat te vervangen. WEG = 0.0; if ((IV == 2)) { EA = VP[0][I - 1]; SA = VP[1][I - 1]; } else { // 10 ISHIFT = 0; KOMP = 0; if ((IV != 0)) { YFDFV(N, I, IV, X0, VP, ref FD, ref FV, ref KOMP); } do { // 20 // 100 klaar = ((FV > 0.0) && (FV < 1.0) && (FD > 0.0)); if (!klaar) { ISHIFT ++; if ((IV == 0)) { if ((FV >= 1.0)) { KOMP = 1; } FD = FormCalculation.CEMIN2; FV = FormCalculation.CEMIN2; klaar = true; } // 110 if (!klaar) { if ((ISHIFT <= 1)) { WEG = VP[1][I - 1]/100.0; if ((KOMP == 1)) { WEG = -WEG; } } // 120 X0 = X0 + WEG; YFDFV(N, I, IV, X0, VP, ref FD, ref FV, ref KOMP); // =regel 20.. } } } while (!(klaar)); // 150 // if (ISHifT <> 0) then PRINT 1000,I, ISHifT; // (*//1000 // //FORMAT (1H ,7HPT: I :=,I2,4X,4HSH :=,I2) (**) // 200 FVS = YNINVF(FV); SA = YNORMD(FVS, 0.0, 1.0)/FD; FVS = SA*FVS; if ((KOMP != 0)) { FVS = -FVS; FV = 1.0 - FV; } EA = X0 - FVS; } } // *************************************************************** // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected void YPFEX(int IT, double TL, double FDI, double FVI, int KOMPI, ref double FD, ref double FV, ref int KOMP) { // (*C // C POISSONFOLGE VON UNABHAENGIGEN ZUFALLSVARIABLEN X // C MIT TL = (ZEITDAUER) * (AUFTRETENSRATE/ZEITEINHEIT) // C FDI = WERT DER DICHTE VON X AN DER STELLE X0 // C FVI = WERT DER VERTEILUNG VON X AN DER STELLE X0 // C = P ( X < X0 ) KOMPI = 0 // C = P ( X > X0 ) KOMPI = 1 // C // C IT STEUERT, OB DIE VERTEILUNG DES MAXIMUMS ODER DES MINIMUMS // C BERECHNET WIRD. // C // C IT = 1 VERTEILUNG DES MAXIMUMS // C // C FV = EXP ( - TL * ( 1 - FVI ) ) // C // C FD = TL * FDI * FV // C // C IT = 2 VERTEILUNG DES MINIMUMS // C // C FV = 1 - EXP ( - TL * FVI ) // C // C FD = TL * FDI * ( 1 - FV ) // C // C BENOETIGTES UP: YREIHE // C(**) int ITK; double F; double H; FD = 0.0; FV = 0.0; KOMP = 0; ITK = IT + KOMPI; F = FVI; if ((ITK != 2)) { F = 1.0 - F; } H = -TL*F; if ((H <= -0.1)) { if ((IT == 2)) { KOMP = 1; } if ((H > -FormCalculation.CEXP2)) { FV = Math.Exp(H); } } else { if ((IT == 1)) { KOMP = 1; } FV = -YREIHE(H); } F = FV; ITK = IT + KOMP; if ((ITK == 2)) { F = 1.0 - F; } FD = TL*FDI*F; } // ************************************************************* // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected double YPHIN(double X0) { double result; // (*C // C STANDARDNORMALVERTEILUNGSINTEGRAL // C // C YPHIN(X0) = P (X < X0) // C NACH ABRAMOWITZ/SEGUN FORMEL 26.2.17 // C(**) double R; double T; R = Math.Abs(X0); if ((R > FormCalculation.CEXP1)) { T = 0.0; } else { R = Math.Exp(-R*R/2.0)/2.506628275; T = 0.5; if ((X0 != 0)) { if ((X0 < 0)) { T = 1.0/(1.0 - 0.2316419*X0); } else { T = 1.0/(1.0 + 0.2316419*X0); } T = ((((1.330274429*T - 1.821255978)*T + 1.781477937)*T - 0.356563782)*T + 0.31938153)*T; T = R*T; } } // 8 R = T; if ((X0 > 0.0)) { R = 1.0 - T; } result = R; return result; } // Hoe F te declareren als parameter? 990920,Faber // F is een functie -> een pointer ofzo? Voorlopig double gemaakt om de compiler te foppen 990925,Faber protected void YREFA(TFFunctiePointerType F, ref double X, double X0, double EPS, int MAXIT, ref int IERror) { // (*C REGULA FALSI (ILLINOIS-ALGORITHMUS) // C // C F function F(X) // C X VARIABLE WIRD SO BESTIMMT, DASS F(X)=0 // C X0 ANFANGSLOESUNG // C EPS GENAUIGKEIT // C MAXIT MAXIMALE ZAHL DER ITERATIONEN // C IER FEHLERINDEX: 0 KEIN FEHLER // C 1 ABLEITUNG = 0 // C 2 KEINE KONVERGENZ // C(**) // EXTERNAL F // C int MI; int Ind; double X1; double X2; double XA; double FX1; double FX2; double FXA; double U; double DF; bool klaar; // om de break's te vervangen in de repeats. klaar = false; IERror = 0; MI = MAXIT; XA = X0; Ind = 0; FXA = F(XA); if ((Math.Abs(FXA) <= EPS)) { X = XA; // --> return } else { // 1 X1 = XA*1.1 + FormCalculation.CFIX2; FX1 = F(X1); if ((Math.Abs(FX1) <= EPS)) { X = X1; // --> return } else { // 2 if ((Math.Abs(FXA) <= Math.Abs(FX1))) { U = XA; XA = X1; X1 = U; U = FXA; FXA = FX1; FX1 = U; } do { // 3 X2 = X1; MI = MI - 1; if ((MI <= 0)) { IERror = 2; // 4 X = X2; klaar = true; } // 5 if (!klaar) { DF = FXA - FX1; while ((Math.Abs(DF) <= FormCalculation.COGE8)) { Ind = Ind + 1; if ((Ind >= 3)) { IERror = 1; X = X2; klaar = true; } if (!klaar) { // (*6*) X1 = X1*1.1 + FormCalculation.CFIX2; FX1 = F(X1); DF = FXA - FX1; } } // 10 if (!klaar) { X2 = X1 - (XA - X1)*FX1/DF; FX2 = F(X2); if ((Math.Abs(FX2) < EPS)) { X = X2; // ipv GOTO 4 klaar = true; } if (!klaar) { // C CONVERGENCE WITH EPS:=COGE27, DEFINED AS INPUT IN 'YREFA'. if (((FX1*FX2 <= 0.0) || (FXA*FX2 >= 0.0))) { XA = X1; FXA = FX1; // 15 } else { FXA = FXA/2.0; } // 20 X1 = X2; FX1 = FX2; } } } } while (!(klaar)); // GOTO 3 } } } // ******************************************************************** // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected double YREIHE(double Y) { double result; // (*C // C EVALUATION OF -(1-EXP(Y)) WITH A POWER-SERIES : // C Y + Y*Y/(1*2) + Y*Y*Y/(1*2*3) + .... // C // C(**) int I; double SE; double SG; bool klaar; // Gebruikt om de break in de until te omzeilen. klaar = false; I = 0; SE = 1.0; SG = Y; // C BELOW MACHINE PRECISION THERE IS NO NEED FOR ITERATION if ((Math.Abs(Y) >= FormCalculation.COGEN)) { SG = 0.0; do { I = I + 1; SE = SE*Y/I; SG = SG + SE; if ((I >= 100)) { klaar = true; } } while (!((Math.Abs(SE/SG) <= FormCalculation.COGE8) || klaar)); } result = SG; return result; } // **************************************************************** // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected void YSORT(int N, ref double[][] VP, ref int[] ISORT) { // (*C // C REIHENFOLGE DER VARIABLEN X[i-1] // C AUFSTEIGEND NACH WIEDERHOLUNGSZAHLEN VP(5,I) // C // C DABEI WERDEN ISORT(1) := INDEX I VON MIN (VP(5,I)) // C UND ISORT(N) = INDEX J VON MAX (VP(5,J)) // C BEI VP(5,I) = VP(5,J>I) WIRD ISORT(K) < ISORT(L>K) // C (**) // DIMENSION VP(10,N), ISORT(N) int I; int J; int Ind; double A; for (I = 1; I <= N; I ++) { VP[4][I - 1] = -VP[4][I - 1]; } Ind = 0; for (I = 1; I <= N; I ++) { A = 0.0; for (J = 1; J <= N; J ++) { if ((VP[4][J - 1] <= A)) { A = VP[4][J - 1]; Ind = J; } } VP[4][Ind - 1] = -VP[4][Ind - 1]; ISORT[N + 1 - I - 1] = Ind; } } // ************************************************************************* // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected void YTXZUY(int N, int[] IV, double[] EX, double[] SX, ref double[] EY, ref double[] SY, ref double[][] VP, double[] X, ref double[][] COV, int NC, int IRHO) { // (*C // C ANFANGSTRANSFORMATIONEN FUER FORM UND SORM // C // C BESTIMMUNG DER PARAMETER DER VARIABLEN // C Y := X X NORMALVERTEILT // C Y := DLOG(X) X LOGNORMALVERTEILT // C Y := INVERSE(PHI(F(X))) SONST // C HIERZU ANFANGSLOESUNG FUER X ERFORDERLICH // C // C EINGANG: // C N DIMENSION // C IV VEKTOR DER VERTEILUNGSTYPEN DER X // C EX VEKTOR DER ERWARTUNGSWERTE DER X // C SX VEKTOR DER STANDARDABWEICHUNGEN DER X // C X ANFANGSLOESUNG (ERFORDELICH, WENN X NICHT // C NORMAL ODER LOGNORMAL VERTEILT IST) // C COV MATRIX DER KORRELATIONSZifFERN DER X[i-1], X[j-1] // C BESETZT MUSS COV[i,j-1] MIT I>:=J SEIN // C NC DIMENSION DER QUADRATISCHEN MATRIX COV // C IRHO KORRELATIONSINDEX // C // C AUSGANG: // C EY VEKTOR DER ERWARTUNGSWERTE DER Y // C SY VEKTOR DER STANDARDABWEICHUNGEN DER Y // C COV MATRIX DER KOVARIANZEN VON Y[i-1],Y[j-1] // C BESETZT WIRD DIE GANZE MATRIX // C VP VERTEILUNGSPARAMETER VON X // C // C ENTSPRECHendE RUECKTRANSFORMATION VON COV DURCH YTYZUX // C // C BENOETIGTE UP: YPARB, YPARTR // C(**) // DIMENSION EX(N), SX(N), EY(N), SY(N), VP(10,N), // 2 X(N), COV(NC,NC), IV(N) int i = 0; int im1 = 0; int j = 0; int IT = 0; int JT = 0; double FD = 0; double FV = 0; double COVIJ = 0; double H = 0; for (i = 1; i <= N; i ++) { IT = IV[i - 1]; EY[i - 1] = EX[i - 1]; SY[i - 1] = 0.0; if ((IT != 0)) { SY[i - 1] = SX[i - 1]; YPARB(N, i, IT, EX[i - 1], SX[i - 1], ref VP); if (((IT == 3) || (IT == 10))) { // (*10*) SY[i - 1] = VP[1][i - 1]; EY[i - 1] = Math.Log(Convert.ToDouble(EX[i - 1] - VP[2][i - 1])) - 0.5*Math.Pow(Convert.ToDouble(SY[i - 1]), 2.0); } else { if ((IT != 2)) { YPARTR(N, i, IT, X[i - 1], VP, ref EY[i - 1], ref SY[i - 1], ref FD, ref FV); } } } } // 20 if ((IRHO != 0)) { // C UMWandLUNG DER KORRELATIONSZifFERN IN KOVARIANZEN COV[0][0] = SY[0]*SY[0]; for (i = 2; i <= N; i ++) { im1 = i - 1; COV[i - 1][i - 1] = SY[i - 1]*SY[i - 1]; for (j = 1; j <= im1; j ++) { if (((IV[i - 1] != 0) && (IV[j - 1] != 0))) { COV[i - 1][j - 1] = COV[i - 1][j - 1]*SX[i - 1]*SX[j - 1]; } } } // 30 // C // C LOGNORMALTRANSFORMATION VON COV for (i = 2; i <= N; i ++) { im1 = i - 1; IT = IV[i - 1]; for (j = 1; j <= im1; j ++) { JT = IV[j - 1]; COVIJ = COV[i - 1][j - 1]; if (((IT == 2) || (IT == 3) || (IT == 10))) { if (((JT == 2) || (JT == 3) || (JT == 10))) { if (((IT != 2) || (JT != 2))) { if (((IT != 2) && (JT != 2))) { // COVIJ := Ln(1.+COVIJ/((EX[i-1]-VP[3,I-1])*(EX[j-1]-VP[3,J-1]))); COVIJ = Math.Log(1.0 + COVIJ/((EX[i - 1] - VP[2][i - 1])*(EX[j - 1] - VP[2][j - 1]))); } else { // (*35*) if ((IT != 2)) { H = EX[i - 1] - VP[2][i - 1]; // (*38*) COVIJ = COVIJ/H; } else { // 37 H = EX[j - 1] - VP[2][j - 1]; // 38 COVIJ = COVIJ/H; } } } } } // 40 COV[i - 1][j - 1] = COVIJ; // 45 COV[j - 1][i - 1] = COVIJ; } // 50 } } } // ************************************************************* // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected void YTYZUX(int N, int[] IV, double[] SX, double[] EY, double[] SY, ref double[][] COV, int NC) { // (*C // C RUECKTRANSFORMATION DER MATRIX COV // C // C SIEHE BESCHREIBUNG DER procedure TDForm4. YTXZUY // C // C EINGANG: N, IV, SX, EY, SY, COV, NC // C // C AUSGANG: COV // C BESETZT WIRD NUR DER UNTERE TEIL MIT COV[i,j-1] (I>:=J) // C(**) // DIMENSION SX(N), EY(N), SY(N), COV(NC,NC), IV(N) // C int I; int IT; int IM1; int J; int JT; double H; double COVIJ; double SXIJ; for (I = 2; I <= N; I ++) { IT = IV[I - 1]; if (((IT == 2) || (IT == 3) || (IT == 10))) { IM1 = I - 1; for (J = 1; J <= IM1; J ++) { JT = IV[J - 1]; if (((JT == 2) || (JT == 3) || (JT == 10))) { if (((IT != 2) || (JT != 2))) { if (((IT != 2) && (JT != 2))) { COV[I - 1][J - 1] = (Math.Exp(Convert.ToDouble(COV[I - 1][J - 1])) - 1.0)*Math.Exp(EY[I - 1] + EY[J - 1] + 0.5*(Math.Pow(Convert.ToDouble(SY[I - 1]), 2.0) + Math.Pow(Convert.ToDouble(SY[J - 1]), 2.0))); } else { H = EY[I - 1] + 0.5*Math.Pow(Convert.ToDouble(SY[I - 1]), 2.0); if ((IT == 2)) { H = EY[J - 1] + 0.5*Math.Pow(Convert.ToDouble(SY[J - 1]), 2.0); } COV[I - 1][J - 1] = COV[I - 1][J - 1]*Math.Exp(H); } } } } } } COV[0][0] = 1.0; for (I = 2; I <= N; I ++) { COV[I - 1][I - 1] = 1.0; IM1 = I - 1; for (J = 1; J <= IM1; J ++) { COVIJ = 0.0; SXIJ = SX[I - 1]*SX[J - 1]; if ((SXIJ > FormCalculation.COGE14)) { COVIJ = COV[I - 1][J - 1]/SXIJ; } COV[I - 1][J - 1] = COVIJ; } } } // ********************************************************************** // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected void YUNI(double X, double A, double B, ref double FD, ref double FV, ref int KOMP) { // (*C // C R E C H T E C K V E R T E I L U N G // C // C X IST ZWISCHEN A UND B RECHTECKVERTEILT // C // C 1 / ( B - A ) A <= X <= B // C FD = // C 0 SONST // C // C 0 A >= X // C FV = ( X - A ) / ( B - A ) A <= X <= B // C 1 X >= B // C // C EINGANG: X, A, B // C AUSGANG: FD, FV, KOMP // C(**) double H; KOMP = 0; FV = 0.0; if ((X > A)) { if ((X >= B)) { KOMP = 1; } else { H = A + (B - A)/2.0; if ((X > H)) { KOMP = 1; FV = (B - X)/(B - A); } else { FV = (X - A)/(B - A); } } } FD = 0.0; if ((X >= A)) { if ((X <= B)) { FD = 1.0/(B - A); } } } // afgeknotte normale verdeling // ************************************************************************ // ********************************************************************** // WAT TE DOEN MET ONDERSTAANDE DUMMY'S? // Voorlopig heb ik ze niet overgenomen naar Pascal,19990923,HenkJanFaber. // ********************************************************************** // (* ->onderstaande procedures in commentaarregels (worden aangeroepen in YFDFV) 19990925,HenkJanFaber // // {======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // =======================================================================================================================} // procedure TDForm4.YV96 (X0, P1, P2, P3, P4, FD, FV, KOMP) // C DUMMY FUER TYPENKENNUNG 96 // * *) // *********************************************************** // * // * PURPOSE: // * ------- // * Bepaal voor een afgekapte normale verdeling de // * verdelings functie en kansdichtheids functie // * // * INPUT: // * ----- // * X0 - Actuele waarde van de parameter // * P1 - Verwachtings waarde van de verdeling // * P2 - Standaard afwijking van de verdeling // * P3 - De waarde warop de verdeling is afgekapt // * P4 - Dummy variabele // * KOMP - Geeft aan of het om overschrijdings // * dan wel onderschrijdings kansen gaat // * // * OUTPUT: // * ------ // * FD - Kans dichtheids functie // * FV - Verdelings functie // * // * NOTE: // * ----- // * // * PROGRAMMER: Han Best // * ---------- // * LAST MODIFIED: 14 maart 1989 P // * ------------- // ********************************************************************* protected void YV96(double X0, double P1, double P2, double P3, double P4, ref double FD, ref double FV, ref int KOMP) { double Hulp = 0; double Dum = 0; if ((X0 <= P3)) { FD = 0.0; FV = 0.0; } else { YNORM(P3, P1, P2, ref Dum, ref Hulp, ref KOMP); if ((KOMP == 1)) { Hulp = 1 - Hulp; } YNORM(X0, P1, P2, ref FD, ref FV, ref KOMP); if ((KOMP == 1)) { FV = 1 - FV; } FD = FD/(1.0 - Hulp); FV = (FV - Hulp)/(1.0 - Hulp); } } // (* procedure TDForm4.YV97 (X0, P1, P2, P3, P4, FD, FV, KOMP) // C DUMMY FUER TYPENKENNUNG 97 // STOP 141 // end // // procedure TDForm4.YV98 (X0, P1, P2, P3, P4, FD, FV, KOMP) // C DUMMY FUER TYPENKENNUNG 98 // STOP 142 // end // // procedure TDForm4.YV99 (X0, P1, P2, P3, P4, FD, FV, KOMP) // C DUMMY FUER TYPENKENNUNG 99 // STOP 143 // end // // // // // *) // ********************************************************************** // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected double YVGAM(double RK) { double result; // (*C // C FUNKTIONSUNTERPROGRAMM ZUR ERMITTLUNG VON VERTEILUNGSPARAMETERN // C IN YPARB BEI VERTEILUNGSTYPEN 8 UND 9 // C WIRD DURCH YREFA AUFGERUFEN // C // C BENOETIGTE UP: GAMMA AUS BIBLIOTHEK IMSL // C (VOLLST. GAMMAFUNKTION) // C(**) // COMMON /FREWEI/ VAR, IVC double H; // Deze functie wordt aangeroepen vanuit YREFA, die op haar beurt // weer wordt aangeroepen door YPARB. Daar wordt ervoor gezorgd // dat IVC en Variantie correcte waarden krijgen; de compiler geeft // natuurlijk wel een warning. 19990928, HenkJanFaber. H = 1.0/RK; if ((IVC == 8)) { H = -H; } result = Variantie*Variantie - Gamma(1.0 + 2.0*H)/Math.Pow(Gamma(1.0 + H), 2.0) + 1.0; return result; } // ********************************************************************** // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected void YWEIB(double X, double U, double RK, double EPS, ref double FD, ref double FV, ref int KOMP) { // (*C // C TYP III W E I B U L L KLEINSTWERTE // C // C FV = 1 - EXP ( - ( (X-EPS) / (U-EPS) ) ** RK ) X > EPS // C // C FD = RK / (U-EPS) * ( (X-EPS) / (U-EPS) ) ** ( RK - 1 ) * ( 1 - FV ) // C // C BENOETIGTES UP: YREIHE // C(**) double H; double V; double Y; KOMP = 0; FD = 0.0; FV = 0.0; Y = X - EPS; V = U - EPS; if ((Y > 0.0)) { KOMP = 1; H = Math.Pow((Y/V), RK); if ((H <= FormCalculation.CEXP2)) { FV = Math.Exp(-H); FD = RK*H*FV/Y; if ((X < U)) { KOMP = 0; FV = -YREIHE(-H); } } } } // procedure YPRUEF (IRHO: Integer); // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected void YPRUEF(int IRHO, int NM, int NC, int[] IV, double[][] VP, double[][] COV, double[] SX, double[] EX, ref int Ierror) { // PRUEFUNG DES STOCHASTISCHEN MODELLS // // AUSGANG IST NUR DER FEHLERPARAMETER IER // BEDEUTUNG VON IER SIEHE HEFT 43 / ANHANG D // Slechte routine. M.i. beter opzetten met case statements // Te generiek opgezet, terwijl de uitbreidingen er niet zijn/ frank int V; int IC; int IVS; int I; int IM1; int J; double VIJ; Ierror = 0; IVS = 0; for (I = 1; I <= NM; I ++) { if (IV[I - 1] > 0) { // TYPENVORRAT : 0 - 10, 96 - 99, // 101 - 110, 196 - 199, // 201 - 210, 296 - 299, // 302, 303, 310, 402, 403, 410, // 501 - 510, 596 - 599, // 601 - 610, 696 - 699 // V = IV[I - 1]; V = V - Convert.ToInt32(Math.Floor(V/100d))*100; // V bevat de tientallen IC = IV[I - 1] - V; // IC bevat de honderdtallen IVS ++; // aantal variabelen met verdeling if ((IC > 0) || (V < 96)) { // IC=0 and 96<=V<=99 ? -> naar 99. Ierror = FormCalculation.YPruefcheck(IC, V); // 1100 // Standard afwijking moet positief zijn if ((SX[I - 1] <= 0.0)) { Ierror = 112000; // 1120 } // WIEDERHOLUNGSZAHL if ((VP[4][I - 1] <= 0.0) || (VP[4][I - 1] > 1.0E+6)) { Ierror = 121000; // 1210 } // AUFTRETENSWAHRSCHEINLICHKEIT if ((VP[5][I - 1] < 0.0) || (VP[5][I - 1] > 1.0)) { Ierror = 122000; // 1220 } // BINOMIALFOLGE : WIEDERHOLUNGSZAHL if ((IC == 100) || (IC == 200)) { if ((VP[6][I - 1] <= 0.0) || (VP[6][I - 1] > 1.0E+6)) { Ierror = 123000; } // BINOMIALFOLGE : AUFTRETENSWAHRSCHEINLICHKEIT if ((VP[7][I - 1] < 0.0) || (VP[7][I - 1] > 1.0)) { Ierror = 124000; // 1240 } } Ierror = FormCalculation.verwachtingCheck(I, EX, VP, V); // 1130 } } } // aanwezige stochastische variabelen if ((IVS == 0)) { Ierror = 131000; } // Korrelatie STuurparameter if ((IRHO > 1) || (IRHO < 0)) { // WriteLn (lferr,'IRHO = ',IRHO,' unzul�ssig.'); // 1410 Ierror = 141000; // -> RETURN } else { if ((IRHO == 1)) { if ((NM > NC)) { // 1420 Ierror = 142000; // -> RETURN } else { // KORRELATIONS-/ KOVARIANZENMATRIX // ERLAUBTE BEZIEHUNGEN UNTEREIN and ER // 500 for (I = 1; I <= NM; I ++) { if ((IV[I - 1] == 2) || (IV[I - 1] == 3) || (IV[I - 1] == 10)) { FIH[I - 1] = 0; } else { FIH[I - 1] = 1; } } for (I = 1; I <= (NM - 1); I ++) { IM1 = I - 1; for (J = 1; J <= IM1; J ++) { VIJ = VP[4][I - 1]/VP[4][J - 1]; if (((Math.Abs(VIJ - 1.0) > FormCalculation.COGE14) || ((FIH[I - 1] != 0) || (FIH[J - 1] != 0)))) { // 530 if ((Math.Abs(COV[I - 1][J - 1]) >= FormCalculation.COGE14)) { // 1440 Ierror = 144000; } } } } } } } } // Error numer // *********************************************************** // * // * PURPOSE: // * ------- // * Deze subroutine berekent de verwachtingswaarde en // * standaard afwijking van de op a afgeknotte normale verdeling // * // * INPUT: // * ----- // * AMU - Verwachtings waarde van de verdeling // * ZIGMA - Standaard afwijking van de verdeling // * A - De waarde waarop de verdeling is afgekapt // * // * OUTPUT: // * ------ // * VERW - Verwachtingswaarde van de afgeknotte verdeling // * STAN - Standaard afwijking van de afgeknotte verdeling // * // * NOTE: // * ----- // * // * PROGRAMMER: Han Best // * ---------- // * LAST MODIFIED: 14 maart 1989 P // written in Pascal 1999-11-05 best // ********************************************************************* // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected double signUnit(double value) { double result; result = 1.0; if (value < 0) { result = -1.0; } return result; } // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected double KdfVariabele(double Average, double Sigma, double Realisatie, int verdeling) { double result; double Value; Value = 1; if (Sigma > 0.001) { if (verdeling == 2) { // Normale verdeling Value = Math.Exp(-0.5*((Realisatie - Average)/Sigma)*((Realisatie - Average)/Sigma))/(Sigma*Math.Sqrt(2.0*Math.PI)); } else if (verdeling == 3) { // lognormale verdeling Value = Math.Exp(-0.5*((Math.Log(Realisatie) - Average)/Sigma)*((Math.Log(Realisatie) - Average)/Sigma))/(Sigma*Realisatie*Math.Sqrt(2.0*Math.PI)); } else if (verdeling == 7) { // Gumbel verdeling Value = Sigma*Math.Exp(-Sigma*(Realisatie - Average) - Math.Exp(-Sigma*(Realisatie - Average))); } } result = Value; return result; } // ======================================================================================================================= // Description: Clculates the combinded density of the variables // // Date ID Modification // 1999-10-01 Best Created // ======================================================================================================================= protected double CombinedKdf(int Nm, double[] Ex, double[] Sx, double[] X, int[] IV) { double result; double Value; double val1; int i; // Var M : wvect; Value = 1; // gezamelijke kansdichtheid k // (* If ( abs(Dk) > 0.001 ) then // begin // With Kl do // M[1] := Ln(realisatie) - average; // With Kr do // M[2] := Ln(realisatie) - average; // Value := 1/ (Sqrt(Dk)*(2 * pi) ) * Exp( -0.5 * MatrixMult(InvCorMatK,2,M)); // end // else // Value := KdfVariabele(Kl); // // { gezamelijke kansdichtheid d } // If ( abs(Dd) >0.001 ) then // begin // With Dl do // M[1] := realisatie - average; // With Dr do // M[2] := realisatie - average; // Val1 := // 1/(Sqrt(Dd)*(2 * pi) ) * Exp( -0.5 * MatrixMult(InvCorMatd,2,M)); // end // else *) for (i = 0; i < Nm; i ++) { val1 = KdfVariabele(Ex[i], Sx[i], X[i], IV[i]); Value = Value*val1; } result = Value; return result; } // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= protected double getAlfa(int Verdeling, double Average, double sigma, double Ontwerp, double Beta) { double result; // kdf gumbel // verdelings functie gumbel // kdf gumbel // verdelings functie gumbel // (x - mu )/sigma genorm gumbel // (x - mu )/sigma genorm Lognorm verdeling // Sigma van genorm. gumbel in ontwerppunt double pfgumbel; double pgumbel; double pfLogNorm; double MuLog; double BetaGum; double BetaLog; double tussig; double MuGumbel; // Mu van genorm. gumbel in ontwerppunt result = 0; if (sigma > 0.001) { if (Verdeling == 2) { // normal result = (Average - Ontwerp)/(Beta*sigma); } else if (Verdeling == 3) { // lognormal pfLogNorm = 1/(Math.Sqrt(2*Math.PI)*sigma*Ontwerp)*Math.Exp(-0.5*((Average - Math.Log(Ontwerp))*(Average - Math.Log(Ontwerp))/(sigma*sigma))); BetaLog = (Math.Log(Ontwerp) - Average)/sigma; tussig = 1/(Math.Sqrt(2*Math.PI)*pfLogNorm)*Math.Exp(-0.5*BetaLog*BetaLog); MuLog = Ontwerp - tussig*BetaLog; result = (MuLog - Ontwerp)/(Beta*tussig); } else if (Verdeling == 7) { pfgumbel = sigma*Math.Exp(-sigma*(Ontwerp - Average) - Math.Exp(-sigma*(Ontwerp - Average))); pgumbel = Math.Exp(-Math.Exp(-sigma*(Ontwerp - Average))); BetaGum = NormalDistrInverse(pgumbel); tussig = 1/(Math.Sqrt(2*Math.PI)*pfgumbel)*Math.Exp(-0.5*BetaGum*BetaGum); MuGumbel = Ontwerp - tussig*BetaGum; result = (MuGumbel - Ontwerp)/(Beta*tussig); } } return result; } // ======================================================================================================================= // Description: Main procedure of a monte carlo simulation // // Date ID Modification // 1999-10-01 Best Created // ======================================================================================================================= } // end TDForm4 } namespace Deltares.Stability.Calculation.Inner { public class FormCalculation { // uses // MSTCalcProgressDlg; // CEMIN = 0.124E-305; // CEMAX = 0.695E293; // CEMAX4 = 0.162E74; // COGEN = 0.111E-15; // COMPI = 3.14159265358979; {=PI} // CFIX2 = 1.E-02; // CFIX4 = 1.E-04; // CEXP1 = 37.5; // CEXP2 = 704.0; // CEXP3 = 6.56; // CEMIN2 = 0.112E-152; // CEMIN3 = 0.108E-101; // CONE1 = 8.44; // COGE27 = 0.276E-4; // COGE3 = 0.481E-5; // COGE23 = 0.231E-10; // COGE8 = 0.172E-12; // COGE14 = 0.461E-22; public const double CEMIN = 0.794E-36; public const double CEMAX = 0.208E37; public const double CEMAX4 = 0.120E10; public const double COGEN = 0.218E-10; public const double COMPI = 3.1415926536; public const double CFIX2 = 1.0E-02; public const double CFIX4 = 1.0E-04; public const double CEXP1 = 12.9; public const double CEXP2 = 83.1; public const double CEXP3 = 4.42; public const double CEMIN2 = 0.819E-18; public const double CEMIN3 = 0.926E-12; public const double CONE1 = 6.85; public const double COGE27 = 0.900E-03; public const double COGE3 = 0.279E-03; public const double COGE23 = 0.781E-7; public const double COGE8 = 0.296E-8; public const double COGE14 = 0.119E-14; // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= public static int YPruefcheck(int Ic, int V) { int result; int IError; IError = 111000; if ((Ic == 300) || (Ic == 400)) { switch (V) { case 2: IError = 0; break; case 3: IError = 0; break; case 10: IError = 0; break; } } else if ((Ic < 700)) { if ((V >= 1) && (V <= 10)) { IError = 0; } } else { switch (V) { case 96: IError = 0; break; case 97: IError = 0; break; case 98: IError = 0; break; case 99: IError = 0; break; } } result = IError; return result; } // ======================================================================================================================= // Date ID Modification // 2008-01-16 Best Created // ======================================================================================================================= public static int verwachtingCheck(int i, double[] ex, double[][] vp, int V) { int result; int ierror; ierror = 113000; switch (V) { case 1: ierror = 0; break; case 2: ierror = 0; break; case 3: if (ex[i - 1] > 0) { ierror = 0; } break; case 4: ierror = 0; break; case 5: if (ex[i - 1] > 0) { ierror = 0; } break; case 6: if ((ex[i - 1] > vp[2][i - 1]) && (ex[i - 1] <= vp[3][i - 1])) { ierror = 0; } break; case 7: ierror = 0; break; case 8: case 9: case 10: if ((ex[i - 1] > vp[2][i - 1])) { ierror = 0; } break; case 96: ierror = 0; break; case 97: ierror = 0; break; case 98: ierror = 0; break; case 99: ierror = 0; break; } result = ierror; return result; } } // end FormCalculation } // **********************************************************************