Thanks for the detailed answer. It kind of surprises me, that the Poisson distr. has the same accuracy to the third decimal place as the binomial, but that shows my naivety with statistics. Thanks for this insight. In my opinion then both approaches should be a reasonable answer.

Right, one should use $N=213$ since it's the possible positive solutions. Also, I could simply use $1-binom(k=0, N=213, p=1/300) since then I don't need to use a sum over a range, couldn't I?

Okay, thank you for your help and advice, that is also exactly what I did with the extra columns.

I finally made a working MICE code that works quite good on my dataset, but I didn't used the indicator variable for this step. However, should I (or can I even) give the model for my ML code, BDT or NN, the indicator variabel or should I neglect it for the fitting of the model (or just let it in as a column of 0 and 1)?

The indicator variable can be the nan in for example a panda dataframe or even the mentioned -999 if I understand that part of your answer correctly? Or can one even give the model the indicator manually? And have you maybe used a method like MICE or KNNImputation yourself or can comment how this could maybe even manipulate the data in a negative way?

Yes, I tried different orders and a few other typical models. The thing is, there ain't a formula which describes this. The purpose is only to find the local minimum. If you have an Idea of a fitting method, please let me know.

Yes ok, so I thought of scipy.odr, Orthogonal distance regression . But nethertheless the minimum stays more or less the same (ca. 154.6) and the errors are stil quite big (speaking of up to 12%). I could also give a minimalistic code example of it, if necessary.