MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn moreOpportunities for recent engineering grads.

Apply Today**New to MATLAB?**

Asked by rbx250
on 21 Aug 2013

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.

*No products are associated with this question.*

Answer by Walter Roberson
on 21 Aug 2013

Accepted answer

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.

Show 4 older comments

rbx250
on 21 Aug 2013

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
on 21 Aug 2013

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

## 1 Comment

## Walter Roberson

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/85263#comment_165301

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