lda.c

LDA程式設計

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_vector.h>
#include <assert.h>


#define ni 3
 
void gsl_matrix_mul_matrix(const gsl_matrix * , const gsl_matrix * , gsl_matrix * );

int main (void)
{
	int i, j;

	//double data[9];	

	double x[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
	double x1[] = { 31, 32, 28, 34, 31, 33, 31, 29, 29, 30, 31, 32 };
	double x2[] = { 33, 30, 33, 33, 31, 28, 32, 32, 29, 34, 35, 31 };
	double x3[] = { 31, 26, 31, 34, 30, 31, 29, 29, 34, 28, 23, 32 };
	double mean[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };

	gsl_vector * v = gsl_vector_alloc (12);
	gsl_vector * v1 = gsl_vector_alloc (12);
	gsl_vector * v2 = gsl_vector_alloc (12);
	gsl_vector * v3 = gsl_vector_alloc (12);
	gsl_vector * mean_ = gsl_vector_alloc (12);
       
       for (i = 0; i < 12; i++)
       {
		gsl_vector_set (v, i, x[i]);
		gsl_vector_set (v1, i, x1[i]);
		gsl_vector_set (v2, i, x2[i]);
		gsl_vector_set (v3, i, x3[i]);
		gsl_vector_set (mean_, i, mean[i]);
       }

/*************************累加與求平均值******************************/
	gsl_vector_add( v, v1 );
	gsl_vector_add( v, v2 );
	gsl_vector_add( v, v3 );
	gsl_vector_memcpy( mean_, v );
	gsl_vector_scale( mean_, 0.33333333 );
/***********************************************************************/

/*************************向量減去平均值******************************/
	gsl_vector_sub( v1, mean_ );
	gsl_vector_sub( v2, mean_ );
	gsl_vector_sub( v3, mean_ );
/***********************************************************************/

/************************** 求E[(x-E[x])(x-E[x])T] *********************/
	gsl_matrix * m = gsl_matrix_alloc ( 3, 12 );
	gsl_matrix_set_zero( m );
	gsl_matrix_set_row( m, 0, v1 );
	gsl_matrix_set_row( m, 1, v2 );
	gsl_matrix_set_row( m, 2, v3 );

	gsl_matrix * m_T = gsl_matrix_alloc ( 12, 3 );
	gsl_matrix * m_mul = gsl_matrix_alloc ( 3, 3 );
	gsl_matrix_transpose_memcpy ( m_T, m);
	gsl_matrix_mul_matrix( m, m_T, m_mul );
	gsl_matrix_scale( m_mul, 0.33333333 );
/***********************************************************************/

	
	for (i = 0; i < 3; i++)
	{
		for (j = 0; j <3; ++j)
		{
			printf("%f  ", gsl_matrix_get (m_mul, i, j));
		}
		printf("\n");
	}
     
	return 0;
}


/************************** 矩陣乘法 *********************/
void gsl_matrix_mul_matrix(const gsl_matrix * a, const gsl_matrix * b, gsl_matrix * c)
{
	assert(a->size2 == b->size1 && a->size1 == c->size1 && b->size2 == c->size2);

	/* set all elements in c to zero */
	gsl_matrix_set_zero(c);

	/* calculate mul */
	int i, j, k;
	for(i = 0; i < a->size1; i++)
	{
   		for(j = 0; j < b->size2; j++)
   		{
       			for(k = 0; k < a->size2; k++)
       			{
           			gsl_matrix_set(c, i, j, gsl_matrix_get(c, i, j) + gsl_matrix_get(a, i, k) * gsl_matrix_get(b, k, j));
       			}
   		}
	}
}