|
| 1 | +#include <math.h> |
| 2 | +#include <stdio.h> |
| 3 | +#include <stdlib.h> |
| 4 | +#include <malloc.h> |
| 5 | +#include <time.h> |
| 6 | +#include <io.h> |
| 7 | +void frequecy_table(int n,double *getvalue1) |
| 8 | +{ |
| 9 | + int i,ni,classno1,count; |
| 10 | + double len1,tmin1,tmax1,width1,*pdfvalue1,*upper1,*lower1,cuumlative_pr=0,class_midpoint; |
| 11 | + FILE *file; |
| 12 | + |
| 13 | + tmin1=getvalue1[0]; tmax1=getvalue1[0]; |
| 14 | + for(ni=1;ni<n;++ni) |
| 15 | + { if (getvalue1[ni]<tmin1) tmin1=getvalue1[ni]; |
| 16 | + if (getvalue1[ni]>tmax1) tmax1=getvalue1[ni]; |
| 17 | + } |
| 18 | + len1=tmax1-tmin1; |
| 19 | + classno1=200; |
| 20 | + width1=len1/(double)classno1; |
| 21 | + pdfvalue1=(double *)malloc(classno1*sizeof(double)); |
| 22 | + upper1=(double *)malloc(classno1*sizeof(double)); |
| 23 | + lower1=(double *)malloc(classno1*sizeof(double)); |
| 24 | + for(ni=0;ni<classno1;++ni) pdfvalue1[ni]=0; |
| 25 | + upper1[0]=tmin1+width1; |
| 26 | + for(ni=1;ni<classno1;++ni) |
| 27 | + { lower1[ni]=upper1[ni-1]; upper1[ni]=lower1[ni]+width1; } |
| 28 | + upper1[classno1-1]=tmax1+1; |
| 29 | + lower1[0]=tmin1-1; |
| 30 | + for(ni=0;ni<n;++ni) |
| 31 | + { |
| 32 | + count=(long)((getvalue1[ni]-tmin1)/width1); |
| 33 | + if (count>=0 && count<classno1) pdfvalue1[count]+=1; |
| 34 | + } |
| 35 | + for(ni=0;ni<classno1;++ni) pdfvalue1[ni]/=(double)(n); |
| 36 | + lower1[0]=tmin1;upper1[classno1-1]=tmax1; |
| 37 | + file=fopen("C:\\C_Bernoulli\\tep\\frequency_table.txt","w"); |
| 38 | + cuumlative_pr=0; |
| 39 | + for(ni=0;ni<classno1;++ni) |
| 40 | + { |
| 41 | + class_midpoint=(lower1[ni]+upper1[ni])/2; |
| 42 | + cuumlative_pr=cuumlative_pr+pdfvalue1[ni]; |
| 43 | + fprintf(file,"%20.10f %20.10f %20.10f %20.10f\n",class_midpoint,pdfvalue1[ni]/width1,class_midpoint,cuumlative_pr); |
| 44 | + } |
| 45 | + fclose(file); |
| 46 | + free(pdfvalue1); free(upper1); free(lower1); |
| 47 | +} |
| 48 | +// the random number |
| 49 | +double simulated_uniform(double alpha,double beta) |
| 50 | +{ |
| 51 | + double tep1,tep=0,dividewangx; |
| 52 | + |
| 53 | + dividewangx=(double)32767*(double)32767*(double)32767; |
| 54 | + again: |
| 55 | + tep1=((double)rand()*(double)32767*(double)32767+(double)rand()*(double)32767+(double)rand())/dividewangx; |
| 56 | + if (tep1<0 || tep1>1) goto again; |
| 57 | + tep=alpha+tep1*(beta-alpha); |
| 58 | + return(tep); |
| 59 | +} |
| 60 | +// The comparison of two strings and the index upper limit is len |
| 61 | +int strcompare_spc(char *s,char *ss,int len) |
| 62 | +{ |
| 63 | + int i=0; |
| 64 | + for(;s[i]!=0 && ss[i]!=0 && i<len;++i) |
| 65 | + if (s[i]!=ss[i]) break; |
| 66 | + return(s[i]-ss[i]); |
| 67 | +// return 0 when two strings are eqully |
| 68 | +} |
| 69 | +// Continuous Bernoulli distribution simulator |
| 70 | +double simulated_continue_bernoulli(double lamda) |
| 71 | +{ |
| 72 | + double tep=0,tep1; |
| 73 | + char ss[1024]; |
| 74 | + int kk=1; |
| 75 | + |
| 76 | + // checking the lamda is 0.5 or not |
| 77 | + if (lamda>=0.499 && lamda<=0.501) |
| 78 | + { |
| 79 | + sprintf(ss,"%12.10f",lamda); |
| 80 | + kk=strcompare_spc(ss,(char *)"0.5000000000",12); |
| 81 | + } |
| 82 | + // kk is 0 when lamda=0.5, X~Uniform(0,1) if lamda=0.5 |
| 83 | + if (kk!=0) |
| 84 | + { |
| 85 | + // the inverse distribution function is the simulated data |
| 86 | + tep1=simulated_uniform(0,1); |
| 87 | + tep1=(2*lamda-1)*tep1-lamda+1; |
| 88 | + tep1=tep1/(1-lamda); |
| 89 | + tep=log(tep1)/log(lamda/(1-lamda)); |
| 90 | + } |
| 91 | + else |
| 92 | + tep=simulated_uniform(0,1); |
| 93 | + return(tep); |
| 94 | +} |
| 95 | +// computing the mean,variance,skewed coefficient and kurtosis coefficient |
| 96 | +void moment(int n,double *datastore) |
| 97 | +{ |
| 98 | + int i; |
| 99 | + double EX,VarX,skew,kurtosis; |
| 100 | + double sumX,sumX2,sumX3,sumX4,sigmaX,tep; |
| 101 | + |
| 102 | + sumX=0; sumX2=0; |
| 103 | + for(i=0;i<n;++i) |
| 104 | + sumX=sumX+datastore[i]; |
| 105 | + sumX=sumX/(double)n; |
| 106 | + for(i=0;i<n;++i) |
| 107 | + sumX2=sumX2+(datastore[i]-sumX)*(datastore[i]-sumX); |
| 108 | + sumX2=sumX2/(double)n; |
| 109 | + if (sumX2>0) |
| 110 | + { |
| 111 | + sigmaX=exp(0.5*log(sumX2)); |
| 112 | + sumX3=0; |
| 113 | + sumX4=0; |
| 114 | + for(i=0;i<n;++i) |
| 115 | + { |
| 116 | + tep=(datastore[i]-sumX)/sigmaX; |
| 117 | + sumX3=sumX3+tep*tep*tep; |
| 118 | + sumX4=sumX4+tep*tep*tep*tep; |
| 119 | + } |
| 120 | + sumX3=sumX3/(double)n; |
| 121 | + sumX4=sumX4/(double)n; |
| 122 | + EX=sumX;VarX=sumX2;skew=sumX3;kurtosis=sumX4; |
| 123 | + printf(" E(X)=%20.10f,\n Var(X)=%20.10f,\n",EX,VarX); |
| 124 | + printf(" skewed coefficient=%20.10f,\n kurtosis coefficient=%20.10f\n",skew,kurtosis); |
| 125 | + } |
| 126 | + else |
| 127 | + printf(" E(X)=%20.10f\n",EX); |
| 128 | +} |
| 129 | +int main() |
| 130 | +{ |
| 131 | + int i,n=100000000; |
| 132 | + double *datastore,lamda; |
| 133 | + char *ss; |
| 134 | + |
| 135 | + ss=(char *)malloc(1024*sizeof(char)); |
| 136 | + // input the parameter lamda |
| 137 | + loop1: |
| 138 | + printf(" The Continuous Bernoulli, lamda="); |
| 139 | + scanf("%s",ss); |
| 140 | + lamda=atof(ss); |
| 141 | + if (lamda<=0 || lamda>=1 ) |
| 142 | + { |
| 143 | + printf(" The Continuous Bernoulli(lamda=%f), it is error, 0<lamda<1,\n please input again.\n",lamda); |
| 144 | + goto loop1; |
| 145 | + } |
| 146 | + datastore=(double *)malloc(n*sizeof(double)); |
| 147 | + printf(" please wait a moment to collect the simulated data,....\n"); |
| 148 | + // getting the simulated data and the amount=100,000,000, |
| 149 | + // the moment of this database is closed to the distribution |
| 150 | + for(i=0;i<n;++i) |
| 151 | + { |
| 152 | + if ((i+1)%10000000==0) printf(" finished %10.8f%%\n",(double)(i+1)/(double)n*100); |
| 153 | + datastore[i]=simulated_continue_bernoulli(lamda); |
| 154 | + } |
| 155 | + printf(" X~Continuous Bernoulli(lamda=%f),\n",lamda); |
| 156 | + moment(n,datastore); |
| 157 | + printf(" the frequency table is c:\\C_Bernoulli\\tep\\frequency_table.txt\n"); |
| 158 | + i=chdir("c:\\C_Bernoulli"); // check this subdirectory is existed |
| 159 | + if (i!=0) mkdir("c:\\C_Bernoulli"); // make this subdirectory |
| 160 | + i=chdir("c:\\C_Bernoulli\\tep"); // check this subdirectory is existed |
| 161 | + if (i!=0) mkdir("c:\\C_Bernoulli\\tep"); // make this subdirectory |
| 162 | + frequecy_table(n,datastore); |
| 163 | + free(ss); free(datastore); |
| 164 | + system("pause"); |
| 165 | + return 0; |
| 166 | +} |
| 167 | + |
| 168 | + |
0 commit comments