Probability estimation in a finite set of probabilities I have the following problem. I have I finite family $\mathcal{F}=\{ P_1,...,P_2 \} $ of discrete probability distributions on the same set $S$. Now I observe experiments from a random variable X which is distributed as one of the probabilities in $\mathcal{F}$. I want to compute the probability that $X \sim P_i$ for each $i$, or an estimate of these probabilities with some thoeretical guarantee.
The real problem I have is that I have to keep a belief on these probabilities and take decisions based on this belief, update it as I observe new samples, and change my decisions accordingly.
The two approaches I thought are with a Dirichlet-categorical Bayesian model, but it is continuous so the probabilities of single probabilities are 0, or with concentration inequalities, but cannot manage to find a solution.
 A: This seems to be a straightforward application of Bayes' Theorem. Your initial belief on the distributions' probabilities yields your priors $P(P_i)$. Given a distribution $P_i$, you can calculate the conditional probability 
$$P(X|P_i) = \prod_{j=1}^n P_i(x_j),\text{ where } X=\{x_1, \dots, x_n\},$$
and then you simply update
$$ P(P_i|X) \propto P(X|P_i)P(P_i). $$

As an illustration, assume your discrete distributions are binary, and that we have three candidates
$$ \mathcal{F} = \{(0.2,0.8), (0.4,0.6), (0.9,0.1)\} $$
with equal prior probabilities. Assume further that the true distribution is $(0.5,0.5)$, so not one of the candidates, although the second one comes closest (but this doesn't matter). Here is how the posterior probabilities develop as we draw 100 samples $x$ and update the probabilities after each observation:

distributions <- rbind(c(0.2,0.8),c(0.4,0.6),c(0.9,0.1))
nn <- 100
probabilities <- matrix(NA,nrow=nn,ncol=nrow(distributions))
probabilities[1,] <- 1/3

for ( ii in 2:nn ) {
    set.seed(ii)    # for reproducibility
    xx <- sample(c(1,2),size=1,prob=c(0.5,0.5))
    probabilities[ii,] <- probabilities[ii-1,]*distributions[,xx]       # update
    probabilities[ii,] <- probabilities[ii,]/sum(probabilities[ii,])    # normalize
}

plot(c(1,nn),c(0,1),type="n",xlab="Step",ylab="Probabilities")
for ( ii in 1:ncol(probabilities) ) lines(1:nn,probabilities[,ii],col=ii)
legend("right",lwd=1,col=1:ncol(probabilities),
    legend=paste0("(",apply(distributions,1,paste,collapse=","),")"))

