公卫人

 找回密码
 立即注册

QQ登录

只需一步,快速开始

查看: 756|回复: 2

[求助] 有序样本聚类???SAS??? wenti

[复制链接]
cherieling 发表于 2019-8-14 16:01:23 | 显示全部楼层 |阅读模式

注册后推荐绑定QQ,之后方才可以使用下方的“用QQ帐号登录”。

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
本人参考胡良平教授的教材《面向问题的统计学3》中的范例,进行有序样本聚类,但是未能得到想要的结果,请大佬指点一二
PS:本例一个分类指标,多个分类指标也是可以的。

Data mxwttjxt10_1;
  input x@@;
  cards;
3.19 5.36 4.58 5.39 4.97 4.93 4.67 2.41 1.64 0.35
  ;
run;
proc iml;
  use mxwttjxt10_1;
  read all into matrix_a;
  close mxwttjxt10_1;
  n=nrow(matrix_a);
  call symputx("n",nrow(matrix_a));
  call  symputx("m",ncol(matrix_a));

/*一、产生直径矩阵*/
  D=j(n,n,0);            /*建立n行n列的矩阵,元素初值均为0,以用来存储类的直径*/
  do i=1 to n;
   do j=1 to n;
    if i<j then D[i,j]=.;   /*用以得到下三角矩阵形式的直径矩阵*/
    else do;        
     temsum=matrix_a[i:j,];             /*选取第i行到第j行作为子矩阵*/
         temsum2=temsum-shape(temsum[:,],i-j+1,&m);          /*用子矩阵中各元素减去本列元素的均值,得到离均差*/
     D[i,j]=round(ssq(temsum2),0.001);          /*计算矩阵中所有元素的平方和,并保留3位小数*/
    end;
   end;
  end;
  create D var("obs1":"obs&n");     *此处&n为实际样本数;
  append from D;
/*产生直径矩阵结束*/

/*二、产生最小损失矩阵*/
  Lb=j(n,n,0);         /*建立n行n列的矩阵,元素均为0,以用来存储最小损失函数值*/
  lastN=j(n,n,0);     /*建立n行n列的矩阵,元素均为0,以用来存储最后一个分类的样品起始号*/
  do j=1 to n;
    do i=1 to n;        
          if i<j | j=1 then do;   *排除无意义的点,以建立下三角矩阵;
           Lb[i,j]=.;
           lastN[i,j]=.;
          end;
          else do;
           temp=j(i-j+1,1)-1;
            do k=j to i;     *设最后一段的起始点为k,由于有j段,所以k至少为j;
             if j=2 then temp[k-j+1]=D[k-1,1]+D[i,k];
             else temp[k-j+1]=Lb[k-1,j-1]+D[i,k];   *当j>2时,即分段数大于3时,用递归公式;
            end;                                
           Lb[i,j]=temp[><];        /*找出最小损失函数值*/
           lastN[i,j]=temp[>:<]+j-1;   /*由当前位置m=k-j+1推出最后一段的起始点k=m+j-1*/
          end;
    end;
  end;
  Lb_lastN=char(Lb,15,3)+"("+char(lastN,2,0)+")";    *将数值矩阵转化为字符矩阵,然后合并两矩阵的对应位置元素的内容,char函数的第三个参数为小数点的位数*;
  create Lb from Lb;
  append from Lb;
  create Lb_lastN from Lb_lastN;
  append from Lb_lastN;
/*产生最小损失矩阵结束*/

/*三、分别找出k分段点*/
  blockNum=j(n,n,.);    *初始化一个元素个数为n*n的矩阵;
  do k=n to 2 by -1;    *至少分2段,最多分n段;
    do j=k to 1 by -1;
     if j=k then blockNum[k,j]=lastN[n,j];    *首先得到最后一分段点,由矩阵lastN中获得,lastN[n,j]即表示n个样品分为j类时最后一类的起始点;
         else if j>1 then blockNum[k,j]=lastN[blockNum[k,j+1]-1,j];   /*欲求中间某段的起始点,也可由lastN中获得,相当于剩余的样品(其个数可用后一类起始点-1计算得到)分为j类时最后一类的起始点*/
         else blockNum[k,j]=1;                        
    end;
  end;
/*分别找出k分段点结束*/
  create blockNum from blockNum;
  append from blockNum;
quit;
ods html;
proc print data=D;             /*输出直径表*/
  title "&n.个有序样品的直径D表";
run;
proc print data=Lb_lastN;     /*输出最小损失表*/
  title "&n.个有序样品的最小损失表";
run;
proc print data=blockNum;    /*输出分类结果*/
  title "&n.个有序样品的聚类结果";
run;
title " ";
ods html close;
data a;                     
  set Lb;
  where col&n^=.;
  drop col1;
run;
proc transpose data=a out=b;
run;
data c;
  set b;
  cluster=_n_+1;
  drop _NAME_;
  rename col1=Lp;
run;
proc gplot data=c;     /*绘制最小损失函数值随分类数变化的趋势图*/
  plot Lp*cluster;
  symbol interpol=join;
run;

结果是只有直径矩阵的结果与书本例子的结果同,目标函数矩阵的结果又相同和不同,趋势图也没有画出来

下面贴对比图;
论坛3.jpg

直径矩阵

直径矩阵


论坛1.jpg

最小损失函数和趋势图

最小损失函数和趋势图

论坛2.jpg


在线等。。。大佬。。。

补充内容 (2019-8-31 16:59):
跟新:今日看到了胡良平的解答过程,这个结果是正确的,是书本上给的有问题。。另外,例子给的是一组数据,有个问题:如果是多组数据,应该怎么输入呢?
风儿dreaming 发表于 2020-1-8 20:59:09 | 显示全部楼层

回帖奖励 +1 分钢镚

请问楼主解决这个问题了么?我目前也遇到了要输入多行多列的数据的问题呢~
回复

使用道具 举报

 楼主| cherieling 发表于 2020-3-1 16:29:59 | 显示全部楼层
风儿dreaming 发表于 2020-1-8 20:59
请问楼主解决这个问题了么?我目前也遇到了要输入多行多列的数据的问题呢~

不好意思,好久没登录了。我后来放弃了该方法
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

会员|至尊|接种|公卫人 ( 沪ICP备06060850号-3 )

GMT+8, 2020-8-10 01:48 , Processed in 0.762320 second(s), 20 queries , Gzip On, MemCached On.

Powered by Discuz! X3.4

Copyright © 2001-2020, Tencent Cloud.

快速回复 返回顶部 返回列表