作者markoo (上行下笑阿..)
看板NCTU-STAT98G
標題0730 程式 EX 4
時間Sun Aug 2 01:42:01 2009
EX4: 由回歸模型生成一組樣本Y = X*beta + error
這邊只介紹最小平方法(MLE解再normal時與LSE相同）
當error是exponential(1/16)時,
程式執行結果可以發現回歸係數估計為(20,3,-6)左右與真實回歸係數(4,3-6)
剛好差了exponential的mean 16.
這邊可以想像當樣本夠大時Exp(1/16) - 16 可以趨近常態分配,
所以截距項的部分就會估計成4+16=20的結果.
#include <stdio.h>
#include "matrix.h"
#include "rand.h"
#include <math.h>
int main()
{
int i,j;
int n=2000,k=3;
//n:sample size, k: model dimension (intercept included)
double *temp;
double seed = 100;
temp = rand_uniform(2*n,1,2,seed);
// 這邊假設解釋變數都是Uniform(1,2)的數值
double **X;
X=matrix_zero(n,k);
for (i=0;i<n;i++)
{
X[i][0]=1;
for (j=1;j<k;j++) X[i][j] = temp[i*2+(j-1)];
}
//這邊把解釋變數排成矩陣X
free(temp);
matrix_print(n,k,X);
double **error1;
double **error2;
temp = rand_normal(n,0,9,seed);
error1 = matrix_zero(n,1);
for (i=0;i<n;i++) error1[i][0] = temp[i];
free(temp);
//生成normal的error 並排成矩陣 然後把後續沒有用到的temp 記憶體釋放掉
temp = rand_exponential(n,1/16.0,seed);
//printf("%lf\n",rand_mean(n,temp));
error2 = matrix_zero(n,1);
for (i=0;i<n;i++) error2[i][0] = temp[i];
free(temp);
//生成exponential的error 並排成矩陣,同樣做free的動作
double **beta;
beta = matrix_zero(k,1);
beta[0][0] = 4;
beta[1][0] = 3;
beta[2][0] = -6;
//設定回歸係數
double **Y;
Y = matrix_product(n,k,X,k,1,beta);
Y = matrix_plus(n,1,Y,error2);
matrix_print(n,1,Y);
//生成樣本Y = X*beta + error
double **hat_beta;
hat_beta = matrix_trans(n,k,X);
hat_beta = matrix_product(k,n,hat_beta,n,k,X);
hat_beta = matrix_inverse(k,k,hat_beta);
hat_beta = matrix_product(k,k,hat_beta,k,n,matrix_trans(n,k,X));
hat_beta = matrix_product(k,n,hat_beta,n,1,Y);
printf("Hat_beta:\n");
matrix_print(k,1,hat_beta);
//這邊做最小平方法
return 0;
}
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 118.166.48.157