%---------------------------------------------------
%   M script and M function files  
%---------------------------------------------------

Example 1
function [s p]=fcircle(r)
% FCIRCLE calculates the area and perimeter of a circle of radii r
% r  圆半径
% s  圆面积
% p  圆周长

%2009.11.30日编写
s=pi*r*r;
p=2*pi*r;


%-END of FUNCTION---------------
%在命令窗口
[s p]=fcircle(10)


Example 2
% tran.m
function [rho theta]=tran(x,y)
rho=sqrt(x^2+y^2);
theta=atan(y/x);
%-END of FUNCTION---------------
%在命令窗口
x=input('Please input x=');
y=input('Please input y=');
[rho,theta]=tran(x,y)

Example 3
% wadd.m
function f=wadd(x,y)
global alpha beta
f=alpha*x+beta*y;
%-END of FUNCTION---------------
%在命令窗口
global alpha beta
alpha=1; beta=2;
wadd(1,2)


Example 4
% anonymous functions
sqr = @(x) x.^2;
a = sqr(5)
quad(sqr, 0, 1)  % integrate from 0 to 1

sumAxBy = @(x, y) (A*x + B*y);
sumAxBy(5, 7)

A = {@(x)x.^2, @(y)y+10, @(x,y)x.^2+y+10} %cell array
A{1}(4) + A{2}(7)


a = 1.3;   b = .2;   c = 30;
parabola = @(x) a*x.^2 + b*x + c;
fplot(parabola, [-25 25])
-------------------------
a = -3.9;   b = 52;   c = 0;
fplot(parabola, [-25 25])
-------------------------
a = -3.9;   b = 52;   c = 0;
parabola = @(x) a*x.^2 + b*x + c;
fplot(parabola, [-25 25])

Example 5 
% Gamma density plot
x = 0:.05:3;
y1 = gampdf(x,1,1/1);
y2 = gampdf(x,2,1/2);
y3 = gampdf(x,3,1/3);
% Plot the functions.
plot(x,y1,'r',x,y2,'g',x,y3,'b')
title('Gamma Distribution')
xlabel('X')
ylabel('f(x)')
text(0.141,0.9,'\leftarrow\lambda=t=1','FontSize',18)
text(1.075,0.7,'\leftarrow\lambda=t=3','FontSize',18)
text(0.68,0.4,'\lambda=t=2\rightarrow','FontSize',18)

%Example Linear regression with matlab
load carsmall
x1 = Weight;
x2 = Horsepower; % Contains NaN data
y = MPG;

%Compute regression coefficients for 
%a linear model with an interaction term:

X = [ones(size(x1)) x1 x2 x1.*x2];
b = regress(y,X) % Removes NaN data
Plot the data and the model:

scatter3(x1,x2,y,'filled')
hold on
x1fit = min(x1):100:max(x1);
x2fit = min(x2):10:max(x2);
[X1FIT,X2FIT] = meshgrid(x1fit,x2fit);
YFIT = b(1) + b(2)*X1FIT + b(3)*X2FIT + b(4)*X1FIT.*X2FIT;
mesh(X1FIT,X2FIT,YFIT)
xlabel('Weight')
ylabel('Horsepower')
zlabel('MPG')
view(50,10)
%----------------------------

Example 6
% Monte Carlo Hypothesis testing 
x=[465 449 468 446 447 456 442 433 447 456 438 449 ...
463 452 435 447 447 450 444 450 454 456 456 454 446];
n=length(x);
M = 1000;
alpha = 0.05;
sigma=7.8; sigxbar=sigma/sqrt(n);
% Get the critical value, using z as test statistic.
cv = norminv(alpha,0,1);
% Start the simulation.
Im = 0;
for i = 1:M
% Generate a random sample under H_0.
xs = sigma*randn(1,n) + 454;
Tm = (mean(xs)-454)/sigxbar;
if Tm <= cv % then reject H_0
Im = Im +1;
end
end
alphahat = Im/M;

% Now check the probability of Type II error.
% Get some alternative hypotheses:
mualt = 445:458;
betahat = zeros(size(mualt));
for j = 1:length(mualt)
Im = 0;
% Get the true mean.
mu = mualt(j);
for i = 1:M
% Generate a sample from H_1.
xs = sigma*randn(1,n) + mu;
Tm = (mean(xs)-454)/sigxbar;
if Tm > cv % Then did not reject H_0.
Im = Im +1;
end
end
betahat(j) = Im/M;
end
% Get the estimated power.
powhat = 1-betahat;
plot(powhat,'b*-')

Example 7
% MCMC
% function scsample.m for below Gibbs sampling
function rs = cssample(x,p,n)
% CSSAMPLE Random sample from an arbitrary discrete/finite distribution.
%
%	R = CSSAMPLE(X,P,N) This function will take an arbitrary discrete.
%	finite distribution and return a random sample from it.
%	The domain of the function is X. These are the values that the random
%	variable can assume. The probability associated with each one is given in the
% 	vector P. The number of variates that will be generated is N.

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

if length(x) ~= length(p)
   error('The size of the input vectors do not match.')
   return
end
% sort just in case they are not in order
[xs,ind]=sort(x);
ps = p(ind);	% sort these in the same order as x
% Get the cdf
F = cumsum(ps);
% Find all of the required variates
for i=1:n
   u = rand(1,1);
   if u<= F(1)
      rs(i) = x(1);
   elseif u > F(end-1)
      rs(i) = x(end);
   else
      ind = find(u <= F);
      rs(i) = xs(ind(1));
   end
end
%-END of FUNCTION---------------
%在命令窗口
% y contains number of disasters.
y=[4,5,4,1,0,4,3,4,0,6,3,3,4,0,2,6,3,3,5,4,5,3,...
1,4,4,1,5,5,3,4,2,5,2,2,3,4,2,1,3,2,2,1,1,1,1,...
3,0,0,1,0,1,1,0,0,3,1,0,3,2,2,0,1,1,1,0,1,0,1,...
0,0,0,2,1,0,0,0,1,1,0,2,3,3,1,1,2,1,1,1,1,2,4, ...
2,0,0,0,1,4,0,0,0,1,0,0,0,0,0,1,0,0,1,0,1];
% year contains the year.
year=1:112+1850;

n = length(y);
m = 1100;	% number in chain
a1 = 0.5; a2 = 0.5; c1 = 0; c2 = 0; d1 = 1; d2 = 1;
theta = zeros(1,m);
lambda = zeros(1,m);
k = zeros(1,n);
% Holds probabilities for k.
like = zeros(1,n);
% Get starting points.
k(1) = unidrnd(n,1,1);
% Note that k will indicate an index to the year
% that corresponds to a hypothesized change-point.
theta(1) = 1; lambda(1) = 1; b1 = 1; b2 = 1;
% Start the Gibbs Sampler.
for i = 2:m
   kk = k(i-1);
   % Get parameters for generating theta
   t = a1 + sum(y(1:kk));
   lam = kk + b1;
   % Generate the variate for theta.
   theta(i) = gamrnd(t,1/lam,1,1);
   % Get parameters for generating lambda.
   t = a2 + sum(y) - sum(y(1:kk));
   lam = n-kk+b2;
   % Generate the variate for lambda.
   lambda(i) = gamrnd(t,1/lam,1,1);
   % Generate the parameters b1 and b2.
   b1 = gamrnd(a1+c1,1/(theta(i)+d1),1,1);
   b2 = gamrnd(a2+c2,1/(lambda(i)+d2),1,1);
   % Now get the probabilities for k.
   for j = 1:n
      like(j) = exp((lambda(i)-theta(i))*j)*(theta(i)/lambda(i))^sum(y(1:j));
   end
   like = like/sum(like);
   % Now sample the variate for k
   k(i) = cssample(1:n,like,1);
end
% convergence plots
subplot(3,1,1)
plot(theta)
ylabel('Theta')
subplot(3,1,2)
plot(lambda)
ylabel('Lambda')
subplot(3,1,3)
plot(k)
ylabel('Change-Point k')


%-----------------------------------------------------------
% Symbolic computation with matlab
%----------------------------------------------------------
x = sym('x')
a = sym('alpha')
rho = sym('(1 + sqrt(5))/2')
f = rho^2 - rho - 1


syms a b c x;
f = sym('a*x^2 + b*x + c');


syms x;
f = sin(5*x);
diff(f)
g = exp(x)*cos(x);
diff(g,2) % take the second derivative of g
c = sym('5');
diff(c)
syms s t  % Derivatives of Expressions with Several Variables
f = sin(s*t);
diff(f,t) 


syms h n x
limit((cos(x+h) - cos(x))/h,h,0)
limit( (1 + x/n)^n,n,inf )
limit(x/abs(x),x,0,'left') %One-Sided Limits


syms a positive;
syms x;
f = exp(-a*x^2);
int(f,x,-inf,inf)

syms x k
s1 = symsum(1/k^2,1,inf)
s2 = symsum(x^k,k,0,inf)


syms a b c x
S = a*x^2 + b*x + c;
solve(S)
b = solve(S,b) % solve for b



%------------------------------------------------
%    Optimization
%------------------------------------------------
% Step 1: Create M file
function f = objfun(x)
f = exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);

% Step 2: 
x0 = [-1,1];    % Starting guess
options = optimset('LargeScale','off');
[x,fval,exitflag,output] = fminunc(@objfun,x0,options)


% Step 1: Create M file
function [c, ceq] = confun(x)
% Nonlinear inequality constraints
c = [1.5 + x(1)*x(2) - x(1) - x(2);
     -x(1)*x(2) - 10];
% Nonlinear equality constraints
ceq = [];
% Step 2:
x0 = [-1,1];     % Make a starting guess at the solution
options = optimset('LargeScale','off');
[x, fval] = ...
fmincon(@objfun,x0,[],[],[],[],[],[],@confun,options) 
