Error using ==> histc; Edge vector must be monotonically non-decreasing.

So I will admit up front that I am terrible at coding, but I have run into an issue that I haven't been able to solve. The error referenced in the title seems to occur about 1/10000 runs of the simulation (which in turn means 1 in a million uses of the algorithm). I believe I know where in the code the error arises (if I comment out the block of code I no longer receive the error), but not what the error means or how to remedy it.

The error occurs in the following block of code:

[Xn,Yn] = selection(X,Y,sI,sB,r,u,N);

%%Expected Genotypes
XAXA = Xn(1,1)*Xn(2,1);
XAXa = Xn(1,2)*Xn(2,1) + Xn(1,1)*Xn(2,2);
XaXa = Xn(1,2)*Xn(2,2);

XAYA = Xn(1,1)*Yn(1);
XAYI = Xn(1,1)*Yn(2);
XAYa = Xn(1,1)*Yn(3);
XaYA = Xn(1,2)*Yn(1);
XaYI = Xn(1,2)*Yn(2);
XaYa = Xn(1,2)*Yn(3);

H1 = [XAXA,XAXa,XaXa];
H2 = [XAYA,XAYI,XAYa,XaYA,XaYI,XaYa];

%%Actual Genotypes

G1 = mnrnd(N/2,H1)/(N/2);
G2 = mnrnd(N/2,H2)/(N/2);
G = [G1,G2];

But, I believe it is a result of something happening in the function "selection". Particularly, in this segment of the code:

%Allele frequencies of gametes produced in the current generation before mutations

Xt(1,1) = (G_now(1)+0.5*G_now(2))/sum(G_now(1:3)); %XA in eggs
Xt(1,2) = (0.5*G_now(2)+G_now(3))/sum(G_now(1:3)); %Xa in eggs
if(Xt(1,2) < 1e-8)
    Xt(1,2) = 0;
end

Xt(2,1) = (G_now(4) + G_now(5) + (1-r)*G_now(6) + r*G_now(7))/sum(G_now(4:9)); %XA in sperm
Xt(2,2) = (r*G_now(6) + (1-r)*G_now(7) + G_now(8) + G_now(9))/sum(G_now(4:9)); %Xa in sperm
if(Xt(2,2) < 1e-8)
    Xt(2,2) = 0;
end

Yt(1) = (G_now(4) + r*G_now(6)+(1-r)*G_now(7))/sum(G_now(4:9)); %YA in sperm
Yt(2) = (G_now(5) + G_now(8))/sum(G_now(4:9)); %YI in sperm
Yt(3) = (G_now(9) + (1-r)*G_now(6) + r*G_now(7))/sum(G_now(4:9)); %Ya in sperm
if(Yt(3) < 1e-8)
    Xt(3) = 0;
end

%Allele frequency of gametes produced in the current generation after mutations

Xn(1,1) = Xt(1,1) + N*u*Xt(1,2);
Xn(1,2) = Xt(1,2) - N*u*Xt(1,2);
Xn(2,1) = Xt(2,1) + N*u*Xt(2,2);
Xn(2,2) = Xt(2,2) - N*u*Xt(2,2);
Yn(1) = Yt(1) + N*u*Yt(3);
Yn(2) = Yt(2);
Yn(3) = Yt(3) - N*u*Yt(3);

Previously, I had frequencies for Xt(1,2), Xt(2,2), and Yt(3) defined as 1 minus their respective complements. For whatever reason, that produced more infrequent errors than the code as it is currently shown. Basically, I have no clue what is wrong.

[Walter Roberson]

http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

Answer 1

[Walter Roberson]

You do not have any obvious use of histc() in the code you show. Is the error being reported as part of one of the calls to mnrnd() ?

The error means that when histc() is called, the bin list given must not have any element that is less than a previous element, but elements equal to the immediately previous element of the bin list are accepted.

[rbx250]

Yes, histc is being called during mnrnd. I suppose I understood what a non-monotonically decreasing edge vector means, what I am unclear on is how the edge vector that defines those bins is constructed and by extension, what in my code is causing the edge vector to not be monotonically decreasing.

I am also unclear why different formulations for the same value result in different frequencies of occurrence of the error.

[Walter Roberson]

There can be round-off error that differs depending on exactly how you calculate the values.

Unfortunately I do not have the stats toolbox that has mnrnd so I cannot test your code. My guess is that you have some fairly small probabilities and that the way the histogram bins are being calculated is slightly wrong.

Could you use "dbstop if error" and when the program stops, "dbup" to find the list of bins being passed to histc() and show those to us (in high precision, at least "format long g"); and then dbup some more to get to the call to mnrnd() and show us the parameters you passed in, again in high precision ?

[rbx250]

Sure. The call to mnrnd is as follows:

G2 = mnrnd(N/2,H2)/(N/2);

N is a scalar set equal to 10000. H2 is as follows (in long g):

H2 =

Columns 1 through 4

         0.946692712763803        0.0531858873091101      0.000121399927087186                         0

Columns 5 through 6

                         0                         0

The edge vector being passed to histc is as follows:

edges =

Columns 1 through 4

                         0         0.946692712763803         0.999878600072913                         1

Columns 5 through 7

                         1                         1                         1

And the the "p" vector that histc is using is:

p =

Columns 1 through 4

         0.946692712763803        0.0531858873091101      0.000121399927087186                         0

Columns 5 through 6

                         0                         0

I don't see how "edge" isn't monotonically non-decreasing as the error implies.

[Walter Roberson]

Hmmmm... Try diff(edges) to see whether the "1" is round-off even to the "format long g" level. And try

fprintf('%.99g\n', edges)

to see the values to more precision

[rbx250]

I think I have a completely crappy work-around to the problem. I played with the value from dbstop (thank you for telling me about it) and mnrnd and what it did not like was the trailing zeros in the probability vector H2. If I only fed mnrnd H2(1:3) it went through fine. So, I set up a switch that truncated the vector if the last elements were zeros and so far it seems to not be getting the same error. I have my fingers crossed that I haven't just made a 1/10000 problem a 1 in a billion problem.

[Walter Roberson]

histc itself is happy with a bin vector of [0 0.946692712763803 1 1 1], so if histc complained then at least one of the "1" must have been slightly different from true 1.0

[rbx250]

Yeah, it is entirely possible that that is the case. I have run close to a million sims with the patch in place and so far it hasn't broken.