How to infer the parameter $p = f(n)$ of different Bernoulli distributions $X_{n}$? I have a dataset corresponding to the results of independant Bernoulli trials. Each trial is associated to a number $n \in ]1;+\infty[$ and follows a Bernoulli distribution with parameter $p=f(n)=b + a/n$. The problem is almost each trial has a different $n$ so follows a different Bernoulli distribution.
Here is a part of my dataset:
 n      result of the trial
 1.47   1
 3.21   0
10.22   0
 1.17   1
 1.17   0
 2.12   1
 1.82   0
...

How can I find the most likely $a$ and $b$?
 A: This seems a problem easily solved by additive risk models. Create the regressor 1/n and estimate the additive risk model relating 1/n to the outcome $y$, you will obtain maximum likelihood estimates of $a$ and $b$. This is one of a general class of "binomial models" which are among a general class of "generalized linear models". 
The most famous example in this class is logistic regression which attributes to the outcome outcome a binomial variance, i.e. that the variance is equal to the mean times one minus the mean, and that the mean response is related to the predictors through a "logit link". $$\log \left( \text{Pr}(Y|X) / (1-\text{Pr}(Y|X)) \right) = X^T\beta$$.
However one can maintain the binomial variance attribute and vary the link. In particular, additive risk uses the identity function as the "link":
$$\text{Pr}(Y|X) = X^T\beta$$
which in 1-dimension is exactly the linear model you express above.
Estimation is done through maximum likelihood and Fisher Scoring in most default software. In R, for instance, the "call" is glm(y ~ x, family=binomial(link="identity"). In fact, there are some antiquated warnings in the ?glm help file about not fitting such a model. But these dismiss the importance of additive risk models and of using alternative fitting strategies (there is ample discussion on that topic through this SE).
A: Question: Do you have any bounds or relation between $a$ and $b$? 
If you don't, suppose $n = 10$, you have many answers for $a$ and $b$ that satisfies $p = 0.5$:


*

*$b = 0$ and $a = 5$

*$b = 1$ and $a = -5$

*$b = 10$ and $a = -95$
If you do have bounds or relationship between your parameters, since you have clear Bernoulli experiment with a custom parameter, I would recommend you to write your likelihood as a function of $a$ and $b$ and use optim function in R.
There is an incomplete sample code that should help you to implement in your experiment, but please notice the parts where you should insert your values based on the bounds and relationship I mentioned before.
loglik = function(pars, n, result){
a = pars[1]
b = pars[2]
my_p = b + a/n

## Your log likelihood
lk = sum(result*log(my_p) + (1-result)*log(1 - my_p))

## Return negative to maximize it
return(-lk)
}

## init_par = <insert your initial value>
my_n = c(1.47, 3.21, 10.22, 1.17, 1.17, 2.12, 1.82)
my_result = c(1,0,0,1,0,1,0)

optim(par = init_par, 
      fun = loglik,  
      method = "L-BFGS-B",
      n = my_n, 
      result = my_result,
##    lower = c(<low_a>, <low_b>),
##    upper = c(<up_a>, <up_b>)
     )

```

A: I did a very rough approximation to $a$ and $b$ via a system of equations. In your data, you sort all of the $n$ values and bin them (in the code below I made bins of same length, it would be probably better to make bins with roughly equal number of $n$'s inside). You count how many $1$'s and $0$'s are in your bins ($\#_1$ and $\#_{1\&0}$), and how many $n$'s fall into this range ($\#_n$). Then you have: $\#_1/\#_{1\&0} = \sum_{n's}(b + a\cdot\frac{1}{n})/\#_n$ for every bin, where only $a$ and $b$ are unknown.
This gives a system of equations. See my example solution in R below:
a <- 0.5
b <- 0.2

B <- 1000000

set.seed(2407)
n <- runif(B, 1, 1000)
bn <- rbinom(B, 1, b+a/n) # result of the trials

bins <- 10
grid <- seq(min(n), max(n), length.out = bins+1)
lhs <- rep(0, bins)
for_a <- rep(0, bins)
for (i in 1:bins) {
  lhs[i] <- sum(bn[n>=grid[i] & n<grid[i+1]])/length(bn[n>=grid[i] & n<grid[i+1]])
  for_a[i] <- sum(1/n[n>=grid[i] & n<grid[i+1]])/length(n[n>=grid[i] & n<grid[i+1]])
}

indices <- !is.na(for_a) & !is.na(lhs)
for_a <- for_a[indices]
lhs <- lhs[indices]
for_b <- rep(1, length(for_a))

lm(lhs ~ -1+for_a+for_b)

The results look okay, but some trials showed that especially in $a$ is a lot of variation:
Call:
lm(formula = lhs ~ -1 + for_a + for_b)

Coefficients:
 for_a   for_b  
0.5227  0.1994

