marple.c

常用算法与数据结构原代码

#include "stdlib.h"
#include "stdio.h"

// 函数说明: 生成符合均匀分布规律的随机数
// 输入参数: double a, 所定区间的下限
//           double b, 所定区间的上限
//           long  *s, 产生随机数的种子
// 返回值  : 生成的随机数
  
double uniform(double a, double b, long *s)
{
	double t;
	
	*s = 2045 * (*s) +1;
	*s = *s - (*s / 1048576) * 1048576;
	t  = (*s) / 1048576.0;
	t  = a + (b - a) * t;
	
	return t;
}

// 函数说明: 生成符合正态分布规律的随机数
// 输入参数: double mean,	均值
//           double sigma,	均方差
//           long  *s,      产生随机数的种子
// 返回值  : 生成的随机数

double gauss(double mean, double sigma, long *s)
{
	int		i;
	double	x, y;
	
	for (x=0, i=0; i<12; i++)
	{
		x += uniform(0.0,1.0, s);
	}
	x -= 6.0;
	y =  mean + x * sigma;

	return y;
}

// 函数说明: 生成符合ARMA(p,q)模型的一组随机数
// 输入参数: double a[],	自回归参数
//           double b[],	滑动平均参数
//           int    p,      自回归部分的阶数
//           int    q,      滑动平均部分的阶数
//           double mean,   均值
//           double sigma,  均方差
//           long   *seed,  产生随机数的种子
//           double x[],    存放产生的随机数
//           int    n,      需要产生的随机数个数
// 返回值  : 无

void arma(double a[], double b[], int p, int q, double mean,
		  double sigma, long *seed, double x[], int n)
{
	int		i, k, m;
	double	s, *w;

	w = malloc(n * sizeof(double));
	for (k=0; k<n; k++)
	{
		w[k] = gauss(mean, sigma, seed);
	}

	x[0] = b[0] * w[0];
	for (k=1; k<=p; k++)
	{
		s = 0.0;
		for (i=1; i<=k; i++)
		{
			s += a[i] * x[k-i];
		}
		s = b[0] * w[k] - s;
		if (q == 0)
		{
			x[k] =s;
			continue;
		}
		m = (k>q) ? q : k;
		for (i=1; i<=m; i++)
		{
			s += b[i] * w[k-i];
		}
		x[k] = s;
	}

	for (k=p+1; k<n; k++)
	{
		s = 0.0;
		for (i=1; i<=p; i++)
		{
			s += a[i] * x[k-i];
		}
		s = b[0] * w[k] -s;
		if (q == 0)
		{
			x[k] =s;
			continue;
		}
		for (i=1; i<=q; i++)
		{
			s += b[i] * w[k-i];
		}
		x[k] = s;
	}
	free(w);
}

// 函数说明: 根据所给的随机数据，拟合AR模型
// 输入参数: double a[],	自回归部分的参数
//           int    *p,     自回归部分的阶数
//           double x[],    随机数据数组
//           int    n,      随机数据的总数
//           int    nmax,   所拟合模型的最大阶数
// 返回值  : 无

void marple(double x[], int n, double a[], int *p, int nmax)
{
	double	w0, g, v, z, u, h, s, DEN, w1, FPL, FPE;
	double	f, b, *c, *d, *c1, *d1, *r, *fi, *fi1;
	double	delta, alpha, beta1, gama1, beta2, gama2, beta3, gama3;
	int		i, j, k, tempj;

	c  = malloc((nmax+1)*sizeof(double));
	d  = malloc((nmax+1)*sizeof(double));
	c1 = malloc((nmax+1)*sizeof(double));
	d1 = malloc((nmax+1)*sizeof(double));
	r  = malloc((n+1)*sizeof(double));
	fi = malloc((nmax+1)*sizeof(double));
	fi1= malloc((nmax+1)*sizeof(double));

	for (w0=0, i=0; i<n; i++)
	{
		w0 += 2 * x[i] * x[i];
	}
	for (r[0]=0, i=0; i<n-1; i++)
	{
		r[0] += 2 * x[i] * x[i+1];
	}
	g     = x[0] * x[0] / w0;
	v     = g;
	z     = x[n-1] * x[n-1] / w0;
	u     = z;
	h     = x[0] * x[n-1] / w0;
	s     = h;
	DEN   = 1 - g - z;
	f     = x[0];
	b     = x[n-1];
	w1    = w0 * DEN;
	c[0]  = x[0] / w1;
	d[0]  = x[n-1] / w1;
	fi[0] = -r[0] / w1;
	w1    = w1 * (1 - fi[0]*fi[0]);
	FPL   = w1 / (2 * (n-1));
	FPL   = (n+1) * FPL / (n-1);

	for (i=0; i<nmax; )
	{
		f = x[i+1];
		b = x[n-i-2];
		for (j=0; j<=i; j++)			// calculate f and b
		{
			f += x[i-j] * fi[j];
			b += x[n-i+j-1] * fi[j];
		}
			
		for(j=0; j<=i;j++)								// 6.4.44
		{
			c1[j+1] = c[j] + f * fi[j] / w1;
			d1[j+1] = d[j] + b * fi[j] / w1;
		}
		for (j=0; j<=i; j++)
		{
			c[j+1] = c1[j+1];
			d[j+1] = d1[j+1];
		}
		c[0] = f / w1;
		d[0] = b / w1;
			
		g += f*f/w1 + (v*v*(1-z) + s*s*(1-g) + 2*s*h*v) / DEN;  // 6.4.47
		z += b*b/w1 + (s*s*(1-z) + u*u*(1-g) + 2*u*h*s) / DEN;
		
		for (s=0, h=0, v=0, u=0, j=0; j<=i+1; j++)		   
		{
			h += x[n-i-2+j] * c[j];
			s += x[n-j-1] * c[j];
			u += x[n-j-1] * d[j];
			v += x[j] * c[j];
		}

		DEN = (1-g)*(1-z) - h*h;
		if (DEN <= 0)
			break;

		alpha = 1 / (1 + (f*f*(1-z) + b*b*(1-g) +2*f*b*h)/w1/DEN);	// 6.4.38
		w1    = w1 * alpha;
		beta1 = (f*(1-z) + b*h) / DEN;
		gama1 = (b*(1-g) + f*h) / DEN;
		beta2 = (s*h + v*(1-z)) / DEN;
		gama2 = (v*h + s*(1-g)) / DEN;
		beta3 = (u*h + s*(1-z)) / DEN;
		gama3 = (s*h + u*(1-g)) / DEN;

		for (j=0; j<=i; j++)
			fi[j] = alpha * (fi[j] + beta1*c[j+1] + gama1*d[j+1]);     // 6.4.37

		k = (i+1) / 2;											 // 6.4.39
		for (j=0; j<=k; j++)
		{
			tempj = i - j + 1;
			c1[j] = c[j] + beta2*c[tempj] + gama2*d[tempj];
			d1[j] = d[j] + beta3*c[tempj] + gama3*d[tempj];
			if (tempj != j)
			{
				c1[tempj] =c[tempj]+ beta2*c[j] + gama2*d[j];
				d1[tempj] =d[tempj]+ beta3*c[j] + gama3*d[j];
			}
		}

		for (j=0; j<=i+1; j++)
		{
			c[j] = c1[j];
			d[j] = d1[j];
		}

		delta = 0.0;
		for (j=i; j>=0; j--)
		{
			r[j+1] = r[j] - x[n-j-2]*x[n-i-1] - x[i]*x[j];			   // 6.4.46
			delta += r[j+1] * fi[j];								   // 6.4.40
		}

		i++;

		for (r[0]=0, j=0; j<n-i-1; j++)
		{
			r[0] += 2 * x[j] * x[j+i+1];
		}
		delta += r[0];

		fi[i] = (-1) * delta / w1;										// 6.4.41

		for (j=0; j<i; j++)
		{
			tempj = i-j-1;
			fi1[j] = fi[j] + fi[i]*fi[tempj];
		}
		for (j=0; j<i; j++)
		{
			fi[j] = fi1[j];
		}

	    w1 *= (1 - fi[i]*fi[i]);                     // 6.4.43
	
		if ( (fi[i]*fi[i]) > 1)
		{
			printf("越界错误!:%d\n",i);
			break;
		}
		FPE = w1 * (n+i+2) / (2*(n-i-1)) / (n-i-2);
		if (FPE < FPL)
		{
			*p   = i+1;
			FPL  = FPE;
			for (j=0; j<=i; j++)
			{
				a[j] = fi[j];
			}
        }

    }

	printf("\nFPE: %f\n\n",FPE);

	free(c);
	free(d);
	free(c1);
	free(d1);
	free(r);
	free(fi);
	free(fi1);
}
   
int main(void)
{
    double a[5] = {1.0, -1.352, 1.338, -0.662, 0.24};
	double b[1] = {1.0};
	int	   p = 4;
	int    q = 0;
	long   seed = 13579;
	double mean = 0.0;
	double var = 1.0;
	double x[1000];
	double c[100];
    int    cc = 0;
	int    nmax = 20;
	int    i  = 0;
	int    n  = 1000;

	printf("所给AR模型的阶数为 %d。\n", p);
	for (i=0; i<=p; i++)
	{
		printf("a(%d) = %10.7lf\n", i, a[i]);
	}

	arma(a, b, p, q, mean, var, &seed, x, n);
	marple(x, n, c, &cc, nmax);

	printf("拟合的AR模型的阶数为 %d。\n", cc);
	printf("a(0) = %10.7lf\n", 1.0);
	for (i=0; i<cc; i++)
	{
		printf("a(%d) = %10.7lf\n", i+1, c[i]);
	}

	return 1;
}