Option Explicit 
 '
 '   Module Description - Triangular Routines
 '
 '   Provides Three double precision VBA (believed generally VB capable)
 '   routines dealing with the Triangular Distribution. Main use
 '   is expected to be in Monte Carlo Simulation
 '
 '   TrigenVB(dMin,dMl,dMax,dP1,dP2,RandNo)
 '   TrigenSizeVB(dMin,dMl,dMax,dP1,dP2,dLower,dUpper)
 '   TriangVB(dMin, dMl, dMax, dPVal)
 '
 '   Coded by UlricTheRed 21/10/2009
 '
 
Function TrigenVB(dMin, dMl, dMax, dPP1, dPP2, dRandNo) 
     '
     '   Returns a random number from the Trigen Distribution
     '   used in @Risk (?and elsewhere)
     '   Does this by calling the TrigenSize Subroutine which
     '   calculates  the dLower and dUpper bounds of a new triangle
     '   with the same most likely value as the original distribution
     '   but dPP1 of the weight to the left of dMin, and dPP2 of the
     '   weight to the left of dMax.
     '
     '
    Dim dLower, dUpper As Double ' Variables to hold "true" Upper and Lower limits _
     ' of triangular Distribution
    Call TrigenSize(dMin, dMl, dMax, dPP1, dPP2, dLower, dUpper) 
    TrigenVB = TriangVB(dLower, dMl, dUpper, dRandNo) 
End Function 
 
 
Sub TrigenSize(dMin, dMl, dMax, dPP1, dPP2, dLower, dUpper) 
    Dim dA, dB, dC, dP1, dP2, dL, dU, dM, dQ, dR, dN, dSign As Double 
    Dim dBracket1, dSqrtBracket1, dBracket2, dF1n, dF2n, dF1nDash, dF2nDash As Double 
    Dim dNplusQ, dOneMinusP1, dCorr, dH, dH1, dH2 As Double 
    Dim iKtr As Integer 
    Const sOne As Single = 1# 
     
     '
     '   Calculates  the lower and upper bounds of a new triangle
     '   with the same most likely value as the original distribution
     '   but pp1 of the weight to the left of dMin, and pp2 of the
     '   weight to the left of dMax.
     
    dA = dMin 
    dB = dMax 
    dC = dMl 
    dP1 = dPP1 / 100# 
    dP2 = sOne - (dPP2 / 100) 
    dSign = sOne 
    dR = dB - dA 
    dQ = dC - dA 
     'If ((p1 + p2) = 0) Then GoTo endhere
     'SImple Data Checks
    If ((dA > dC) Or (dC > dB)) Then 
        MsgBox ("Check Your Entry data it is inconsistent!") 
        GoTo endhere 
    End If 
     '
     ' If p1 = 0 then m becomes the origin
     '
    If (dP1 = 0) Then 
        dM = 0# 
        If (dP2 = 0) Then 
            dN = dR 
            GoTo closeenough 
        End If 
        dBracket1 = (2 * dR - dQ * dP2) 
        dN = (dBracket1 + Sqr(dBracket1 * dBracket1 - 4# * (sOne - dP2) * dR * dR)) / (2# * (sOne - dP2)) 
        GoTo closeenough 
    End If 
     '
     ' If p2 = 0 then n=r
     '
    If (dP2 = 0) Then 
        dN = dR 
        dBracket1 = dP1 * (dQ + dR) 
        dM = (dBracket1 + Sqr((dBracket1 * dBracket1) + 4 * (sOne - dP1) * dQ * dP1 * dR)) / (2# * (sOne - dP1)) 
        GoTo closeenough 
    End If 
     '
     '
     'Carry out transformation - move dMin estimate to origin
     'i.e X=x-a
     '
     ' Generate an initial guess fo n before using
     ' Newton-Raphson iterative technique as in detailed description
     '
    dN = dR * (sOne + 4 * dP2) 
     '
     '   Begin Iteration
     '
    For iKtr = 1 To 50 Step 1 
        dNplusQ = dN + dQ 
        dOneMinusP1 = sOne - dP1 
        dBracket1 = dP1 * dP1 * dNplusQ * dNplusQ + 4 * dOneMinusP1 * dP1 * dN * dQ 
        If (dBracket1 < 0) Then 
            MsgBox ("Too Much Weight in the Tails") 
            dM = 0 
            dN = dR 
            GoTo closeenough 
        End If 
        dSqrtBracket1 = Sqr(dBracket1) 
        dF2n = ((dP1 * dNplusQ) + dSign * dSqrtBracket1) / (2# * dOneMinusP1) 
        dBracket2 = (2 * dR + dP2 * dF2n - dP2 * dQ) 
        dF1n = dN * dN * (sOne - dP2) - dN * dBracket2 + (dR * dR + dP2 * dQ * dF2n) 
        dF2nDash = (dP1 + (0.5 * dSign * (2 * dP1 * dP1 * (dN + dQ) + 4# * dOneMinusP1 * dP1 * dQ) / dSqrtBracket1)) / (2# * dOneMinusP1) 
        dF1nDash = (2 * dN * (sOne - dP2)) - (dN * dP2 * dF2nDash) - (dBracket2) + (dP2 * dQ * dF2nDash) 
        dCorr = dF1n / dF1nDash 
        dN = dN - dCorr 
        dM = dF2n 
         '
         '   If correction is small enough stop iterating
         '
        If (Abs(dCorr / dN) < 0.00000001) Then GoTo closeenough 
    Next iKtr 
     '
     '   If we get here iteration hasn't converged
     '
    MsgBox ("No convergence") 
closeenough: 
    dU = dN + dA 
    dL = -dM + dA 
     '
     '   Return to original co-ordinates
     '
    dH = 2 / (dU - dL) 
    dH1 = dH * (dA - dL) / (dC - dL) 
    dH2 = dH * (dU - dB) / (dU - dC) 
    If (dH1 < 0 Or dH2 < 0) Then 
        MsgBox ("Too Much Weight in the Tails") 
        dU = dB 
        dL = dA 
    End If 
    dLower = dL 
    dUpper = dU 
endhere: 
End Sub 
 
 
Function TriangVB(dMin, dMl, dMax, dPval) As Double 
    Dim dA, dB, dC, dP, dBminA, dCminA, dBminC, dPminl As Double 
     '
     '  This function returns thd Inverse CDF - i.e
     '  the "x" value associated with a cumulative
     '  distribution probability of p. If p is
     '  chosen at random from a uniform [0,1]
     '  then it returns a random number from the Triangular Distribution
     '
    dA = dMin 
    dB = dMax 
    dC = dMl 
    dP = dPval 
     '
     '  Check i/p in order
     '
    If ((dB < dC) Or (dC < dA) Or (dP > 1#) Or (dP < 0#)) Then 
        MsgBox ("Error1 in i/p data") 
        TriangVB = 0# 
        GoTo endroutine 
    End If 
    dCminA = dC - dA 
    dBminA = dB - dA 
    dBminC = dB - dC 
    dPminl = dCminA / dBminA 
     '
     ' Now check which side of b we fall into
     ' pml is the Cumulative Probability at ml
     ' value (knowing overall area =1)
     '
    If (dP = 0) Then 
        TriangVB = dA 
    ElseIf (dP = 1#) Then TriangVB = dB 
    ElseIf (dP < dPminl) Then TriangVB = dA + ((dP * dBminA * dCminA) ^ 0.5) 
    Else 
        TriangVB = dB - (((1# - dP) * dBminA * dBminC) ^ 0.5) 
    End If 
endroutine: 
End Function