//-----------------------------------------------------------------
// list_66.c:
//     ２項分布のエミュレート
//-----------------------------------------------------------------
#include <stdio.h>
#include <math.h>

#define N 35
#define P 0.2
//-----------------------------------------------------------------
double binomial_p(int n,double p,int bk);
unsigned long long calc_casen(int n,int r);
//-----------------------------------------------------------------

int main()
{
    int i,pi;
    //--------------------------------------------------
    // 二項分布のエミュレート(確率の変化)：
    // 生起確率が大きくなると分布が回数大の方へ移動
    //-------------------------------------------------
    printf("**確率の変化\n");
    printf("生起回数\t生起確率\n");
    printf("\t");
    for(pi=1;pi<=5;pi++) printf("%.1lf\t",P*pi);
    printf("\n");
    for(i=0;i<=N;i++){
        printf("%d:\t",i);
        for(pi=1;pi<=5;pi++){
            printf("%.1lf\t",binomial_p(N,P*pi,i)*100);
        }
        printf("\n");
    }
    //--------------------------------------------------
    // 二項分布のエミュレート(総事象数の変化)：
    // 総事象数が増すと、分布が左右対照に近くなる
    //--------------------------------------------------
    printf("\n**事象数の変化\n");
    printf("生起回数\tその確率\n");
    printf("\t");
    for(pi=1;pi<=5;pi++) printf("%d\t",5*pi);
    printf("\n");
    for(i=1;i<=N;i++){
        printf("%d:\t",i);
        for(pi=1;pi<=5;pi++){
            printf("%.2lf\t",binomial_p(5*pi,P,i)*100);
        }
        printf("\n");
    }
    return 0;
}
//-----------------------------------------------------------------
// ２項分布(B(n,p))の確率計算
//-----------------------------------------------------------------
// 確率pで毎回独立に起こる事象がn回のうちk回起こる確率:
//   P{k}=nCk*p^k*(1-p)^(n-k)
//-----------------------------------------------------------------
double binomial_p(int n,double p,int bk)
{
    double q=(1-p),bp;
    bp = calc_casen(n,bk) * pow(p,bk) * pow(q,n-bk);
    return bp;
}

//-----------------------------------------------------------------
// 組合せ数 nCr の計算[Pascalの3角形による方法]
//-----------------------------------------------------------------
unsigned long long calc_casen(int n,int r)
{
    int i,j,nn=n;
    unsigned long long* a=
        (unsigned long long*)calloc(n,sizeof(unsigned long long));
    unsigned long long retv;
    if(nn-r<r) r=nn-r;
    if(r==0) return 1; if(r==1) return nn;
    for(i=1;i<r;i++) a[i]=i+2;
    for(i=3;i<=nn-r+1;i++){
        a[0]=i;
        for(j=1;j<r;j++) a[j]+=a[j-1];
    }
    retv = a[r-1];
    free(a);
    return retv;
}