语法

plot(X,Y)

plot(X,Y,LineSpec)

plot(X1,Y1,...,Xn,Yn)

plot(X1,Y1,LineSpec1,...,Xn,Yn,LineSpecn)

plot(Y)

plot(Y,LineSpec)

plot(___,Name,Value)

plot(ax,___)

h = plot(___)

说明

示例

plot(X,Y)创建Y中数据对X二维线图中对应值。

语法

scatter3(X,Y,Z)

scatter3(X,Y,Z,S)

scatter3(X,Y,Z,S,C)

scatter3(___,'filled')

scatter3(___,markertype)

scatter3(___,Name,Value)

scatter3(ax,___)

h =scatter3(___)

说明

示例

scatter3(X,Y,Z)在向量X、Y和Z圆圈显示在指定位置。

示例

scatter3(X,Y,Z,S)使用S画出指定大小的每个圆圈。画一个大小相等的圆圈，请画S指定为标量。为了绘制每个具有特定大小的圆，请将其绘制S指定为向量。

示例

scatter3(X,Y,Z,S,C)使用C画出每个圆圈的指定颜色。

• 如果C是 RGB 三元组，或包含颜色名称的字符向量或字符串，用指定的颜色绘制所有圆圈。

• 如果C是三列矩阵之一C行数等于X、Y和Z的长度，则C指定相应圆圈的每行 RGB 颜色值。

• 如果C是长度与X、Y和Z相同长度的向量，C值线性映射到当前颜色图中的颜色。

示例

scatter3(___,'filled')在前面的语法中使用任何输入参数组合来填充这些圆。

示例

scatter3(___,markertype)指定标记类型。

示例

scatter3(___,Name,Value)修改散点图，使用一个或多个名称-值对组参数。

示例

scatter3(ax,___)绘制图形ax在指定的坐标区域而不是当前坐标区 (gca) 中。选项 ax 可以位于前面的语法中的任何输入参数组合之前。

示例

h = scatter3(___) 返回 Scatter 对象。在创建散点图后，可使用 h 修改其属性。

双箱模型 优化 PCA 箱式 pcolor

箱式 不优化

clc;clear all; close all; V1 = 3e16   ; % volume of upper ocean (m^3) V2 = 100e16 ; % volume of deeper ocean (m^3) Fr = 3e13   ; % river inputs (m^3/yr) Fo = 6e14   ; % turnover rate (m^3/yr) xr = 1.3    ; % concentration of P in rivers (mmol/m3) Tau = 50  ; % 100 year residence time of P in the surface; Pu = 0 ; % initial concentration in the upper ocean  Pd = 0 ; % initial concentration in the deeper ocean  Prodtvy = Pu*V1/Tau ; % Production rate  ReminEff = 0.97 ;     % remineralization ratio 

dPudt = (Fr*xr - Fo*Pu + Fo*Pd - Prodtvy)/V1 ; dPddt = (Fo*Pu - Fo*Pd + ReminEff*Prodtvy)/V2 ;

n  = 1  ; dt = 10 ; figure(1) % start time stepping  for ji = 1:dt:1e7     n = n + 1;     set(0,'CurrentFigure',1);     Pu = Pu + dPudt * dt ;     Pd = Pd + dPddt * dt ;     Prodtvy = Pu*V1/Tau ; % Production     dPudt = (Fr*xr - Fo*Pu + Fo*Pd - Prodtvy)/V1;     dPddt = (Fo*Pu - Fo*Pd + ReminEff*Prodtvy)/V2;     % plot a dot every 5000 iterations.     if mod(n, 5000) == 0         hp = plot(ji,Pu,'r.:',ji,Pd,'G.:');         set(hp,'MarkerSize',16)         set(gca,'FontSize',16)         xlabel('time (years)')         ylabel('concentration of P (mmol/m3)')         legend('Pu','Pd','Location','best')         drawnow; hold on;     end  end

箱式 优化

addpath('~/Dropbox/myfunc') par.Tau = 100 ;  % 100 year residence time of P in the surface; par.V1 = 3e16 ;  % volume of upper ocean (m^3) par.V2 = 100e16; % volume of deeper ocean (m^3) par.Fr = 3e13;   % river inputs (m^3/yr) par.xr = 1.5 ;   % concentration of P in rivers (mmol/m3) par.Fo = 6e14;   % turnover rate (m^3/yr) par.ReminEff = 0.99 ; % remineralization ratio  par.Puobs = 0.5; par.Pdobs = 2.5; x0 = log([par.xr, par.Tau, par.ReminEff]) ;% 防止算出负值 最后要换算回去 myfun = @(x) neglogpost2(x, par);

options = optimoptions(@fminunc                  , ...                        'Algorithm','quasi-newton', ...%拟牛顿法                        'GradObj','off'            , ...                        'Hessian','off'            , ...                        'Display','iter'          , ...                        'MaxFunEvals',2000        , ...                        'MaxIter',2000            , ...                        'TolX',1e-7               , ...                        'TolFun',1e-7             , ...                        'DerivativeCheck','off'   , ...                        'FinDiffType','central'   , ...                        'PrecondBandWidth',Inf)   ; % nip = length(x0); [xhat,fval,exitflag] = fminunc(myfun,x0,options);%求无约束多变量函数最小值 按option中的选项 [f,mod] = neglogpost2(xhat,par); % save results  save('pho_xhat', 'xhat','f') fprintf('-------------- end! ---------------\n');  

function [f,mod] = neglogpost2(x, par)     [Pu,Pd] = fpho(par,x);     mod = [Pu; Pd] ;     obs = [par.Puobs; par.Pdobs];     V = [par.V1; par.V2];     W = d0(V./sum(V));%求权矩阵 且存储为稀疏矩阵     err = (mod - obs) ;     f = err'*W*err ;%最小二乘法 end

function [Pu,Pd] = fpho(par,x)     V1 = par.V1 ;     V2 = par.V2 ;      Fr = par.Fr ; % m^3 per year     Fo = par.Fo ; % m^3 per year

    xr = exp(x(1)) ; % mmol/m3     Tau = exp(x(2));     ReminEff = exp(x(3)) ; % ratio to sediment

    D = [Fo+ReminEff*V1/Tau, -Fo; -Fo-V1/Tau, Fo]; %方程组左边     rhs = [0; -Fr*xr];%方程组右边          X = inv(D)*rhs;%解方程组     Pu = X(1);     Pd = X(2);     Prodtvy = Pu*V1/Tau ; % Production           rhs = [0; -Fr*xr];          X = inv(D)*rhs;     Pu = X(1);     Pd = X(2); end 

function X=d0(v) X=spdiags(v(:),0,length(v(:)),length(v(:))); end

内置函数PCA

% using built-in function; [data,txt] = xlsread('scores'); [nx,ny] = size(data); ave  = mean(data); Dstd = std(data) ; demean = data - repmat(ave,[nx,1]); normD  = demean./repmat(Dstd,[nx,1]); % do the PCA using Matlab built-in function PCA; [coeff, score, latent, tsquared, explaned, mu] = pca(norm, 'Algorithm', 'eig'); pc1 = coeff(:,1); pc2 = coeff(:,2); pc3 = coeff(:,3); figure(1) yline(0,'r-') hold on h1 = plot(pc1,'ro-','linewidth',1.5,'MarkerSize',3); hold on h2 = plot(pc2,'k*-.','linewidth',1.5,'MarkerSize',3); hold on h3 = plot(pc3,'b^:','linewidth',1.5,'MarkerSize',3); hold off grid on legend([h1,h2,h3],{'pc1','pc2','pc3'}) xticks([1:13]) header = txt(1,:); kind  = txt(2:end,end); xticklabels({'Chinese','Math','English','Physi.','Chem.','Bio.',...              'Politics','History','Geography','Music','Physical','Arts','Tech'})

% the original data are grouped based on PCA score; figure(2) PC1 = score(:,1); PC2 = score(:,2); PC3 = score(:,3); labels = kind; scatter3(PC1,PC2,PC3,'ro') text(PC1,PC2,PC3,labels) xlabel('PC1') ylabel('PC2') zlabel('PC3')

% demonstrate how the score is calculated; figure(3) scores = demean*coeff; score1 = scores(:,1); score2 = scores(:,2); score3 = scores(:,3);

scatter3(score1,score2,score3,'ro') text(score1,score2,score3,labels) xlabel('PC1') ylabel('PC2') zlabel('PC3')

自己写PCA

[data,txt] = xlsread('scores'); [nx,ny] = size(data); ave  = mean(data); demean = data - repmat(ave,[nx,1]); % raw data using corrcoef ; R = corrcoef(data);%相关系数 [V, Lambda] = eigsort(R);%特征值分解+ 把特征值矩阵的diag按从大到小排列，且返回对应的列号+把Lambda变成一个对角矩阵 PoV = diag(Lambda)/trace(Lambda);%求每个主成分的贡献 [PoV, cumsum(PoV)] pc1 = V(:,1); pc2 = V(:,2); pc3 = V(:,3); figure(1) yline(0,'r-') hold on h1 = plot(pc1,'ro-','linewidth',1.5,'MarkerSize',3); hold on h2 = plot(pc2,'k*-.','linewidth',1.5,'MarkerSize',3); hold on h3 = plot(pc3,'b^:','linewidth',1.5,'MarkerSize',3); hold off legend([h1,h2,h3],{'pc1','pc2','pc3'}) xticks([1:13]) header = txt(1,:); kind  = txt(2:end,end); xticklabels({'Chinese','Math','English','Physi.','Chem.','Bio.',...              'Politics','History','Geography','Music','Physical','Arts','Tech'})

figure(2) scores = demean*V; score1 = scores(:,1); score2 = scores(:,2); score3 = scores(:,3); labels = kind; scatter3(score1,score2,score3,'ro') text(score1,score2,score3,labels) xlabel('PC1') ylabel('PC2') zlabel('PC3') title('PCA based on data*V')

% to demonstrate the without normalization, the seperation is worse. figure(3) R = cov(data); [V, Lambda] = eigsort(R); PoV = diag(Lambda)./trace(Lambda); [PoV, cumsum(PoV)] scores = data*V; score1 = scores(:,1); score2 = scores(:,2); score3 = scores(:,3); labels = kind; scatter3(score1,score2,score3,'ro') text(score1,score2,score3,labels) xlabel('PC1') ylabel('PC2') zlabel('PC3') title('PCA based on cov(data)') [coeff,sc,latent,tsquared,explained,mu] = pca(demean);

function [V, Lambda] = eigsort(R)     [V, Lambda] = eig(R);     [lambda, ilambda] = sort(diag(Lambda),'descend');     Lambda = diag(lambda);     V = V(:,ilambda); end  

箱式 demo

clear all; close all; clc % A = rand(500,1); A = normrnd(5,1,100,1);%正态随机数 B = normrnd(5,1,100,1); Y = quantile(A,[0.25,0.5,0.75]); %quantile 分位数

w = 1.0 ;  % boxplot([A,B], 'Whisker',w) boxplot(A, 'Whisker',w)% （认为是）合理值的范围 w越大则两个横线向上下延伸越多 q1 = Y(1); q2 = Y(2); q3 = Y(3);

text(1.1,q1,'25% percentile') text(1.1,q2,'median') text(1.1,q3,'75% percentile')

uwhickler = q3 + w * (q3 - q1); dwhickler = q1 - w * (q3 - q1); text(1.1, uwhickler, 'q3 + w * (q3 - q1)') text(1.1, dwhickler, 'q1 - w * (q3 - q1)')

drawArrow = @(x,y) quiver( x(1),y(1),x(2)-x(1),y(2)-y(1),0 );

箱式图和阴影

load MLD1850.mat MLD lat lon [Lat,Lon] = meshgrid(lat,lon) ; ibig = find(MLD(:) >250); MLD(ibig) = nan; data = [Lat(:), MLD(:)] ; ikp = find(~isnan(sum(data,2)) & ~isinf(sum(data,2)));

lat_kp = Lat(ikp); mld_kp = MLD(ikp); figure(1) edges = linspace(min(lat_kp), max(lat_kp), 18); % edges = min(lat_kp):10:max(lat_kp)+5; [N,ibin] = histc(lat_kp,edges); boxplot(mld_kp,ibin) xticks(1:18)

xticklabels({'89.5S', '78.97S', '68.44S', '57.91S', '47.38S', '36.85S', ...              '26.32S', '15.79S', '5.26S', '5.26N', '15.79N', '26.32N', ...              '36.85N', '47.38N', '57.91N', '68.44N', '78.97N', '89.50N'});

for ji = 1:length(N)      ikp = find(ibin == ji);     ave(ji) = mean(mld_kp(ikp));     stdn(ji) = std(mld_kp(ikp)); end

figure(2) inbetweenx = [1:length(edges), fliplr(1:length(edges))]; inbetweeny = [ave + stdn, fliplr(ave - stdn)]; fill(inbetweenx, inbetweeny,[0.9,0.9,0.9]); hold on plot(ave,':r^') hold off xticks(1:18) xticklabels({'89.5S', '78.97S', '68.44S', '57.91S', '47.38S', '36.85S', ...              '26.32S', '15.79S', '5.26S', '5.26N', '15.79N', '26.32N', ...              '36.85N', '47.38N', '57.91N', '68.44N', '78.97N', '89.50N'});  

mk Anomaly

addpath('~/Dropbox/myfunc')

fname = 'mlotst_Omon_CESM2-WACCM_historical_r2i1p1f1_gr_185001-201412.nc';

lat = ncread(fname,'lat'); lon = ncread(fname,'lon'); time = ncread(fname,'time'); MLD = ncread(fname,'mlotst',[1,1,1],[Inf,Inf,12]);

MLD_Jan = MLD(:,:,1); MLD_Jul = MLD(:,:,7);

anom = permute(MLD_Jul - MLD_Jan,[2,1]) ; pcolor(anom);colorbar;shading flat; set(gca,'color','black')  %colormap(darkb2r(-500,500)); colorbar 

pcolor作图和动画

clc; clear; fname='mlotst_Omon_CESM2-WACCM_historical_r2i1p1f1_gr_185001-201412.nc'; ncdisp(fname); %ocean data regridded from native gx1v7 displaced pole grid (384x320 %latxlon) 如果把极点放在真正南北极点模型会产生很大的误差 lat=ncread(fname,'lat');%lat=ncread(fname,'lat',[1,1,1],[Inf,Inf,12]) 可以取出第一年的Data  [long,lat,month] lon=ncread(fname,'lon'); time=ncread(fname,'time'); MLD=ncread(fname,'mlotst'); MLD1=MLD(:,:,1); figure(1); pcolor(MLD1);colorbar;shading flat;%pcolor: pixel color %shading faceted：默认模式，在曲面或图形对象上叠加黑色的网格线； % shading flat：是在shading faceted的基础上去掉图上的网格线； % shading interp：对曲面或图形对象的颜色着色进行色彩的插值处理，使色彩平滑过渡 ； caxis([0 400]);

MLDp=permute(MLD1,[2,1,3]);%permute 指的是把矩阵中不同维度的数据进行交换 这里是交换1,2维度 figure(2); pcolor(MLDp);colorbar;shading flat; caxis([0 400]);% 把colorbar的一部分归0

MLD=permute(MLD,[2,1,3]); [x,y,z]=size(MLD); figure(3); % vidfile=VideoWriter('MLD_movie.mp4','MPEG-4'); % open(vidfile); % for ji=1:20 %     iMLD=MLD(:,:,ji); %     pcolor(iMLD);colorbar;shading flat; %     caxis([0 400]);%否则colorbar的数值会上下横跳 %     drawnow; % %     pause(0.1);%去掉括号里的数字可以手动切图(用回车） %     F(ji)=getframe(gcf);%捕获坐标区或图窗作为影片帧 %     writeVideo(vidfile,F(ji)); % end % close(vidfile); %save MLD1850 MLD lat lon;% MLD1850是文件名 MLD=MLD(:,:,1); [Lat,Lon]=meshgrid(lat,lon); ibig=find(MLD(:)>250); MLD(ibig)=nan;

data=[Lat(:),MLD(:)]; ikp=find(~isnan(sum(data,2))&~isinf(sum(data,2)));%找到不是nan和inf的值

lat_kp=Lat(ikp); mld_kp=MLD(ikp); figure (1) edges=linspace(min(lat_kp),max(lat_kp),18); [N,ibin]=histc(lat_kp,edges); get(gca) boxplot(mld_kp,ibin)

xticks(1:18) xticklabels({'89.50S','78.97S','68.44S','57.91S','47.38S','36.85S','26.32S',...     '15.79S','5.26S','5.25N','15.79N','26.32N','36.85N','47.38N','57.91N',...     '68.44N','78.97N','89.50N'}) for ji=1:length(N)     ikp=find(ibin==ji);     ave(ji)=mean(mld_kp(ikp));     stdn(ji)=std(mld_kp(ikp)); end figure(2) inbetweenx = [1:length(edges), fliplr(1:length(edges))]; inbetweeny = [ave + stdn, fliplr(ave - stdn)]; fill(inbetweenx, inbetweeny,[0.9,0.9,0.9]); hold on plot(ave,':r^') hold off xticks(1:18) xticklabels({'89.5S', '78.97S', '68.44S', '57.91S', '47.38S', '36.85S', ...              '26.32S', '15.79S', '5.26S', '5.26N', '15.79N', '26.32N', ...              '36.85N', '47.38N', '57.91N', '68.44N', '78.97N', '89.50N'}); %load时同理