13

図書館内のさまざまな本の数の信頼区間を作成するコードを作成しようとしています (また、有益なプロットを作成します)。

私のいとこは小学校に通っており、毎週彼の先生から本が渡されます。それから彼はそれを読み、次の週に別のものを手に入れるのに間に合うようにそれを返します. しばらくすると、彼が以前に読んだ本を入手していることに気付き始め、これは時間の経過とともに徐々に一般的になりました.

図書館にある実際の本の数が N で、教師が一様に無作為に (置き換えて) 1 冊を選び、毎週あなたに渡すとします。週 t に、読んだ本を受け取った回数が x である場合、 https://math.stackexchange.com/questions/に従って図書館の本の数の最尤推定値を作成できます。 615464/how-many-books-are-in-a-library .

例: A、B、C、D、および E の 5 冊の本がある図書館を考えてみましょう。7 週間連続して [A、B、A、C、B、B、D] の本を受け取った場合、x の値 (重複の数) は、それらの各週の後に [0, 0, 1, 1, 2, 3, 3] になります。つまり、7 週間後に、既に 3 回読んだ本を受け取ったことを意味します。

尤度関数を視覚化するために (何が正しいかを理解していると仮定して)、尤度関数をプロットすると思われる次のコードを作成しました。最大値は約 135 です。これは、上記の MSE リンクによると、実際には最尤推定値です。 

from __future__ import division
import random
import matplotlib.pyplot as plt
import numpy as np

#N is the true number of books. t is the number of weeks.unk is the true number of repeats found 
t = 30
unk = 3
def numberrepeats(N, t):
    return t - len(set([random.randint(0,N) for i in xrange(t)]))

iters = 1000
ydata = []
for N in xrange(10,500):
    sampledunk = [numberrepeats(N,t) for i in xrange(iters)].count(unk)
    ydata.append(sampledunk/iters)

print "MLE is", np.argmax(ydata)
xdata = range(10, 500)
print len(xdata), len(ydata)
plt.plot(xdata,ydata)
plt.show()

出力は次のようになります

ここに画像の説明を入力

私の質問は次のとおりです。

  • 95% 信頼区間を取得して図にプロットする簡単な方法はありますか?
  • 平滑化された曲線をプロットに重ねるにはどうすればよいですか?
  • 私のコードを書くべきだったより良い方法はありますか? あまりエレガントではなく、かなり遅いです。

95% 信頼区間を見つけるということは、サンプリングによって得られる経験的最尤推定値 (この例では理論的には 135 になるはずです) が 95% の確率でその範囲内に収まるように、x 軸の範囲を見つけることを意味します。@mbatchkarov が与えた答えは、現在これを正しく行っていません。

https://math.stackexchange.com/questions/656101/how-to-find-a-confidence-interval-for-a-maximum-likelihood-estimateに数学的な答えがあります。

4

3 に答える 3

8

最初の部分は大丈夫そうなので、2 番目と 3 番目のポイントに取り組みます。

scipy.interpolateとスプライン、またはscipy.optimize.curve_fitを使用して、滑らかな曲線に合わせる方法はたくさんあります。個人的にはcurve_fit、独自の関数を指定してパラメーターに適合させることができるので、私は を好みます。

または、パラメトリック関数を学習したくない場合は、numpy.convolveを使用して単純なローリング ウィンドウ スムージングを行うことができます。

コードの品質に関しては、純粋な python で作業しているため、numpy の速度を利用していません。私はあなたの（既存の）コードを次のように書きます：

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt

# N is the true number of books.
# t is the number of weeks.
# unk is the true number of repeats found 
t = 30
unk = 3
def numberrepeats(N, t, iters):
    rand = np.random.randint(0, N, size=(t, iters))
    return t - np.array([len(set(r)) for r in rand])

iters = 1000
ydata = np.empty(500-10)
for N in xrange(10,500):
    sampledunk = np.count_nonzero(numberrepeats(N,t,iters) == unk)
    ydata[N-10] = sampledunk/iters

print "MLE is", np.argmax(ydata)
xdata = range(10, 500)
print len(xdata), len(ydata)
plt.plot(xdata,ydata)
plt.show()

これをさらに最適化することはおそらく可能ですが、この変更により、私のマシンではコードの実行時間が最大 30 秒から最大 2 秒になりました。

于 2014-02-01T16:56:31.613 に答える
6

信頼区間を取得する簡単な (数値的な) 方法は、単純にスクリプトを何度も実行して、推定値がどれだけ変化するかを確認することです。その標準偏差を使用して、信頼区間を計算できます。

時間の都合上、別のオプションとして、N の各値 (私は 2000 を使用) で一連の試行を実行し、それらの試行のランダムなサブサンプリングを使用して、推定量の標準偏差の推定値を取得します。基本的に、これには、試行のサブセットを選択し、そのサブセットを使用して尤度曲線を生成し、その曲線の最大値を見つけて推定量を取得することが含まれます。これを多くのサブセットに対して行うと、推定量の信頼区間を見つけるために使用できる推定量が得られます。私の完全なスクリプトは次のとおりです。

import numpy as np

t = 30
k = 3
def trial(N):
    return t - len(np.unique(np.random.randint(0, N, size=t)))

def trials(N, n_trials):
    return np.asarray([trial(N) for i in xrange(n_trials)])

n_trials = 2000
Ns = np.arange(1, 501)
results = np.asarray([trials(N, n_trials=n_trials) for N in Ns])

def likelihood(results):
    L = (results == 3).mean(-1)

    # boxcar filtering
    n = 10
    L = np.convolve(L, np.ones(n) / float(n), mode='same')

    return L

def max_likelihood_estimate(Ns, results):
    i = np.argmax(likelihood(results))
    return Ns[i]

def max_likelihood(Ns, results):
    # calculate mean from all trials
    mean = max_likelihood_estimate(Ns, results)

    # randomly subsample results to estimate std
    n_samples = 100
    sample_frac = 0.25
    estimates = np.zeros(n_samples)
    for i in xrange(n_samples):
        mask = np.random.uniform(size=results.shape[1]) < sample_frac
        estimates[i] = max_likelihood_estimate(Ns, results[:,mask])

    std = estimates.std()
    sterr = std * np.sqrt(sample_frac) # is this mathematically sound?
    ci = (mean - 1.96*sterr, mean + 1.96*sterr)
    return mean, std, sterr, ci

mean, std, sterr, ci = max_likelihood(Ns, results)
print "Max likelihood estimate: ", mean
print "Max likelihood 95% ci: ", ci

この方法には 2 つの欠点があります。1 つは、同じ試行セットから多くのサブサンプルを取得しているため、推定値が独立していないことです。この影響を抑えるために、各サブセットの結果の 25% のみを使用しました。もう 1 つの欠点は、各サブサンプルがデータの一部にすぎないことです。そのため、これらのサブセットから得られた見積もりは、完全なスクリプトを何度も実行して得られた見積もりよりも分散が大きくなります。これを説明するために、標準偏差を 4 の平方根で割った値として標準誤差を計算しました。これは、サブサンプルの 1 つよりもフル データ セットに 4 倍のデータがあるためです。ただし、これが数学的に正しいかどうかを知るには、モンテカルロ理論に精通していません。スクリプトを何度も実行したところ、結果が妥当であることがわかりました。

最後に、尤度曲線にボックスカー フィルターを使用して、それらを少し滑らかにしました。理想的には、これにより結果が改善されるはずですが、フィルタリングを行っても結果にはかなりのばらつきがありました。全体的な推定量の値を計算するとき、すべての結果から 1 つの尤度曲線を計算し、その最大値を使用する方がよいか (これが私が最終的に行ったことです)、それともすべての平均を使用する方がよいか確信が持てませんでした。サブセット推定量。サブセット推定量の平均を使用すると、フィルタリング後に残る曲線の粗さの一部を相殺するのに役立つ場合がありますが、これについてはわかりません。

于 2014-02-01T19:37:30.710 に答える