Programming Issue

[greentangerine]

Hi All,

Have programmed this code for a water security issue.
The prob.mod was supplied to use along with our own program.

nrepmax needs to run for 500 000 replicates. It works for around 100 000 replicates, but after that there is an error from the prod.mod code. I have highlighted the error message that comes up, not sure if the highlight works.

Thanks in Advance!

Imports System.Math
Module warragamba

    Sub main()
        ' declare variables
        Dim nrep, nyear, nrepmax As Integer      ' loop variables
        Dim demand, volume, capacity, area, evap, VolumeNew As Double
        Dim restrict, desal, empty As Integer         ' counters
        Dim Zt, Z_old, lambda, mean, stdev, correction As Double
        Dim percentfull As Double
        Dim Prob_restrict, Prob_desal, Prob_empty As Double

        ' initial values
        demand = 400000
        capacity = 1886000
        nrepmax = [COLOR="red"]100[/COLOR]
        'random variable values
        lambda = 0.1
        mean = 29.277
        stdev = 3.589
        correction = 0.2184
        ' replicate loop
        For nrep = 1 To nrepmax
            ' initial loop values
            volume = 0.5 * capacity
            restrict = 0
            desal = 0
            empty = 0
            ' year loop
            For nyear = 1 To 40
                ' increase demand each year (population increase)
                demand = demand ^ (1.02)
                ' restrictions check
                If volume < 0.45 * capacity Then
                    demand = 0.75 * 400000
                    restrict = restrict + 1
                Else
                    demand = 400000
                End If
                ' Find surface area
                area = 17.1 + 57.9 * (volume / 1886000) ^ 0.7555    'm2
                ' Find evaporation volume
                evap = (1530 / 1000) * 0.7 * area       ' m3
                ' Find annual stream inflow into warragamba
                Zt = mean + correction * (Z_old - mean) + normal(0.0, stdev)
                Qinflow = (Zt * lambda + 1) ^ (1 / lambda)
                Zold = Zt ' for next loop
                ' Calculate volume of reservoir
                VolumeNew = volume - demand - evap + Qinflow
                ' Desalination plant check
                If VolumeNew <= 0.3 * capacity Then
                    'desalination plant built in 2 years
                    For nyear = nyear + 2 To 40
                        VolumeNew = VolumeNew + (250 * 365)
                    Next nyear
                    desal = desal + 1
                End If
                ' Make sure vol < cap
                If VolumeNew > capacity Then
                    VolumeNew = capacity
                ElseIf VolumeNew <= capacity Then
                    VolumeNew = VolumeNew
                End If
                ' Find % full
                percentfull = (VolumeNew / capacity) * 100
                ' count # times empty
                If percentfull = 0 Then
                    empty = empty + 1
                End If
                ' update volume for next loop
                volume = VolumeNew
            Next nyear
        Next nrep
        ' calculate probabilities
        ' restrictions
        Prob_restrict = (restrict / CDbl(nrepmax)) * 100
        ' desalination
        Prob_desal = (desal / CDbl(nrepmax)) * 100
        ' empty
        Prob_empty = (empty / CDbl(nrepmax)) * 100
        ' output
        ' probability of empty
        If Prob_empty > 0.1 Then
            Console.WriteLine("The probability of Warragamba Dam being empty is: " & Prob_empty)
            Console.WriteLine("Therefore, it will be necessary to augment the system with a new low cost source of water")
        Else
            Console.WriteLine("The probability of Warragamba Dam being empty is: " & Prob_empty)
            Console.WriteLine("The water supply is sufficient")
        End If

        'probability of restrictions
        Console.WriteLine("The probability of Sydney having water restrictions is: " & Prob_restrict)
        'probability of desal
        Console.WriteLine("The probability of Sydney needing a Desalination Plant to be built is: " & Prob_restrict)
        Console.ReadKey()
    End Sub

End Module

Module probMod
        '
        ' Module for sampling from probability distributions
        ' The following module procedures are public:
        '   uniRan          --> uniRan ()
        '                       Function returns a random number sampled from a uniform distribution with limits 0 to 1
        '   normal          --> normal (mean as double, stdDev as double) as double
        '                       Function returns a random number or deviate sampled from a normal
        '                       distribution with Expected value mean and standard devaition stdDev
        '   poidev          --> poidev (mean as double) as long
        '                       Function returns a random number or deviate sampled from a Possion
        '                       distribution with expected value equal to the variable mean
        '   createHistTable --> CreateHistTable(x() As Double, nX As Long, nBins As Long, firstBin As Double, _
        '                                       incrBin As Double, densHist As Boolean, histBins() As String, _
        '                                       binCounts() As Double, errorFlag As Boolean, errorMessage As String)
        '                       Subroutine generates a histogram table which may be printed in the calling
        '                       code. The sub supports frequency histograms and density histograms. See agrument
        '                       list below:
        '                        INPUTS:
        '                         x = samples to construct histogram from
        '                             Double type array ()
        '                         nX = number of samples to use in histogram calcs
        '                             Long type scalar
        '                         nBins = number of histogram bins (columns/bars in plot)
        '                             Long type scalar
        '                             **** NOTE: histogram bin and bin count lists are of length nBins+1,
        '                                         i.e. have max index nBins and min index 0
        '                         firstBin = smallest value in historagm bin list
        '                             Double type scalar
        '                         incrBin = bin interval or width
        '                             Double type scalar
        '                         densHist = create density historgram if true
        '                             Boolean type
        '                        OUTPUTS:
        '                         histBins = list of histogram bin labels
        '                             String type array
        '                         binCounts = returns frequency or density of each histogram bin
        '                             Double type array
        '                         errorFlag = Boolean type, TRUE if error has occured
        '                         errorMessage = string giving warning or error message

        Function normal(mean As Double, stdDev As Double) As Double
            '----
            '
            If stdDev <= 0.0 Then
                Call fatalError("f-normal/Negative standard deviation")
            Else
                normal = mean + stdDev * ppnd(uniRan())
            End If
        End Function

        Function ppnd(p As Double) As Double
            '-----
            ' Algorithm AS 111 Applied Statistics (1977) vol 26, no.1
            ' "The percentage points of the normal distribution" by J.D. Beasley and S.G. Springer

            Dim q, r As Double
            Const zero As Double = 0.0, split = 0.42, half = 0.5, one = 1.0,
                  a0 = 2.50662823884, a1 = -18.61500062529,
                  a2 = 41.39119773534, a3 = -25.44106049637,
                  b1 = -8.4735109309, b2 = 23.08336743743,
                  b3 = -21.06224101826, b4 = 3.13082909833,
                  c0 = -2.78718931138, c1 = -2.29796479134,
                  c2 = 4.85014127135, c3 = 2.32121276858,
                  d1 = 3.54388924762, d2 = 1.63706781897
            '-----
            q = p - half
            If Abs(q) <= split Then
                r = q * q
                ppnd = q * (((a3 * r + a2) * r + a1) * r + a0) / ((((b4 * r + b3) * r + b2) * r + b1) * r + one)
            Else
                r = p
                If q > zero Then
                    r = one - p
                End If
                If r <= zero Then
                    Call fatalError[COLOR="Red"]("f-ppnd/probability argument p is not in range 0 to 1")[/COLOR]
                End If
                r = Sqrt(-Log(r))
                ppnd = (((c3 * r + c2) * r + c1) * r + c0) / ((d2 * r + d1) * r + one)
                If q < zero Then
                    ppnd = -ppnd
                End If
            End If
        End Function

        Function uniRan() As Double
            uniRan = Rnd()
        End Function
        Function poidev(xm As Double) As Long
            '------
            ' Poisson random number generator
            '
            ' Based on code in 'Numerical Recipes' by Press et al. (1986)
            '----

            Dim pi As Double = 3.141592654
            Dim em, t, y As Double
            Dim alxm, g, sq As Double
            '------
            If xm < 12.0 Then
                g = Exp(-xm)
                em = -1
                t = 1.0
                Do
                    em = em + 1.0
                    t = t * uniRan()
                    If t <= g Then
                        Exit Do
                    End If
                Loop
            Else
                sq = Sqrt(2.0 * xm)
                alxm = Log(xm)
                g = xm * alxm - gammln(xm + 1.0)
                Do
                    y = Tan(pi * uniRan())
                    em = sq * y + xm
                    If em < 0.0 Then
                        Continue Do
                    End If
                    em = Int(em)
                    t = 0.9 * (1.0 + y ^ 2) * Exp(em * alxm - gammln(em + 1.0) - g)
                    If uniRan() <= t Then
                        Exit Do
                    End If
                Loop
            End If
            poidev = Round(em)
        End Function

        Sub fatalError(text As String)
            Console.WriteLine(text)
            Stop
        End Sub

        Function gammln(xx As Double) As Double
            Dim j As Integer
            Dim ser, tmp, x, y As Double
            Dim cof() As Double = {76.180091729471457, -86.505320329416776,
                                   24.014098240830911, -1.231739572450155,
                                    0.001208650973866179, -0.000005395239384953}
            Const stp As Double = 2.5066282746310007
            '-----
            x = xx
            y = x
            tmp = x + 5.5
            tmp = (x + 0.5) * Log(tmp) - tmp
            ser = 1.0000000001900149
            For j = 0 To 5
                y = y + 1.0
                ser = ser + cof(j) / y
            Next j
            gammln = tmp + Log(stp * ser / x)
        End Function

        ' This subroutine generates a histogram table which may be printed in the calling
        ' code. The sub supports frequency histograms and density histograms. See agrument
        ' list below:
        '  INPUTS:
        '   x = samples to construct histogram from
        '       Double type array ()
        '   nX = number of samples to use in histogram calcs
        '       Long type scalar
        '   nBins = number of histogram bins (columns/bars in plot)
        '       Long type scalar
        '       **** NOTE: histogram bin and bin count lists are of length nBins+1,
        '                   i.e. have max index nBins and min index 0
        '   firstBin = smallest value in historagm bin list
        '       Double type scalar
        '   incrBin = bin interval or width
        '       Double type scalar
        '   densHist = create density historgram if true
        '       Boolean type
        '  OUTPUTS:
        '   histBins = list of histogram bin labels
        '       String type array
        '   binCounts = returns frequency or density of each histogram bin
        '       Double type array
        '   errorFlag = Boolean type, TRUE if error has occured
        '   errorMessage = string giving warning or error message
        Sub CreateHistTable(x() As Double, nX As Long, nBins As Long, firstBin As Double,
                        incrBin As Double, densHist As Boolean, histBins() As String,
                        binCounts() As Double, errorFlag As Boolean, errorMessage As String)
            ' declare local args
            Dim minX As Double
            Dim tempBins() As Double
            Dim i As Long, j As Long
            ' initialise
            errorFlag = False
            ' check nBins > 0
            If (nBins <= 0) Then
                errorFlag = True
                errorMessage = "HISTTABLE: ERROR nBins must be > 0"
                Exit Sub
            End If
            ' check if minX <= firstBin-incrBin and warn if this is the case
            minX = x.Min
            If minX <= (firstBin - incrBin) Then
                errorMessage = "HISTTABLE: WARNING min(x) is outside lowest bin; histogram may be distorted"
            End If
            ' set up bin values and labels
            ReDim tempBins(nBins - 1)
            tempBins(0) = firstBin
            histBins(0) = CStr(firstBin)
            For i = 1 To nBins - 1
                tempBins(i) = tempBins(i - 1) + incrBin
                histBins(i) = CStr(tempBins(i - 1) + incrBin)
            Next i
            'histBins(nBins) = "> " & Format(tempBins(nBins - 1), "######0.0#")
            histBins(nBins) = "> " & CStr(tempBins(nBins - 1))
            For i = 0 To nBins
                binCounts(i) = 0
            Next i
            '--------- histogram calcs
            ' bin counts
            For i = 0 To nX - 1
                For j = 0 To nBins - 1
                    If (x(i) <= tempBins(j)) Then   ' update bin count
                        binCounts(j) = binCounts(j) + 1
                        Exit For
                    ElseIf (j = (nBins - 1)) Then   ' must belong in > bin
                        binCounts(nBins) = binCounts(nBins) + 1
                    End If
                Next j
            Next i
            errorMessage = "HISTTABLE: FrequencyCompletedOK"
            ' determine density if needed
            If densHist Then
                For j = 0 To nBins
                    binCounts(j) = binCounts(j) / (CDbl(nX) * incrBin)
                Next j
                errorMessage = "HISTTABLE: DensityCompletedOK"
            End If
        End Sub

End Module

 

[CJ_London]

Wecome to the forum.

If you want help ou will need to address a couple of things,

Your code is long, difficult to read and no explanation of what it is trying to do. Suggest you repost and use the code tags (the # button) to preserve indenting and also document the code so potential responders understand the flow.

At the moment you are highlighting a bit of text and saying it is an error. But it is an error of your own making in that it gets to that part of the code because r<=zero. This implies you have an issue with the data, not an error in the code. So you will need to identify how r is calculated and work back from that.

Are you aware of debug.print and debug.assert? If not, google them and use them to track variable values

 

[jdraw]

Welcome from NSW.

That's a lot of code with so few comments.
What exactly do you expect/want readers to do?
What is your working environment? ?Console.WriteLine??

This does not do what you think (in vba)
--Dim nrep, nyear, nrepmax As Integer ' loop variables

nrep and nyear will be Variant data type

only nrepmax will be integer

 

[pbaldy]

Suggest you repost and use the code tags (the # button) to preserve indenting

Code tags added.

 

[Cronk]

"nrepmax needs to run for 500 000 replicates"

One problem for a start is that nrepmax is decared an integer which won't handle values above 32767

Declare it as a long variable

 

[gemma-the-husky]

As others noted,

put each variable, and variable type on a separate line (easiest)
use "as long" rather than "as integer". Generally there is nothing to be gained by using a smaller variable type

put "option explicit" at the top of the module - where it says "option compare database", just in case you have any undeclared variables

I can't see the error you are getting, but the problem is most likely just an overflow error caused by the loop counter exceeding the max value for the datatype

 

[Solo712]

Hi All,

Have programmed this code for a water security issue.
The prob.mod was supplied to use along with our own program.

nrepmax needs to run for 500 000 replicates. It works for around 100 000 replicates, but after that there is an error from the prod.mod code. I have highlighted the error message that comes up, not sure if the highlight works.

Thanks in Advance!

        Function uniRan() As Double
             [COLOR="Red"]Randomize[/COLOR]             
             uniRan = Rnd()
        End Function

End Module

The first things that jumps at me is that you are not using the Randomizing function correctly. You must use the Randomize keyword to reseed the random number generator. This may not not explain your problem but you are generating the same random number (p argument in the ppnd function) on each call to the "percentage points of normal distribution" function. So this is one problem that you need to fix regardless.

 

[The_Doc_Man]

Your highlighted error indicates that you are getting results for a probability computation that are not in the double-ended range 0-1. (Single-ended would of course be 0-0.5)

Looking at the code, it is hard to be sure. I think you are building a probability calculator for a "normal" distribution, but you are using what appears to be what is called a "telescoped" function to do this. If you get a large enough sample, you are accumulating a LOT of data and you are inviting all sorts of numeric errors to creep in due to a sad but true fact: Even DOUBLE only gives you a limited number of bits to work with.

Your coefficients have positive and negative signs intermixed, which tells me you are potentially dealing with what is called "catastrophic cancellation" - where you are subtracting two numbers that are so close to each other that most of the bits of precision in either number cancel each other out and you get a VERY rough approximation of the real value that you should get.

The last time I dealt with numbers like that, the problem kicked in when I went a little bit past the six-sigma range in probabilities. I started to get numbers that just barely went beyond a probability of 1, maybe in the 1.0000001 range - even though as sigma approaches infinity the double-sided probability equation as derived symbolically still only approaches one.. And I was using a function that was numerically exact (i.e. not a telescoped function) with all DOUBLE variables as intermediates. When the numbers start to add up, truncation errors creep in and you get these little variations due to various degrees of catastrophic cancellation.

I had the opportunity to try the same function on a machine that could use IEEE "X" format, which was something rude like 128-bit floating and over 100 bits of binary mantissa. Damned probability still started acting up at the 10-sigma level. I didn't actually have anything that precise, I was generating a table for someone. But when you are dealing with that many bits and a polynomial that starts hitting 10, 15, or 20 terms before it loses significance, all sorts of issues can crop up.

My thought is that you need to be more tolerant of these errors and put in a final test so that your function just returns 1 when it gives you that kind of round-off or truncation induced error. OR you could take the alternative approach of running up a bunch of numbers testing the magnitude of your inputs until you start getting the anomaly. At that point, put a test in your function so that if it reaches that number where your answer gets overwhelmed by catastrophic cancellation, just short-circuit the function and return 1.