diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region.h index 4851f92def..3a9148d0cd 100644 --- a/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region.h +++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region.h @@ -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 diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region_dev.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region_dev.h index d883956b82..e685d5b800 100644 --- a/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region_dev.h +++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region_dev.h @@ -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::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::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( diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/dcm_table.h b/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/dcm_table.h index a493f0f055..69b80dab73 100644 --- a/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/dcm_table.h +++ b/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/dcm_table.h @@ -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; diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_arguments.h b/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_arguments.h index 617bf0ad6e..283e7b4c3e 100644 --- a/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_arguments.h +++ b/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_arguments.h @@ -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 class DCMTable; +/** @brief The argument list for the mixed logit discrete choice + * model. + */ template 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 diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_dev.h b/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_dev.h index 4157cadc7d..833a3b3c91 100644 --- a/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_dev.h +++ b/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_dev.h @@ -22,20 +22,29 @@ double MixedLogitDCM::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::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( 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::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::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::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::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( integration_samples.size() * (integration_samples.size() - 1)); @@ -164,10 +199,18 @@ double MixedLogitDCM::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::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(core::math::Pow<4, 1>(table_->num_people())); + static_cast(core::math::Pow<4, 1>(table_.num_people())); return gradient_error; } @@ -263,30 +306,34 @@ void MixedLogitDCM::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( table_.num_people() * arguments_in.initial_dataset_sample_rate_); - int num_integration_samples = + int initial_num_integration_samples = std::max( static_cast( 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::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::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::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::Compute( template bool MixedLogitDCM::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::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::ConstructBoostVariableMap_( boost::program_options::value()->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()->default_value(10.0), + "OPTIONAL The maximum trust region radius used in the trust region " + "search." )( "trust_region_search_method", boost::program_options::value()->default_value("cauchy"), @@ -437,6 +539,12 @@ bool MixedLogitDCM::ConstructBoostVariableMap_( std::cerr << "Missing required --num_discrete_choices_per_person_in.\n"; exit(0); } + if((*vm)[ "trust_region_search_method" ].as() != "cauchy" && + (*vm)[ "trust_region_search_method" ].as() != "dogleg" && + (*vm)[ "trust_region_search_method" ].as() != "steihaug") { + std::cerr << "Invalid option specified for --trust_region_search_method.\n"; + exit(0); + } return false; } @@ -490,9 +598,23 @@ void MixedLogitDCM::ParseArguments( // Parse where to output the discrete choice model predictions. arguments_out->predictions_out_ = vm[ "predictions_out" ].as(); - // Parse how to perform the trust region search. - arguments_out->trust_region_search_method_ = - vm[ "trust_region_search_method" ].as(); + // Parse the maximum trust region radius. + arguments_out->max_trust_region_radius_ = + vm[ "max_trust_region_radius" ].as(); + + // Parse how to perform the trust region search: cauchy, dogleg, steihaug + if(vm[ "trust_region_search_method" ].as() == "cauchy") { + arguments_out->trust_region_search_method_ = + core::optimization::TrustRegionSearchMethod::CAUCHY; + } + else if(vm[ "trust_region_search_method" ].as() == "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. diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_distribution.h b/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_distribution.h index 14e687db73..0a56a9f127 100644 --- a/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_distribution.h +++ b/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_distribution.h @@ -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 ¶meters_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 ¶meters_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)$. diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_sampling.h b/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_sampling.h index 8dba43441d..3e5b4d61e4 100644 --- a/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_sampling.h +++ b/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/mixed_logit_dcm_sampling.h @@ -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 =