编程习惯的重要性—

最近学习GSL，顺便把手上的一个非参估计的东西重写。大家知道，非参估计里面有个很大的问题就是窗宽选择，Cross-Validation是个很不错的想法，但是因为要最优化，很慢，所以我打算用C重写一遍。

先是非参数估计的程序，代码如下：

免费会员网

void KernelSelfRegs(double *yhat, const double *y, const double *x, const double *h, int N, int K, int LOO)
{
	double sum_k[N],sum_y[N];
	int i,j;
	double ker;
	double sk,sy;
	for (i=1;i<N;++i)
	{
	  sk=0.0;sy=0.0;
	  #pragma omp parallel for reduction(+:sk,sy) private(ker)
	  for (j=0;j<i;++j)
	  {
	    ker=MultiKernel(x+i*K,x+j*K,h,K);
	    *(sum_k+j)+=ker;
	    *(sum_y+j)+=(ker*(*(y+i)));
	    sk=sk+ker;
	    sy=sy+ker*(*(y+j));
	  }

  

   免费会员网
  

	  *(sum_k+i)+=sk;
	  *(sum_y+i)+=sy;
	}
	#pragma omp parallel for
	for (i=0;i<N;++i)
	{
		if (LOO==1){
			*(yhat+i)=*(sum_y+i)/(*(sum_k+i));
		}
		else{
			*(yhat+i)=(*(sum_y+i)+pow(Kernel(0),K)*(*(y+i)))/(*(sum_k+i)+pow(Kernel(0),K));
		}
	}
}

免费会员网

然后写了一个测试的程序：

void main(){
	int N=1000,K=2;
	double y[N],u[N],x[N][K];
	double h[]={1.807767,1.713608};
	int i;
	double yhat[N];
	for (i=0;i<N;i++){
		u[i]=(double)rand()/RAND_MAX-0.5;//random number in (-1,1)
		x[i][0]=10*((double)rand()/RAND_MAX);//random number in (1,10)
		x[i][1]=5*((double)rand()/RAND_MAX);//random number in (1,5)
		y[i]=sin(x[i][0])+x[i][1]+u[i];
	}
	time_t start=time(NULL);
	// KernelRegs(yhat,y,*x,h,*x,N,K,N);
	KernelSelfRegs(yhat,y,*x,h,N,K,1);
	time_t end=time(NULL);
	for (i=0;i<N;i++)
		printf("y=%f, The true value is %f, and the est-value is %f\n",y[i],sin(x[i][0])+x[i][1],yhat[i]);
	printf("Time: %f\n",(double)(end-start));