Late, but as promised

[memnus]

I promised a matlab puzzle.

So I'm trying to get the probability of an event from a multinomial distribution. I currently have a function that takes in a vector x and returns the number of permutations of that many elements:

d = factorial(sum(x));
for y = x
    d = d / factorial(y);
end
return d

What I'd like to be able to do is accept a 2-d matrix, and have this function act on each row and return a column vector, without using any for loops. I can write factorial any of the following ways for a single value:

prod(1:x)
prod(cumsum(ones(1, x)))

But again, those don't generalize so well. A factorial that generalizes to act on each element of an array would be spiffy.

Edit: The exact function I'm looking for exists in Matlab 7. In fact, the factorial() function in matlab 7 works exactly how I want it to. I'm using matlab 6.5. If anyone out there actually has matlab 7, and can open their factorial.m and tell me how it works, I'd be much obliged.

click

Comments

[mackenzie]

Can we have something where you need to correct the grammar of a sentence or something?

[squirrelloid.livejournal.com]

i think i need to know what x is before i can answer that for Matlab 6.5... Its obvious x is a vector of integers, so lets see if i cant figure something out...

The solution involves creating y as an nxk matrix where k is size x and n is max(x). The rows of y should be [1:max(x)]. Then, make X a matrix of the same size as y where the columns are all x. This is much easier if k is constant, and certainly if max(x) is bounded (because it doesnt really need to be max(x), but just greater than or equal to max(x)).

Ok, then you perform the operation y = y.*(y<=X), {retain only the relevant numbers for each row}

then perform y = y + y(==0) {which replaces all the 0s with 1s, which are invisible to products}

and then perform prod(y,2). {returns a vector where each element is the product of each row of y}

That will give you all your factorials. Even if you need a loop to make X and y, its probably faster than looping factorials. And if you can bound max(x), and/or k is constant, then you can do things like X = [x,x,x,x...,x] (bound max(x) times) and y = [1:max(x);1:max(x);...] (k times, max(x) meaning the bound of max(x)), which avoids loops entirely. actually, as long as k is bounded, then you can use the bound of k, and refer to y(1:k,:) instead of just y. So, as long as you can bound those quantities, no loops are needed.

(There's got to be an easier way of turning a vector into a matrix where each row/column is that vector, but i'm not thinking of it right now).

[squirrelloid.livejournal.com]

oh, and thats the hard part. obviously the rest of the computation is trivial once you have that.