gum_least.c

计算gumbul极值的一段程序,没有附带算例,因为太简单了

		const double Gumbel_Least(const double OrginalData[],const int iSize,const int iYear,double& gumbelA,double& gumbelB) 
		//OrginalData	:orginal input data vector of size iSize;
    //iSize		:size of OrginalData
    //gumbelA,B :parameter of a,b,c in gumbel extreme statistic
     {
			double   x_ave, y_ave, sigma_x, sigma_y, rxy, p, xtr_value, *y;
			int     i, isize;

			y=new double[iSize];
			x_ave = 0.0;
			y_ave = 0.0;
			sigma_x = 0.0;
			sigma_y = 0.0;
			rxy = 0.0;
			isize = iSize - 1;
            
			for ( i = 0; i <= isize; i++ )  {
			    
				y[i] = (double) (i + 1) / ( iSize + 1 );
				y[i] = -log ( -log ( 1 - y[i] ) );
				
				x_ave += OrginalData[i];
				y_ave += y[i];
				
				sigma_x += OrginalData[i] * OrginalData[i];
				sigma_y += y[i] * y[i];
				rxy += OrginalData[i] * y[i];
			}

			x_ave /= iSize;
			y_ave /= iSize;
			
			sigma_x = sigma_x / iSize - x_ave * x_ave;
			sigma_y = sigma_y / iSize - y_ave * y_ave;
			rxy     = rxy / iSize -x_ave * y_ave;
			
			sigma_x = sqrt ( sigma_x );
			sigma_y = sqrt ( sigma_y );
			rxy     = rxy / ( sigma_x * sigma_y );
			
			gumbelB = rxy * sigma_x / sigma_y;
	    gumbelA = x_ave - gumbelB * y_ave;
	    
	    p = 1 - 1.0 / (int) (1==iYear?1.01:iYear);
		  xtr_value = gumbelA - gumbelB * log( -log(p) );
	   			
	  	return xtr_value;
		}