分类目录归档:SAS应用

批量拟合并提取结果参数的SAS程序

提取某个文件夹下.fit文件,然后批量拟合,并提取拟合后的参数到一个数据集中。
效率很高。
 
filename f_list pipe "dir D:\three\dat /b";
data mylist;
infile f_list;
input file :$14.;
run;
data mylist;
  set mylist;
where file ? ‘fit’;
call symput(‘n’,_n_);
put file _n_;
run;
%macro getn;
   data _null_;
      set mylist;
   %do i=1 %to &n.;
          if _n_=&i. then call symput("fn&i.",file);
 %end;
%mend;
%getn;
%macro inverse;
    
  %do i=1 %to &n.;
      data ex;
       infile "D:\three\dat\%CMPRES(&&fn&i.)" firstobs=2 ;
    input x y;
    run;

    ods listing close;
    ods output ParameterEstimates = paras ;
     proc nlin data=ex ;
     parms a1=100  a2=100 a3=100  ta1=10 ta2=10 ta3=10   e=10;
     model y=a1*exp(-x/ta1)+a2*exp(-x/ta2)+a3*exp(-x/ta3)+e;
    run;
    ods listing;

          data inver (keep= name parameter estimate);
               set paras;
      name=scan("&&fn&i.",1);
      run;
                      
    proc append base=inversx data=inver ;                                           
    run;
       %end;
%mend;
%inverse
proc sort data=inversx;
by  name;
run;
 proc transpose data=inversx out=final (drop=_name_);
  var  estimate ;
  id parameter;
  by name;
run;

data ex (drop=aa1-aa3 xa1-xa3 i j);
   set final;
        aa1=ta1;aa2=ta2;aa3=ta3; 
xa1=a1;xa2=a2;xa3=a3;
   call sortn(of ta:);
   array arr_a{3} a1-a3;
   array arr_aa{3} aa:;
   array arr_xa{3} xa:;
    do i=1 to 3;
      do j= 1 to 3; 
       if arr_aa{j}=ordinal(i,of ta:) then arr_a{i}=arr_xa{j};
       end;
    end; 
run;
proc print data=ex;
run;

SAS函数&CALL例程介绍与实例精选

SAS函数&CALL例程介绍与实例精选

AN INTRODUCTION OF SAS FUNCTION & CALL ROUTINE AND SELECTED SAMPLE

Bruce Shao

摘要:
SAS提供了大量的内置函数,并且在数量和内容上不断的发展改进(注:本文所述SAS函数指SAS/BASE模块内的函数)。函数数量的变化非常显著,数量从SAS6.12的二百多个增加到SAS9.2的近五百个。这篇文章前半部分将简单介绍SAS函数的形式和种类,以及SAS系统版本更新过程中数量的变化,对于从SAS9.13升级到SAS9.2后内置函数的数量功能和类型上的变化做详细的介绍,特别会介绍到SAS9.2系统中用户可以自定义函数的情况;后半部分将精选典型类型的常用函数来做实例,你将会在这些例子中得到启示。一旦你学会使用合适的函数,将会使你的工作效率加倍。希望本篇文章能引起你对SAS系统提供的有用工具——SAS函数的关注。
正文:
SAS/BASE提供的函数有两种形式:一种是“FUNCTION”(后面用“函数”),可以进行计算或者系统操作并且会返回值;另一种是“CALL routine”(后面用“CALL例程”),用来改变变量的值或者执行其他的系统函数,但不返回值。SAS函数可以非常方便的用于DATA步中,WHERE子句和SQL查询语句中;CALL例程一般用于DATA步中。SAS提供了种类繁多的函数和CALL例程,并且一直在变化发展。
表1中列出了从SAS6.12版本到最新的SAS9.2版本语言参考:字典中函数和CALL例程的个数,图1中描述了函数和CALL例程总数量的变化趋势图。

表1 各SAS版本号中函数和例程的个数

图1 SAS函数和Call例程总数量的变化趋势

SAS/BASE有几类常用的函数:字符函数、字符串匹配函数、数学函数、日期和时间函数、统计函数、随机函数、金融函数、分位数函数等。最新的SAS9.2/BASE模块的语言参考:字典中所列的内置函数发生很多的变化。

在种类方面,SAS9.2语言参考:字典中内置函数,与有26类内置函数的SAS9.13比较,SAS9.2减少了程序响应度量函数(ARM)、货币折算函数(Currency Conversion)、双字节函数(DBCS)等函数种类,新增了算术函数(Arithmetic)、组合函数(Combinatorial)、距离函数(Distance)、财务函数(Finance)、数字函数(Numeric)、搜索函数(Search)、排序函数(Sort),共有30类函数。

在数量及功能方面,SAS9.2语言参考:字典中内置函数比SAS9.13增加了44个函数(包括Call例程),并且有些类型的函数如外部文件函数(External Files)中已有的函数如DOPEN、FOPEN、FEXIST,FILENAME,FILEREF,MOPEN,PATHNAME等函数参数说明被改进;有些函数是从其他模块转过来的,如SAS/ETS的INTCINDEX,INTCYCLE等等; SAS High-Performance Forecasting的HOLIDAY和NWKDOM函数等等; 有些已存在的函数和CALL例程得到增强,如函数CALL POKE, CALL PLKELONG, CALL SCAN, DATDIF, FSEP, INDEX, LAG, SCAN,ZIPSTATE等;PRX函数和例程取代了RX函数和例程,函数RXMATCH, RXPARSE,例程RXCHANGE CALL,RXFREE CALL,REXSUBSTR CALL将不再存在于SAS9.2语言参考:字典中。另外功能强大的SCAN函数和CALL SCAN例程取代了SCANQ函数和CALL SCANQ例程。

SAS9.2系统中函数使用发生的巨大变革是用户可以自定义函数。这样对于用户来说用自定义函数可以取代部分宏语言,模块化编程更容易。与函数比较起来,SAS宏有着不足的地方: 如1,写宏要非常小心谨慎,特别是符号”&”;2,宏可读性差,3自定义函数可以对某些变量其保护作用,不会乱窜现象 4自定义函数是程序更容易模块化,自定义函数可以保存起来,同部门直接可以共享。
上面简单的介绍了SAS函数的发展变化情况,和新版本SAS9.2中函数的比较详细的变化,下面将通过一些实例讨论一下SAS9版本中的一些常用的函数的使用比较和SAS9.2自定义函数的方法。
例一:字符处理函数
首字母大写的代码实现
这是字符处理经常会遇到的情况,下面给出三种通过函数的实现方法:
方法一:使用函数CAT,LOWCASE,UPCASE,SUBSTR

DATA ex;
INFORMAT first last $30.;
input first last;
first=CAT(UPCASE(SUBSTR(first,1,1)), LOWCASE(SUBSTR(first,2)));
last =CAT(UPCASE(SUBSTR(last,1,1)), LOWCASE(SUBSTR(last,2)));
cards;
ronald cODy
THomaS eDISON
albert einstein
;
run;

方法二:使用函数LOWCASE, UPCASE, SUBSTR

DATA ex;
INFORMAT first last $30.;
INPUT first last;
first = LOWCASE(first);
LAST = LOWCASE(last);
SUBSTR(first,1,1) = UPCASE(SUBSTR(first,1,1));
SUBSTR(last,1,1) = UPCASE(SUBSTR(last,1,1));
DATALINES;
ronald cODy
THomaS eDISON
albert einstein
;

方法三:使用函数PROPCASE

DATA ex;
INFORMAT first last $30.;
input first last;
first=PROPCASE(first);
last =PROPCASE(last);
cards;
ronald cODy
THomaS eDISON
albert einstein
;
run;

前面两个例子使用几个函数的组合的代码可以实现这个要求。第三个例子使用函数PROPCASE可以实现首字母大写转变的功能,另外有 ”delimiter” 可选项,可以满足字符串中出现其他特殊符号后要求字母大写的转变功能。如name=PROPCASE(“jeans*west”,’*’),可得到name的值为“Jeans*West”。
例二 数值处理函数
比赛计分计算,计算去掉一个最高分和一个最低分后的总分。
方法一:使用函数MIN,MAX,SUM

DATA ex;
INPUT Name:$ q1-q10 ;
ARRAY QQ q1-q10;
TOTAL=SUM (of q1-q10)-MIN(of q1-q10)-MAX(of q1-q10);
CARDS;
Tom 90 80 81 79 78 76 84 82 84 50
Lucy 89 78 79 77 76 75 81 79 83 74
Jack 47 78 80 81 82 82 78 85 77 85
;
RUN;

方法二:使用函数ORDINAL,SUM

DATA ex;
INPUT Name:$ q1-q10 ;
ARRAY QQ q1-q10;
TOTAL=SUM (of q1-q10)-ORDINAL(1,of q1-q10)-ORDINAL(10,of q1-q10);
CARDS;
Tom 90 80 81 79 78 76 84 82 84 50
Lucy 89 78 79 77 76 75 81 79 83 74
Jack 47 78 80 81 82 82 78 85 77 85
;
RUN;

方法三:使用函数SMALLEST, LARGEST,SUM

DATA ex;
INPUT Name:$ q1-q10 ;
ARRAY QQ q1-q10;
TOTAL=SUM (of q1-q10)-SMALLEST (1, of q1-q10)-LARGEST (1,of q1-q10);
CARDS;
Tom 90 80 81 79 78 76 84 82 84 50
Lucy 89 78 79 77 76 75 81 79 83 74
Jack 47 78 80 81 82 82 78 85 77 85
;
RUN;

方法四:使用例程CALL SORTN,SUM

DATA ex;
INPUT Name:$ q1-q10 ;
ARRAY QQ q1-q10;
CALL SORTN (of qq [*]) ;
TOTAL=SUM (of q2-q9);
CARDS;
Tom 90 80 81 79 78 76 84 82 84 50
Lucy 89 78 79 77 76 75 81 79 83 74
Jack 47 78 80 81 82 82 78 85 77 85
;
RUN;

上面给出了四种不同的函数方法来实现这个计算,但是这几种函数用法和功能并不是完全一样的,如对缺失值的处理不一样,当观察值中出现缺失值时,MIN,MAX将忽略缺失值后计算、CALL SORTN则会把缺失值排为最小值的位置后计算,而SMALLEST, LARGEST和ORDINAL则直接生成缺失值。当计算要求改变如去掉最大两个值和最小两个值后的总分,函数MIN,MAX就难于实现,而SMALLEST, LARGEST和ORDINAL可以实现,但是不如例程CALL SORTN的代码简洁。
例三:日期处理函数
SAS提供了很多日期处理的函数如SECOND,MINUTE,HOUR,DAY,MONTH,YEAR;DATEPART,TIMEPART;INTNX,INTCK,YRDIF等等,但是实际应用中还是不够用。SAS9.2/BASE模块的语言参考:字典中新增了2个内置函数HOLIDAY和NWKDOM,分别用来返回某一年指定假期的具体日期和返回指定年月某一周第几天的具体日期。
例3:2009年的母亲节是哪一天?

data _null_;
MOTHERS = holiday('MOTHERS', 2009);
format MOTHERS date9.;
put MOTHERS=;
run;

SAS9.2的新函数HOLIDAY目前暂时提供22个美国和加拿大的节假日选项。
例四:自定义函数
在SAS9.2系统中,我们可以用C,C++或者SAS语言来自定义函数,分别提供了两种编译语句:用PROC PROTO语句来编译C或C++语言定义的函数;用PROC FCMP语句来编译SAS语言定义的函数,这样就可以建立自己的函数库。SAS本身已经提供了一些专业的函数库如统计、金融函数库,以后来自不同行业的用户可以根据自己的需要建立各自专业的函数库共享使用。
例4:在临床医学中,计算一个study­_ day是个常规任务。

/*自定义函数*/
proc fcmp outlib=sasuser.funcs.trial;
function study_day(start, event);
n = event - start;
if n >= 0 then
return(n + 1);
else
return(n);
endsub;
/*调用函数*/
options cmplib=sasuser.funcs;
data _null_;
start = '15Feb2006'd;
today = '27Mar2006'd;
sd = study_day(start, today);
put sd=;
run;

参考资料6陈述了实现相同目的的宏在使用中将会出现意外问题的讨论,具体见6。
小结:函数是SAS系统提供的重要工具之一,数量、内容和功能上不断的得到完善,并且在新版本9.2提供了一个重要的自定义函数功能。如果你是对函数不熟悉,建议你开始学习使用SAS函数,一旦你学会使用合适的函数,将会使你的工作效率加倍;如果你是一个经验丰富的用户,建议多查看新版本中的函数,也许可以找到更简便的方法来解决以前同样的工作难题。
参考:

  1. SAS Institute. 1996. SAS system for windows. V6.12 TS020. Carey, NC: SAS Institute, Inc.
  2. SAS Institute. 2000. SAS Language Reference: Dictionary, Version8. Carey, NC: SAS Institute, Inc.
  3. Deb Cassidy. 2004. ”An Introduction to SAS Functionality” SUGI29 Montréal, Canada Paper 255-29
  4. SAS Institute Inc. 2008. SAS®9.2 Language Reference dictionary Cary, NC: SAS Institute Inc.
  5. SAS Institute Inc. 2008. What’s New in SAS® 9.2. Carey, NC: SAS Institute, Inc.
  6. Jason Secosky. 2007. A Sampler of What’s New in Base SAS® 9.2 RTSUG Cary, NC: SAS Institute Inc.
  7. Web Sites: SAS爱好者论坛 http://www.sasor.com

联系方式CONTACT INFORMATION
资料整理,难免出错,希望不吝赐教;并欢迎建议和交流,联系方式如下:
Your comments and questions are valued and encouraged. I can be contacted at:
电子信箱:sxlion8 (at) gmail.com

简单线性规划问题SAS编程之二

钢管长度定制问题:

假设实际生产需要303mm,251mm,202mm,151mm,107mm的钢管各100根【理论上需要一根(303+251+202+151+107)*100=101400的长钢管,这样才不会多余的边角料浪费,实际上这么长的一根钢管供货商是不便于提供的】。现供货商只能提供长度相同且长度范围为1300mm至1650mm的钢管,请问需要向订货商定制多少根和长度为多少的钢管才能保证浪费最少(即边角料的总长度最小)?

引用:http://www.mysas.net/forum/viewtopic.php?f=4&t=4768&start=10   ID:tianlai888

ps:原文中的表达不够清晰,且10根算起来不过瘾,故修改。

代码:

options nomprint;

data ex;
input  _type_$ _rhs_;
cards;
ge 100
ge 100
ge 100
ge 100
ge 100
min .
upperbd .
integer .
;
run;
%macro sxl(out=out);

proc sql;
create table &out.
  (
    numm num format = 12.       label = ‘number’,
    i  num format = 12. label = ‘length’

  );
quit;

%do i = 1300 %to 1650;

data _aa (keep=c d e f g i );
   a1=int(&i./303);
   a2=int(&i./251);
   a3=int(&i./202);
   a4=int(&i./151);
   a5=int(&i./107);
   w=&i.-107;
    do c=0 to a1;
          do d=0 to a2;
              do e=0 to a3;
                  do f=0 to a4;
                     do g=0 to a5;
                        x=303*c+251*d+202*e+151*f+107*g;
                        y=&i.-x;
                         if x>w and x<=&i. then do;
                                                     i+1;
                                                     output;
                                                  end;
                      end;
                  end;
             end;
           end;
     end;
run;

proc sort data=_aa out=_aa(drop=i);
by descending i ;
run;

data _null_;
set _aa(drop=_all_) nobs=rows;
call symput(‘nobs’,rows);
stop;
run;

data _collect(drop=i);
LENGTH _VAR_ $8.;
set _aa;
i+1;
_VAR_=catt(‘COL’,i);
run;

proc transpose data=_aa out=_aa1;
run;
data _bb(drop=i);
do i =1 to &nobs.;
    a=1;
    b=100;
    c+1;
    output;
    end;
    run;
    proc transpose data=_bb out=_bb1;
    run;

data ex1(drop=_NAME_);
set _aa1 _bb1;
run;

data ex2;
merge ex1 ex ;
run;

ods listing close;
ods output solutionsummary=_p;
proc lp data=ex2  IMAXIT=500;
run;
ods _all_ close;
ods listing ;

proc sql noprint;
select CVALUE1  
   into :numm from _p;

    data _tem;
      numm=&numm.;
      i  = &i.;
    run;

    proc sql;       
      insert into
        &out.
      select
        *
      from
        _tem;
    quit;

%end;
data &out.;
   set &out.;
   length=numm*i;
   if numm ne lag(numm) then output;
   run;

proc sort data=&out. out=&out.;
   by length ;
   run;
proc print data=&out.;
run;
%mend sxl;
%sxl;

result:

Obs            numm               i    length

  1              69            1471    101499
  2              67            1515    101505
  3              65            1562    101530
  4              72            1411    101592
  5              63            1613    101619
  6              77            1320    101640
  7              75            1360    102000
  8              78            1310    102180
  9              70            1460    102200
10              68            1509    102612
11              66            1557    102762
12              73            1408    102784
13              64            1608    102912
14              76            1357    103132
15              79            1307    103253
16              71            1457    103447
17              80            1300    104000
18              74            1406    104044

从结果中看出,当定制1471mm的钢管时,浪费最少。需要定制69根,浪费101499-101400=99mm。

小结:实际程序跑了100s,理论上猜测应该少于30s,代码还有待优化。

简单线性规划问题SAS编程

问题:一批铁丝,数量不限,每一根铁丝长 1650 mm ,现需要截取下列品种、长度、数量的铁丝:
品种 长度/mm 需要数量/根
a1 303 10
a2 251 10
a3 202 10
a4 151 10
a5 107 10
注:铁丝的长度、品种的个数、品种的长度、需要数量都是变量。
怎样计算,在达到需要截取的品种、长度、数量的前提下,使用的铁丝根数最少,浪费的最小?
这里有两个要求:1,使用1650mm铁丝的根数最少
                        2,有一根剩下的最长
代码1:
 
data sxl;
input _ROW_$ x1 x2 x3 x4 x5 x6 _TYPE_$ _RHS_;
cards;
object 0 0 0 0 0 1 MIN .
proc1 303 251 202 151 107 -1650 lE 0
proc2 1 0 0 0 0 0 GE 10
proc3 0 1 0 0 0 0 GE 10
proc4 0 0 1 0 0 0 GE 10
proc5 0 0 0 1 0 0 GE 10
proc6 0 0 0 0 1 0 GE 10
bound 100 100 100 100 100 100 UPPERBD .
inbd 1 2 3 4 5 6 INTEGER .
;
proc print;
proc lp data=sxl;
run;

RESULT:

Variable Reduced
ol Name Status Type Price Activity Cost

1 x1 ALTER INTEGER 0 11 0
2 x2 ALTER INTEGER 0 11 0
3 x3 ALTER INTEGER 0 12 0
4 x4 ALTER INTEGER 0 11 0
5 x5 BASIC INTEGER 0 10 0
6 x6 INTEGER 1 7 1

 
代码2:
data aa (keep=c d e f g i );
   a1=int(1650/303);
   a2=int(1650/251);
   a3=int(1650/202);
   a4=int(1650/151);
   a5=int(1650/107);
   w=1650-107;
    do c=0 to a1;
          do d=0 to a2;
              do e=0 to a3;
                  do f=0 to a4;
                     do g=0 to a5;
                        x=303*c+251*d+202*e+151*f+107*g;
                        y=1650-x;
      
                         if x>w and x<=1650 then do;
                                                     i+1;
                                                     output;
                                end;
                      end;
                  end;
             end;
           end;
     end;
run;
proc print;
run;
proc sort data=aa out=aa(drop=i);
by descending i ;
run;
data collect(drop=i);
LENGTH _VAR_ $8.;
set aa;
i+1;
 _VAR_=catt(‘COL’,i);
run;
proc print;
run;
proc transpose data=aa out=aa1;
run;
proc print;
run;
data bb(drop=i);
do i =1 to 334;
    a=1;
 b=100;
 c+1;
 output;
 end;
 run;
 proc transpose data=bb out=bb1;
 run;
data ex1(keep=col1-col334);
set aa1 bb1;
run;
data ex;
input  _type_$ _rhs_;
cards;
ge 10
ge 10
ge 10
ge 10
ge 10
min .
upperbd .
integer .
;
run;
data ex2;
merge ex1 ex ;
run;
proc print;
run;
proc lp data=ex2 primalout=p;
run;
data sxl;
set p;
where _VALUE_ NE 0;
run;
proc sort data=collect out=collect;
by _VAR_;
run;
proc print;
run;
proc sort data=sxl out=sxl;
by _VAR_;
run;
proc print;
run;
data res(keep= _VAR_    _LBOUND_  _VALUE_   c  d  e  f   g);
merge sxl collect;
by _VAR_;
run;
proc print;
WHERE _VALUE_ GE 0;
run;
 

分析比较:代码1中给出的目标不准确,要求截断后的短铁丝都为10根时,结果一致,都为7根。

不过如果要求短铁丝为100根时,就会出现代码1的结果小于代码2. 因为代码1中的目标不严格。代码2还给出了具体的截断方法。

但是截断方法是有很多种的。可以根据调整某种截断方法的次数来修改重新求解。

当然上面仅仅满足了要求1.

为了满足要求2,我们需要枚举出截断最后一根铁丝的可能性方法。

代码3:

data ex;
do  a1=1 to int(1650/303);
      x1+1;  x2=0;   x3=0;   x4=0;   x5=0;
    do a2=1 to int(1650/251);
      x2+1;  x3=0;   x4=0;   x5=0;
          do a3=1 to int(1650/202);
        x3+1;   x4=0;   x5=0;
             do a4=1 to int(1650/151);
       x4+1;    x5=0;
                  do a5=1 to int(1650/107);
                       x5+1;
                       y=1650-(303*(x1-1)+251*(x2-1)+202*(x3-1)+151*(x4-1)+107*(x5-1));
                     if y>=0 then do ;
                                      r1=x1-1;
                    r2=x2-1;
                    r3=x3-1;
                    r4=x4-1;
                    r5=x5-1;
                             output;
           end;
                   end; 
     end;
   end;
  end;
end;
run;

对于最后一根(如果是10根的话对应了第7根,如果是100根对应了第90根)可能存在1492种截法。

按上面枚举的结果修改一下代码2约束要求即可。

如第一次尝试:

eq 10
eq 10
eq 10
eq 10
eq 9

对于1492次,我们可以用宏来完成这项工作(宏代码略),遇到有合理解即停止循环,得到符合要求的结果(一般来说是唯一的,某些情况下有几种)。

小结:1,这个是简单的整数线性规划问题。由于约束问题的选择的不严格,代码1出现了误差,并且当要求的量(短铁丝根数)变大时容易从结果中看出来。

         2,由于约束不够所以要求1会出现很多种结果,代码2只给出了其中的一种(按照默认的先多后少原则)。问题中的要求2,则增加了一个约束,这时结果的唯一性才有可能。对于多次的修改代码工作可以通过宏来完成。但是由于这个约束部分基于要求1的结果,所以还需要人工来识别。

       3,实际中生产还存在很多问题要求来约束解,如截断方式个数的要求,可以通过修改代码来完成这些要求。

       4,关于对1650mm的原料铁丝的长度要求的优化问题,则大大的增加了问题的难度,不过可以根据代码1的想法及非线性整数规划得到接近于最优值的解(不过不一定是正确解)。

参考:
巧用SAS/OR软件求解多解整数线性规划 <<五邑大学学报(自然科学版)>>2004年 第18卷 第03期 邹祥福
倪勤. S AS最优化软件速成 .北京:科学出版社, 1 9 9 8 . 2 5 — 7 0 .