-1

I've already posted a similar question under this link but it hasn't worked maybe due to the fact that I've expressed myself not clearly enough or I was not able to comprehend the answer. However this question is more precise since if I've realised my project using R (R Studio). Unfortunately I've recognised that R Studio calculates my Poisson-Matrix as slow as Excel that's why I am forced to realise my project using Matlab.

Code in R:

a <- Table$ValueA
b <- Table1$ValueB
lapply(a,function(a) {lapply(b,function(b) {outer(dpois(0:20,a), dpois(0:20,b))})})

Calculation in Excel:

= POISSON(x;a;FALSE) * POISSON(x;b;FALSE)

, where a and b are variables to which all pairs of input values are assigned to, one after another. X is a variable that runs from 0 to 20. So for each record of input-value pairs (a,b) (located in two columns) I want to create a Matrix that is based on the Poisson calculations above, with X running from 0 to 20.

So e.g. the very first result in the matrix, in cell(1,1) is based on the following calculation:

= POISSON(0;a;FALSE) * POISSON(0;b;FALSE)

The following graphic shows the matrix in Excel:

enter image description here

Finally I want to sum up three parts of that matrix and present these three values next to the two columns (containing the input values a and b) for each record of input values:

1) Sum of: The diagonal from the top left corner to the bottom right corner
2) Sum of: The the upper rest of the matrix
3) Sum of: The lower rest of the matrix

It would be so great if there is anybody to help since I really hope to get a performance boost using Matlab!

EDIT: Since there was some confusion I wanted to clarify the following: a and b are column vectors (double 306x1). There are two columns full of different values of a and b. I want to run through all records. In each record there is an "input value pair a,b" (my two lambdas) that should lead to the calculation to the corresponding matrix (21x21) and the presentation of the three results (s1 = sum(diag(M)); s2 = sum(triu(M),'all'); s3 = sum(tril(M),'all');) in three columns next to the columns with the input values (for each specific record).

blowbuh
  • 37
  • 6
  • It would be useful to [edit] the question to clarify what the dimensions of `a` and `b` are and what these values mean conceptually. Look at your [other question](https://stackoverflow.com/q/58932515/8239061), they appear to be the *rate parameter* (lambda) for the Poisson distribution. – SecretAgentMan Nov 22 '19 at 23:01
  • I've edited the question. Your 3rd approach is (almost) working since it delivers a 6426x6426 matrix with all the results. But I need to have a separate matrix for each record of the input values a,b. to sum up the three parts of each matrix and finally have three columns with all the results for each corresponding record of input values. – blowbuh Nov 25 '19 at 20:07

1 Answers1

2

Approach 1: Requires Statistics toolbox
You can use the Probability Distribution Objects if you have the Statistics toolbox.

By lower rest and upper rest, I think you mean the upper and lower triangular matrices. That means the "summing stuff up" part is made easier with diag, triu, and tril. 

% MATLAB R2019a
a = 1;   % Rate for Poisson dist A
b = 2;   % Rate for Poisson dist B

X = 0:20;

pd_a = makedist('Poisson',a);
pd_b = makedist('Poisson',b);

A = pdf(pd_a,X);
B = pdf(pd_b,X);

M = A(:)*B(:).';

% Sum stuff up
s1 = sum(diag(M));
s2 = sum(triu(M),'all');
s3 = sum(tril(M),'all');

Approach 2: Requires Statistics toolbox
This is also achievable with the poispdf function. 

A2 = poisspdf(X,a);
B2 = poisspdf(X,b);
M2 = A2(:)*B2(:).';

Approach 3: No toolbox required
If you don't want to use any toolbox at all, then you can hard code this without too much trouble based on the Poisson distribution and an anonymous function. 

% Poisson pmf
fh =@(k,lam) ((lam.^k).*exp(-lam))./factorial(k);    % k must be integer

A3 = fh(X,a);
B3 = fh(X,b);
M3 = A3(:)*B3(:).';

Note: A low tech way
The above approaches are just using the fact that a column multiplied by a row (of appropriate length) produces a matrix. You can verify this with the code below. 

M4 = zeros(numel(X));
for ii = 1:numel(X)
    for jj = 1:numel(X)
        % use any approach for the calculation itself
        M4(ii,jj) = poisspdf(X(ii),a)*poisspdf(X(jj),b);  % approach chosen arbitrarily

        % M4(ii,jj) = fh(X(ii),a)*fh(X(jj),b);            % alternate
        % M4(ii,jj) = pdf(pd_a,X(ii))*pdf(pd_b,X(jj));    % alternate
    end
end

It is easy to verify that M == M2 == M3 == M4 to an acceptable precision.

SecretAgentMan
  • 2,856
  • 7
  • 21
  • 41
  • This is really amazing, thank you for your detailed answers. Nevertheless I need to say that it is not working, probably by virtue of my non existent Matlab skills. **Approach 2** leads to error "Error using poisspdf (line 29) Requires non-scalar arguments to match in size." **Approach 3** shows a M3 6426x6426 matrix when calculating with 306 records of input values a,b **Approach 4** shows error Error using * Incorrect dimensions for matrix multiplication.. – blowbuh Nov 22 '19 at 22:42
  • **Approach 1** leads to Error using prob.PoissonDistribution>checkargs (line 213) The parameter LAMBDA must be a nonnegative finite numeric scalar. Error in prob.PoissonDistribution (line 91) checkargs(lambda) Error in makedist (line 88) pd = feval(theDefiningClassName, varargin{:}); – blowbuh Nov 22 '19 at 22:43
  • What version of MATLAB are you using? I cleared my workspace and re-verified each code block it all runs. For **#1**, indeed *lambda* must be nonnegative and a scalar value. This imples `a` and `b` are nonnegative scalars (not vectors). – SecretAgentMan Nov 22 '19 at 22:45
  • For **#2**, it looks like you passed `a` and `b` as vectors. Those are the rates. Are they changing? I think this is part of [my confusion from your original question](https://stackoverflow.com/questions/58932515/how-to-calculate-with-the-poisson-distribution-in-matlab/58934162?noredirect=1#comment104166213_58932515) as well. – SecretAgentMan Nov 22 '19 at 22:50
  • Sorry for my late answer, they are changing. So a and b are column vectors (double 306x1). There are two columns full of different values. I want to run through all records. In each record there is an "input value pair a,b" that leads to a matrix (21x21) and the presentation of the three results (s1 = sum(diag(M)); s2 = sum(triu(M),'all'); s3 = sum(tril(M),'all');) in three columns next to the columns with input values (for each specific record). So right now Matlab expects not column vectors but scalars? So I would need to change the values of all records to scalars? – blowbuh Nov 24 '19 at 17:57