话不多说,直接上代码了,c代码
// 2.1、使用中位数法,而非算术中位法
qsort(tmp_v_res, ns, sizeof(double), cmpres);
median_v_res =median(tmp_v_res, ns);// 观测值中数 Eq.14
// 2.1.1 求取偏离量,中心化
for (j = 0; j < ns; j++) {
tmp_v_res[j] -= median_v_res;// Eq.15
/*v_res[j] -= median_v_res;*/
if (tmp_v_res[j] < 0) tmp_v_res[j] *= -1;//取绝对值
}
// 2.1.2 使用中位数计算方差因子
qsort(tmp_v_res, ns, sizeof(double), cmpres);
median_v_res = median(tmp_v_res, ns);//观测值中数
sigma1 = median_v_res / 0.6745;// Eq.17
printf("抗差单位权中误差(weighted)sigma1= %f\n", sigma1);
/* 残差比较 */
static int cmpres(const void *p1, const void *p2){
double *q1 = (double *)p1, *q2 = (double *)p2;
double tt = (*q1) - (*q2);
return tt<-0.0 ? -1 : (tt>0.0 ? 1 : 0);
}
/* 求中位数 */
static double median(const double* A,unsigned int n){
if (n>0){
return (n & 1 == 1) ? A[(n - 1) / 2] : (A[n / 2] + A[n / 2 - 1]) / 2;
}
return 0;
}