Matlab version for the Poisson-Binomial distribution?

jscargle commented 6 years ago

Hi ... I want to use the Poisson-Binomial distribution for statistical assessment of some results I am getting from LIGO data. It would be easiest for me if there were a MatLab version, although I may know enough rudimentary Python to translate or even run it.

Has anyone done comparisons with randomized simulations to compare with the theoretical formula?

Any advice, comments, volunteering to translate to MatLab, etc. are welcome!

tsakim commented 6 years ago

Hello @jscargle,

unfortunately I am not an expert in MatLab and I don't know of any implementation. However, there seems to be a workaround to execute Python code in MatLab:

There is also an implementation in R from the author of the original article himself which may be useful for you:

https://cran.r-project.org/web/packages/poibin/index.html

Regarding the randomized simulations, the code was tested using an accuracy check proposed in the original article by Hong and some probability results were compared with the R package. Randomized tests would be a good idea though.

I am happy to hear that you are working on the LIGO project, hopefully this helps!

jscargle commented 6 years ago

Hi ... thanks very much for these comments and suggestions. It is intriguing that one can execute Python from MatLab.

Cheers, Jeff

Jeffrey D. Scargle Space Science Division Mail Stop 245-3 NASA Ames Research Center Moffett Field, CA 94035-0001

Phone: 650 604 6330 Cell: 415 385 4297 Fax: 650 604 6779

From: tsakim [notifications@github.com] Sent: Monday, August 06, 2018 1:54 PM To: tsakim/poibin Cc: Scargle, Jeffrey D. (ARC-SST); Mention Subject: Re: [tsakim/poibin] Matlab version for the Poisson-Binomial distribution? (#5)

Hello @jscarglehttps://github.com/jscargle,

unfortunately I am not an expert in MatLab and I don't know of any implementation. However, there seems to be a workaround to execute Python code in MatLab:

There is also an implementation in R from the author of the original article himself which may be useful for you:

https://cran.r-project.org/web/packages/poibin/index.html

Regarding the randomized simulations, the code was tested using an accuracy check proposed in the original article by Hong and some probability results were compared with the R package. Randomized tests would be a good idea though.

I am happy to hear that you are working on the LIGO project, hopefully this helps!

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/tsakim/poibin/issues/5#issuecomment-410849872, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ABo4AmsJV3DqklWxq4H4bAwR5tPVPPigks5uOK0YgaJpZM4VotTl.

aravinat21 commented 1 year ago

In case someone finds it useful:

function [pmf, cdf] = poibin(p)
% Created on Tue Mar 29, 2016
% Author:
%     Mika Straka
% 
% Description:
%     Implementation of the Poisson Binomial distribution for the sum of
%     independent and not identically distributed random variables as described
%     in the reference [Hong2013]_.
% 
%     Implemented method:
% 
%         * ``pmf``: probability mass function
%         * ``cdf``: cumulative distribution function
% References:
% .. [Hong2013] Yili Hong, On computing the distribution function for the Poisson
%     binomial distribution,
%     Computational Statistics & Data Analysis, Volume 59, March 2013,
%     Pages 41-51, ISSN 0167-9473,
%     http://dx.doi.org/10.1016/j.csda.2012.10.006.
sz = size(p);n = sz(2);
omega = 2 * pi / (n + 1);
chi = ones([n + 1, 1,  sz(3:end)]);
nh = floor(n / 2 + mod(n, 2));
exp_value = exp(omega * repmat(1: nh, [1, 1, sz(3:end)]) * 1j);
xy = 1 - p + p .* pagetranspose(exp_value);
argz_sum = sum(atan2(imag(xy), real(xy)), 2);
exparg = sum(log(abs(xy)), 2);
d_value = exp(exparg);
chi(2:nh + 1, :) = reshape(pagetranspose(d_value) .* pagetranspose(exp(argz_sum * 1j)), nh, []);
chi(nh + 2:n + 1, :) = flip(conj(chi(2:n - nh + 1, :)), 1);
chi = chi / (n + 1);

pmf = fft(chi, [], 1);
if any(imag(pmf) > eps)
    error 'pmf values have to be real.'
end
pmf = real(pmf) + eps;
cdf = cumsum(pmf);
end

Or the following:

function [P, Q]=poibin2(p)
p(p==0) = [];
N = length(p);
alpha = prod(p);
s = -(1-p) ./ p;
S = poly(s);S=S(N+1:-1:1);
P = alpha * S;
Q = cumsum(P);
end

as mentioned here.

tsakim / poibin

Matlab version for the Poisson-Binomial distribution? #5