1

多項式 (または条件付き) ロジット モデルの OxMetrics コードを開発しようとしています ( http://data.princeton.edu/wws509/notes/c6s2.htmlで説明されています)。私は OxMetrics を初めて使用し、最適化ルーチン (MaxBFGS) がどのように機能するかを理解するのにまだ苦労しています。

どんな助けでも大歓迎です！

// FUNCTION FOR MULTINOMIAL LOGIT MODELLING
//X = Independent variables (e.g. product features)
//Y = Dependent variable (i.e. discrete choices (0/1))
//N = Total number of individuals (N=500)
//T = Number of observations (choice tasks) per respondent (T=20)
//J = Number of products per task (J=2 => choice between two options)
//sv => starting values
//llik => model log-likelihood

LOGIT(const sv, const llik, X, Y, N, T, J)
{
decl num, den, prbj, prbi, maty;
num = reshape(exp(X*sv[0:6]), N*T, J);
den = sumr(num);
prbj = num ./ den;
maty = reshape(Y .== 1, N*T, J);
prbi = sumr(prbj .* maty);
llik[0] = -sumc(log(prbi));
return llik[0];
}

main()
{
decl data, N, T, J, X, Y, sv, dFunc, fit;
data = loadmat("C:/Users/.../Data/data.csv");
X = data[][33] ~ data[][5:10];
Y = data[][12];
N = 500;
T = 20;
J = 2;
sv = <0;0;0;0;0;0;0>;
println ("\nEstimating using MaxBFGS");
fit = MaxBFGS(LOGIT, X=X, Y=Y, N=N, T=T, J=J, &sv, &dFunc, 0, TRUE);
println (MaxConvergenceMsg(fit), " at parameters ", sv', "with function value ", double(dFunc));
}
4

2 に答える 2

1

どうもありがとう、あなたの答えはとても役に立ちました（どういうわけか私の脳はRからOxに切り替えるのに苦労しています）これが混合ロジット（またはランダムパラメータロジット、またはカーネル密度）のOxを投稿するのに適切な場所かどうかはわかりませんが、他の人の役に立つ。もちろん、次のコードを改善するためのアイデア/考え/トリックは大歓迎です!

データセットの宣言形式

decl N = 433;     // Number of respondents
decl R = 100;     // Number of draws
decl T = 8;  // Number of tasks per respondent (!!! same for all respondents)
decl J = 2;   // Number of options per task (!!! same for all tasks)
decl Y = <0>;      // choice variable
decl X = <1;2;3;4;5;6>; // 6 param, including both fixed effects + mean of random effects
decl RC = <2;3;4;5;6>; // 5 param corresponding to SD of random effects
decl H, mY, mX, mRC;

混合ロジット確率を計算する関数

fn_MXL(const beta, const obj, const grad, const hess) 
{
decl i, seqprb, num, den, sllik;
seqprb = zeros(N,R);
for(i = 0; i < R; ++i)
{
num = reshape(exp(mX * beta[0:5][] + mRC .* (H[N*T*J*i:(N*T*J*(i+1))-1][]) * beta[6:][]), N*T, J); // !!! Edit "beta" accordingly
den = sumr(num);
seqprb[][i] = prodr(reshape(sumr(num ./ den .* mY), N, T));
}
sllik = sumc(log(meanr(seqprb)));
obj[0] = double(sllik);
return 1;
}

モデルのコーディング/推定

main()
{
decl data, i, id1, id2, Indiv, prime, Halton, sv, ir, obj;
// 1. Select variables
data = loadmat("C:/Users/.../choice_data.xls");
mY = reshape(data[][Y], N*T, J);
mX = data[][X];
mRC = data[][RC];
delete data;
// 2. Generate halton draws
id1 = id2 = zeros(N,R); // Create ID variables for both respondents and draws to re-order the draws
for(i = 0; i < N; ++i)
{
id1[i][] = range(1,R); // ID for draws
id2[i][] = i + 1; // ID for respondents
}
Indiv = reshape(id1, N*R, 1) ~ reshape(id2, N*R, 1);
Halton = zeros(N*R, sizeof(RC));
prime = <2;3;5;7;11;13;17;19;23;29;31;37;41;43;47;53;59;61;67;71;73;79;83;89;97;101;103;107;109;113>; // List of prime numbers
for(i = 0; i < sizeof(RC); ++i)
{
Halton[][i] = haltonsequence(prime[i], N*R);    // prime number ; number of draws
}
Halton = Indiv ~ Halton;
H = zeros(rows(Halton)*T*J,columns(Halton));     // Duplicate draws (T*J) times to match structure of "mX" (!!! possible to run out of memory)
for(i = 0; i < (T*J); ++i)
{
H[rows(Halton)*i:(rows(Halton)*(i+1)-1)][] = Halton;
}
H = sortbyc(H, <0;1>)[][2:]; // Draws are organised as following: Draw > Indiv (e.g. first 5 rows would correspond to 1st draw for first 5 indiv)
H = quann(H); // Transform halton draws into normal draws
// 3. Estimation
sv = <0;0;0;0;0;0;0;0;0;0;0>;      // Starting values for X and RC
MaxControl(-1, 1, 1);
ir = MaxBFGS(fn_MXL, &sv, &obj, 0, TRUE);
print("\n", MaxConvergenceMsg(ir), " using numerical derivatives", "\nFunction value = ", obj, "; parameters:", sv);
}
于 2016-05-19T00:18:06.527 に答える
0

公式の MaxBFGS ヘルプ (関数で「F1」を押す) またはここを読む必要があります。

この関数は次のように宣言されます。

MaxBFGS(const func, const avP, const adFunc, const amInvHess, const fNumDer);
  • 引数func: 最大化する関数
  • 引数 avP: 入力 --> 開始値の行列、出力 --> 最終係数。したがって、最初に最適化する必要があるすべての引数を単一の行列 (sv変数) に格納する必要があります。
  • 引数amInvHess: 0 に設定すると単純に推定されます。
  • 引数fNumDer: 0 に設定: 解析的な 1 次導関数の場合、1 に設定すると数値的な 1 次導関数を使用します。

最適化する必要がある関数には、次のプロトタイプが必要です。

LOGIT(const vP, const adFunc, const avScore, const amHessian)

簡単な紹介として、次の例 ( OxMetrics7\ox\samples\maximize\probit1.ox) を参照して、MaxBFGS を介してプロビット モデルを推定することができます - これは、公式ドキュメント " Introduction to Ox -Jurgen A. Doornik and Marius Ooms -2006 " に記載されています)。

// file : ...\OxMetrics7\ox\samples\maximize\probit1.ox
#include <oxstd.oxh>
#import <maximize>

decl g_mY;                                  // global data
decl g_mX;                                  // global data


///////////////////////////////////////////////////////////////////
// Possible improvements:
// 1) Use log(tailn(g_mX * vP)) instead of log(1-prob)
//    in fProbit. This is slower, but a bit more accurate.
// 2) Use analytical first derivatives. This is a bit more robust.
//    Add numerical computation of standard errors.
// 3) Use analytical second derivatives and Newton instead
//    of Quasi-Newton (BFGS). Because the log-likelihood is concave,
//    we don't really need the robustness of BFGS.
// 4) Encapsulate in a class to avoid global variables.
// 5) General starting value routine, etc. etc.
//
// probit2.ox implements (2)
// probit3.ox implements (4)
// probit4.ox implements (4) in a more general way
///////////////////////////////////////////////////////////////////

fProbit(const vP, const adFunc, const avScore,
    const amHessian)
{
    decl prob = probn(g_mX * vP);   // vP is column vector

    adFunc[0] = double(
        meanc(g_mY .* log(prob) + (1-g_mY) .* log(1-prob)));

return 1;                           // 1 indicates success
}

main()
{
    decl vp, dfunc, ir;

    print("Probit example 1, run on ", date(), ".\n\n");

    decl mx = loadmat("data/finney.in7");

    g_mY = mx[][0];       // dependent variable: 0,1 dummy
    g_mX = 1 ~ mx[][3:4]; // regressors: 1, Lrate, Lvolume
    delete mx;

    vp = <-0.465; 0.842; 1.439>;        // starting values

    MaxControl(-1, 1, 1);          // print each iteration
                                               // maximize
    ir = MaxBFGS(fProbit, &vp, &dfunc, 0, TRUE);

    print("\n", MaxConvergenceMsg(ir),
        " using numerical derivatives",
        "\nFunction value = ", dfunc * rows(g_mY),
        "; parameters:", vp);
}

Ps: 変数にグローバル変数を使用してX,Y,N,T,J..、コードを改善できます (probit2.ox、probit3.ox を参照...)。

于 2016-05-16T14:20:48.243 に答える