generate_phi.c

一个通用的隐性马尔可夫C代码库 开发环境:C语言 简要说明:这是一个通用的隐性马尔可夫C代码库

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_histogram.h>
#include <ghmm/rng.h>
#include <ghmm/randvar.h>
#include <stdio.h>
#include <stdlib.h>


#ifndef GSL_HISTOGRAM_SET_RANGES_UNIFORM
int gsl_histogram_set_ranges_uniform(gsl_histogram* h, double xmin, double xmax )
{
  size_t i;

  if (xmin >= xmax)
    {
      GSL_ERROR_VAL ("xmin must be less than xmax", GSL_EINVAL, 0);
    }

  for (i = 0; i < h->n + 1; i++)
    {
      h->range[i] = xmin + ((double) i / (double) h->n) * (xmax - xmin);
    }
  return h->n;
}
#endif /* HAVE_GSL_HISTOGRAM_SET_RANGES_UNIFORM */

/* bisection algorithms to explore precision*/
int greatest_thats_different_from_1(void)
{
  double low,up,half;
  low=0;
  up=10;
  /* find start interval */
  while (randvar_get_PHI(up)<1.0)
    {
      low=up;
      up*=2;
    }
  /* shrink it */
  while (up-low>0.001)
    {
      half=(low+up)/2.0;
      if (randvar_get_PHI(half)<1.0)
	low=half;
      else
	up=half;
    }
  printf("%f: greatest int, for which PHI is different from 1.0\n",low);
  return 0;
}

int least_thats_different_from_0(void)
{
  double low,up,half;
  low=-10;
  up=0;
  /* find start interval */
  while (randvar_get_PHI(low)>0.0)
    {
      up=low;
      low*=2;
    }
  /* shrink it */
  while (up-low>0.001)
    {
      half=(low+up)/2.0;
      if (randvar_get_PHI(half)>0.0)
	up=half;
      else
	low=half;
    }
  printf("%f: least int, for which PHI is different from 0.0\n",low);
  return 0;
}

int print_PHI_table(void)
{
  long i;
  double a;
  gsl_sf_result r;
  i=0;
  for(a=0.0;a>-20.0;a-=0.001)
    {
      /* with gsl special function */
      (void)gsl_sf_erf_e(a,&r);
      /* with randvar function */
      printf("%f,%e,%e\n",a,r.val/2.0+0.5,r.err/2.0);
      i++;
    }
  printf("Erzeugte %ld Werte\n",i);
  return 0;
}

/* create a gnuplot file with pictures from randomnumbers*/
void gnuplot_show_distributions(void)
{
  /* where output goes */
  FILE* file;
  const int generate_count=100000;
  /* count how often generator was used*/
  int generate_counter;
  /* row of bins */
  gsl_histogram* bins;

  /* now to stdout*/
  file=stdout;
  /* init bin_row */
  bins=gsl_histogram_calloc_uniform(100,-10.0,10.0);
  /* initialise generator */
  gsl_rng_init();

  /* prepare output for gnuplot */
  fprintf(file,"set terminal postscript eps\n");

  /* generate a normal distribution picture */
  fprintf(file,"set output 'std_normal_dist.eps'; plot '-' using 1:3 title 'standard normal distribution'\n");
  for (generate_counter=0; generate_counter<generate_count; generate_counter++)
    gsl_histogram_increment(bins,randvar_std_normal(0));
  (void)gsl_histogram_fprintf (file, bins, "%f","%f");
  fprintf(file,"e\n");

  /* generate  three truncated normal distribution in one picture */
  gsl_histogram_set_ranges_uniform(bins,-2,10);
  fprintf(file,"set output 'trunc_normal_dist.eps';\n\
plot '-' using 1:3 title 'mue=-4','-' using 1:3 title 'mue=0', '-' using 1:3 title 'mue=4'\n");

  /* generate histogram with mue=-4 */
  gsl_histogram_reset(bins);
  for (generate_counter=0; generate_counter<generate_count; generate_counter++)
    gsl_histogram_increment(bins,randvar_normal_pos(-4,1,0));
  (void)gsl_histogram_fprintf (file, bins, "%f","%f");
  fprintf(file,"e\n");

  /* generate histogram with mue=0 */
  gsl_histogram_reset(bins);
  for (generate_counter=0; generate_counter<generate_count; generate_counter++)
    gsl_histogram_increment(bins,randvar_normal_pos(0,1,0));
  (void)gsl_histogram_fprintf (file, bins, "%f","%f");
  fprintf(file,"e\n");

  /* generate histogram with mue=4 */
  gsl_histogram_reset(bins);
  for (generate_counter=0; generate_counter<generate_count; generate_counter++)
    gsl_histogram_increment(bins,randvar_normal_pos(4,1,0));
  (void)gsl_histogram_fprintf (file, bins, "%f","%f");
  fprintf(file,"e\n");

  /* terminate gnuplot */
  fprintf(file,"set output; quit\n");

  /* cleanup */
  gsl_histogram_free(bins);
}

void gnuplot_show_distributions_pdf(void)
{
  /* where output goes */
  FILE* file;
  /* ranges */
  double xmin=-2;
  double xmax=10;
  double now;
  double samples;
  double step;

  file=stdout;
  samples=100;
  step=(xmax-xmin)/samples;

  /* prepare output for gnuplot */
  fprintf(file,"set terminal postscript eps\n");

  fprintf(file,"set output 'trunc_normal_pdf.eps';\n\
plot '-' title 'mue=-4','-' title 'mue=0', '-' title 'mue=4'\n");
  
  for (now=xmin; now<=xmax;now+=step)
    {
      fprintf(file,"%f\t%f\n",now,randvar_normal_density_pos(now,-4.0,1.0));
    }
  fprintf(file,"e\n");
  for (now=xmin; now<=xmax;now+=step)
    {
      fprintf(file,"%f\t%f\n",now,randvar_normal_density_pos(now,0,1.0));
    }
  fprintf(file,"e\n");
  for (now=xmin; now<=xmax;now+=step)
    {
      fprintf(file,"%f\t%f\n",now,randvar_normal_density_pos(now,4.0,1.0));
    }
  fprintf(file,"e\n");

  /* terminate gnuplot */
  fprintf(file,"set output; quit\n");
}

int main()
{
#if 0
  (void)greatest_thats_different_from_1();
  (void)least_thats_different_from_0();
#endif

  gnuplot_show_distributions_pdf();
  return 0;
}