OxMetrics - 条件逻辑模型
OxMetrics - Conditional logit model
我正在尝试为多项式(或条件)logit 模型开发 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));
}
您应该阅读 MaxBFGS 官方帮助(在函数上按 "F1")或 here。
此函数声明为:
MaxBFGS(const func, const avP, const adFunc, const amInvHess, const fNumDer);
- 参数
func
:函数最大化
- 参数
avP
:输入--> 起始值矩阵,输出 -> 最终系数。因此,您应该首先将所有必须优化的参数存储在一个矩阵中(您的 sv
变量)。
- 参数
amInvHess
: 设置为 0 以进行简单估计。
- 参数
fNumDer
: 设置为 0:for 解析一阶导数,设置为 1 以使用数值一阶导数。
需要优化的函数必须具有以下原型:
LOGIT(const vP, const adFunc, const avScore, const amHessian)
作为简单介绍,您可以查看以下示例 (OxMetrics7\ox\samples\maximize\probit1.ox
) 通过 MaxBFGS 估计概率模型 - 它记录在官方文档“Introduction to Ox -Jurgen A. Dornik 和 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...)
非常感谢,您的回答很有帮助(不知何故我的大脑难以从 R 切换到 Ox)
不确定这是否是 post 用于混合 logit(或随机参数 logit,或核密度)的 Ox 的正确位置,但它可能对其他人有用。当然ideas/thoughts/tricks非常欢迎改进以下代码!
声明数据集的格式
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;
}
型号CODING/ESTIMATION
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);
}
我正在尝试为多项式(或条件)logit 模型开发 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));
}
您应该阅读 MaxBFGS 官方帮助(在函数上按 "F1")或 here。
此函数声明为:
MaxBFGS(const func, const avP, const adFunc, const amInvHess, const fNumDer);
- 参数
func
:函数最大化 - 参数
avP
:输入--> 起始值矩阵,输出 -> 最终系数。因此,您应该首先将所有必须优化的参数存储在一个矩阵中(您的sv
变量)。 - 参数
amInvHess
: 设置为 0 以进行简单估计。 - 参数
fNumDer
: 设置为 0:for 解析一阶导数,设置为 1 以使用数值一阶导数。
需要优化的函数必须具有以下原型:
LOGIT(const vP, const adFunc, const avScore, const amHessian)
作为简单介绍,您可以查看以下示例 (OxMetrics7\ox\samples\maximize\probit1.ox
) 通过 MaxBFGS 估计概率模型 - 它记录在官方文档“Introduction to Ox -Jurgen A. Dornik 和 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...)
非常感谢,您的回答很有帮助(不知何故我的大脑难以从 R 切换到 Ox) 不确定这是否是 post 用于混合 logit(或随机参数 logit,或核密度)的 Ox 的正确位置,但它可能对其他人有用。当然ideas/thoughts/tricks非常欢迎改进以下代码!
声明数据集的格式
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;
}
型号CODING/ESTIMATION
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);
}