C++

Forum

 
P-value calculation is not exact via F distribution in C++

yo2xia (20)
New code in C++, no precision error. share it with you!

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
 
 //*---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*
//*---*---*---*P-value for comparison of Full and Reduced models--*---*
//*---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*

//---get P value from F value and df1 df2---------------------------start
double c[11] = { 0.0000677106, -0.0003442342, 0.0015397681, -0.0024467480,
		0.0109736958, -0.0002109075, 0.0742379071, 0.0815782188, 0.4118402518,
		0.4227843370, 1.0000000000 };

double gammaln(double xx) {
	double x, y, tmp, ser;
	static const double cof[6] = {
		76.18009172947146,
		-86.50532032941677,
		24.0140982408091,
		-1.231739572460155,
		0.1208650973866179e-2,
		-0.5395239384953e-5
	};
	y = x = xx;
	tmp = (x + 0.5) * log(x + 5.5) - (x + 5.5);
	ser = 1.000000000190015;
	for (int j = 0; j < 6; j ++) {
		ser += cof[j] / (y + 1);
		y = y + 1;
	}
	return tmp + log(2.5066282746310005 * ser / x);
}

double beta(double x, double y) {
	if (x <= 0 || y <= 0) {
		return 0;
	}
	return exp(gammaln(x)+gammaln(y)-gammaln(x + y));
}

double abslt(double x) {
	if (x < 0.0)
		return -x;
	else
		return x;
}

double fi(int N, double x, double a, double b) {
	int n = N / 2;
	double f = 0.0, f1, s1, s2, tmpU, tmpV;
	int i;
	for (i = n; i >= 1; i--) {
		tmpU = (a + 2.0 * i - 1.0) * (a + 2.0 * i);
		s2 = i * (b - i) * x / tmpU;
		f1 = s2 / (1.0 + f);
		tmpV = (a + 2.0 * i - 2.0) * (a + 2.0 * i - 1.0);
		s1 = -(a + i - 1.0) * (b + a + i - 1.0) * x / tmpV;
		f = s1 / (1.0 + f1);
	}
	return 1.0 / (1.0 + f);
}

double incomBeta(double x, double a, double b) {
	if (a <= 0.0 || b <= 0.0) {
		return 0.0;
	}
	if (abslt(x - 0.0) < 1.0e-30 || abslt(x - 1.0) < 1.0e-30) {
		return 0.0;
	}

	double c1, c2, c3, f1, f2;
	int n;
	c1 = pow(x, a);
	c2 = pow(1.0 - x, b);
	c3 = beta(a, b);
	if (x < (a + 1.0) / (a + b + 2.0)) {
		n = 1;
		while (1) {
			f1 = fi(2 * n, x, a, b);
			f2 = fi(2 * n + 2, x, a, b);
			if (abslt(f2 - f1) < 1.0e-30)
				return f2 * c1 * c2 / a / c3;
			else
				n++;
		}
	} else {
		if (abslt(x - 0.5) < 1.0e-30 && abslt(a - b) < 1.0e-30)
			return 0.5;
		else {
			n = 1;
			while (1) {
				f1 = fi(2 * n, 1.0 - x, b, a);
				f2 = fi(2 * n + 2, 1.0 - x, b, a);
				if (abslt(f2 - f1) < 1.0e-30)
					return 1.0 - f2 * c1 * c2 / b / c3;
				else
					n++;
			}
		}
	}
	return 0;
}

double getPvalue(double f, double n1, double n2) {
	if (f < 0.0)
		f = -f;
	return incomBeta(n2 / (n2 + n1 * f), n2 / 2.0, n1 / 2.0);
}
//---get P value from F value and df1 df2---------------------------end
     
    int main()  
    {  
      double y;  
      y=FDist(304.94,17,18);  
      std::cout << "P-value =" << std::endl;
      std::cout << y << std::endl; 
      system("pause");  
    }



Sorry, I am late to reply all of you. And this code is right. share it with you!
Last edited on
MiiNiPaa (8886)
1) This might be precision issue. Try to step with a debugger and find out when results are not expected.
2) Comparing floats with == and != is dangerous as (because of aformentioned precision issues) even values which logically should be equal is not. (for example 0.1 + 0.1 +0.1 +0.1 +0.1 == 0.5 might yield false)
TheIdeasMan (6817)
Hi,

std::cout defaults to 6 dp, so on line 213, it probably will print 0. Look up how to std::cout in scientific form. It might just be that.

Does that level of precision really matter? Is the original data that accurate? I don't know anything about the subject - just putting it out there.

Should there be some precision value where anything less effectively is zero? I am sure there should be if you want to compare properly like MiiNiPaa is saying - you need to check if a value is in the range of x +/- precision for it to be equal. Maybe you already done that - not sure what all the variables mean. I read somewhere that exact decimals are coming in C++17.

The other thing that bugs me is all the un-initialised variables at declaration - you probably have set them all, but it is asking for trouble IMO. For instance, I can never remember whether local variables are initialised to 0 automatically, so I always set everything to something at declaration. What is the value of p on lines 17 & 19? p is set on line 23, which makes the assignment on line 19 irrelevant. One can set warnings for un-initialised variables.

Hope this helped - even a little :+)

I am supposed to be in bed !! Cheers
jlb (4973)
This looks problematic to me.
1
2
3
4
5
6
7
8
9
10
11
12
 
    double FDist(double F,double m,double n)  
    {  
        double xx,p;  
      
        if(m<=0 || n<=0) p=-1;  
        else if(F>0)  
        {  
            xx=F/(F+n/m);  
            p=betainc(xx,m/2,n/2);  
        }  
        return(1-p);  
    }


It seems that it might be possible to use p before it is initialized. What happens is m > 0 and n > 0 and F < 0?
yo2xia (20)
Sorry, I am late to reply you. And this is the code with no problem, share it with you. Good luck to all of you.
Topic archived. No new replies allowed.