/*
指数分布のシミュレーション
<Time-stamp: <2013/09/17 17:54:05 JST nakata> 

・λ>0について、X 〜 Exp(λ) <==> P(X≦x) = 1-exp(-λx), (x>0)
・U〜Unif([0,1])<==> P(U≦x)=x, (0≦x≦1)

定理: U〜Unif([0,1])⇒ - (log(1-U))/λ 〜 Exp(λ)

を用いてパラメータλの指数分布に従う乱数を発生する。
チェックのために平均1/λと分散1/λ^2 を表示する。
*/

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#define N 100000 /* 乱数の出力数 */

int main(void)
{
  int i;
  double ex1, lambda=2.0;
  double sum1=0.0, sum2=0.0;
  double ave, var; 
  /* 乱数の種をtime からとるおまじない */
  srand((unsigned)time(NULL));

  /* パラメータ λ>0 に従う指数分布に従う乱数列の出力 */
  printf("\nパラメータ %f>0 に従う指数分布に従う乱数列の出力\n",lambda);
  for(i = 0; i < N; i++){
    ex1 = -log(rand()/(1.0+RAND_MAX))/lambda; // U〜Unif([0,1])⇒「1-U 〜 U」
    sum1 = sum1 + ex1;
    sum2 = sum2 + ex1*ex1;
  }
  ave = sum1/N; var = sum2/N-ave*ave; // 2乗の平均 - 平均の2乗
  printf("Kurikaeshi: %d\nHeikin %f (riron %f)\nBunsan %f (riron %f)\n", 
	 N,  ave, 1.0/lambda, var, 1.0/(lambda*lambda));
  return(0);
}