#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define ITER pow(10, 4)
double* rnorm(double mu, double sigma, int n) {
double *samples = (double *)malloc(n * sizeof(double));
for (int i = 0; i < n; i++) {
double U1, U2, W, mult;
static double X1, X2;
static int call = 0;
if (call == 1) {
call = !call;
samples[i] = mu + sigma * X2;
continue;
}
do {
U1 = -1 + ((double) rand () / RAND_MAX) * 2;
U2 = -1 + ((double) rand () / RAND_MAX) * 2;
W = pow (U1, 2) + pow (U2, 2);
} while (W >= 1 || W == 0);
mult = sqrt ((-2 * log (W)) / W);
X1 = U1 * mult;
X2 = U2 * mult;
call = !call;
samples[i] = mu + sigma * X1;
}
return samples;
}
double var(double *arr, int size) {
double mean_val = 0;
double sum_squares = 0;
double var_sum = 0;
// Calculate the mean and sum of squares
for (int i = 0; i < size; i++) {
mean_val += arr[i];
sum_squares += arr[i] * arr[i];
}
mean_val /= size;
// Calculate the variance
for (int i = 0; i < size; i++) {
var_sum += (arr[i] - mean_val) * (arr[i] - mean_val);
}
return (sum_squares - (size * mean_val * mean_val)) / (size - 1);
}
int main() {
srand(time(NULL));
int iter = ITER;
int n = 5;
double mu0 = 0.0, sig0 = 1.0;
double an = -0.8969, bn = 2.3647, cn = 0.5979, muT = 0.00748, sigT = 0.9670;
double b = bn;
double c = cn * pow(sig0, 2);
double a = an - 2 * bn * log(sig0);
double L = 3.262, k = 0.45;
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};
int gm_length = sizeof(gm) / sizeof(gm[0]);
double *ARL = (double *)malloc(gm_length * sizeof(double));
double *rl = (double *)malloc(iter * sizeof(double));
double *SCn = (double *)malloc(iter * sizeof(double));
double *SCp = (double *)malloc(iter * sizeof(double));
//double LL;
//double UL;
double *Cp = (double *)malloc(iter * sizeof(double));
double *Cn = (double *)malloc(iter * sizeof(double));
double *dp = (double *)malloc(iter * sizeof(double));
double *dn = (double *)malloc(iter * sizeof(double));
printf("Gamma\tARL\n");
printf("--------------\n");
for (int d = 0; d < gm_length; d++) {
double sig = gm[d] / sig0;
double sum_rl = 0;
for (int j = 0; j < iter; j++) {
for (int i = 0; i < iter; i++) {
double *x = rnorm(mu0, sig, n);
double s_sq = var(x, n);
double t = a + b * log(s_sq + c);
double z = (t - muT) / sigT;
if (i == 0) {
Cp[i] = fmax(0, z - k);
Cn[i] = fmax(0, -z - k);
} else {
Cp[i] = fmax(0, z - k + Cp[i - 1]);
Cn[i] = fmax(0, -z - k + Cn[i - 1]);
}
if (i == 0) {
dp[i] = 1;
dn[i] = 1;
} else {
dp[i] = (Cp[i - 1] == 0) ? 1 : dp[i - 1] + 1;
dn[i] = (Cn[i - 1] == 0) ? 1 : dn[i - 1] + 1;
}
if (i == 0) {
SCn[i] = (1 / dn[i]) * z + (1 - (1 / dn[i]));
SCp[i] = (1 / dp[i]) * z + (1 - (1 / dp[i]));
} else {
SCn[i] = (1 / dn[i]) * z + (1 - (1 / dn[i])) * SCn[i - 1];
SCp[i] = (1 / dp[i]) * z + (1 - (1 / dp[i])) * SCp[i - 1];
}
double LL = -L / sqrt(dn[i]);
double UL = L / sqrt(dp[i]);
if (SCn[i] < LL || SCp[i] > UL) {
rl[j] = i + 1;
break;
} else {
rl[j] = 0;
}
free(x);
}
}
for (int i = 0; i < iter; i++) {
sum_rl += rl[i];
}
ARL[d] = sum_rl / iter;
}
for (int i = 0; i < gm_length; i++) {
printf("%.2f\t%.2f\n", gm[i], ARL[i]);
}
free(ARL);
free(rl);
free(SCn);
free(SCp);
free(Cp);
free(Cn);
free(dp);
free(dn);
return 0;
}