業務でやりたいなーと思ったことをできる方法が見つけられなかったので自分で考えました。
そんなに珍しい対象でも、方法論として最新でもなんでもないので、「それこんな名前が付いてるよ」とか「こんな解き方がメジャーだよ」いった情報があれば教えて頂けると助かります。
問題設定
簡単に言うと、目的変数は、〇〇ratio系です。
説明変数xから目的変数yを求めるような問題を回帰というわけですが、今回は目的変数が、「$N$回コインを投げた際、$k$回表が出た」といった結果が渡されます。
コインは複数種類あり、どのコインを渡されたかは、説明変数$x$に依存して決まることがわかっています。
また、各コインを投げる回数は$x$に依存しないものとします。
試行回数が多いコインは信頼できるけど、試行回数が少ないコインは信頼できません。
さて、いずれのコインも無限回振れば、真の確率pは得られるわけですが、今回のケースではそうはいきません。
対処法の一つは、十分大きい回数のサンプルだけを残して回帰することです。
では、十分大きい回数のサンプルが少ないといった状況で、サンプルをできるだけ無駄にしないで、良い（と思える）回帰ができるでしょうか。
今回は、一番簡単なケースで、真の確率pが説明変数xに比例する（$p = wx$, $x>0$）ケースを考えてみます。
面倒なので、値域の話は考えないことにしましょう。（アルゴリズムでフォローします。）
定式化
定義
観測データは、
(X,Y)=(x_n,y_n)_{n=1}^{N}\\
x_n \in R^+\\
y_n = (y_{n+}, y_{n-})\\
y_{n+}, y_{n-} \in Z\\
N_n = y_{n+}+y_{n-}
確率モデル
表の確率$p_n=wx_n$でのコインを$N_n$回試行して$y_{n+}$回表が出たサンプルの生成確率（尤度）は二項分布で表すことができ、
p(y|w;x)\equiv Binom(y_+;N, wx)\\
=\left(\begin{matrix}N\\y_+\end{matrix}\right)(wx)^{y_+}(1-wx)^{y_-}\\
p(Y|w;X)=\prod_n p(y_n|w,x_n)
となります。最尤推定の場合、上記をwについて最大化すれば良いですね。
最尤推定
\ln p(Y|w,X) = \sum_n y_{n+}\ln w+y_{n-}\ln(1-wx)+const.\\
\frac{\partial}{\partial w}\ln p(Y|w,X)
= \sum_n \frac{y_{n+}}{w}-\frac{y_{n-}x_n}{1-wx_n}\\
\frac{\partial^2}{\partial w\partial w}\ln p(Y|w,X)
= \sum_n -\frac{y_{n+}}{w^2}-\frac{y_{n-}x_n^2}{(1-wx_n)^2}< 0
二次の勾配が負なので、凸最適化になることが証明できました。
パラメータ一つなので、二分法でも使いましょう。
ここで、一つ注意したいのが、wの取り得る範囲です。
定義から、$1 - wx_n$は全ての$n$について正でなくてはいけないので、$w<1/x_n$となるように最大値を決めます。（0になると発散してしまうので注意が必要です。）
実装
def binom_regression(yp, yn, x):
if yp.sum() + yn.sum() == 0:
return 0, 0 * x
def binom_grad(w):
return sum(yp) / w - (yn * x / (1 - w * x)).sum()
min_w = 1e-5
min_grad = binom_grad(min_w)
if min_grad <= 0:
return 0, 0 * x
max_w = (1 - 1e-5) / x.max()
max_grad = binom_grad(max_w)
if max_grad >= 0:
return max_w, max_w * x
c = 0
while (max_w - min_w > 1e-5) & (c < 100):
mid_w = (max_w + min_w) / 2
mid_grad = binom_grad(mid_w)
if mid_grad > 0:
min_w, min_grad = mid_w, mid_grad
else:
max_w, max_grad = mid_w, mid_grad
c += 1
return mid_w, mid_w * x
数値実験（をやりたい）
いつかやります
解釈
今回定義した確率モデルで意図したかったことは、信頼度が高い（$N$が大きい）サンプルはなるべく正確に、信頼度が低い（$N$が小さい）サンプルはざっくりと再現するような挙動です。
$N$が大きいサンプルでは、最頻値$pN$に大きな確率を割り振るため、$k=pN$となるような$p=wx$を取りやすくなります。$N$が小さいサンプルでは、比較的なだらかな分布になるため、多少外しても構わないような形になります。
イメージとしては、各サンプルが持つ分布の分散に反比例する重みを付与したガウス回帰のような挙動になること思います。ただし、上振れと下振れに対する感度が変わり得るといったところでしょうか。
まとめ
各サンプルが二項変数の場合の確率モデルとその最尤推定法を定義しました。
今回は一番シンプルな形で定義しましたが、例えば$p=wx$の部分をシグモイド関数にすると、もう少し自由度の高いモデルになると思います。（単調増加関数なので、凸性も失われません。）
また、$x$側にも同様に二項変数を使ったモデルもできそうな気がします。$x$が確率変数になって未知パラメータ数がサンプル数を上回るので、追加の仮定や正則化の導入が必須になりそうです。
同じようなモデルやライブラリが見つからなかったので自分で作りましたが、どこかにあるんじゃないかと思っています。