作者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