fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include <time.h>
  5.  
  6.  
  7. #define ITER pow(10, 4)
  8.  
  9. double* rnorm(double mu, double sigma, int n) {
  10.     double *samples = (double *)malloc(n * sizeof(double));
  11.     for (int i = 0; i < n; i++) {
  12.         double U1, U2, W, mult;
  13.         static double X1, X2;
  14.         static int call = 0;
  15.  
  16.         if (call == 1) {
  17.             call = !call;
  18.             samples[i] = mu + sigma * X2;
  19.             continue;
  20.         }
  21.  
  22.         do {
  23.             U1 = -1 + ((double) rand () / RAND_MAX) * 2;
  24.             U2 = -1 + ((double) rand () / RAND_MAX) * 2;
  25.             W = pow (U1, 2) + pow (U2, 2);
  26.         } while (W >= 1 || W == 0);
  27.  
  28.         mult = sqrt ((-2 * log (W)) / W);
  29.         X1 = U1 * mult;
  30.         X2 = U2 * mult;
  31.  
  32.         call = !call;
  33.  
  34.         samples[i] = mu + sigma * X1;
  35.     }
  36.     return samples;
  37. }
  38.  
  39.  
  40.  
  41. double var(double *arr, int size) {
  42.     double mean_val = 0;
  43.     double sum_squares = 0;
  44.     double var_sum = 0;
  45.  
  46.     // Calculate the mean and sum of squares
  47.     for (int i = 0; i < size; i++) {
  48.         mean_val += arr[i];
  49.         sum_squares += arr[i] * arr[i];
  50.     }
  51.     mean_val /= size;
  52.  
  53.     // Calculate the variance
  54.     for (int i = 0; i < size; i++) {
  55.         var_sum += (arr[i] - mean_val) * (arr[i] - mean_val);
  56.     }
  57.  
  58.     return (sum_squares - (size * mean_val * mean_val)) / (size - 1);
  59. }
  60.  
  61. int main() {
  62.     srand(time(NULL));
  63.     int iter = ITER;
  64.  
  65.     int n = 5;
  66.     double mu0 = 0.0, sig0 = 1.0;
  67.     double an = -0.8969, bn = 2.3647, cn =  0.5979, muT =  0.00748, sigT =  0.9670;
  68.  
  69.     double b = bn;
  70.     double c = cn * pow(sig0, 2);
  71.     double a = an - 2 * bn * log(sig0);
  72.     double L = 3.262, k = 0.45;
  73.  
  74.     double gm[] = {1};//0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 1, 1.05, 1.1, 1.2, 1.3, 1.4, 1.5, 2, 3};
  75.     int gm_length = sizeof(gm) / sizeof(gm[0]);
  76.  
  77.     double *ARL = (double *)malloc(gm_length * sizeof(double));
  78.     double *rl = (double *)malloc(iter * sizeof(double));
  79.     double *SCn = (double *)malloc(iter * sizeof(double));
  80.     double *SCp = (double *)malloc(iter * sizeof(double));
  81.     //double LL;
  82.     //double UL;
  83.     double *Cp = (double *)malloc(iter * sizeof(double));
  84.     double *Cn = (double *)malloc(iter * sizeof(double));
  85.     double *dp = (double *)malloc(iter * sizeof(double));
  86.     double *dn = (double *)malloc(iter * sizeof(double));
  87.  
  88.  
  89.     printf("Gamma\tARL\n");
  90.     printf("--------------\n");
  91.  
  92.     for (int d = 0; d < gm_length; d++) {
  93.         double sig = gm[d] / sig0;
  94.         double sum_rl = 0;
  95.  
  96.         for (int j = 0; j < iter; j++) {
  97.  
  98.             for (int i = 0; i < iter; i++) {
  99.                 double *x = rnorm(mu0, sig, n);
  100.                 double s_sq = var(x, n);
  101.                 double t = a + b * log(s_sq + c);
  102.                 double z = (t - muT) / sigT;
  103.  
  104.                 if (i == 0) {
  105.                     Cp[i] = fmax(0, z - k);
  106.                     Cn[i] = fmax(0, -z - k);
  107.                 } else {
  108.                     Cp[i] = fmax(0, z - k + Cp[i - 1]);
  109.                     Cn[i] = fmax(0, -z - k + Cn[i - 1]);
  110.                 }
  111.  
  112.                 if (i == 0) {
  113.                     dp[i] = 1;
  114.                     dn[i] = 1;
  115.                 } else {
  116.                     dp[i] = (Cp[i - 1] == 0) ? 1 : dp[i - 1] + 1;
  117.                     dn[i] = (Cn[i - 1] == 0) ? 1 : dn[i - 1] + 1;
  118.                 }
  119.  
  120.                 if (i == 0) {
  121.                     SCn[i] = (1 / dn[i]) * z + (1 - (1 / dn[i]));
  122.                     SCp[i] = (1 / dp[i]) * z + (1 - (1 / dp[i]));
  123.                 } else {
  124.                     SCn[i] = (1 / dn[i]) * z + (1 - (1 / dn[i])) * SCn[i - 1];
  125.                     SCp[i] = (1 / dp[i]) * z + (1 - (1 / dp[i])) * SCp[i - 1];
  126.                 }
  127.  
  128.                 double LL = -L / sqrt(dn[i]);
  129.                 double UL = L / sqrt(dp[i]);
  130.  
  131.                 if (SCn[i] < LL || SCp[i] > UL) {
  132.                     rl[j] = i + 1;
  133.                     break;
  134.                 } else {
  135.                     rl[j] = 0;
  136.                 }
  137.  
  138.                 free(x);
  139.             }
  140.  
  141.  
  142.         }
  143.  
  144.         for (int i = 0; i < iter; i++) {
  145.             sum_rl += rl[i];
  146.         }
  147.  
  148.         ARL[d] = sum_rl / iter;
  149.     }
  150.  
  151.  
  152.     for (int i = 0; i < gm_length; i++) {
  153.         printf("%.2f\t%.2f\n", gm[i], ARL[i]);
  154.     }
  155.  
  156.     free(ARL);
  157.     free(rl);
  158.     free(SCn);
  159.     free(SCp);
  160.     free(Cp);
  161.     free(Cn);
  162.     free(dp);
  163.     free(dn);
  164.  
  165.  
  166.     return 0;
  167. }
  168.  
Success #stdin #stdout 1.08s 5260KB
comments (?)
stdin
copy
stdout
copy 
Gamma	ARL
--------------
1.00	334.44
https://ideone.com/mKL2S6
language:
C++14 (gcc 8.3)
created:
12 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!