プロジェクトの場合、現在Stata出力ファイル（.dta）に存在し、古いStataスクリプトから計算されたいくつかの結果を複製する必要があります。プロジェクトの新しいバージョンはPythonで作成する必要があります。
私が苦労している特定の部分は、Stataのxtileコマンドの加重バージョンに基づいて分位ブレークポイントの計算を一致させることです。データポイント間の結びつきは重みとは関係がなく、使用している重みは連続量からのものであるため、結びつきはほとんどありません（テストデータセットには結びつきがありません）。ですから、同点による誤分類はそうではありません。
加重パーセンタイルに関するウィキペディアの記事と、Rのタイプ7分位数を複製する必要がある代替アルゴリズムについて説明しているこの相互検証済みの投稿を読みました。
両方の重み付けされたアルゴリズム（下部のコード）を実装しましたが、Stata出力で計算された分位数とまだ十分に一致していません。
Stataルーチンで使用される特定のアルゴリズムを知っている人はいますか？ドキュメントはこれを明確に説明していませんでした。CDFの平坦な部分で平均をとってそれを反転することについて何かを述べていますが、これは実際のアルゴリズムをほとんど説明しておらず、他の補間を行っているかどうかについてあいまいです。
とは重みを受け入れず、重み付き分位数を実行できないことに注意してくださいnumpy.percentile。通常の等しい重み付き分位数だけです。scipy.stats.mstats.mquantiles私の問題の核心は、ウェイトを使用する必要性にあります。
注：以下の両方の方法をかなり多くデバッグしましたが、バグが見つかった場合は、コメントでバグを提案してください。私はより小さなデータセットで両方の方法をテストしましたが、結果は良好であり、Rが使用している方法を保証できる場合のRの出力とも一致します。コードはまだそれほどエレガントではなく、2つのタイプ間でコピーされすぎていますが、出力が必要なものであると確信したときに、後で修正される予定です。
問題は、Stataが使用する方法がわからないことです。同じデータセットで実行した場合に、xtile以下のコードとStataの間の不一致を減らしたいと考えています。xtile
私が試したアルゴリズム：
import numpy as np
def mark_weighted_percentiles(a, labels, weights, type):
# a is an input array of values.
# weights is an input array of weights, so weights[i] goes with a[i]
# labels are the names you want to give to the xtiles
# type refers to which weighted algorithm.
# 1 for wikipedia, 2 for the stackexchange post.
# The code outputs an array the same shape as 'a', but with
# labels[i] inserted into spot j if a[j] falls in x-tile i.
# The number of xtiles requested is inferred from the length of 'labels'.
# First type, "vanilla" weights from Wikipedia article.
if type == 1:
# Sort the values and apply the same sort to the weights.
N = len(a)
sort_indx = np.argsort(a)
tmp_a = a[sort_indx].copy()
tmp_weights = weights[sort_indx].copy()
# 'labels' stores the name of the x-tiles the user wants,
# and it is assumed to be linearly spaced between 0 and 1
# so 5 labels implies quintiles, for example.
num_categories = len(labels)
breaks = np.linspace(0, 1, num_categories+1)
# Compute the percentile values at each explicit data point in a.
cu_weights = np.cumsum(tmp_weights)
p_vals = (1.0/cu_weights[-1])*(cu_weights - 0.5*tmp_weights)
# Set up the output array.
ret = np.repeat(0, len(a))
if(len(a)<num_categories):
return ret
# Set up the array for the values at the breakpoints.
quantiles = []
# Find the two indices that bracket the breakpoint percentiles.
# then do interpolation on the two a_vals for those indices, using
# interp-weights that involve the cumulative sum of weights.
for brk in breaks:
if brk <= p_vals[0]:
i_low = 0; i_high = 0;
elif brk >= p_vals[-1]:
i_low = N-1; i_high = N-1;
else:
for ii in range(N-1):
if (p_vals[ii] <= brk) and (brk < p_vals[ii+1]):
i_low = ii
i_high = ii + 1
if i_low == i_high:
v = tmp_a[i_low]
else:
# If there are two brackets, then apply the formula as per Wikipedia.
v = tmp_a[i_low] + ((brk-p_vals[i_low])/(p_vals[i_high]-p_vals[i_low]))*(tmp_a[i_high]-tmp_a[i_low])
# Append the result.
quantiles.append(v)
# Now that the weighted breakpoints are set, just categorize
# the elements of a with logical indexing.
for i in range(0, len(quantiles)-1):
lower = quantiles[i]
upper = quantiles[i+1]
ret[ np.logical_and(a>=lower, a<upper) ] = labels[i]
#make sure upper and lower indices are marked
ret[a<=quantiles[0]] = labels[0]
ret[a>=quantiles[-1]] = labels[-1]
return ret
# The stats.stackexchange suggestion.
elif type == 2:
N = len(a)
sort_indx = np.argsort(a)
tmp_a = a[sort_indx].copy()
tmp_weights = weights[sort_indx].copy()
num_categories = len(labels)
breaks = np.linspace(0, 1, num_categories+1)
cu_weights = np.cumsum(tmp_weights)
# Formula from stats.stackexchange.com post.
s_vals = [0.0];
for ii in range(1,N):
s_vals.append( ii*tmp_weights[ii] + (N-1)*cu_weights[ii-1])
s_vals = np.asarray(s_vals)
# Normalized s_vals for comapring with the breakpoint.
norm_s_vals = (1.0/s_vals[-1])*s_vals
# Set up the output variable.
ret = np.repeat(0, N)
if(N < num_categories):
return ret
# Set up space for the values at the breakpoints.
quantiles = []
# Find the two indices that bracket the breakpoint percentiles.
# then do interpolation on the two a_vals for those indices, using
# interp-weights that involve the cumulative sum of weights.
for brk in breaks:
if brk <= norm_s_vals[0]:
i_low = 0; i_high = 0;
elif brk >= norm_s_vals[-1]:
i_low = N-1; i_high = N-1;
else:
for ii in range(N-1):
if (norm_s_vals[ii] <= brk) and (brk < norm_s_vals[ii+1]):
i_low = ii
i_high = ii + 1
if i_low == i_high:
v = tmp_a[i_low]
else:
# Interpolate as in the type 1 method, but using the s_vals instead.
v = tmp_a[i_low] + (( (brk*s_vals[-1])-s_vals[i_low])/(s_vals[i_high]-s_vals[i_low]))*(tmp_a[i_high]-tmp_a[i_low])
quantiles.append(v)
# Now that the weighted breakpoints are set, just categorize
# the elements of a as usual.
for i in range(0, len(quantiles)-1):
lower = quantiles[i]
upper = quantiles[i+1]
ret[ np.logical_and( a >= lower, a < upper ) ] = labels[i]
#make sure upper and lower indices are marked
ret[a<=quantiles[0]] = labels[0]
ret[a>=quantiles[-1]] = labels[-1]
return ret
