apaam/mlpack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Done for tonight.

Browse Source
This commit is contained in:
(no author)
2011-01-28 04:57:21 +00:00
parent bdc0713f65
commit dd4e3d0bd4
7 changed files with 249 additions and 96 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region.h
+2 -1
View File
@@ -68,7 +68,8 @@ class TrustRegionUtil {
static double ReductionRatio(
const arma::vec &step,
double iterate_function_value, double next_iterate_function_value,
const arma::vec &gradient, const arma::mat &hessian);
const arma::vec &gradient, const arma::mat &hessian,
double *decrease_predicted_by_model);
};
template<typename FunctionType>
fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region_dev.h
+8 -5
View File
@@ -16,7 +16,8 @@ namespace optimization {
double TrustRegionUtil::ReductionRatio(
const arma::vec &step,
double iterate_function_value, double next_iterate_function_value,
const arma::vec &gradient, const arma::mat &hessian) {
const arma::vec &gradient, const arma::mat &hessian,
double *decrease_predicted_by_model) {
// The actual function value decrease.
double function_value_decrease =
@@ -24,10 +25,10 @@ double TrustRegionUtil::ReductionRatio(
// Predicted objective value decrease by the trust region model:
// -g'p-0.5*p'Hp
double decrease_predicted_by_model =
*decrease_predicted_by_model =
- arma::dot(gradient, step) -
0.5 * arma::as_scalar(arma::trans(step) * hessian * step);
return function_value_decrease / decrease_predicted_by_model;
return function_value_decrease / (* decrease_predicted_by_model);
}
void TrustRegionUtil::ComputeSteihaugDirection(
@@ -327,8 +328,9 @@ void TrustRegion<FunctionType>::Optimize(
// Whether to optimize until convergence.
bool optimize_until_convergence = (num_iterations <= 0);
// The initial radius is set to the 10 % of the maximum radius.
double p_norm = 0;
double current_radius = 0.1;
double current_radius = 0.1 * max_radius_;
int it_num;
// Step direction.
@@ -361,10 +363,11 @@ void TrustRegion<FunctionType>::Optimize(
arma::vec next_iterate = (*iterate) + p;
double iterate_function_value = this->Evaluate_(*iterate);
double next_iterate_function_value = this->Evaluate_(next_iterate);
double decrease_predicted_by_model;
double rho =
core::optimization::TrustRegionUtil::ReductionRatio(
p, iterate_function_value, next_iterate_function_value,
gradient, hessian);
gradient, hessian, &decrease_predicted_by_model);
// Update the trust region radius.
core::optimization::TrustRegionUtil::TrustRadiusUpdate(
fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/dcm_table.h
+2 -2
View File
@@ -67,7 +67,7 @@ class DCMTable {
*/
void choice_probabilities(
int person_index, const arma::vec &beta_vector,
arma::vec *choice_probabilities) {
arma::vec *choice_probabilities) const {
int num_discrete_choices = this->num_discrete_choices(person_index);
choice_probabilities->set_size(num_discrete_choices);
@@ -221,7 +221,7 @@ class DCMTable {
*/
void get_attribute_vector(
int person_index, int discrete_choice_index,
arma::vec *attribute_for_discrete_choice_out) {
arma::vec *attribute_for_discrete_choice_out) const {
int index = cumulative_num_discrete_choices_[person_index] +
discrete_choice_index;
fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_arguments.h
+21 -3
View File
@@ -1,4 +1,6 @@
/** @file mixed_logit_dcm_arguments.h
*
* The arguments used for the mixed logit discrete choice model.
*
* @author Dongryeol Lee (dongryel@cc.gatech.edu)
*/
@@ -7,14 +9,20 @@
#define MLPACK_MIXED_LOGIT_DCM_MIXED_LOGIT_DCM_ARGUMENTS_H
#include "core/table/table.h"
#include "core/optimization/trust_region.h"
#include "mlpack/mixed_logit_dcm/mixed_logit_dcm_distribution.h"
namespace mlpack {
namespace mixed_logit_dcm {
/** @brief The discrete choice model table type.
*/
template<typename TableType>
class DCMTable;
/** @brief The argument list for the mixed logit discrete choice
* model.
*/
template<typename TableType>
class MixedLogitDCMArguments {
public:
@@ -56,7 +64,8 @@ class MixedLogitDCMArguments {
/** @brief The trust region search method to be used.
*/
std::string trust_region_search_method_;
enum core::optimization::TrustRegionSearchMethod::SearchType
trust_region_search_method_;
/** @brief Used for determining the stopping condition based on
* the gradient norm.
@@ -72,8 +81,14 @@ class MixedLogitDCMArguments {
*/
double integration_sample_error_threshold_;
/** @brief The maximum trust region radius.
*/
double max_trust_region_radius_;
public:
/** @brief The default constructor.
*/
MixedLogitDCMArguments() {
attribute_table_ = NULL;
num_discrete_choices_per_person_ = NULL;
@@ -83,8 +98,11 @@ class MixedLogitDCMArguments {
gradient_norm_threshold_ = 0;
max_num_integration_samples_per_person_ = 0;
integration_sample_error_threshold_ = 0;
max_trust_region_radius_ = 0;
}
/** @brief The destructor.
*/
~MixedLogitDCMArguments() {
delete attribute_table_;
attribute_table_ = NULL;
@@ -94,7 +112,7 @@ class MixedLogitDCMArguments {
distribution_ = NULL;
}
};
};
};
}
}
#endif
fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_dev.h
+172 -50
View File
@@ -22,20 +22,29 @@ double MixedLogitDCM<TableType>::GradientErrorSecondPart_(
// gradient for the second error component. Again, careful aliasing
// is done here.
arma::vec second_tmp_vector;
second_tmp_vector.set_size(2 *(table_->num_parameters() + 1));
second_tmp_vector.set_size(2 *(table_.num_parameters() + 1));
arma::vec second_tmp_choice_probability_gradient_outer(
second_tmp_vector.memptr() + 1, table_->num_parameters(), false);
second_tmp_vector.memptr() + 1, table_.num_parameters(), false);
arma::vec second_tmp_choice_probability_gradient_inner(
second_tmp_vector.memptr() + table_->num_parameters() + 2,
table_->num_parameters(), false);
second_tmp_vector.memptr() + table_.num_parameters() + 2,
table_.num_parameters(), false);
// Temporary vector.
arma::vec outer_choice_probabilities;
arma::vec inner_choice_probabilities;
arma::vec outer_choice_prob_weighted_attribute_vector;
arma::vec inner_choice_prob_weighted_attribute_vector;
// The quantity to be eventually returned.
double second_part = 0;
// The outer sum loop.
for(int i = 0; i < table_->num_people(); i++) {
for(int i = 0; i < table_.num_people(); i++) {
// Get the outer person index.
int outer_person_index = table_->shuffled_indices_for_person(i);
// Get the outer person index and its discrete choice index.
int outer_person_index = table_.shuffled_indices_for_person(i);
int outer_discrete_choice_index =
table_.get_discrete_choice_index(outer_person_index);
// Get the simulated choice probability and the simulated choice
// probability gradient for the given outer person.
@@ -47,16 +56,18 @@ double MixedLogitDCM<TableType>::GradientErrorSecondPart_(
// The integration samples for the given outer person.
const std::vector< arma::vec > &integration_samples =
table_->integration_samples(outer_person_index);
sample.integration_samples(outer_person_index);
double normalization_factor =
1.0 / static_cast<double>(
integration_samples.size() * (integration_samples.size() - 1));
// The inner sum loop.
for(unsigned int k = i + 1; k < table_->num_people(); k++) {
for(unsigned int k = i + 1; k < table_.num_people(); k++) {
// Get the inner person index.
int inner_person_index = table_->shuffled_indices_for_person(k);
// Get the inner person index and its discrete choice.
int inner_person_index = table_.shuffled_indices_for_person(k);
int inner_discrete_choice_index =
table_.get_discrete_choice_index(inner_person_index);
// Get the simulated choice probability and the simulated choice
// probability gradient for the given inner person.
@@ -68,7 +79,7 @@ double MixedLogitDCM<TableType>::GradientErrorSecondPart_(
// Compute delta_hik (Equation 2.3).
arma::vec delta_hik;
delta_hik.set_size(2 *(table_->num_parameters() + 1));
delta_hik.set_size(2 *(table_.num_parameters() + 1));
double first_factor =
-1.0 / (
core::math::Sqr(outer_simulated_choice_probability) *
@@ -86,24 +97,42 @@ double MixedLogitDCM<TableType>::GradientErrorSecondPart_(
outer_simulated_choice_probability_gradient,
inner_simulated_choice_probability_gradient);
delta_hik[0] = first_factor * dot_product;
delta_hik.submat(1, table_->num_parameters(), 0, 0) =
delta_hik.submat(1, table_.num_parameters(), 0, 0) =
second_factor * outer_simulated_choice_probability_gradient;
delta_hik[table_->num_parameters() + 1] = third_factor * dot_product;
delta_hik[table_.num_parameters() + 1] = third_factor * dot_product;
delta_hik.submat(
table_->num_parameters() + 2, table_->num_parameters(), 0, 0) =
table_.num_parameters() + 2, table_.num_parameters(), 0, 0) =
second_factor * inner_simulated_choice_probability_gradient;
// Loop through each integration sample.
for(unsigned int j = 0; j < integration_samples.size(); j++) {
const arma::vec &integration_sample = integration_samples[j];
table_.choice_probabilities(
outer_person_index, integration_sample, &outer_choice_probabilities);
second_tmp_vector[0] =
table_->choice_probability(outer_person_index, integration_sample);
table_->choice_probability_gradient(
outer_person_index, &second_tmp_choice_probability_gradient_outer);
second_tmp_vector[ table_->num_parameters() + 1] =
table_->choice_probability(inner_person_index, integration_sample);
table_->choice_probability_gradient(
inner_person_index, &second_tmp_choice_probability_gradient_inner);
outer_choice_probabilities[ outer_discrete_choice_index ];
table_.distribution()->ChoiceProbabilityWeightedAttributeVector(
table_, outer_person_index, outer_choice_probabilities,
&outer_choice_prob_weighted_attribute_vector);
table_.distribution()->ChoiceProbabilityGradientWithRespectToParameter(
sample.parameters(), table_,
outer_person_index, outer_choice_probabilities,
outer_choice_prob_weighted_attribute_vector,
&second_tmp_choice_probability_gradient_outer);
table_.choice_probabilities(
inner_person_index, integration_sample, &inner_choice_probabilities);
second_tmp_vector[ table_.num_parameters() + 1] =
inner_choice_probabilities[ inner_discrete_choice_index ];
table_.distribution()->ChoiceProbabilityWeightedAttributeVector(
table_, inner_person_index, inner_choice_probabilities,
&inner_choice_prob_weighted_attribute_vector);
table_.distribution()->ChoiceProbabilityGradientWithRespectToParameter(
sample.parameters(), table_,
inner_person_index, inner_choice_probabilities,
inner_choice_prob_weighted_attribute_vector,
&second_tmp_choice_probability_gradient_inner);
// Take the dot product between the two vectors and square it.
second_part +=
@@ -125,15 +154,21 @@ double MixedLogitDCM<TableType>::GradientErrorFirstPart_(
// The temporary vector for extracting choice probability and its
// gradient. Careful aliasing is done here.
arma::vec first_tmp_vector;
first_tmp_vector.set_size(table_->num_parameters() + 1);
first_tmp_vector.set_size(table_.num_parameters() + 1);
arma::vec first_tmp_choice_probability_gradient(
first_tmp_vector.memptr() + 1, table_->num_parameters(), false);
first_tmp_vector.memptr() + 1, table_.num_parameters(), false);
// The choice probability vector (temporary space) and the attribute
// vector weighted by it.
arma::vec choice_probabilities;
arma::vec choice_prob_weighted_attribute_vector;
double first_part = 0;
for(int i = 0; i < table_->num_people(); i++) {
for(int i = 0; i < table_.num_people(); i++) {
// Get the person index.
int person_index = table_->shuffled_indices_for_person(i);
// Get the person index and its discrete choice.
int person_index = table_.shuffled_indices_for_person(i);
int discrete_choice_index = table_.get_discrete_choice_index(person_index);
// Get the simulated choice probability and the simulated choice
// probability gradient for the given/ person.
@@ -156,7 +191,7 @@ double MixedLogitDCM<TableType>::GradientErrorFirstPart_(
// The integration samples for the given person.
const std::vector< arma::vec > &integration_samples =
table_->integration_samples(person_index);
sample.integration_samples(person_index);
double normalization_factor =
1.0 / static_cast<double>(
integration_samples.size() * (integration_samples.size() - 1));
@@ -164,10 +199,18 @@ double MixedLogitDCM<TableType>::GradientErrorFirstPart_(
// Loop through each integration sample.
for(unsigned int j = 0; j < integration_samples.size(); j++) {
const arma::vec &integration_sample = integration_samples[j];
first_tmp_vector[0] =
table_->choice_probability(person_index, integration_sample);
table_->choice_probability_gradient(
person_index, &first_tmp_choice_probability_gradient);
// Get the choice probability vector and its weighted version.
table_.choice_probabilities(
person_index, integration_sample, &choice_probabilities);
table_.distribution()->ChoiceProbabilityWeightedAttributeVector(
table_, person_index, choice_probabilities,
&choice_prob_weighted_attribute_vector);
first_tmp_vector[0] = choice_probabilities[discrete_choice_index];
table_.distribution()->ChoiceProbabilityGradientWithRespectToParameter(
sample.parameters(), table_, person_index,
choice_probabilities, choice_prob_weighted_attribute_vector,
&first_tmp_choice_probability_gradient);
// Take the dot product between the two vectors and square it.
first_part +=
@@ -192,7 +235,7 @@ double MixedLogitDCM<TableType>::GradientError_(
// Divide by the normalization term, which is the total number of
// people in the dataset raised to the 4-th power.
gradient_error /=
static_cast<double>(core::math::Pow<4, 1>(table_->num_people()));
static_cast<double>(core::math::Pow<4, 1>(table_.num_people()));
return gradient_error;
}
@@ -263,30 +306,34 @@ void MixedLogitDCM<TableType>::Compute(
mlpack::mixed_logit_dcm::MixedLogitDCMResult *result_out) {
// Here is the main entry of the algorithm.
int num_data_samples =
int initial_num_data_samples =
static_cast<int>(
table_.num_people() *
arguments_in.initial_dataset_sample_rate_);
int num_integration_samples =
int initial_num_integration_samples =
std::max(
static_cast<int>(
arguments_in.initial_integration_sample_rate_ *
arguments.max_num_integration_samples_per_person_), 36);
arguments_in.max_num_integration_samples_per_person_), 36);
// Initialize the starting optimization parameter $\theta_0$ and its
// associated sampling information.
arma::vec gradient;
arma::mat hessian;
SamplingType iterate;
parameters_for_init.zeros(table_.num_parameters());
iterate.Init(
&table_, num_data_samples, num_integration_samples);
&table_, initial_num_data_samples, initial_num_integration_samples);
iterate.parameters().zeros(table_.num_parameters());
iterate.NegativeSimulatedLogLikelihoodGradient(&gradient);
iterate.NegativeSimulatedLogLikelihoodHessian(&hessian);
// The step direction.
// The step direction and its 2-norm.
arma::vec p;
double p_norm;
// The initial trust region radius is set to the 10 % of the maximum
// trust region radius.
double current_radius = 0.1 * arguments_in.max_trust_region_radius_;
// Enter the trust region loop.
do {
@@ -294,7 +341,8 @@ void MixedLogitDCM<TableType>::Compute(
// Obtain the step direction by solving Equation 4.3
// approximately.
core::optimization::TrustRegionUtil::ObtainStepDirection(
search_method_, current_radius, gradient, hessian, &p, &p_norm);
arguments_in.trust_region_search_method_, current_radius, gradient,
hessian, &p, &p_norm);
// Get the reduction ratio rho (Equation 4.4)
SamplingType next_iterate;
@@ -302,10 +350,11 @@ void MixedLogitDCM<TableType>::Compute(
double iterate_function_value = iterate.NegativeSimulatedLogLikelihood();
double next_iterate_function_value =
next_iterate.NegativeSimulatedLogLikelihood();
double decrease_predicted_by_model;
double model_reduction_ratio =
core::optimization::TrustRegionUtil::ReductionRatio(
p, iterate_function_value, next_iterate_function_value,
gradient, hessian);
gradient, hessian, &decrease_predicted_by_model);
// Compute the data sample error and the integration sample error.
double data_sample_error = this->DataSampleError_(iterate, next_iterate);
@@ -314,11 +363,54 @@ void MixedLogitDCM<TableType>::Compute(
// Determine whether the termination condition has been reached.
if(TerminationConditionReached_(
arguments_in, model_reduction_ratio,
arguments_in, decrease_predicted_by_model,
data_sample_error, integration_sample_error,
theta_sampling, gradient)) {
iterate, gradient)) {
break;
}
// If we are not ready to terminate yet, then consider one of the
// three options. (1) Increase the data sample size; (2) Increase
// the integration sample size; (3) Do a step to the new iterate.
// Increase the data sample size.
if(iterate.num_active_people() < table_.num_people() &&
sqrt(data_sample_error) >= sqrt(integration_sample_error) &&
fabs(decrease_predicted_by_model) <
arguments_in.gradient_norm_threshold_ * sqrt(data_sample_error)) {
int max_allowable_people = table_.num_people() -
iterate.num_active_people();
int predicted_increase =
core::math::Sqr(
arguments_in.gradient_norm_threshold_ *
sqrt(data_sample_error) / decrease_predicted_by_model) *
iterate.num_active_people() - iterate.num_active_people();
int num_additional_people =
std::min(std::max(1, predicted_increase), max_allowable_people);
iterate.AddActivePeople(
num_additional_people, initial_num_integration_samples);
// Update the gradient and the hessian.
iterate.NegativeSimulatedLogLikelihoodGradient(&gradient);
iterate.NegativeSimulatedLogLikelihoodHessian(&hessian);
continue;
}
// Increase the integration sample size.
if(sqrt(integration_sample_error) > sqrt(data_sample_error) &&
fabs(decrease_predicted_by_model) <
arguments_in.gradient_norm_threshold_ *
sqrt(integration_sample_error)) {
// Update the gradient and the hessian.
iterate.NegativeSimulatedLogLikelihoodGradient(&gradient);
iterate.NegativeSimulatedLogLikelihoodHessian(&hessian);
continue;
}
// Now the third case: adjust the trust region radius and make the
// next iterate based on the trust region optimization.
}
while(true);
}
@@ -326,7 +418,7 @@ void MixedLogitDCM<TableType>::Compute(
template<typename TableType>
bool MixedLogitDCM<TableType>::TerminationConditionReached_(
const ArgumentType &arguments_in,
double model_reduction_ratio, double data_sample_error,
double predicted_objective_value_improvement, double data_sample_error,
double integration_sample_error, const SamplingType &sampling,
const arma::vec &gradient) const {
@@ -334,15 +426,20 @@ bool MixedLogitDCM<TableType>::TerminationConditionReached_(
// people in the sampling (outer term consists of all people).
if(sampling.num_active_people() == table_.num_people()) {
// If the predictived improvement in the objective value is less
// If the predicted improvement in the objective value is less
// than the integration sample error and the integration sample
// error is small,
if(model_reduction_ratio < c_factor *) {
if(predicted_objective_value_improvement <
arguments_in.gradient_norm_threshold_ * sqrt(integration_sample_error) &&
sqrt(integration_sample_error) <
arguments_in.integration_sample_error_threshold_) {
// Compute the gradient error.
double gradient_error = this->GradientError_(sampling);
if() {
if(arma::dot(gradient, gradient) +
arguments_in.gradient_norm_threshold_ *
sqrt(gradient_error) <= 0.001) {
return true;
}
}
@@ -395,6 +492,11 @@ bool MixedLogitDCM<TableType>::ConstructBoostVariableMap_(
boost::program_options::value<double>()->default_value(0.01),
"OPTIONAL The threshold for determining whether the integration sample "
"error is small or not."
)(
"max_trust_region_radius",
boost::program_options::value<double>()->default_value(10.0),
"OPTIONAL The maximum trust region radius used in the trust region "
"search."
)(
"trust_region_search_method",
boost::program_options::value<std::string>()->default_value("cauchy"),
@@ -437,6 +539,12 @@ bool MixedLogitDCM<TableType>::ConstructBoostVariableMap_(
std::cerr << "Missing required --num_discrete_choices_per_person_in.\n";
exit(0);
}
if((*vm)[ "trust_region_search_method" ].as<std::string>() != "cauchy" &&
(*vm)[ "trust_region_search_method" ].as<std::string>() != "dogleg" &&
(*vm)[ "trust_region_search_method" ].as<std::string>() != "steihaug") {
std::cerr << "Invalid option specified for --trust_region_search_method.\n";
exit(0);
}
return false;
}
@@ -490,9 +598,23 @@ void MixedLogitDCM<TableType>::ParseArguments(
// Parse where to output the discrete choice model predictions.
arguments_out->predictions_out_ = vm[ "predictions_out" ].as<std::string>();
// Parse how to perform the trust region search.
arguments_out->trust_region_search_method_ =
vm[ "trust_region_search_method" ].as<std::string>();
// Parse the maximum trust region radius.
arguments_out->max_trust_region_radius_ =
vm[ "max_trust_region_radius" ].as<double>();
// Parse how to perform the trust region search: cauchy, dogleg, steihaug
if(vm[ "trust_region_search_method" ].as<std::string>() == "cauchy") {
arguments_out->trust_region_search_method_ =
core::optimization::TrustRegionSearchMethod::CAUCHY;
}
else if(vm[ "trust_region_search_method" ].as<std::string>() == "dogleg") {
arguments_out->trust_region_search_method_ =
core::optimization::TrustRegionSearchMethod::DOGLEG;
}
else {
arguments_out->trust_region_search_method_ =
core::optimization::TrustRegionSearchMethod::STEIHAUG;
}
// The number of parameters that generate each $\beta$ is fixed now
// as the Gaussian example in Appendix.
fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_distribution.h
+41 -29
View File
@@ -55,19 +55,19 @@ class MixedLogitDCMDistribution {
virtual void Init(const std::string &file_name) const = 0;
void ChoiceProbabilityWeightedAttributeVector(
DCMTableType *dcm_table_in, int person_index,
const DCMTableType &dcm_table_in, int person_index,
const arma::vec &choice_probabilities,
arma::vec *choice_prob_weighted_attribute_vector) const {
// Get the number of discrete choices for the given person.
int num_discrete_choices = dcm_table_in->num_discrete_choices(
int num_discrete_choices = dcm_table_in.num_discrete_choices(
person_index);
choice_prob_weighted_attribute_vector->set_size(
dcm_table_in->num_attributes());
dcm_table_in.num_attributes());
choice_prob_weighted_attribute_vector->zeros();
for(int i = 0; i < choice_probabilities.n_elem; i++) {
arma::vec attribute_vector;
dcm_table_in->get_attribute_vector(person_index, i, &attribute_vector);
dcm_table_in.get_attribute_vector(person_index, i, &attribute_vector);
(*choice_prob_weighted_attribute_vector) +=
choice_probabilities[i] * attribute_vector;
}
@@ -78,15 +78,19 @@ class MixedLogitDCMDistribution {
* \beta^{\nu}(\theta))$ (Equation 8.5)
*/
double ChoiceProbabilityHessianWithRespectToAttribute(
DCMTableType *dcm_table_in,
int person_index, int discrete_choice_index,
const DCMTableType &dcm_table_in, int person_index,
int row_index, int col_index,
const arma::vec &choice_probabilities,
const arma::vec &choice_prob_weighted_attribute_vector) const {
// Get the discrete choice index.
int discrete_choice_index =
dcm_table_in.get_discrete_choice_index(person_index);
// Get the choice probability for the chosen discrete choice.
double choice_probability = choice_probabilities[discrete_choice_index];
arma::vec discrete_choice_attribute_vector;
dcm_table_in->get_attribute_vector(
dcm_table_in.get_attribute_vector(
person_index, discrete_choice_index, &discrete_choice_attribute_vector);
// The (row, col)-th entry of $\bar{X}_i res_{i,j_i^*}(\beta) (
@@ -106,11 +110,11 @@ class MixedLogitDCMDistribution {
// The (row, col)-th entry of $\bar{X}_i \tilde{P}_i(\beta)
// \bar{X}_i'$
double third_part = 0;
int num_discrete_choices = dcm_table_in->num_discrete_choices(
int num_discrete_choices = dcm_table_in.num_discrete_choices(
person_index);
for(int i = 0; i < num_discrete_choices; i++) {
arma::vec attribute_vector;
dcm_table_in->get_attribute_vector(person_index, i, &attribute_vector);
dcm_table_in.get_attribute_vector(person_index, i, &attribute_vector);
third_part += attribute_vector[row_index] * choice_probabilities[i] *
attribute_vector[col_index];
}
@@ -120,20 +124,19 @@ class MixedLogitDCMDistribution {
}
void ChoiceProbabilityHessianWithRespectToAttribute(
DCMTableType *dcm_table_in,
int person_index, int discrete_choice_index,
const DCMTableType &dcm_table_in, int person_index,
const arma::vec &choice_probabilities,
const arma::vec &choice_prob_weighted_attribute_vector,
arma::mat *hessian_out) const {
hessian_out->set_size(
dcm_table_in->num_attributes(), dcm_table_in->num_attributes());
for(int j = 0; j < dcm_table_in->num_attributes(); j++) {
for(int i = 0; i < dcm_table_in->num_attributes(); i++) {
dcm_table_in.num_attributes(), dcm_table_in.num_attributes());
for(int j = 0; j < dcm_table_in.num_attributes(); j++) {
for(int i = 0; i < dcm_table_in.num_attributes(); i++) {
hessian_out->at(
i, j) =
this->ChoiceProbabilityHessianWithRespectToAttribute(
dcm_table_in, person_index, discrete_choice_index, i, j,
dcm_table_in, person_index, i, j,
choice_probabilities, choice_prob_weighted_attribute_vector);
}
}
@@ -145,13 +148,16 @@ class MixedLogitDCMDistribution {
*/
void HessianProducts(
const arma::vec &parameters_in,
DCMTableType *dcm_table_in,
int person_index, int discrete_choice_index,
const DCMTableType &dcm_table_in, int person_index,
const arma::vec &choice_probabilities,
const arma::vec &choice_prob_weighted_attribute_vector,
arma::mat *hessian_first_part,
arma::vec *hessian_second_part) const {
// Get the discrete choice index.
int discrete_choice_index =
dcm_table_in.get_discrete_choice_index(person_index);
// Compute $\frac{\partial}{\partial \theta} \beta^{\nu}(\theta)
// \frac{\partial^2}{\partial \beta^2} P_{i j_i^*} (
// \beta^{\nu}(\theta)) ( \frac{\partial}{\partial \theta}
@@ -159,10 +165,10 @@ class MixedLogitDCMDistribution {
arma::mat first;
arma::mat second;
this->AttributeGradientWithRespectToParameter(
parameters_in, dcm_table_in->num_attributes(), &first);
parameters_in, dcm_table_in.num_attributes(), &first);
this->ChoiceProbabilityHessianWithRespectToAttribute(
dcm_table_in, person_index, discrete_choice_index,
choice_probabilities, choice_prob_weighted_attribute_vector, &second);
dcm_table_in, person_index, choice_probabilities,
choice_prob_weighted_attribute_vector, &second);
(*hessian_first_part) = first * second * arma::trans(first);
// Compute $\frac{\partial}{\partial \theta}
@@ -170,8 +176,8 @@ class MixedLogitDCMDistribution {
// \beta^{\nu}(\theta))$.
arma::vec third;
this->ChoiceProbabilityGradientWithRespectToAttribute(
dcm_table_in, person_index, discrete_choice_index,
choice_probabilities, choice_prob_weighted_attribute_vector, &third);
dcm_table_in, person_index, choice_probabilities,
choice_prob_weighted_attribute_vector, &third);
(*hessian_second_part) = first * third;
}
@@ -179,15 +185,18 @@ class MixedLogitDCMDistribution {
* P_{i,j}(\beta)$ (Equation 8.2).
*/
void ChoiceProbabilityGradientWithRespectToAttribute(
DCMTableType *dcm_table_in,
int person_index, int discrete_choice_index,
const DCMTableType &dcm_table_in, int person_index,
const arma::vec &choice_probabilities,
const arma::vec &choice_prob_weighted_attribute_vector,
arma::vec *gradient_out) const {
// Get the discrete choice inde.
int discrete_choice_index =
dcm_table_in.get_discrete_choice_index(person_index);
// Get the discrete choice attribute vector.
arma::vec discrete_choice_attribute;
dcm_table_in->get_attribute_vector(
dcm_table_in.get_attribute_vector(
person_index, discrete_choice_index, &discrete_choice_attribute);
(*gradient_out) = discrete_choice_attribute -
choice_prob_weighted_attribute_vector;
@@ -205,19 +214,22 @@ class MixedLogitDCMDistribution {
*/
void ChoiceProbabilityGradientWithRespectToParameter(
const arma::vec &parameters_in,
DCMTableType *dcm_table_in,
int person_index, int discrete_choice_index,
const arma::vec &choice_probabilities,
const DCMTableType &dcm_table_in,
int person_index, const arma::vec &choice_probabilities,
const arma::vec &choice_prob_weighted_attribute_vector,
arma::vec *product_out) const {
// Find the person's discrete choice index.
int discrete_choice_index =
dcm_table_in.get_discrete_choice_index(person_index);
// Initialize the product.
product_out->set_size(this->num_parameters());
// Compute the gradient matrix times vector for the person's
// discrete choice.
arma::vec attribute_vector;
dcm_table_in->get_attribute_vector(
dcm_table_in.get_attribute_vector(
person_index, discrete_choice_index, &attribute_vector);
// Compute $\bar{X}_i res_{i,j_i^*}(\beta)$.
fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_sampling.h
+3 -6
View File
@@ -146,7 +146,7 @@ class MixedLogitDCMSampling {
dcm_table_->get_discrete_choice_index(person_index);
dcm_table_->distribution()->
ChoiceProbabilityGradientWithRespectToParameter(
parameters_, dcm_table_, person_index, discrete_choice_index,
parameters_, *dcm_table_, person_index,
beta_vector, choice_probabilities, &beta_gradient_product);
// Update the simulated choice probabilities
@@ -166,8 +166,7 @@ class MixedLogitDCMSampling {
arma::mat hessian_first_part;
arma::vec hessian_second_part;
dcm_table_->distribution()->HessianProducts(
parameters_, dcm_table_, person_index,
discrete_choice_index, beta_vector,
parameters_, *dcm_table_, person_index, beta_vector,
choice_probabilities, &hessian_first_part, &hessian_second_part);
simulated_loglikelihood_hessians_[person_index].first.push_back(
hessian_first_part);
@@ -188,7 +187,7 @@ class MixedLogitDCMSampling {
// Get the index of the active person.
int person_index = dcm_table_->shuffled_indices_for_person(i);
for(int j = simulated_choice_probabilities_[j].num_samples();
for(int j = simulated_choice_probabilities_[i].num_samples();
j < num_integration_samples_[person_index]; j++) {
// Draw a beta from the parameter theta and add it to the
@@ -300,8 +299,6 @@ class MixedLogitDCMSampling {
int person_index = dcm_table_->shuffled_indices_for_person(i);
// Get the simulated choice probability for the given person.
int discrete_choice_index =
dcm_table_->get_discrete_choice_index(person_index);
double simulated_choice_probability =
this->simulated_choice_probability(person_index);
double inverse_simulated_choice_probability =