【求助】hill图 平均超额图 用matlab怎么做啊

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

请问 matlab 里 mean excess 和hill图怎么画呢  哪位高手指点一下我这只菜鸟吧~

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：MATLAB atlab matla Hill Hil matlab

1. function [c,uband,lband]=hillplot(data,param,opt,ci,plotlen)
   
2. %Plots the Hill estimate of the tail index of heavy-tailed data
   
3. %
   
4. %  USAGE: [out,uband,lband]=hillplot(data,param,opt,ci,plotlen)
   
5. %
   
6. %   data: Data vector
   
7. %  param: whether 'alpha' or 'xi' (1/alpha) should be plotted
   
8. %    opt: Option for plotting with respect to threshold ('t') or number of exceedances  
   
9. %         ('n').
   
10. %     ci: Confidence interval
    
11. %plotlen: Optional argument that determines plotting range of hill estimate starting from
    
12. %         first order statistics. It can be entered as an integer (no quotes) or as a  
    
13. %         percentage as '%5' including the quotes.
    
14. %
    
15. %    out: Parameter plotted
    
16. %  uband: Upper confidence interval
    
17. %  lband: Lower confidence interval
    
18.          
    
19. 
    
20. data=surecol(data);
    
21. ordered=flipud(sort(data));
    
22. ordered=ordered(ordered>0);
    
23. n=length(ordered);
    
24. loggs=log(ordered);
    
25. avesumlog=cumsum(loggs)./(1:n)';
    
26. diffi=avesumlog-loggs;
    
27. diffs=[NaN; diffi(2:n)];
    
28. if strcmp(param,'alpha'),
    
29.     diffs=1./diffs;
    
30. elseif strcmp(param,'xi'),
    
31.     diffs=diffs;
    
32. else disp('Parameter should be one of ''xi'' or ''alpha'' ');
    
33.     return
    
34. end
    
35.   
    
36. ses=diffs./sqrt(1:n)';
    
37. x=1:n;
    
38. y=diffs;
    
39. last=n;
    
40. if nargin==5,
    
41.    if isstr(plotlen)==1,
    
42.                 perc=str2num(plotlen(2:end))/100;
    
43.                 last=floor(perc*n);
    
44.    else
    
45.       last=plotlen;
    
46.    end
    
47. end
    
48. if strcmp(opt,'n')  
    
49. plot(x(1:last),y(1:last))
    
50. c=y;
    
51. qq=norminv(1-(1-ci)/2);
    
52. uband=y+ses(x)*qq;
    
53. lband=y-ses(x)*qq;
    
54. hold on
    
55. plot(x(1:last),uband(1:last),'r:');
    
56. plot(x(1:last),lband(1:last),'r:');
    
57. xlabel('Order Statistics')
    
58. ul=[uband;lband];
    
59. scx=range(x(1:last))/20;
    
60. scy=range(ul(1:last,:))/20;
    
61. axis([min(x(1:last))-scx max(x(1:last))+scx min(ul)-scy max(ul(1:last,:))+scy]);
    
62. 
    
63. 
    
64. elseif strcmp(opt,'t'),
    
65.    
    
66.    threshs=findthresh(data,1:n);
    
67.    plot(threshs(1:last),y(1:last));
    
68.         c=y;
    
69.         qq=norminv(1-(1-ci)/2);
    
70.         uband=y+ses(x)*qq;
    
71.         lband=y-ses(x)*qq;
    
72.         hold on
    
73.         plot(threshs(1:last),uband(1:last),'r:');
    
74.         plot(threshs(1:last),lband(1:last),'r:');
    
75. 
    
76.         xlabel('Thresholds')
    
77.    ul=[uband;lband];
    
78.    scx=range(threshs(1:last))/20;
    
79.    scy=range(ul(1:last,:))/20;
    
80.         axis([min(threshs(1:last))-scx max(threshs(1:last))+scx min(ul)-scy max(ul(1:last,:))+scy]);
    
81.    set(gca,'XDir','Reverse');
    
82. else  
    
83.    error('Last argument should be ''n'' or ''t''')
    
84. end
    
85. 
    
86. if strcmp(param,'alpha'),
    
87.     ylabel('alpha');
    
88. end
    
89. if strcmp(param,'xi'),
    
90.     ylabel('xi');
    
91. end
    
92. title('Hillplot')
    
93. hold off

复制代码

两个子函数

1. function c=surecol(data)
   
2. %Assures that data is column vector
   
3. [m,n]=size(data);
   
4. if n>m,
   
5.    c=data';
   
6. else
   
7.    c=data;
   
8. end
   
9. 
   
10. function out=findthresh(data,ne)
    
11. %Finds a threshold so that a given number of extremes lie above  
    
12. %
    
13. %USAGE: out=findthresh(data,ne)
    
14. %
    
15. % data: Data vector
    
16. %   ne: Number of extremes to lie above.
    
17. %       It can be supplied as an integer or as a percentage. If as a percentage it should
    
18. %       be entered between quotes i.e. '%5' (5% of the data lie above the threshold found).  
    
19. %
    
20. %  out: threshold calculated
    
21. data=surecol(data);
    
22. if isstr(ne)==1,
    
23.         len=length(data);
    
24.         perc=str2num(ne(2:end))/100;
    
25.         ne=floor(perc*len);
    
26. end
    
27. data=flipud(sort(data));
    
28. thresholds=flipud(unique(data));
    
29. lne=length(ne);
    
30. for i=1:lne,
    
31. index=data(ne(i))==thresholds;
    
32. index=find(index);
    
33. ind(i)=min(index+1,length(thresholds));
    
34. end
    
35. out=thresholds(ind);

复制代码