sample from a fixed population

[ajmac]

hey all

I need to sample from a fixed population (in this case animals) along a negative exponential distribution of group sizes (on the x axis).

eg. 20,000 animals in total in group sizes that follow a negative exp. distn. (i.e. large number of small groups tailing off exponentially to small number of large groups). The total of all the groups should be the fixed number, 20,000.

Is there any way I can find out, if I sample 1000 times, how many of those samples would produce zeros? i.e. x=0

Please ask if anything is unclear.

Any help would be much appreciated,

ajmac

[triggernum5]

Tell me if I'm not getting your jist.. A group could be considered a family unit in any population of 20000?
Eg:
f(0) = 20000 Loners
f(50)= 10000 Pairs
f(100) = one 20000 member family
???
(Except it runs exponentially not linearly)
OR
Integralx(f(0)->f(100)) = 20000??

Where I'm lost is that traditionally x is your independant variable.. The Only place where x=0 is at f(0).. What does x represent physically if anything at all?
In anycase, I don't see how you intend you derive anykind of probability model.. Rather samples would be 100% predictable and never zero unless there are eigenvalues.. Maybe whip up a rough sketch of your graph to clarify what you mean..

[Rassis]

How many groups do you have? Or, alternatively, what is the group mean size?
What is the rationale behind a sample producing a zero? Aren´t all groups supposed to contain some animals?

Will you please explain a little better what you mean?

[ajmac]

The model is essentially a spatial model. Each sample represents the number of animals in 1km2: therefore, if I sample 200 times, I have defined the area as 200km2. zeros represent sqaure kilometres where no animals are present. This I am doing by sampling from the cumulative probability, which starts at x=0.

However, I need to be able to tell the model how many animals are in the landscape, without changing the size of the area.

My question is: how do I tell the model to splice 20,000 animals into in group sizes that it gets from my cumulative probability, without changing the number of samples. i.e. squeezing 20,000 animals into 200km, in group sizes that it gets from the cumulative distribution (i.e. more small groups than bigger ones).

Hope this makes it a little clearer. Again, give me a shout if not, and I'll think it through again.

cheers

[IMG]modelgraph[/IMG]

[triggernum5]

So in layman terms you're asking:

If I have 20000 animals over 200km^2, what is the probable number of animals in any given 1km^2 area?
This might be the opposite of what you need but I don't think so.. It will give you a typical result, ignoring all physical factors.. For physicality purposes, first apply it to subsets that share particular charicteristics..

'This procedure calculates (in fact estimates) the Mean Value and the Standard Deviation of an
'arbitrary function f(x,y) of two independent normally distributed random variables x and y.
'The program uses the Monte Carlo Simulation Method. The function f(x,y) is specified in
'run-time and its values are calculated using the MS ScriptControl (msscript.ocx).
'Add Component: Microsoft script control 1.0 (msscript.ocx)
 
'Example:
'--------
'Given that random variable x follows the normal distribution (5,1) - with Mean value = 5 and
'Standard Deviation = 1 and that random variable y follows the normal distribution (4,2) - with
'Mean value = 4 and Standard Deviation = 2, What is the Mean value and the Standard Deviation
'of the function f(x,y) = x + x*y - y ?
 
'The Standard Deviation is the most difficult part to calculate here,
'yet note that not always even Mean_f = f(Mean_x, Mean_y) !
 
'Written by Vagelis Plevris, Greece
'mail to: [email][email protected][/email]
 
Private Sub cmdCALCULATE_Click()
 
Dim X_Values(), Y_Values(), F_Values() 'Allocatable arrays that will hold the values of x, y and f(x,y)
Const Pi = 3.14159265358979 'Pi value
 
Mean_X = Val(txtMeanX.Text) 'Mean value of Random Variable x (Given)
StDev_X = Val(txtStDevX.Text) 'Standard deviation of Random Variable x (Given)
 
Mean_Y = Val(txtMeanY.Text) 'Mean value of Random Variable y (Given)
StDev_Y = Val(txtStDevY.Text) 'Standard deviation of Random Variable y (Given)
 
N_MCS = Val(txtN.Text) 'Number of Monte Carlo Simulations to be done
 
'Redim arrays to hold the N_MCS data
ReDim X_Values(1 To N_MCS)
ReDim Y_Values(1 To N_MCS)
ReDim F_Values(1 To N_MCS)
 
'Randomize in order to obtain different results every time the program is run
Randomize
 
For j = 1 To Int(N_MCS / 2)
'Int(N_MCS / 2) because the random numbers are calculated in pairs, so we obtain 2 numbers each time (for every loop)
 
    'Two Random numbers uniformly distributed in the [0,1) space
    '(will be used later to generate the values of x)
    u1x = Rnd
    u2x = Rnd
    
    'Another two Random numbers uniformly distributed in the [0,1) space
    '(will be used later to generate the values of y)
    u1y = Rnd
    u2y = Rnd
    
    'A pair of random numbers that follow the normal distribution (0,1):
    '(0,1) means: Mean value = 0, Standard Deviation = 1
    '(will be used later to generate the values of x)
    Z1x = Sqr(-2 * Log(u1x)) * Sin(2 * Pi * u2x)
    Z2x = Sqr(-2 * Log(u1x)) * Cos(2 * Pi * u2x)
    
    'Another pair of random numbers that follow the normal distribution (0,1):
    '(0,1) means: Mean value = 0, Standard Deviation = 1
    '(will be used later to generate the values of y)
    Z1y = Sqr(-2 * Log(u1y)) * Sin(2 * Pi * u2y)
    Z2y = Sqr(-2 * Log(u1y)) * Cos(2 * Pi * u2y)
 
    'A pair of random numbers that follow the given normal distribution for x (Mean_X, StDev_X):
    '(Mean_X, StDev_X) means: Mean value = Mean_X, Standard Deviation = StDev_X
    '(will be used as values for x)
    X1 = Mean_X + Z1x * StDev_X
    X2 = Mean_X + Z2x * StDev_X
 
    'Assign these values to the corresponding array X_Values
    X_Values(2 * j - 1) = X1
    X_Values(2 * j) = X2
 
    'A pair of random numbers that follow the given normal distribution for y (Mean_Y, StDev_Y):
    '(Mean_Y, StDev_Y) means: Mean value = Mean_Y, Standard Deviation = StDev_Y
    '(will be used as values for y)
    Y1 = Mean_Y + Z1y * StDev_Y
    Y2 = Mean_Y + Z2y * StDev_Y
 
    'Assign these values to the corresponding array Y_Values
    Y_Values(2 * j - 1) = Y1
    Y_Values(2 * j) = Y2
 
    'Now calculate the values of the function f(x,y) and assign these values to the corresponding array F_Values
    F_Values(2 * j - 1) = Func_f_xy(X1, Y1)
    F_Values(2 * j) = Func_f_xy(X2, Y2)
Next j
 
'The arrays are ready, so we can now calculate the Mean value of f(x,y) and its Standard Deviation
 
'For verification purposes, we will also calculate the Mean value and the Standard Deviation of both x and y variables
'The obtained values for the Mean and Standard Deviation of  x and y variables
'should be close to the initially given values Mean_X, StDev_X, Mean_Y, StDev_Y
 
'Calculate the sum of the values of every array
sum_x = 0: sum_y = 0: sum_f = 0 'Set sum to zero
For n = 1 To N_MCS
    sum_x = sum_x + X_Values(n)
    sum_y = sum_y + Y_Values(n)
    sum_f = sum_f + F_Values(n)
Next n
 
'Calculate the mean value (average) of every array
ave_x = sum_x / N_MCS 'Should be close to the given Mean value Mean_X, but not exactly equal to it
ave_y = sum_y / N_MCS 'Should be close to the given Mean value Mean_Y, but not exactly equal to it
ave_f = sum_f / N_MCS
 
'Calculate the Variance and the Standard Deviation of every array
adevx = 0: Varx = 0: epx = 0 'Set values to zero
adevy = 0: Vary = 0: epy = 0 'Set values to zero
adevf = 0: Varf = 0: epf = 0 'Set values to zero
For j = 1 To N_MCS
    Sx = X_Values(j) - ave_x: epx = epx + Sx: adevx = adevx + Abs(Sx): px = Sx * Sx: Varx = Varx + px
    Sy = Y_Values(j) - ave_y: epy = epy + Sy: adevy = adevy + Abs(Sy): py = Sy * Sy: Vary = Vary + py
    Sf = F_Values(j) - ave_f: epf = epf + Sf: adevf = adevf + Abs(Sf): pf = Sf * Sf: Varf = Varf + pf
Next j
Varx = (Varx - epx ^ 2 / N_MCS) / (N_MCS - 1): Sdevx = Sqr(Varx)
'Sdevx should be close to the given Standard Deviation StDev_X, but not exactly equal to it
Vary = (Vary - epy ^ 2 / N_MCS) / (N_MCS - 1): Sdevy = Sqr(Vary)
'Sdevy should be close to the given Standard Deviation StDev_Y, but not exactly equal to it
Varf = (Varf - epf ^ 2 / N_MCS) / (N_MCS - 1): sdevf = Sqr(Varf)
 
'Write the result to the textbox
txtRESULT.Text = ""
txtRESULT.Text = txtRESULT.Text + "Mean of x = " + CStr(ave_x) + vbCrLf + "StDev of x = " + CStr(Sdevx) + vbCrLf + vbCrLf
txtRESULT.Text = txtRESULT.Text + "Mean of y = " + CStr(ave_y) + vbCrLf + "StDev of y = " + CStr(Sdevy) + vbCrLf + vbCrLf
txtRESULT.Text = txtRESULT.Text + "Mean of f(x,y) = " + CStr(ave_f) + vbCrLf + "StDev of f(x,y) = " + CStr(sdevf) + vbCrLf
 
End Sub

[Rassis]

Ajmac,

Your problem must be solved by now…Anyway, here’s what I think if I correctly understood the case.

You have a population of 20,000 animals which are spread over an area of 200 Km2. The area is divided into 200 squares of 1 Km per 1 Km. Every animal laying in any one square at a certain moment form a group. The groups are not randomly spread over the whole area but rather according to some decay function (the larger the group the smaller the number of squares occupied). Empty squares are possible to occur. The function is a discrete function of the Pareto modified type, which you already have in your model and whose parameters are made to vary over time in such a way, that the total number of animals and the total number of squares remain constant. Now you want to sample at random from these 200 squares.

If all that I wrote above is right, then, it is quite obvious how to proceed. Just give numbers to the squares (from 1 to 200) and then sample at random from the equation: 1 + rand()*(200 – 1) and make sure that the same number doesn’t come out more than once, by neglecting the output when it doesn’t comply with this condition.

“rand()” means any number between 0 and 1 sampled from an uniform probability distribution.

If you are still interested and my understanding of the facts is uncorrect, please let me know and I will be glad to try once more.

[ajmac]

Hello all,
the problem is not yet solved, more due to time constraints. Both your interpretations seem right. The parameters of the distribution function (possibly Pareto or negative binomial) I can change depending on whether I want to see greater clumping or more animals spread out, it doesn't necessarily change over time.

Do you know how you would distribute 20,000 animals over 200km2? I didn't quite understand all the code described by Triggernum5, despite his very good efforts at clarity....

[Rassis]

I attach an EXCEL file with a function used to represent learning curves in operations management which, I think, can fairly fit your case.

Because it is a three parameters function, it is a little bit difficult to tackle:

1. The parameter A (cell C3) means that the number of animals diminish by A% whenever the number of squares double (1, 2, 4, 8 and so forth);
2. The parameter Ymin (cell C7) must be such that Yn = 0, that is, we accept to have a few empty squares. You don´t have to change it as I already coded the cell in order it returns a value that satisfies this condition;
3. The parameter Y1 (cell C6) is the only one that you have to adjust by trial and error until you get 20,000 animals in total in cell F7. This you can achieve faster by using the Goal-seek EXCEL capability (Set cell: F7; To value: 20,000; By changing cell: C6).

If you find this function not suitable, visit Wikipedia and try the Pareto distribution modified.

http://en.wikipedia.org/wiki/Pareto_...o_distribution

I hope this helps.