randomvariable.cc.svn-base

Probabilistic graphical models in matlab.

}


unsigned int DiscreteRandomVariable::sample()
{
	unsigned int i = internal_sample();
	
	++prior_samples[i];
	
	return i;
}

unsigned int DiscreteRandomVariable::sample_without_recording ()
{
	return internal_sample();
}

unsigned int DiscreteRandomVariable::sample_without_recording (const vector <double> & sp_vec)
{
	assert ( sp_vec.size() == sampling_probabilities.size() );
	
	if ( is_conditionned () )
	{
		return last_sampled_value;
	}
	
	copy(sp_vec.begin(), sp_vec.end(), sampling_probabilities.begin() );
	
	return internal_sample();
}

unsigned int DiscreteRandomVariable::sample_from_particles (const vector <unsigned int> & positions, const vector <double> & weights)
{
	assert ( positions.size() == weights.size());
	
	if ( is_conditionned () )
	{
		return last_sampled_value;
	}
	
	unsigned int i;
	
	double a = 0;
	double random_number = uniform_distribution_0_1();
	
	for (i = 0; i < weights.size(); ++i )
	{
		a+= weights[i];
		
		if (random_number < a)
		{
			break;
		}
		
	}
	
	assert ( i < weights.size());
	
	++prior_samples[positions[i]];
	
	last_sampled_value = positions[i];
	
	return positions[i];	
}

// Private methods

void DiscreteRandomVariable::set_value_at_random()
{	
	if (conditionned) // we don't want to change the value if it is already observed!
	{
		return;
	}
	
	uniform_int <int> uni_int_dist(0, number_values-1);
	
	variate_generator <boost::lagged_fibonacci607 &, uniform_int <int> > vg (fibonnacci_number_generator, uni_int_dist);
	
	last_sampled_value = vg();	
}

unsigned int DiscreteRandomVariable::internal_sample()
{
	if ( is_conditionned () )
	{
		return last_sampled_value;
	}
	
	unsigned int i;
	
	double a = 0;
	double random_number = uniform_distribution_0_1();
	
	for (i = 0; i < number_values; ++i )
	{
		a+= sampling_probabilities[i];
		
		if (random_number < a)
		{
			break;
		}
		
	}
	
	assert ( i < number_values);
	
	last_sampled_value = i;
	
	return i;
}

// Unknown methods

vector < DiscreteRandomVariable *> DiscreteRandomVariable::uncombine()
{
	vector < DiscreteRandomVariable *> a;
	cout << "ERROR (CRITICAL): we called uncombine() from the DiscreteRandomVariable class" << endl;
	return a;
}

void DiscreteRandomVariable::compute_mean()
{
	if (particle_mean)
	{
		//cout << "particle_mean" << endl;
		
		double a(0.0);
		
		for ( unsigned int i = 0; i <  particle_means_vec.size();++i)
		{		
			// cout << "Particle " << i << " " << particle_vec[i] << endl;
			
			a += particle_means_vec[i];
		}
		
		mean = a / particle_means_vec.size();
		
		return;
	}
	
	
	double a(0.0);
	unsigned int b(0);
	
	for ( unsigned int i = 0; i <  prior_samples.size();++i)
	{
		b += prior_samples[i];;
	}
	
	for ( unsigned int i = 0; i <  prior_samples.size();++i)
	{
		a += (double) (i * prior_samples[i])/b;
	}
	
	mean = a;
}

void DiscreteRandomVariable::compute_true_mean()
{
	double a(0.0);
	
	for ( unsigned int i = 0; i <  sampling_probabilities.size();++i)
	{
		a += i * sampling_probabilities[i];
	}
	
	mean = a;
}

// Testing purposes

void DiscreteRandomVariable::set_mean_particles( vector <unsigned int> & particle_vec )
{
	double a(0.0);
	
	for ( unsigned int i = 0; i <  particle_vec.size();++i)
	{		
		// cout << "Particle " << i << " " << particle_vec[i] << endl;
		
		a += particle_vec[i];
	}
	
	//mean = a / particle_vec.size();
		
	particle_mean = true;

	//cout << "Mean (" << index << " ) " << mean << endl;
	
	particle_means_vec.push_back(a / particle_vec.size());
	
}

void DiscreteRandomVariable::set_mean_bw_particles( vector <unsigned int> & particle_vec, vector <double> & weights_vec )
{
	double a(0.0);
	
	for ( unsigned int i = 0; i <  particle_vec.size();++i)
	{
		//cout << "Particle " << i << " " << particle_vec[i] << " " << weights_vec[i] << endl;
		
		if (index == 55)
		{
			
			//cout << particle_vec[i] << " " << index << " ) " << weights_vec[i] << endl;
			
		}
		
		a += (double) (particle_vec[i]) * weights_vec[i];
	}
	
	//mean = a / particle_vec.size();
	
	particle_mean = true;
	
	if (index == 55)
	{
	
	//cout << "******************Mean (" << index << " ) " << a<< endl;
	
	}
	
	particle_means_vec.push_back(a);
	
}

// Meant to be virtual later

bool DiscreteRandomVariable::loop_over ()
{	
	if ( end_loop)
	{	
		end_loop = false;
		return false;
	}
	
	if ( conditionned)
	{
		//cout << "I am conditionned" << endl;
		end_loop = true;
		return true;
	}
	
	if (start_loop)
	{
		start_loop = false;
		last_sampled_value = 0;
	}
	else
	{
		++last_sampled_value;
	}
	
	if (last_sampled_value +1 == number_values)
	{
		start_loop = true;
		end_loop = true;
	}
	
	return true;
}

void DiscreteRandomVariable::serialize_marginals (std::ofstream & output_file) const
{
	output_file << "( " << index << " ) ( " ;
	
	for (vector<double>::const_iterator it = sampling_probabilities.begin(); it != sampling_probabilities.end(); ++it)
	{
		output_file << *it << " ";
	}
	
	output_file << " ) ";
}

#ifndef NDEBUG

void DiscreteRandomVariable::debug_display_probabilities() const
{
	cout << "( " << index << " ) ( " ;
	
	for (vector<double>::const_iterator it = sampling_probabilities.begin(); it != sampling_probabilities.end(); ++it)
	{
		cout << *it << " ";
	}
	
	cout << " ) " << endl;
}

#endif


//****************** End of DiscreteRandomVariable Implementation ********************* //

//****************** CombinedDiscreteRandomVariable Implementation ********************* //

CombinedDiscreteRandomVariable::CombinedDiscreteRandomVariable( list<DiscreteRandomVariable *> l, unsigned int b): DiscreteRandomVariable(b) // constructor, setting the index AND the possible number of values
{
	unsigned int a = 1;
	
	for (list<DiscreteRandomVariable *>::iterator it = l.begin(); it != l.end(); ++it)
	{
		a = a * (*it)->get_number_values();
		
		correspondance_table[ (*it)->get_index() ] = number_of_variable_values.size() ;
		
		if (number_of_variable_values.empty())
		{
			help_vector.push_back(1);
		}
		else
		{
			help_vector.push_back( help_vector.back() * number_of_variable_values.back() );
		}
		
		number_of_variable_values.push_back((*it)->get_number_values());
		variable_values.push_back(0);
		
	}
	
	sampling_probabilities.resize(a,  1.0 / a);
	prior_samples.resize(a);
	number_values = a;
	
}

unsigned int CombinedDiscreteRandomVariable::get_combined_position() const
{
	return inner_product(variable_values.begin(), variable_values.end(), help_vector.begin(), 0);
}

vector < DiscreteRandomVariable * > CombinedDiscreteRandomVariable::uncombine()
{
	vector < DiscreteRandomVariable * > result_vec;
	
	//list <DiscreteRandomVariable *> ert;
	//CombinedDiscreteRandomVariable * er = new CombinedDiscreteRandomVariable( ert , 0);
	
	//DiscreteRandomVariable * yt = new DiscreteRandomVariable( 5 , 5);
	
	//er->sampling_probabilities [0] = 2.0;
	//yt->sampling_probabilities [0] = 2.0;
	
	for ( map <unsigned int, unsigned int>::const_iterator it = correspondance_table.begin(); it != correspondance_table.end(); ++it)
	{
		DiscreteRandomVariable * new_variable = new DiscreteRandomVariable( it->first, number_of_variable_values[ it->second] );
		
		result_vec.push_back( new_variable );
		
		unsigned int a = 0;
		
		for ( vector <double>::iterator current_probability = (new_variable->sampling_probabilities).begin(); current_probability != new_variable->sampling_probabilities.end(); ++current_probability)
		{
			
			set_variable_value_by_index( it->first, a );
			
			*current_probability = internal_sum_over_variable(it->second);
			
			//cout << *current_probability << endl;
			
			++a;
		}
		
	}
	
	return result_vec;
}

double CombinedDiscreteRandomVariable::internal_sum_over_variable(unsigned int index)
{
	
	unsigned int left_block_size,right_block_size,right_blocks, fixed_value;
	
	
	// Right_blocks is the number of blocks we have to "the right", it is in fact the number of times we must sum in the outer loop
	// Right_block_size is the size of the blocks to the right... we add that once we are over in the inner loop
	
	if (index == number_of_variable_values.size()-1)
	{
		right_blocks = 1;
		right_block_size = 0; // Not used in this case
	}
	else
	{
		right_blocks = number_of_variable_values[number_of_variable_values.size()-1]*help_vector[number_of_variable_values.size()-1]/help_vector[index+1];
		right_block_size = help_vector[index+1];
	}
	
	// The fixed value represents how much we must add constantly to read the point where "the sum is frozen"
	// Left block size represents the size of the blocks to the left, ie, 1 if we are not summing over the first variable
	
	left_block_size = help_vector[index];
	fixed_value = variable_values[index]*left_block_size;
	
	//cout << "fixed value " << fixed_value << endl;
	
	double sum = 0;
	unsigned int a = fixed_value;
	
	for (unsigned int i = 0; i < right_blocks; ++i)
	{
		
		for (unsigned int j = 0; j < left_block_size; ++j)
		{
			sum += sampling_probabilities[a+j];
		}
		
		a += right_block_size;
		
	}
	
	return sum;
	
}

//****************** End of CombinedDiscreteRandomVariable Implementation ********************* //