I am not sure if I understand your question correctly. I guess once you get P, then you would like to choose its elements randomly, right? If so, you can use the following code:
alpha = 1; beta = 1; % replace by your desired values
for k = 0:7
fun = @(lambda) (lambda.^k).*(exp(-1.*lambda)).* gampdf(lambda,alpha,beta)./factorial(k);
P(k+1) = integral(fun,0,Inf);
end
out = datasample(P,1) % 1 observation, replace 1 by higher value for more observations
This will be uniformly random. If you are interested then you may use:
out = P(randperm(numel(P),1))
