Path: news.mathworks.com!not-for-mail
From: "Richard " <rfrank@dominionsw.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Random sequence with descending probabilities
Date: Fri, 1 Oct 2010 13:38:08 +0000 (UTC)
Organization: Dominion Software, Inc.
Lines: 110
Message-ID: <i84o80$4cg$1@fred.mathworks.com>
References: <i83oc0$ann$1@fred.mathworks.com> <i83sok$l9n$1@fred.mathworks.com>
Reply-To: "Richard " <rfrank@dominionsw.com>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1285940288 4496 172.30.248.35 (1 Oct 2010 13:38:08 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Fri, 1 Oct 2010 13:38:08 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1953011
Xref: news.mathworks.com comp.soft-sys.matlab:675010

Hi Roger,

Thanks for the help.

I'm not sure that round() isn't ok here...? It seems to be working as I would expect....

There's a couple of things that I think need some tweaks

y = floor(L*((U+1)/L).^rand(1,n));

will always return numbers between 1 and L (because n^0 is 1) - rand here is returning a 1 by n matrix of numbers 0....1

and also L cannot be 0 - I failed to mention that 0 - n is a valid range....
so........ I tweaked things and came up with the code below - one of my concerns
was that there be no numbers in the range that end up (because of rounding or something) with a 0 probability. All integers in the range must have > 0 probability.....

I have to port this code to another language so I can rely on any matlab particular functions (except for just error checking)...

This seems to work in my test cases so far.....let me know what you think...

Thanks!

Rick
clear;

% change for test cases....
lower = 5;
upper = 10;

power = 1;
exponential = 2;

which = exponential;

n = 10000;
x = rand(1,n);
 
if which == exponential
   y = x .^ (1 + rand(1,n));
   y = lower + (y * ((upper - lower))); 
   y = round(y);
   
else
   y = x .^2;
   y = lower + (y * ((upper - lower))); 
   y = round(y);

end;
figure(1),plot(sort(y));
mn = min(y);
mx = max(y);

figure(2),hist(y,(upper - lower) + 1);

x = lower:1:upper
counts = histc(y,lower:1:upper);
if any(counts <= 0.5)
    msg = 'Error = histogram gap'
else
    msg = 'OK  no histogram gap'
end

 csum = cumsum(counts);
 figure(3);
bar(x,csum,'BarWidth',1)



"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i83sok$l9n$1@fred.mathworks.com>...
> "Richard " <rfrank@dominionsw.com> wrote in message <i83oc0$ann$1@fred.mathworks.com>...
> > Hi
> > 
> > I want to create small sets of random integers - say 2- 23, and or 20 - 39.
> > 
> > I would like to have the lower number have a higher probability than the higher number,
> > with some sort of exponentially decreasing curve.
> > 
> > Here's one idea that I tried:
> > 
> > x = rand(1,10000);
> > 
> > y = x .^2;
> > 
> > lower = 5;
> > upper = 20;
> > 
> > y = lower + (y * (upper - lower));
> > 
> > y = round(y);
> > 
> > mn = min(y);
> > mx = max(y);
> > 
> > figure,hist(y,(upper - lower) + 1);
> > 
> > I wonder if this approach is robust? Will there be gaps in this histogram for some ranges?
> > 
> > I'm not sure why it seems to work, which is why I guess I'm unsure of it.
> > 
> > Thanks,
> > 
> > Rick
> - - - - - - - -
>   I would not use the 'round' function here since that would give an anomalous distribution at the two endpoints.  Also what you used was not an exponential curve but a power curve - in this case, the power of 2.  For an exponential distribution curve I would do the following.  Let U and L be your two integer-valued range extremes.
> 
>  y = floor(L*((U+1)/L).^rand(1,n));
> 
> The values of y should all be integers and range from L to U.
> 
> Roger Stafford