From: Roger Stafford on
Shonkhor <will.allen(a)hotmail.com> wrote in message
<27978954.1209894665341.JavaMail.jakarta(a)nitrogen.mathforum.org>...
> Thanks very much to all who responded.

You are quite welcome.

> I have used Roger Stafford's - can I credit your effort in a report I have for
submission?

Yes, you may do that.

Roger Stafford

From: alessandro mura on
Dear Roger,
thank you. Yes, both methods require massive
optimizations. For example, I usually need to
generate N (millions or more) of test particles in my
montecarlo models, and I prefer the first method,
picking N/w particles per iteraction, having more than w iteractions
until I get n> N particles. This leads to a number of "good" values
slighlty higher than N,
but allows some vectorialization in the code.
The second method: I would go for an analytical (or polynomial best fit)
expression for ipdf to avoid "find". I was posting a code like that, then I
preferred the "find" option
for better explain the phylosophy of the method.

Regards

Ale
--
Alessandro Mura
Istituto Nazionale di Astrofisica - IFSI
http://pptt4.ifsi-roma.inaf.it/~mura/index.html
http://www.alessandromura.it


From: Roger Stafford on
"alessandro mura" <fake(a)address.com> wrote in message <fvl8um$1np
$1(a)aioe.org>...
> Dear Roger,
> thank you. Yes, both methods require massive
> optimizations. For example, I usually need to
> generate N (millions or more) of test particles in my
> montecarlo models, and I prefer the first method,
> picking N/w particles per iteraction, having more than w iteractions
> until I get n> N particles. This leads to a number of "good" values
> slighlty higher than N,
> but allows some vectorialization in the code.
> The second method: I would go for an analytical (or polynomial best fit)
> expression for ipdf to avoid "find". I was posting a code like that, then I
> preferred the "find" option
> for better explain the phylosophy of the method.
>
> Regards
>
> Ale
> --
> Alessandro Mura
---------
Yes, some kind of reasonably easy to calculate, analytical fit to the inverse
cumulative probability distribution function would be fine if it can be found.
Even if it is somewhat inaccurate, it can be used in conjunction with a
subsequent search, hopefully requiring only a very few comparisons, to find
the exact value, starting from the estimate and referring to the known cdf.
However, I think you would have trouble doing polynomial fits to inverse cdf's
for distributions with unbounded ranges.

Roger Stafford


From: Peter Perkins on
Shonkhor wrote:
> Can anyone help with this one?
>
> I want to use a probability matching rule in an ideal observer model. From a non-uniform discrete probability distrubution (in this case p of angles between 1 and 360 degreed in 1 degree integers) I want to select a response angle with a probability equal to its probability in this distribution.
>
> I hope this makes sense but any questions just ask. Does anyone have any idea how to implement this in Matlab? It has me stumped.

Will, unless I've misunderstood your description, you have

1) the values 1:360, and
2) a probability for each of those angles.

If you have access to the Statistics Toolbox, all you need is

randsample(1:360,p,n,true)

to generate a sample of size n, drawing values from 1:360 independently with
replacement.
From: Roger Stafford on
Peter Perkins <Peter.PerkinsRemoveThis(a)mathworks.com> wrote in message
<fvnakj$17t$1(a)fred.mathworks.com>...
> Will, unless I've misunderstood your description, you have
>
> 1) the values 1:360, and
> 2) a probability for each of those angles.
>
> If you have access to the Statistics Toolbox, all you need is
>
> randsample(1:360,p,n,true)
>
> to generate a sample of size n, drawing values from 1:360 independently
> with replacement.
-----------
That will teach me not to go diving into a problem without checking to see if
someone else has already done it. I have a weakness in that direction!
Anyway, Shonkhor got three versions that don't require the Statistics Toolbox.

According to the documentation for 'randsample', the call should be made
as:

randsample(1:360,n,true,p);

Does your argument sequence really work also?

I do have a question about 'randsample' if it doesn't violate Mathworks'
proprietary interests. The four-argument case is markedly more difficult to
implement than the others. Does it use something like the binary search I
used to effect the inverse cumulative distribution of p, Simon's 'histc' method,
or perhaps something akin to Alessandro's acceptance-rejection method
directly on the p values? In other words, if m = 360, is it an order(n*log(m))
algorithm, order ((n+m)*log(n+m)) algorithm, or something altogether
different? Users would probably like to know this, but the documentation
doesn't say.

Roger Stafford