一维GMM算法C语言实现

    在模式识别课程中，学习了GMM算法。并实现1维GMM算法。

     代码如下：

    

/*Author : Li PanTime : 2012/9/20Version : 1.0Using this file to test GMM algorithm.Input : input the the gaussian distribution number, the real u1, u2, ...,        the real xigama(1), sigma(2), sigma(3), and input the expected gaussian        distribution number... */ #include <stdio.h> #include <math.h> #include <stdlib.h> #include <assert.h> #define SAMPLE_NUMBER 2000 #define MAX_K 10 #define EPS 0.000001 #define PI 3.141593 typedef struct Element {float numerator;float denominator;float result; }Element;  Element elements[SAMPLE_NUMBER][MAX_K];  int gaussian_distribution_number;  float *means; float *variances; float *probabilities;  float samples[SAMPLE_NUMBER]; void generate_samples(); void free_mvp(); void calculate_arguments(int K); float _fill_in_elements_and_return_q(int K); float _sum(int total_sample,int col_number); float _sum_with_vector(int total_sample,int col_number,float * multiply_vector); float _gaussian_distribution(float x,float mean,float variance); void print_mvp(int expected_gaussian_distribution_number);  int main() {int i;int expected_gaussian_distribution_number;fprintf(stdout, "Please enter the gaussian distribution number : \n");fscanf(stdin, "%d", &gaussian_distribution_number);means = (float *)malloc(gaussian_distribution_number * sizeof(float));assert(NULL != means);variances = (float *)malloc(gaussian_distribution_number * sizeof(float));assert(NULL != variances);probabilities = (float *)malloc(gaussian_distribution_number * sizeof(float));assert(NULL != probabilities);fprintf(stdout, "Please enter the means,variances and probability for each gaussian distribution with format [mean variance probability]\n");for (i = 0; i < gaussian_distribution_number; ++i) {fscanf(stdin, "%f %f %f", &means[i], &variances[i], &probabilities[i]);}fprintf(stdout, "Enter the expected gaussian distribution number :\n");fscanf(stdin, "%d", &expected_gaussian_distribution_number);generate_samples();calculate_arguments(expected_gaussian_distribution_number);   /*Fill in means, variances and probabilities*/print_mvp(expected_gaussian_distribution_number);free_mvp();return 0; }void generate_samples() {int i;int j;int selected_gaussian_distribution;int n = 20; /*When building gaussian distribution, use n.*/float x = 0;float pro;/*Here we change the elements in probabilities*/for (i = 1; i < gaussian_distribution_number; ++i) {probabilities[i] += probabilities[i - 1];}srand(time(0));for (i = 0; i < SAMPLE_NUMBER; ++i) {pro = (float)(rand() % 1001) * 0.001f;for (j = 0; j < gaussian_distribution_number; ++j) {if (pro <= probabilities[j]) {selected_gaussian_distribution = j;break;}}x = 0.0;for (j = 0; j < n; ++j) {x+= (float)((rand() % 1001) * 0.001f);}x = (x - (float)n / 2) / sqrt((float)n / 12); /*This formula is used to build gaussian distribution!*/samples[i] = sqrt(variances[selected_gaussian_distribution]) * x + means[selected_gaussian_distribution];}free_mvp();}void free_mvp() {free(means);free(variances);free(probabilities);}/*It is the core part of this program!*/void calculate_arguments(int K) {int i;int t;float temp[SAMPLE_NUMBER];float f;float previous_q;float current_q = 0.0;/*_allocate_and_initialize_mvp(K);*/    _allocate_and_initialize_mvp(K);    previous_q = _fill_in_elements_and_return_q(K);for (; ;) {/*Calculating probabilities*/for (i = 0; i < K; ++i) {probabilities[i] = 1.0f * _sum(SAMPLE_NUMBER, i)/ SAMPLE_NUMBER;}/*Calculating means*/for (i = 0; i < K; ++i) {means[i] = _sum_with_vector(SAMPLE_NUMBER, i, samples) / (_sum(SAMPLE_NUMBER, i));}/*Calculating variances*/for (i = 0; i < K; ++i) {for (t = 0; t < SAMPLE_NUMBER; ++t) {temp[t] = (samples[t] - means[i]) * (samples[t] - means[i]);}variances[i] = _sum_with_vector(SAMPLE_NUMBER,i,temp)/ (_sum(SAMPLE_NUMBER, i));}current_q = _fill_in_elements_and_return_q(K);if (fabs(current_q - previous_q) < EPS) {fprintf(stdout, "previous_q : %f\n", previous_q);fprintf(stdout, "current_q : %f\n", current_q);break;} else {previous_q = current_q;}}}float _sum(int total_sample, int col_number) {float result = 0.0f;int i;for (i = 0; i < total_sample; ++i) {result += elements[i][col_number].result;}return result;}float _sum_with_vector(int total_sample, int col_number, float *multiply_vector) {float result = 0.0f;int i;for (i = 0; i < total_sample; ++i) {result += elements[i][col_number].result * multiply_vector[i];}return result;}void _allocate_and_initialize_mvp(int K) {int i;float total = 0.0f;float mean;float variance;float step1;float step2;means = (float *)malloc(K * sizeof(float));variances = (float *)malloc(K * sizeof(float));probabilities = (float*)malloc(K * sizeof(float));for (i = 0; i < SAMPLE_NUMBER; ++i) {total += samples[i];}mean = total / SAMPLE_NUMBER;step1 = mean / K;for (i = 0; i < SAMPLE_NUMBER; ++i) {total += pow(samples[i] - mean, 2);}variance = total / SAMPLE_NUMBER;step2 = variance / K;for (i = 0; i < K; ++i) {//means[i] = (float)(rand() % ((int)ceil(mean))) + (float)(rand() % ((int)ceil(mean)));means[i] = mean + pow(-1, i) * step1 * i;fprintf(stdout, "%d : %f\n", i, means[i]);variances[i] = variance + pow(-1, i) * step2 * i;//1.0f;probabilities[i] = 1.0f / K;}}float _fill_in_elements_and_return_q(int K) {int i;int t;float f;float current_q = 0.0f;for (t = 0; t < SAMPLE_NUMBER; ++t) {f = 0.0f;for (i = 0; i < K; ++i) {f+= probabilities[i] * _gaussian_distribution(samples[t],means[i], variances[i]);}        /*f may be 0!    Return value from _gaussian_distribution may be 0! *///fprintf(stdout, "_fill_in... f : %f\n", f);for (i = 0; i < K; ++i) {elements[t][i].numerator = probabilities[i] * _gaussian_distribution(samples[t],means[i], variances[i]);elements[t][i].denominator = f;elements[t][i].result = elements[t][i].numerator / elements[t][i].denominator;}}current_q = 0.0;for (t = 0; t < SAMPLE_NUMBER; ++t) {current_q += log(elements[t][0].denominator);}return current_q;}float _gaussian_distribution(float x, float mean, float variance) {float result = 0.0f;result = 1.0f / (sqrt(2 * PI) * sqrt(variance)) * exp((x - mean) * (x - mean) / (-2 * variance)); return result;}void print_mvp(int expected_gaussian_distribution_number) {int i;fprintf(stdout, "\nidprobabilitymeanvariance");for (i = 0; i < expected_gaussian_distribution_number; ++i) {fprintf(stdout, "\n%d%f%f%f\n", (i + 1), probabilities[i], means[i], variances[i]);}}/*  Finished time: 2012/9/21 12:33   */

     该算法对某些样本可以很好滴估计。