function [db0,da0,da1,da2,da3,pb0,pa0,pa1,pa2,pa3]=shuffER(x1,x2,y1,y2,g01,g02,N)
inob0=0;
inoa0=0;
inoa1=0;
inoa2=0;
inoa3=0;


modelfunER1= @(a,x1)( exp(x1(:,1).*a(1)).* (a(2) + a(3).*x1(:,2) + a(4).*x1(:,3) + a(5).*x1(:,4)) );
modelfunER2= @(b,x2)( exp(x2(:,1).*b(1)).* (b(2) + b(3).*x2(:,2) + b(4).*x2(:,3) + b(5).*x2(:,4)) );
modelfunER1r= @(ar,arn1)( exp(arn1(:,1).*ar(1)).* (ar(2) + ar(3).*arn1(:,2) + ar(4).*arn1(:,3) + ar(5).*arn1(:,4)) );
modelfunER2r= @(br,arn2)( exp(arn2(:,1).*br(1)).* (br(2) + br(3).*arn2(:,2) + br(4).*arn2(:,3) + br(5).*arn2(:,4)) );

fER1=nlinfit(x1,y1,modelfunER1,g01);
b01=fER1(1,1);
a01=fER1(1,2);
a11=fER1(1,3);
a21=fER1(1,4);
a31=fER1(1,5);
fER2=nlinfit(x2,y2,modelfunER2,g02);
b02=fER2(1,1);
a02=fER2(1,2);
a12=fER2(1,3);
a22=fER2(1,4);
a32=fER2(1,5);
db0=b01-b02;
da0=a01-a02;
da1=a11-a12;
da2=a21-a22;
da3=a31-a32;

Nd1=size(x1,1);
Nd2=size(x2,1);
count=1;

 while count<=N
  arn1=[];arn2=[];y1r=[];y2r=[];
  
  ip=randperm(Nd1+Nd2);
  xx=[x1;x2];
  yy=[y1;y2];
  arn1=xx(ip(1:Nd1),:);
  arn2=xx(ip(Nd1+1:Nd1+Nd2),:);
  y1r=yy(ip(1:Nd1));
  y2r=yy(ip(Nd1+1:Nd1+Nd2));

  fER1r=nlinfit(arn1,y1r,modelfunER1r,g01);
  b01r=fER1r(1,1);
  a01r=fER1r(1,2);
  a11r=fER1r(1,3);
  a21r=fER1r(1,4);
  a31r=fER1r(1,5);
  fER2r=nlinfit(arn2,y2r,modelfunER2r,g02);
  b02r=fER2r(1,1);
  a02r=fER2r(1,2);
  a12r=fER2r(1,3);
  a22r=fER2r(1,4);
  a32r=fER2r(1,5);
  
  if(max(abs(fER1r))>100 || max(abs(fER2r))>100 || min(abs(fER1r))<1.e-5 || min(abs(fER2r))<1.e-5)
      continue
  end
  
  db0rn=b01r-b02r;
  da0rn=a01r-a02r;
  da1rn=a11r-a12r;
  da2rn=a21r-a22r;
  da3rn=a31r-a32r;
  
  if(abs(db0rn)>=abs(db0))
      inob0=inob0+1;
  end
  if(abs(da0rn)>=abs(da0))
      inoa0=inoa0+1;
  end
  if(abs(da1rn)>=abs(da1))
      inoa1=inoa1+1;
  end
  if(abs(da2rn)>=abs(da2))
      inoa2=inoa2+1;
  end
  if(abs(da3rn)>=abs(da3))
      inoa3=inoa3+1;
  end
  
  count=count+1;
  
 end

pb0=inob0/N;
pa0=inoa0/N;
pa1=inoa1/N;
pa2=inoa2/N;
pa3=inoa3/N;