# 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.

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=","),")"))