다이어리/다이어리

Syntax Highlighter 테스트 for MATLAB & R

Jae-seong Yoo 2014. 2. 1. 06:56
매트랩 테스트
% Example 2.11
% Computational Statistics Handbook with MATLAB, 2nd Edition
% Wendy L. and Angel R. Martinez

% get the mean and covariance
mu = zeros(1,2);
cov_mat = eye(2);% Identity matrix
% get the domain
% should range (-4,4) in both directions
[x,y]=meshgrid(-4:.2:4,-4:.2:4);
% reshape into the proper format for the function
X=[x(:),y(:)];
Z = csevalnorm(X,mu,cov_mat);
% now reshape the matrix for plotting
z=reshape(Z,size(x));
subplot(1,2,1) % plot the surface
surf(x,y,z),axis square, axis tight
title('BIVARIATE STANDARD NORMAL')
subplot(1,2,2) % look down on the surface
pcolor(x,y,z),axis square
title('BIVARIATE STANDARD NORMAL')

% now do the same thing for a covariance matrix
% with non-zero off-diagonal elements
cov_mat = [1 0.7 ; 0.7 1];
Z=csevalnorm(X,mu,cov_mat);
z=reshape(Z,size(x));
subplot(1,2,1)
surf(x,y,z),axis square, axis tight
title('BIVARIATE NORMAL')
subplot(1,2,2)
pcolor(x,y,z),axis square
title('BIVARIATE NORMAL')
 
 


R 테스트
# Example 2.11
# Computational Statistics Handbook with MATLAB, 2nd Edition
# Wendy L. and Angel R. Martinez
# Edited by Jae-seong Yoo at 20140127

csevalnorm = function(x,mu,cov_mat)
{
	# CSEVALNORM Multivariate normal probability density function.
	#
	#	Y = CSEVALNORM(X,MU,COVM) Returns the value of the multivariate
	#	probability density function at the locations given in X.
	#
	#	INPUTS:		X is an n x d matrix of domain locations.
	#				MU is a 1 x d vector
	#				COVM is a d x d covariance matrix
	#
	# 	See also CSNORMP, CSNORMC, CSNORMQ

	#   W. L. and A. R. Martinez, 9/15/01
	#   Computational Statistics Toolbox 

	n = dim(x)[1];
	d = dim(x)[2];
	prob = matrix(0,n,1)
	a=(2*pi)^(d/2)*sqrt(det(cov_mat));
	covi = solve(cov_mat);
	for (i in 1:n)
	{
		xc = x[i,]-mu;
		arg=xc%*%covi%*%t(xc);
		prob[i]=exp((-0.5)*arg);
	}
	prob=prob/a;
	return(prob)
}



# get the mean and covariance
mu = matrix(0,1,2);
cov_mat = diag(1, 2);	# Identity matrix

# get the domain
# should range (-4,4) in both directions
n <- seq(-4, 4, by=0.2)
x = t(matrix(rep(n, 41), nrow=41))
y = matrix(rep(n, 41), nrow=41)

# reshape into the proper format for the function
X = cbind(as.vector(x), as.vector(y))

Z = csevalnorm(X,mu,cov_mat);
# now reshape the matrix for plotting
z=matrix(Z, dim(x)[1], dim(x)[2])

op <- par(mfrow = c(1, 3))

persp(z)
title("BIVARIATE # NORMAL")

plot(x, z)
title('BIVARIATE # NORMAL')

plot(z)



# now do the same thing for a covariance matrix
# with non-zero off-diagonal elements
cov_mat = matrix(c(1, 0.7, 0.7, 1), nrow=2);
Z=csevalnorm(X,mu,cov_mat);
z=matrix(Z, dim(x)[1], dim(x)[2])

op <- par(mfrow = c(1, 3))

persp(z)
title("BIVARIATE NORMAL")

plot(x, z)
title('BIVARIATE NORMAL')

plot(z)
 
 



잘 되는 것 같기는 하다.

썩 깔끔하지는 않은 것 같지만...