A direct approach is to use a wheel of fortune with $k$ equidistant pointers at angular positions
$$\zeta_j~=~\alpha+ 2\,\pi\,j/k,~~~~j=0,\ldots,k-1,$$ with equi-distributed $\alpha \in (0,2\pi/k]$ and randomly permuted sectors of angular size $2\,\pi\, p_i/k$. Choose the indices whose sectors are indicated by the pointers.
To show that for each $i=1,\ldots,d$ one has $P(i\in I) ~=~p_i$ consider the angular sector belonging to the index $i$ can at most be reached by two pointers $\zeta_j, \,\zeta_{j+1 \mod k}$, and that the sector can be divided in two sub-sectors each of which can be reached by only one of these pointers, but not at the same time. From that it is easily calculated that the combined probability to reach the angular sector belonging to the index $i$ by choosing an equi-distributed $\alpha \in (0,2\pi/k]$ indeed is $p_i$.
A slightly simplified version (using intervals instead of sectors) is implemented in the following Mathematica function that expects a list indprob
of $d$ pairs $(i, p_i)$ and the value of $k$ fulfilling the condition $\sum p_i =k$. It then returns a list of $k$ indices having the corresponding probabilities.
RandomIndices[indprob_, k_] := Module[{ransam = RandomSample[indprob], acc, pointer, ind,
alpha = RandomReal[{0, 1}]},
acc = Accumulate[ransam[[;; , 2]]];
pointer = Range[k] + alpha - 1;
ind = With[{p = #}, FirstPosition[acc, _?((# > p) &)]] & /@ pointer;
ransam[[Flatten[ind]]][[;; , 1]] // Sort]
To destroy potentially remaining correlations between the indices one could further subdivide the sectors of each index into several pieces and permute the then much larger number of sectors.