一部分だけのデータからガウシアンフィッティングを行いたい

Question

tl; dr

下図のように正規分布に従うデータの一部から全体の正規分布を推定したいです。(青点が与えられるデータ、緑が推定したい曲線)

詳細

発生確率がxに対する正規分布に従うような事象があります。その事象のうち、x<aを満たすようなデータのみからガウシアンフィッティングを行いたいです。
ユースケースとしては、時間tに対する正規分布に従うデータがあり、全体のデータが得られる前にそれの最尤値を推定したい、といった形です。

考えたこと

正規分布の累積密度関数はロジスティック曲線と形が似ていることが知られているため、データの累積和をラベルに、経過時間を特徴量に取ったものをロジスティック回帰すればよいのではないかと考えました。
しかし、0/1をラベルにとった二値でのロジスティック回帰しか資料が見つからなかったため、行えませんでした。(多項ロジスティック回帰で行なえこともないが、明らかに非効率であると考えました(ラベルの種類数が特徴量の種類数になってしまうので))

ご教示いただきたいこと

• ラベルを0/1に限定せず、0-1の間で連続的な値をとれるロジスティック回帰を考えることはできるのか。できる場合、使用すれば良いライブラリや実装を添えていただけるとありがたいです。
• ロジスティック回帰によるアプローチ以外で、一般にこのような問題を解決するアプローチがあれば教えていただきたいです。

上の一つでも分かる場合、ご回答をお願いします。

Answer 1

記載のようなケースは、「打ち切りモデル」と呼ばれるモデルで解決できるように思います。打ち切りモデルというのは、本来は存在するデータ（つまり質問の緑線部分）が打ち切られて観測データとして存在しないケースにおいて本来の確率分布のモデルパラメータ（平均・分散）を推計しようというものです。
打ち切りモデルの考え方は、打ち切られてしまった緑線部分を欠損値とみなして欠損値部分のモデルパラメータが存在することを前提に、観測値の尤度が最大になるようにモデルパラメータをフィッティングさせていくというものです。

このモデルを構築するには、確率モデル専用のソフトウェアであるstanが必要です。（もしかするとpymcでもできるかもしれませんが、私は触ったことがないのでわかりません）stanは、pythonとは別もののモデリングソフトウェアですが、pythonからpystanモジュールを介してキックできるようになっているのでstanの使い方をマスターすれば、何とかなるはずです。

それなりに使いこなそうとすると、ベイズ推計・MCMCなどの前提知識が必要となるのであまり手軽とはいえませんが、ほかに方法がなさそうであれば、試してみる価値があるかと思います。

私自身は実装経験がないので、ただちにサンプルコードを示すことができませんが、「stan 打ち切りモデル」で検索すればstanのサンプルコードが見つかるはずです。

stanのチュートリアルに該当するものがあったので、実験しました。

sample = stats.norm.rvs(loc=0,scale=0.5,size=10000)  # 推計したいデータ
Choice = sample[sample<-0.5]   # -0.5以上は観察されないことにした実験用データ

# stan本体のコード
stancode = """
data {
  int<lower=0> N_obs;  // 観測値の件数
  int<lower=0> N_cens; // 打ち切られたデータの件数
  real y_obs[N_obs];   // 観測値
  real<lower=max(y_obs)> U; // 打ち切りポイント
}
parameters {
  real<lower=U> y_cens[N_cens]; // 打ち切られたデータ　必ず打ち切りポイント以上になると設定
  real mu;  // 推計したい平均
  real<lower=0> sigma;  // 推計したい分散
}
model {
  y_obs ~ normal(mu,sigma);  //　観測値が正規分布に従うことを示すモデル
  y_cens ~ normal(mu,sigma); // 欠損値が正規分布に従うことを示すモデル
}
"""
model = pystan.StanModel(model_code=stancode)  # コードのコンパイル

standata = {'N_obs':len(Choice),   # stanに与えるデータの生成
            'N_cens':10000,        # 打ち切られたデータは10000件とした
            'y_obs':Choice,
            'U':np.max(Choice)}

result = model.sampling(data=standata)  # samplingの実行

結果として平均は0.030とそれなりのものが得られました。（観測値の件数が1630件と本来のデータの20％弱からの推計なので、いいほうだということもできそうです）

Answer 2

0-1の連続値（確率や割合など）でフィッティングさせるロジスティック回帰は問題なく成立します。ただし分類向きの実装では0 or 1のバイナリの目的変数しか取らないようにしている場合もありますし（scikit-learnなど）、紹介事例もだいたい分類をやっているでしょう。

私のブログですがこちらとか、

【python】ロジスティック回帰で確率値で学習させる - 静かなる名辞

その中からもリンクしていますがこちらとか、

データ解析・マイニングとR言語

を見てください。

Comments

[hayataka2049]

……ロジスティック回帰ができたとして、どうやって生データから累積分布を求めるんだろう？　累積「度数」なら出せなくはなさそうですが。
累積「度数」を出して「ロジスティック関数」（ロジスティック回帰で使うのとは別、0-1でなくても構わない）に当てはめるか？

[keymoon]

回答いただきありがとうございます。
ちゃんと二値以外でも存在しているんですね。記事まで提示して頂きありがとうございます。(回帰で検索してもほぼ全てがbinaryで困っていました。試してみようと思います。)
累積分布についてですが、今自分が適用したいケースでは全体の起こりうる事象数が確定しているので、途中までで打ち切ったとしても確率分布を出すことができます。(返信が遅くなってしまい申し訳ありません。)