?? linear.java
字號:
/*
* 創建日期 2006-7-12
*
* TODO 要更改此生成的文件的模板，請轉至
* 窗口 － 首選項 － Java － 代碼樣式 － 代碼模板
*/
/**
* @author new
*
* TODO 要更改此生成的類型注釋的模板，請轉至
* 窗口 － 首選項 － Java － 代碼樣式 － 代碼模板
*/
import java.lang.Math;
import java.io.*;
import java.util.Vector;
public class linear{
private double[] dataX; //存放自變量數據;
private double[] dataY; //存放因變量數據;
private double interc; //存放截距;
private double coeffi; //存放回歸系數;
private double F; //存放F檢驗值;
private double t; //存放t檢驗值;
int N;
public linear(double[] X,double[] Y){
dataX=X;
dataY=Y;
N=dataX.length;
}
public linear( ){
}
//一元線性回歸
//datax(1 To n)：自變量，n為觀測次數，已知
//datay(1 To n)：因變量，n為觀測次數，已知
//interc：截距，計算結果
//coeffi：回歸系數，計算結果
//F：F檢驗值，計算結果
//t：t檢驗值，計算結果
public void LinearMethod(){
double Xa=0,Ya=0,Sxx=0,Sxy=0,Syy=0,SSR=0,SSE=0,Ur,Ue,Syx2=0,Sb,Sb2,Sx=0;
int n,i;
n=dataX.length;
for(i=0;i<n;i++){
Xa=Xa + dataX[i];
Ya=Ya + dataY[i];
}
//平均值
Xa=Xa/n;
Ya=Ya/n;
for(i=0;i<n;i++){
Sxx=Sxx + (dataX[i] - Xa)*(dataX[i] - Xa);
Sxy=Sxy + (dataX[i] - Xa)*(dataY[i] - Ya);
}
coeffi = Sxy / Sxx ; //回歸系數
interc = Ya - coeffi * Xa ; //截距
//由回歸所導致的方差
for(i=0;i<n;i++){
SSR = SSR + (interc + coeffi * dataX[i] - Ya)*(interc + coeffi * dataX[i] - Ya);
}
Ur = 1; //回歸方差自由度
SSR = SSR/Ur;
//由剩余所導致的方差
for(i=0;i<n;i++){
SSE = SSE + (dataY[i] - interc-coeffi * dataX[i])*(dataY[i] - interc-coeffi * dataX[i]);
}
//求F值
Ue = n - 2; //剩余方差自由度
SSE = SSE / Ue;
F = SSR / SSE;
//求t值
for(i=0;i<n;i++){
Syx2 = Syx2 + (dataY[i] - (interc+coeffi * dataX[i])) * (dataY[i] - (interc+coeffi * dataX[i]));
}
Syx2 = Syx2 / (n - 2);
for(i=0;i<n;i++){
Sx = Sx + (dataX[i] - Xa) * (dataX[i] - Xa);
}
Sb2 = Syx2 / Sx;
Sb = Math.sqrt(Sb2);
t = Math.abs(coeffi / Sb);
}
//求正態分布的分位數
//Q：上側概率
//x：分位數
public double PNorm(double Q){
double x;
double p,y,z;
double b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10;
double b;
b0 = 1.570796288;
b1 = 0.03706987906;
b2 = -0.0008364353589;
b3 = -0.0002250947176;
b4 = 0.000006841218299;
b5 = 0.000005824238515;
b6 = -0.00000104527497;
b7 = 8.360937017E-08;
b8 = -3.231081277E-09;
b9 = 3.657763036E-11;
b10 = 6.936233982E-13;
if(Q==0.5)
return 0;
if(Q>0.5)
p=1-Q;
else
p=Q;
y = -Math.log(4 * p * (1 - p));
b = y * (b9 + y * b10);
b = y * (b8 + b);
b = y * (b7 + b);
b = y * (b6 + b);
b = y * (b5 + b);
b = y * (b4 + b);
b = y * (b3 + b);
b = y * (b2 + b);
b = y * (b1 + b);
z = y * (b0 + b);
x = Math.sqrt(z);
if (Q > 0.5 )
x = -x;
return x;
}
//計算F分布的分布函數
//n1：自由度，已知
//n2：自由度，已知
//F：F值，已知
//p：下側概率，所求
//d：概率密度，所求
public double[] F_DIST(int n1, int n2,double F){
double p,d;
double[] ret = new double[2];
double x,u,Lu;
int IAI,IBI,nn1,nn2;
int I;
double PI = 3.14159265359;
if(F == 0)
{
p=0;d=0;
ret[0]=p;
ret[1]=d;
return ret;
}
x=n1*F/(n2+n1*F);
if(n1%2==0)
if(n2%2==0)
{u=x*(1-x);p=x;IAI=2;IBI=2;}
else
{
u=x*(double)Math.sqrt(1-x)/2;
p=1-(double)Math.sqrt(1-x);
IAI=2;IBI=1;
}
else
if(n2%2==0){
p=(double)Math.sqrt(x);
u=p*(1-x)/2;
IAI=1;IBI=2;
}
else{
u=(double)Math.sqrt(x*(1-x))/PI;
p=1-2*(double)Math.atan((double)Math.sqrt((1-x)/x))/PI;
IAI=1;IBI=1;
}
nn1=n1-2;
nn2=n2-2;
if(u==0){
d=(double)u/F;
ret[0]=p;
ret[1]=d;
return ret;
}
else{
Lu=(double)Math.log(u);
}
if(IAI!=n1){
for(I=IAI;I<=nn1;I=I+2){
p=p-2*(double)u/I;
Lu=Lu+(double)Math.log((1+IBI/I)*x);
u=(double)Math.exp(Lu);
}
}
if(IBI==n2){
d=(double)u/F;
ret[0]=p;
ret[1]=d;
return ret;
}
for(I=IBI;I<=nn2;I=I+2){
p=p+2*(double)u/I;
Lu=Lu+Math.log((1+(double)n1/I)*(1-x));
u=(double)Math.exp(Lu);
}
d=(double)u/F;
ret[0]=p;
ret[1]=d;
return ret;
}
//計算F分布的分位數
//n1：自由度，已知
//n2：自由度，已知
//Q：上側概率，已知
//F：分位數，所求
public double PF_DIST(int n1, int n2, double Q){
double F=0.0;
double[] re = new double[2];
double DF12,DF22,A,b;
double A1,b1,p,YQ;
double E,FO,pp,d;
double GA1,GA2,GA3;
int K;
DF12=(double)n1/2;
DF22=(double)n2/2;
A=(double)2/(9*n1);
A1=1-A;
b=(double)2/(9*n2);
b1=1-b;
p=1-Q;
YQ=PNorm(Q);////////////////////////////////////////////
E=b1*b1-b*YQ*YQ;
if(E>0.8)
{ FO=(double)Math.pow(((A1*b1+YQ*(double)Math.sqrt(A1*A1*b+A*E))/E),3);
}
else{
GA1=lnGamma(DF12+DF22);
GA2=lnGamma(DF12);
GA3=lnGamma(DF22);
FO=(double)((double)2/n2)*(GA1-GA2-GA3+0.69315+(DF22-1)*(double)Math.log(n2)-DF22*(double)Math.log(n1)-(double)Math.log(Q));
FO=(double)Math.exp(FO);
}
for(K=1;K<=30;K++){
re= F_DIST(n1,n2,FO);
pp=re[0];
d=re[1];
if(d==0){F=FO; return F;}
F = FO - (double)(pp - p) / d;
if ((double)Math.abs(FO - F) < 0.000001 * (double)Math.abs(F))
{ return F;}
else FO = F;
}
return F;
}
//計算GAMMA函數
//x：自變量
//z：GAMMA函數值
public double GAMMA(double x){
double z;
double H,y,y1 ;
H = 1;
y = x;
double temp;
do{
temp=y;
if(y==2) {z=H;return z;}
else if(y<2) {H = H / y; y = y + 1;}
else if(y>=3) {y = y - 1;H = H * y; }
}while(temp<2||temp>=3);
y = y - 2;
y1 = y * (0.005159 + y * 0.001606);
y1 = y * (0.004451 + y1);
y1 = y * (0.07211 + y1);
y1 = y * (0.082112 + y1);
y1 = y * (0.41174 + y1);
y1 = y * (0.422787 + y1);
H = H * (0.999999 + y1);
z = H;
return z;
}
//求Gamma函數的對數LogGamma(x)
//x：自變量
//G：Gamma函數的對數
public double lnGamma(double x){
double G;
double y,z,A;
double b,b1;
int n,I;
if(x<8)
{y=x+8;n=-1;}
else
{y=x;n=1;}
z=1/(y*y);
A=(y-0.5)*(double)Math.log(y)-y+0.9189385;
b1 = (0.0007663452 * z - 0.0005940956) * z;
b1 = (b1 + 0.0007936431) * z;
b1 = (b1 - 0.002777778) * z;
b = (double)(b1 + 0.0833333) / y;
G = A + b;
if (n >= 0 ) return G;
y = y - 1;
A = y;
for(I = 1;I<=7;I++)
A = A * (y - I);
G = G - (double)Math.log(A);
return G;
}
//計算t分布的分布函數
//n：自由度，已知
//T：t值，已知
//pp：下側概率，所求
//dd：概率密度，所求
public double[] T_Dist(int n, double t){
double pp,dd;
double[] ret=new double[2];
int Sign;
double TT,x;
double p,u,GA1,GA2;
int IBI,N2,I;
double PI=3.14159265359;
if(t==0){
GA1=GAMMA(n / 2);
GA2=GAMMA(n/2+0.5);
pp=0.5;
dd=GA2/(Math.sqrt(n*PI)*GA1);
ret[0]=pp;
ret[1]=dd;
return ret;
}
if(t<0)
Sign=-1;
else
Sign=1;
TT = t * t;
x = TT / (n + TT);
if(n%2==0) { //n為偶數
p = Math.sqrt(x);
u = p * (1 - x) / 2;
IBI = 2;
}
else { //n為奇數
u = Math.sqrt(x * (1 - x)) / PI;
p = 1 - 2 * Math.atan(Math.sqrt((1 - x) / x)) / PI;
IBI = 1;
}
if(IBI!=n){
N2=n-2;
for(I=IBI;I<=N2;I=I+2){
p = p + 2 * u / I;
u = u * (1 + I) / I * (1 - x);
}
}
dd=u/Math.abs(t);
pp=0.5+Sign*p/2;
ret[0]=pp;
ret[1]=dd;
return ret;
}
// 求t分布的分位數
//n：自由度，已知
//Q：上側概率（<=0.5），已知
//T：分位數，所求
public double PT_DIST(int n,double Q){
double t=0;
double PIS, DFR2, C;
double Q2, p, YQ=0, E;
double GA1, GA2, GA3;
double T0, pp, d;
int K;
double PI=3.14159265359;
PIS = Math.sqrt(PI);
DFR2 = n / 2;
if (n == 1 ){
t = Math.tan(PI * (0.5 - Q));
return t;}
if (n == 2 ){
if (Q > 0.5 ) C = -1;
else C = 1;
Q2 = (1 - 2 * Q) *(1 - 2 * Q);
t = Math.sqrt(2 * Q2 / (1 - Q2)) * C;
return t;
}
p = 1 - Q;
YQ=PNorm (Q ); //正態分布分位數
E = (1 - 1 / (4 * n)) *(1 - 1 / (4 * n)) - YQ * YQ / (2 * n);
if( E > 0.5 )
T0 = YQ / Math.sqrt(E);
else{
GA1=lnGamma( DFR2);
GA2=lnGamma (DFR2 + 0.5 );
GA3 = Math.exp((GA1 - GA2) / n);
T0 = Math.sqrt(n) / Math.pow((PIS * Q * n) , (1 / n)) / GA3;
}
double[] re=new double[2];
for (K = 1 ;K<=30;K++){
re=T_Dist( n, T0);
pp=re[0];
d=re[1];
if (d == 0) {
t = T0;
return t;}
t = T0 - (pp - p) / d;
if (Math.abs(T0 - t) < 0.000001 * Math.abs(t) )
return t;
else T0 = t;
}
return t;
}
//從文件中讀取要分析的數據
public void ReadText(String fileName){
int num=0;
Vector v1=new Vector();
Vector v2=new Vector();
try {
FileInputStream fstream = new FileInputStream(fileName);
//FileInputStream fstream = new FileInputStream("data.txt");
DataInputStream in = new DataInputStream(fstream);
String currentLine;
char[] c;
int index1,index2;
while (in.available() !=0) {
num++;
currentLine=in.readLine();
System.out.println (currentLine);
c=currentLine.toCharArray();
index1=currentLine.indexOf(' ');
index2=currentLine.indexOf(' ');
String sub1=currentLine.substring(0,index1);
v1.add(sub1);
String sub2=currentLine.substring(index1+1);
v2.add(sub2);
}
in.close();
} catch (Exception e) {
System.err.println("File input error");
}
N=num;
dataX=new double[num];
dataY=new double[num];
for(int i=0;i<num;i++)
{ dataX[i]=i+1;
dataY[i]=Double.parseDouble((String)v2.elementAt(i));
}
}
public static void main(String args[]){
System.out.print("請輸入要分析的文件路徑：");
linear ln = new linear();
String name=null;
byte [] b=null;
try{
BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
name=in.readLine();
} catch (Exception e) {
System.err.println("IO exception!");
}
System.out.println("要分析的數據如下：");
ln.ReadText(name);
ln.LinearMethod();
System.out.println(" 回歸方程為： "+ "Y = "+ln.coeffi+" * X"+" + "+ln.interc);
double F005A=ln.PF_DIST(1,ln.N-2,0.05);
double F001A=ln.PF_DIST(1,ln.N-2,0.01);
System.out.println(" F檢驗值得："+ln.F);
System.out.println(" 顯著水平為0.05時的F臨界值為："+F005A);
System.out.println(" 顯著水平為0.01時的F臨界值為："+F001A);
System.out.println(" F檢驗結論：");
if(ln.F<=F005A)
System.out.println(" 回歸方程線性關系 不顯著 ");
if(ln.F>F005A&&ln.F<=F001A)
System.out.println(" 回歸方程線性關系 顯著 ");
if(ln.F>F001A)
System.out.println(" 回歸方程線性關系 特別顯著 ");
double F005B=ln.PT_DIST(ln.N-2,0.05/2);
double F001B=ln.PT_DIST(ln.N-2,0.01/2);
System.out.println(" t檢驗值得："+ln.t);
System.out.println(" 顯著水平為0.05時的t臨界值為："+F005B);
System.out.println(" 顯著水平為0.01時的t臨界值為："+F001B);
System.out.println(" t檢驗結論：");
if(ln.F<=F005B)
System.out.println(" 回歸系數 不顯著大于0");
if(ln.F>F005B&&ln.F<=F001B)
System.out.println(" 回歸系數 顯著大于0 ");
if(ln.F>F001B)
System.out.println(" 回歸系數 特別顯著大于0 ");
}
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -