Данный раздел предназначен для тех, кто интересуется статистической обработкой данных и хотел бы подробнее узнать о методах и алгоритмах анализа.
Введение в анализ данных
- знакомство с основами обработки данных на
примере программы ПАНДА.
Популярные вычислительные
алгоритмы - совершенно готовые к использованию
и, при этом, оптимизированные!
Пакет ПАНДА 6.хх - средство многомерного статистического анализа данных. Работает под DOS. Обрабатываются прямоугольные матрицы типа ОБЪЕКТЫ х ПРИЗНАКИ. Объектов может быть до 25,000, признаков - до 1,000. Допускаются пропущенные значения. Тpебуемый для запуска минимальный объем свободной опеpативной памяти - 390К. На жестком диске пакет занимает менее 1 Мб, для нормальной работы требуется еще около 1 Мб свободного места.
Основные виды процедур:
|
|
Частотное описание |
|
|
Корреляционный анализ (для всех типов шкал) |
|
|
Кросс-табуляция |
|
|
Регрессионный анализ (линейный, шаговый, нелинейный) |
|
|
Кластерный анализ (5 методов, в т.ч. ISODATA, Кpаб, и дp.) |
|
|
Бинарная иерархия пpизнаков |
|
|
Факторный анализ (Метод главных компонент) |
|
|
Дискриминантный анализ (линейный, шаговый) |
|
|
Дисперсионный анализ (ANOVA, MANOVA) |
|
|
Проверка гипотез |
|
|
Непараметрические критерии |
|
|
Визуализация данных (графики, гистограммы, диаграммы) |
Для нормальной работы пакета необходимо иметь следующие файлы:
|
|
PAS.EXE - основной программный модуль / ОБЯЗАТЕЛЕН |
|
|
HELPMESS.PND - файл ситуационной помощи / желателен |
|
|
DICTMESS.PND - словарь специальных терминов / желателен |
|
|
TESTS.EXE - тестовые набоpы данных / необязателен |
|
|
TUTOR.EXE - обучающий модуль / желателен |
Попробовать программу можно, скачав её отсюда.
Пакет разработан частной фирмой Группа ПАНДА.
Алгоритмы обработки данных
Ниже приведены выдержки из программы ПАНДА, на основе которых вы легко сможете написать свои собственные процедуры на привычном для вас языке. Здесь использовался dBase (Clipper 5).
Базовый метод факторного анализа - Метод Главных Компонент.
Процедура PAS_W65
k1=m
corr=array(k1,k1)
decl x1[k1], mco[k1], s1[k1]
supex=0
supervise=2
* - вычисление корреляций
pas_corr()
set devi to scre
@ 13,3 clea to 22,54
@ 12,5 say "Вычисление собственных чисел ... "
colpo=col()
priv b[k1,k1], m[k1], y1[k1], r[k1]
* ---------------- метод Якоби -------------------
for p=1 to k1
for j=1 to p-1
b[p,j]=0
b[j,p]=0
next
b[p,p]=1
m[p]:=r[p]:=corr[p,p]
next
afill(y1,0)
rp=0
for i=1 to 50 && количество итераций, которое я могу
вытерпеть
@ 13,3 say subs(done,1,i)
s=0
for p=1 to k1-1
for q=p+1 to k1
s+=abs(corr[p,q])
next
next
if s=0.0000000
exit
endi
t1=0
if i<4
t1=0.20*s/k1^2
endi
for p=1 to k1-1
for q=p+1 to k1
g=100*abs(corr[p,q])
if i>4 .and. abs(r[p])+g=abs(r[p]) .and. abs(r[q])+g=abs(r[q])
corr[p,q]=0
loop
endi
if abs(corr[p,q])>t1
h=r[q]-r[p]
if abs(h)+g=abs(h)
t=corr[p,q]/h
else
t0=.5*h/corr[p,q]
t=1/(abs(t0)+sqrt(1+t0*t0))
if t0<0
t=-t
endi
endi
c=1/sqrt(1+t*t)
s=t*c
t2=s/(1+c)
h=t*corr[p,q]
y1[p]-=h
y1[q]+=h
r[p]-=h
r[q]+=h
corr[p,q]=0
for j=1 to p-1
g=corr[j,p]
h=corr[j,q]
corr[j,p]=g-s*(h+g*t2)
corr[j,q]=h+s*(g-h*t2)
next
for j=p+1 to q-1
g=corr[p,j]
h=corr[j,q]
corr[p,j]=g-s*(h+g*t2)
corr[j,q]=h+s*(g-h*t2)
next
for j=q+1 to k1
g=corr[p,j]
h=corr[q,j]
corr[p,j]=g-s*(h+g*t2)
corr[q,j]=h+s*(g-h*t2)
next
for j=1 to k1
g=b[j,p]
h=b[j,q]
b[j,p]=g-s*(h+g*t2)
b[j,q]=h+s*(g-h*t2)
next
@ 12,colpo say ltrim(str(++rp,6))
endi
next
next
for p=1 to k1
r[p]:=m[p]:=m[p]+y1[p]
y1[p]=0
next
show_time()
next
for j=1 to k1
m[j]=j
next
for j=1 to k1-1
for l=j+1 to k1
if r[m[j]]<r[m[l]]
a0=m[j]
m[j]=m[l]
m[l]=a0
endi
next
next
corr=aclone(b)
acopy(r,y1)
for j=1 to k1
r[j]=y1[m[j]]
for i=1 to k1
b[i,j]=corr[i,m[j]]
next
next
set devi to prin
@ prow()+3,1 say "Дисперсии компонент и проценты от
общей дисперсии"
kf=min(10,k1)
@ prow()+2,0 say ""
for i=1 to kf
@ prow(),9*i-1 say tran(r[i],'9999.999')
next
@ prow()+2,1 say "Доля, %"
for i=1 to kf
@ prow(),9*i say tran(100*r[i]/k1,'9999.99')
next
@ prow()+1,1 say "Накопл."
cum=0
for i=1 to kf
cum+=r[i]
@ prow(),9*i say tran(100*cum/k1,'9999.99')
next
@ prow()+3,1 say "Критерий Кеттелла выбора числа
компонент"
mx=19/r[1]
@ prow()+2,4 say "0"
@ prow(),22 say tran(r[1],'99.999')
@ prow()+1,4 say "+"
ct=repl('-',21)
for j=1 to k1
ct=stuff(ct,int(r[j]*mx+1),1,"+")
next
@ prow(),5 say ct
for j=1 to k1
@ prow()+1,1 say tran(j,'99')
if r[j]*mx+5<5
@ prow(),4 say "+"
else
@ prow(),4 say "¦"
@ prow(),r[j]*mx+5 say '.'
endi
next
@ prow()+3,1 say "Нагрузки признаков на главные
компоненты "+if(k1>10,"(первые 10)","")
@ prow()+1,0 say ""
for i=1 to k1
@ prow()+1,1 say prnames[realch[i]]
for j=1 to kf
@ prow(),9*j say tran(b[i,j],'99.9999')
next
next
* матрицу факторных нагрузок записать в файл
* (для регрессии на главные компоненты)
set devi to scre
release corr, m, y1
@ 12,5 say "Вычисляются значения компонент ... "
colpo=col()-2
@ 13,3 say repl('_',50)
@ 15,3 say vrzavoc_
colpo1=col()
use (fname+".pdd")
for t=1 to kf
appe blan
ppt="КМП"+str(t,3)
repl prnam with ppt, high with 998, low with -998, missi with 9999, scala with 3
next
use (fname+".pmf")
copy stru extended to tmpstru
use tmpstru
for t=1 to kf
appe blan
ppp="p"+ltrim(str(k+t,3))
repl Field_name with ppp, Field_type with "N", Field_len with 7, Field_dec with
4
next
use
crea tmpbody from tmpstru
erase tmpstru.dbf
use tmpbody
appe from (fname+".pmf")
use
erase (fname+".pmf")
rename tmpbody.dbf to (fname+".pmf")
use (fname+".pmf")
zamer=seconds()
for l=1 to n
@ 13,3 say subs(done,1,l*50/n)
@ 12,colpo say tran(n-l,'99999')
go l
for j=1 to kf
s=0
for i=1 to k1
ppp="p"+ltrim(str(realch[i],3))
xx=&ppp
s+=b[i,j]*(xx-x1[i])/s1[i]
next
if r[j]#0
s/=r[j]
endi
if abs(s)>=10
s=9.9999*s/abs(s)
endi
ppp="p"+ltrim(str(j+k,3))
repl &ppp with s
next
if l=1
elapsed(n,15,colpo1)
endi
next
clos data
setcolor("BG/")
@ 14,3 clea to 17,54
@ 15,3 say "В файл "+alltrim(fname)+" добавлены
признаки,"
@ 16,3 say "содержащие значения первых компонент
("+ltrim(str(kf,4))+"),"
@ 17,3 say "которые можно проанализировать
отдельно."
Конец процедуры PAS_65
Процедура PAS_CORR
set colo to W*/
@ 12,3 say ""
set colo to W/
@ 12,5 say "Подсчет сумм и квадратов"
colpo=col()
@ 13,3 say repl('_',50)
@ 15,3 say vrzavoc_
colpo1=col()
* если SUPERVISE=1, то:
* поддиагональные элементы CORR содержат
корреляции, а
* наддиагональные - вероятности хвостов (для корр.
анализа)
* если SUPERVISE=2, то:
* наддиагональные элементы CORR содержат
корреляции,
* с проверкой на ненормальность данных (для
факторного)
* если SUPERVISE=3, то:
* наддиагональные элементы CORR содержат
корреляции
* (для бинарной иерархии)
afill(x1,0) && суммы
afill(s1,0) && суммы квадратов
afill(mco,0) && кол-вА нормальных значений
set colo to /W
@ 24,0 clea to 24,79
set colo to W/
zamer=seconds()
for i=1 to n
@ 13,3 say subs(done,1,i*50/n)
go i
for j=1 to k1
ppp="p"+ltrim(str(realch[j],3))
z=&ppp
if z#misscode[realch[j]] .and. z#-999.7
mco[j]=mco[j]+1
x1[j]=x1[j]+z
s1[j]=s1[j]+z*z
endi
next
do show_time
a=inke()
if a=-9 .or. a=301
do way_out
endi
if i=1
do elapsed with n,15,colpo1
endi
next
for j=1 to k1 && вычисление средних и ст. отклонений
if mco[j]>0
s1[j]=sqrt((mco[j]*s1[j]-x1[j]^2)/(mco[j]*(mco[j]-1)))
x1[j]=x1[j]/mco[j]
else
s1[j]=0
x1[j]=0
endi
if s1[j]=0 .and. supervise=2
tone(140,10)
setcolor("W+/")
@ 17,3 say "Признак "+prnames[realch[j]]+" весь пропущен
или константа !"
@ 18,3 say "Вычисления прекращаются ..."
setcolor("W/")
supex=1
a=inke(8)
endi
corr[j,j]=1
next
if supex#1
@ 12,5 say "Подсчет корреляций, осталось "
colpo=col()
@ 13,3 say repl('_',50)
@ 15,3 say vrzavoc_
colpo1=col()
@ 15,colpo1 say " "
arrcnt=0
zamer=seconds()
for j=1 to k1-1
do show_time
ppp1="p"+ltrim(str(realch[j],3))
mcod1=misscode[realch[j]]
for i=j+1 to k1
mcod2=misscode[realch[i]]
ppp2="p"+ltrim(str(realch[i],3))
arrcnt=arrcnt+1
@ 12,colpo say k1*(k1-1)/2-arrcnt+1 pict '9999'
covar=0
if mco[j]>0 .and. mco[i]>0 .and. s1[i]>0 .and. s1[j]>0
for l=1 to n
go l
z1=&ppp1
if z1=mcod1 .or. z1=-999.7
z1=x1[j]
endi
z2=&ppp2
if z2=mcod2 .or. z2=-999.7
z2=x1[i]
endi
covar=covar+(z1-x1[j])*(z2-x1[i])
next
r=covar/((n-1)*s1[j]*s1[i])
do case
case supervise=1
corr[i,j]=r && поддиагональ
prb=0
if abs(r)<1
tstat=r*sqrt((n-2)/(1-r*r))
do student with (tstat),(n-2),prb
endi
corr[j,i]=prb && НАДдиагональ
case supervise>1
corr[j,i]=r && НАДдиагональ
endc
else && если один из признаков ЦЕЛИКОМ пропущен
do case
case supervise=1
corr[j,i]=0
corr[i,j]=7
case supervise=3
corr[j,i]=7
case supervise=2
corr[j,i]=0
endc
endi
if j=1 .and. i=2
do elapsed with k1*(k1-1)/2,15,colpo1
endi
@ 13,3 say subs(done,1,arrcnt*100/(k1*(k1-1)))
next
a=inke()
if a=-9 .or. a=301
way_out()
endi
next
endi
Конец процедуры PAS_CORR
Процедура CHISTEST
set fixed off
set deci to 4
clea
? "Тестирование качества аппроксимации ф-ции
распределения ХИ-квадрат"
? " "
prb=0
a=1
hi=0
v=0
do whil a#27
@ 5,5 say "Введите ХИ-кв." get hi pict '999.9999'
@ 6,5 say "Введите ст.св." get v pict '999'
read
do chi_square with hi,v,prb
? hi,v," аппр.:",prb
a=inke(0)
endd
retu
***************
proc chi_square
***************
parameters hi,v,prb
* Аппроксимация Пизера и Пратта функции
распределения ХИ-квадрат
* см. Дж.Мэйндоналд. Вычисл. алгоритмы в прикл.
стат. -М.: ФиС, 1988, стр. 284
if v<11
z=hi/2
z2=z*z
c=1
g=1
d=v/2
a=d
d3=d+2
yes=.T.
do whil yes
a=a+1
c=c*z/a
g=g+c
if c/g>0.0000005
yes=.T.
else
yes=.F.
endi
endd
a2=d3*d3
fu=1/(12*a2)*(1-1/a2*(1/30-1/a2*(1/105-1/(140*a2))))
g=g*exp(d*log(z)-d3*fu-(d3-.5)*log(d3)+d3-z)*(d+1)
prb=1.00000-g/sqrt(2*3.14159263)
else
d1_=v-1
t2=d1_/hi
d3=hi-v+2/3-.8/v
* Maindonald 287. 5710-5755
g=1
if t2>0
if abs(1.0000-t2)>0.1
g=(1-t2*t2+2*t2*log(t2))/(1-t2)^2
else
g=0
for j=1 to 5
g=g+2*(1-t2)^j/((j+1)*(j+2))
next
endi
endi
* Maindonald 284. 5135
z1=d3*sqrt((1+g)/(2*hi))
* Maindonald 282. 9.4
z1=abs(z1)
z2=z1*z1
if z1>1.9
m1=1-1/(z2+3-1/(.22*(z2+3.2)))
prb=m1/z1/10^(.21714724*z2+.39909)
else
t9=1/(1+.33267*z1)
a1=.4361836
a2=-.1201676
a3=.937298
prb=t9*(a1+a2*t9+a3*t9*t9)/exp(z2/2)/sqrt(2*3.141593)
endi
* добавление для использования в ХИ-кв.
if d3<0
prb=1-prb
endi
endi
retu
Конец процедуры CHISTEST
Процедура NORMAL
set fixed off
set deci to 4
use normtest
clea
? "Тестирование качества аппроксимации
нормальной функции распределения"
? " "
do whil !eof()
dev=deviate
prb=0
do normal with dev,prb
? dev,probab," аппр.:",prb
skip
endd
retu
***********
proc normal
***********
parameters dev,prb
* Аппроксимация Пизера и Пратта нормальной
функции распределения
z1=abs(dev)
z2=z1*z1
if z1>1.9
m1=1-1/(z2+3-1/(.22*(z2+3.2)))
prb=m1/z1/10^(.21714724*z2+.39909)
else
t9=1/(1+.33267*z1)
a1=.4361836
a2=-.1201676
a3=.937298
prb=t9*(a1+a2*t9+a3*t9*t9)/exp(z2/2)/sqrt(2*3.141593)
endi
retu
Конец процедуры NORMAL
Процедура PAS_CRSS
if uup=1 && если с управляющим
признаком (3-мер сопряж)
cgf0=contr[gf0]
else
cgf0=0
endi
use (fname+".pmf")
n=lastrec()
priv f_[40], yy[40], y[n]
g1_p=-1
g2_p=-1
srt='Сортировка значений признака '
pds='Подсчет количества различных значений'
for gf1=1 to k1
g1=realch1[gf1]
for gf2=1 to k2
use (fname+".pmf")
n=lastrec()
@ 4,3 clea to 22,54
@ 4,10 say "Считается сопряженность признаков"
@ 6,27 say "и"
g2=realch2[gf2]
@ 6,15 say tran(g1,'999')+": "+prnames[g1] color "/W"
@ 6,29 say tran(g2,'999')+": "+prnames[g2] color "/W"
helper=33
show_time()
sw_F10(0)
sw_Esc(0)
@ 10,3 clea to 12,52
if g1#g2
if g1#g1_p && если первый признак - следующий в
серии
* т.е., произошла смена первого признака
g1_p=g1
afill(yy,0,1,mx1)
afill(f_,0,1,mx1)
@ 10,8 say '' color "W+*/"
@ 10,10 say srt+prnames[g1]
pp1='p'+ltrim(str(g1,3))
pp2='p'+ltrim(str(g2,3))
h71=0
afill(y,0)
for t=1 to n
go t
if uup=1
z0=&pp0
else
z0=cgf0
endi
if z0=cgf0
z1=&pp1
z2=&pp2
if z2#misscode[g2] .and. z2#-999.7 .and. z1#misscode[g1] .and. z1#-999.7
h71++
y[h71]=z1
endi
endi
next
asort(y,1,h71)
d=1
yy[1]=1
f_[1]=y[1]
show_time()
@ 10,10 clea to 10,54
@ 10,10 say pds
for i=2 to h71
if y[i]#y[i-1]
if d<mx1 && по первому признаку (вертикаль) - до
mx1 значений
d++
f_[d]=y[i] && различные значения
endi
endi
yy[d]++ && сколько d-тых значений
next
show_time()
d1=d
tol='Только одно значение у признака '
endi
if d1>1
if g2_p#g2 && если сменился второй признак
afill(yy,0,mx1+1)
afill(f_,0,mx1+1)
g2_p=g2
@ 10,10 say srt+prnames[g2]+' '
pp2='p'+ltrim(str(g2,3))
h72=0
afill(y,0)
for t=1 to n
go t
if uup=1
z0=&pp0
else
z0=cgf0
endi
if z0=cgf0
z1=&pp1
z2=&pp2
if z2#misscode[g2] .and. z2#-999.7 .and. z1#misscode[g1] .and. z1#-999.7
h72++
y[h72]=z2
endi
endi
next
asort(y,1,h72)
d=1
yy[mx1+1]=1
f_[mx1+1]=y[1]
do show_time
@ 10,10 clea to 10,54
@ 10,10 say pds
for i=2 to h72
if y[i]#y[i-1]
if d<mx2 && по второму признаку (гориз) - до mx2
значений
d++
f_[mx1+d]=y[i] && различные значения
endi
endi
yy[mx1+d]++ && сколько d-тых значений
next
d2=d
endi
if d2>1
@ 10,10 clea to 10,54
@ 10,10 say 'Вычисление таблицы частот : %'
decl b[mx1*mx2]
afill(b,0)
w=0
alw=100/n
for l=1 to n
do show_time
go l
w++
@ 10,38 say w*alw pict '999'
if uup=1
z0=&pp0
else
z0=cgf0
endi
if z0=cgf0
z1=&pp1
z2=&pp2
checkt=z1#misscode[g1] .and. z1#-999.7 .and. z2#misscode[g2] .and. z2#-999.7
if checkt
ii=ascan(f_,z1,1,d1)
ii=iif(ii>0,ii,d1)
jj=ascan(f_,z2,mx1+1,mx1+d2)-mx1
jj=iif(jj>0,jj,d2)
s=(ii-1)*mx2+jj
b[s]++
endi
endi
next
old_n=n
n=0
for j=1 to d2
n+=yy[mx1+j]
next
toscr=(d1<=14 .and. d2<=5)
if toscr
@ 1,0,23,79 box '________ ' color "+/"
nn=n
n=old_n
@ 2,3 say fldat()
n=nn
@ 3,3 say "Таблица сопряженности признаков "
@ 3,col() say prnames[g1] color "W+/"
@ 3,col() say " и "
@ 3,col() say prnames[g2] color "W+/"
set colo to /W
for i=1 to 6
@ 6+i,2 say subs(prnames[g1],i,1)
next
set colo to W/
@ 3,50 say '. Показаны частоты'
else
@ 10,3 clea to 10,54
@ 10,4 say "Таблица слишком большая для экрана -"
@ 11,4 say "выводится только в файл ... "
@ 11,32 say "" color "W+*/"
endi
set devi to prin
nn=n
n=old_n
if file(fname+'.pnv')
use (fname+'.pnv')
bor=.T.
else
bor=.F.
endi
vlna=spac(10)
ifv=.F.
if bor .and. uup=1
loca for varnum=upr .and. varval=cgf0
if !eof()
vlna=alltrim(valname)
ifv=.T.
endi
endi
if uup=1
@ prow()+3,1 say "ТАБЛИЦА для значения
"+alltrim(str(contr[gf0],11,4))+if(ifv,"
("+vlna+")","")+" управляющего признака
"+prnames[upr]
@ prow()+1,1 say prnames[upr]+": "+prlabels[upr]
endi
@ prow()+2,1 say "По вертикали отложен признак
"+prnames[g1]+", по горизонтали - "+prnames[g2]
@ prow()+1,1 say prnames[g1]+": "+prlabels[g1]
@ prow()+1,1 say prnames[g2]+": "+prlabels[g2]
n=nn
plu=d2*5
@ prow()+2,0 say ''
for u=1 to 6
@ prow(),16+plu+u+u say subs(prnames[g2],u,1)
next
@ prow()+2,2 say " Частота |"
if pat=1
for i=1 to d2
vlna=spac(11)+"|"
if bor
loca for varnum=g2 .and. varval=f_[mx1+i]
if !eof()
vlna=valname+' |'
endi
endi
@ prow(),pcol() say vlna
next
@ prow(),d2*wid1+ofs1+12 say 'Итого по'
endi
@ prow()+1,2 say " Процент |"
if pat=2
@ prow()+1,2 say "Ожидаемое значение |"
@ prow()+1,2 say " Процент по строке |"
for i=1 to d2
vlna=spac(11)+"|"
if bor
loca for varnum=g2 .and. varval=f_[mx1+i]
if !eof()
vlna=valname+' |'
endi
endi
@ prow(),pcol() say vlna
next
@ prow(),d2*wid1+ofs1+12 say 'Итого по'
@ prow()+1,2 say "Процент по столбцу |"
endi
set devi to scre
for i=1 to d2
if toscr
@ beg,i*wid+ofs say f_[mx1+i] pict '999999.9999'
endi
set devi to prin
@ prow(),i*wid1+ofs1-3 say f_[mx1+i] pict '999999.9999'
@ prow(),pcol() say "|"
set devi to scre
next
set devi to prin
@ prow(),pcol()+5 say "строке"
set devi to scre
hi=0
output=0
for i=1 to d1
if toscr
@ beg+1+i,ofs say f_[i] pict '999999.9999'
endi
set devi to prin
if i=1
@ prow()+1,2 say
iif(pat=2,'--------',spac(11)+'--------|')+repl('-----------|',d2+pat-1)+'-----------'
else
if output<6
output++
@ prow()+1,1 say subs(prnames[g1],output,1)
else
@ prow()+1,0 say " "
endi
@ prow(),ofs1 say '--------|'+repl('-----------|',d2)+'-----------'
endi
if output<6
output++
@ prow()+1,1 say subs(prnames[g1],output,1)
else
@ prow()+1,0 say " "
endi
@ prow(),ofs1+1 say tran(f_[i],'9999.9')+" |"
for j=1 to d2
s=b[(i-1)*mx2+j]
if s>0
if toscr
set devi to scre
@ row(),wid*j+ofs+2 say s pict '9999'
set devi to prin
endi
@ prow(),wid1*j+ofs1-1 say s pict '9999'
endi
@ prow(),wid1*j+ofs1+8 say "|"
next
if toscr
set devi to scre
@ row(),wid*j+ofs+4 say yy[i] pict '9999'
set devi to prin
endi
@ prow(),wid1*j+ofs1+2 say yy[i] pict '9999'
if output<6
output++
@ prow()+1,1 say subs(prnames[g1],output,1)
else
@ prow()+1,0 say " "
endi
vlna=' '
if bor
loca for varnum=g1 .and. varval=f_[i]
if !eof()
vlna=valname
endi
endi
@ prow(),10 say vlna
prw=prow() && полный процент
@ prw,ofs1+8 say "|"
for j=1 to d2
s=100*b[(i-1)*mx2+j]/n
if s>0
@ prw,wid1*j+ofs1+1 say s pict '99.9'
endi
@ prw,wid1*j+ofs1+8 say "|"
next
@ prow(),wid1*j+ofs1+4 say 100*yy[i]/n pict '99.9'
if pat=2
if output<6
output++
@ prow()+1,1 say subs(prnames[g1],output,1)
else
@ prow()+1,0 say " "
endi
prw=prow() && ожидаемые значения
@ prw,ofs1+8 say "|"
endi
for j=1 to d2
s_=yy[i]*yy[mx1+j]/n
if pat=2
s=int(10*s_)/10
if s>0
@ prw,wid1*j+ofs1 say s pict '999.9'
endi
@ prw,wid1*j+ofs1+8 say "|"
endi
* подсчет статистики ХИ-квадрат
hi=hi+(b[(i-1)*mx2+j]-s_)^2/s_
next
if pat=2
if output<6
output++
@ prow()+1,1 say subs(prnames[g1],output,1)
else
@ prow()+1,0 say " "
endi
prw=prow() && процент по строке
@ prw,ofs1+8 say "|"
for j=1 to d2
s=100*b[(i-1)*mx2+j]/yy[i]
if s>0
@ prw,wid1*j+ofs1 say s pict '999.9'
endi
@ prw,wid1*j+ofs1+8 say "|"
next
if output<6
output++
@ prow()+1,1 say subs(prnames[g1],output,1)
else
@ prow()+1,0 say " "
endi
prw=prow() && процент по столбцу
@ prw,ofs1+8 say "|"
for j=1 to d2
s=100*b[(i-1)*mx2+j]/yy[mx1+j]
if s>0
@ prw,wid1*j+ofs1 say s pict '999.9'
endi
@ prw,wid1*j+ofs1+8 say "|"
next
endi
set devi to scre
next
set devi to prin
if output<6
output++
@ prw+1,1 say subs(prnames[g1],output,1)
else
@ prw+1,0 say " "
endi
@ prw+1,ofs1 say '--------|'+repl('-----------|',d2)+'-----------'
@ prw+2,ofs1+8 say '|'+repl(' |',d2)
@ prw+3,ofs1-1 say "Итого по |"
set devi to scre
do show_time
for i=1 to d2
if toscr
@ beg+3+d1,i*wid+ofs+2 say yy[mx1+i] pict '9999'
endi
set devi to prin
@ prw+3,i*wid1+ofs1-1 say yy[mx1+i] pict '9999'
@ prw+3,pcol()+5 say "|"
set devi to scre
next
if toscr
@ beg+3+d1,i*wid+ofs+4 say n pict '9999'
endi
if toscr
set colo to */W
@ 24,3 say ""
set colo to /W
@ 24,5 say "Подождите, пока рассчитываются
Р-значения ..."
set colo to W/
endi
do show_time
set devi to prin
@ prow(),i*wid1+ofs1+2 say n pict '9999'
@ prow()+1,ofs1 say "столбцу |"
for i=1 to d2
@ prow(),i*wid1+ofs1+1 say 100*yy[mx1+i]/n pict '99.9'
@ prow(),pcol()+3 say "|"
next
@ prow(),i*wid1+ofs1+3 say '100.0'
oldn=old_n-n
if oldn>0
@ prow()+2,1 say "В таблицу не включены объекты с
пропусками, их обнаружено "+ltrim(tran(oldn,'9999'))+" ,
или "+ltrim(tran(100*oldn/old_n,'999.9'))+" %"
endi
* печать мер и статистик
@ prow()+3,1 say "Статистика Значение АСО Уровень (Р)
Границы"
v=(d1-1)*(d2-1)
@ prow()+2,1 say "ХИ-квадрат "+tran(hi,'99999.999')+"
("+ltrim(tran(v,'999'))+")"
* вычисление Р-значения (prb) (уровень значимости) ,
распределение ХИ-квадрат
prb=0
do chi_square with hi,v,prb
@ prow(),57 say tran(prb,'9.99999')
e=min(d1,d2)
ee=iif(e>1,e-1,1)
@ prow()+1,1 say "Коэф. сопряженности Пирсона
"+tran(sqrt(hi/(hi+n)),'9.999')
@ prow(),57 say tran(prb,'9.99999')
@ prow(),72 say "0 , "+tran(sqrt(ee/e),'9.999')
aaso=n*ee
* вычисление Р-значения (prb) , стандартное
нормальное распределение
do normal with (sqrt(hi)),prb
@ prow()+1,1 say "Мера Крамера "+tran(sqrt(hi/aaso),'9.999')
@ prow(),45 say tran(1/sqrt(aaso),'999.999')+' '+tran(prb,'9.99999')+' 0 , 1'
if d1=2 .and. d2=2
f1122=b[1]*b[mx2+2]
f1221=b[2]*b[mx2+1]
if f1221>0
ot=f1122/f1221
@ prow()+2,1 say "Отношение шансов "+tran(ot,'999.999')
if b[1]*b[2]*b[mx2+1]*b[mx2+2]>0
chl=sqrt(1/b[1]+1/b[2]+1/b[mx2+1]+1/b[mx2+2])
as=ot*chl
* вычисление Р-значения (prb) , стандартное
нормальное распределение
do normal with (1/chl),prb
@ prow(),45 say tran(as,'999.999')+' '+tran(prb,'9.99999')
endi
endi
endi
if scale[g1]<3 .or. scale[g2]<3
s=0
t1=0
t2=0
for i=1 to d1
for j=1 to d2
s1=0
s2=0
if i>1
for l=1 to i-1
s1=s1+yy[l]
next
endi
if j>1
for l=1 to j-1
s2=s2+yy[mx1+l]
next
endi
s=s+b[(i-1)*mx2+j]*(s1+(yy[i]-n)/2)*(s2+(yy[mx1+j]-n)/2)
next
t1=t1+yy[i]^3-yy[i]
next
for j=1 to d2
t2=t2+yy[mx1+j]^3-yy[mx1+j]
next
r=12*s/sqrt((n^3-n-t1)*(n^3-n-t2))
@ prow()+2,1 say "Коэф. ранг.корр. Спирмена
"+tran(r,'99.999')
if n>2
as=sqrt((1-r*r)/(n-2))
* вычисление Р-значения (prb) , стандартное
нормальное распределение
if as#0
do normal with abs(r/as),prb
else
prb=0
endi
@ prow(),45 say tran(as,'999.999')+' '+tran(prb,'9.99999')+' -1 , 1'
endi
endi
p=0
q=0
sfaabb=0
for i=1 to d1
for j=1 to d2
* > >
s0=0
for l=i+1 to d1
for l1=j+1 to d2
s0=s0+b[(l-1)*mx2+l1]
next
next
p=p+s0*b[(i-1)*mx2+j]
* < <
s1=0
for l=1 to i-1
for l1=1 to j-1
s1=s1+b[(l-1)*mx2+l1]
next
next
* > <
s2=0
for l=i+1 to d1
for x=1 to j-1
s2=s2+b[(l-1)*mx2+x]
next
next
q=q+s2*b[(i-1)*mx2+j]
* < >
s3=0
for l=1 to i-1
for x=j+1 to d2
s3=s3+b[(l-1)*mx2+x]
next
next
aa=s0+s1
bb=s2+s3
sfaabb=sfaabb+b[(i-1)*mx2+j]*(aa-bb)^2
next
next
f=p-q
t4=2*e*f/(n*n*ee)
@ prow()+1,1 say "Тау-с Стьюарта "+tran(t4,'99.999')
as=2*e/(n^3*ee)*sqrt(n*n*sfaabb-4*n*f*f)
* вычисление Р-значения (prb) , стандартное
нормальное распределение
if as#0
do normal with abs(t4/as),prb
else
prb=0
endi
@ prow(),45 say tran(as,'999.999')+' '+tran(prb,'9.99999')+' -1 , 1'
h=p+q
a1=0
a2=0
for i=1 to d1
for j=1 to d2-1
a3=0
for l=j+1 to d2
a3=a3+b[(i-1)*mx2+l]
next
a1=a1+b[(i-1)*mx2+j]*a3
next
next
for j=1 to d2
for i=1 to d1-1
a3=0
for l=i+1 to d1
a3=a3+b[(l-1)*mx2+j]
next
a2=a2+b[(i-1)*mx2+j]*a3
next
next
as=sqrt( (4*n+10)/(9*n*(n-1)) )
r=f/sqrt((h+a1)*(h+a2))
* вычисление Р-значения (prb) , стандартное
нормальное распределение
if as#0
do normal with abs(r/as),prb
else
prb=0
endi
@ prow()+1,1 say "Тау-b Кендалла "+tran(r,'99.999')+"
"+tran(as,'999.999')+' '+tran(prb,'9.99999')+' -1 , 1'
@ prow()+1,1 say "Гамма Гудмена-Крускала
"+tran(f/h,'99.999')
@ prow(),71 say '-1 , 1'
@ prow()+1,1 say "Меры Сомерса : d(1->2) "+tran(f/(h+a2),'99.999')
@ prow(),71 say '-1 , 1'
@ prow()+1,16 say "d(2->1) "+tran(f/(h+a1),'99.999')
@ prow(),71 say '-1 , 1'
@ prow()+1,0 say ''
if d1=2 .and. d2=2
if f1221>0
@ prow()+2,1 say "Если уровни 1-го признака А1 и А2, а 2-го
- Б1 и Б2, то вероятность"
@ prow()+1,1 say "попадения на Б1 у об'екта из А1 в
(отношение шансов) раз больше, чем"
@ prow()+1,1 say "у об'екта из А2."
endi
endi
set devi to scre
if toscr
do sw_F10 with 1
set colo to +*/W
@ 24,3 say ""
if k1=1 .and. k2=1
set colo to /W
@ 24,5 say "Н"+ajmite_+"любую "+klavishu+" для
вывода процентных таблиц"
keyb ''
a=inke(0)
@ 24,0 say repl(' ',80)
if a=-9 .or. a=301
do way_out
endi
endi
show_time()
set colo to W/
@ 3,50 say '. Показаны проценты'
for i=1 to d1
for j=1 to d2
s=100*b[(i-1)*mx2+j]/n
if s>0
@ beg+1+i,wid*j+ofs+4 say s pict '99.9'
endi
next
@ row(),wid*j+ofs+4 say 100*yy[i]/n pict '99.9'
next
for i=1 to d2
@ beg+3+d1,i*wid+ofs+4 say 100*yy[mx1+i]/n pict '99.9'
next
@ beg+3+d1,i*wid+ofs+3 say '100.0'
if k1=1 .and. k2=1
set colo to /W
@ 24,5 say "Н"+ajmite_+"любую "+klavishu+" для
продолжения"
set colo to W/
keyb ''
a=inke(0)
if a=-9 .or. a=301
way_out()
endi
endi
show_time()
if pat=2
@ 3,50 say '. Проценты по строкам'
for i=1 to d1
for j=1 to d2
s=100*b[(i-1)*mx2+j]/yy[i]
if s>0
@ beg+1+i,wid*j+ofs+3 say s pict '999.9'
endi
next
@ row(),wid*j+ofs+4 say 100*yy[i]/n pict '99.9'
next
for i=1 to d2
@ beg+3+d1,i*wid+ofs+4 say 100*yy[mx1+i]/n pict '99.9'
next
@ beg+3+d1,i*wid+ofs+3 say '100.0'
if k1=1 .and. k2=1
keyb ''
a=inke(0)
show_time()
if a=-9 .or. a=301
do way_out
endi
endi
@ 3,50 say '. Проценты по столбцам'
for i=1 to d1
for j=1 to d2
s=100*b[(i-1)*mx2+j]/yy[mx1+j]
if s>0
@ beg+1+i,wid*j+ofs+3 say s pict '999.9'
endi
next
@ row(),wid*j+ofs+4 say 100*yy[i]/n pict '99.9'
next
for i=1 to d2
@ beg+3+d1,i*wid+ofs+4 say 100*yy[mx1+i]/n pict '99.9'
next
@ beg+3+d1,i*wid+ofs+3 say '100.0'
if k1=1 .and. k2=1
keyb ''
a=inke(0)
do show_time
if a=-9 .or. a=301
do way_out
endi
endi
endi
rest scre from crosst
else
@ 11,32 say " "
if psound
tone(800,1)
endi
endi
else
set colo to W+/
@ 10,10 say tol+prnames[g2]+" ! "
set colo to W/
keyb ''
a=inke(h3)
if a=-9 .or. a=301
way_out()
endi
endi
else
set colo to W+/
@ 10,10 say tol+prnames[g1]+" ! "
set colo to W/
keyb ''
a=inke(h3)
if a=-9 .or. a=301
way_out()
endi
endi
endi
next
next
Конец процедуры PAS_CRSS
Процедура PAS_DIST
set colo to W*/
@ 12,3 say ""
set colo to W/
@ 12,5 say "Подсчет сумм и квадратов"
colpo=col()
@ 13,3 say repl('_',50)
@ 15,3 say vrzavoc_
colpo1=col()
* если SUPERVISE=3, то:
* наддиагональные элементы CORR содержат
расстояния
* (для бинарной иерархии)
afill(x1,0) && суммы
afill(s1,0) && суммы квадратов
afill(mco,0) && кол-вА нормальных значений
set colo to /W
@ 24,0 clea to 24,79
set colo to W/
zamer=seconds()
for i=1 to n
@ 13,3 say subs(done,1,i*50/n)
go i
for j=1 to k1
ppp="p"+ltrim(str(realch[j],3))
z=&ppp
if z#misscode[realch[j]] .and. z#-999.7
mco[j]=mco[j]+1
x1[j]=x1[j]+z
s1[j]=s1[j]+z*z
endi
next
do show_time
a=inke()
if a=-9 .or. a=301
do way_out
endi
if i=1
do elapsed with n,15,colpo1
endi
next
for j=1 to k1 && вычисление средних и ст. отклонений
if mco[j]>0
s1[j]=sqrt((mco[j]*s1[j]-x1[j]^2)/(mco[j]*(mco[j]-1)))
x1[j]=x1[j]/mco[j]
else
s1[j]=0
x1[j]=0
endi
if s1[j]=0 .and. supervise=2
tone(140,10)
setcolor("W+/")
@ 17,3 say "Признак "+prnames[realch[j]]+" весь пропущен
или константа !"
@ 18,3 say "Вычисления прекращаются ..."
setcolor("W/")
supex=1
a=inke(8)
endi
corr[j,j]=0
next
if supex#1
@ 12,5 say "Подсчет расстояний, осталось "
colpo=col()
@ 13,3 say repl('_',50)
@ 15,3 say vrzavoc_
colpo1=col()
@ 15,colpo1 say " "
arrcnt=0
zamer=seconds()
for j=1 to k1-1
do show_time
ppp1="p"+ltrim(str(realch[j],3))
mcod1=misscode[realch[j]]
for i=j+1 to k1
mcod2=misscode[realch[i]]
ppp2="p"+ltrim(str(realch[i],3))
arrcnt=arrcnt+1
@ 12,colpo say k1*(k1-1)/2-arrcnt+1 pict '9999'
covar=0
if mco[j]>0 .and. mco[i]>0 .and. s1[i]>0 .and. s1[j]>0
for l=1 to n
go l
z1=&ppp1
if z1=mcod1 .or. z1=-999.7
z1=x1[j]
endi
z2=&ppp2
if z2=mcod2 .or. z2=-999.7
z2=x1[i]
endi
covar+=((z1-x1[j])/s1[j]-(z2-x1[i])/s1[i])^2
next
r=sqrt(covar)
do case
case supervise=1
corr[i,j]=r && поддиагональ
prb=0
if abs(r)<1
tstat=r*sqrt((n-2)/(1-r*r))
do student with (tstat),(n-2),prb
endi
corr[j,i]=prb && НАДдиагональ
case supervise>1
corr[j,i]=r && НАДдиагональ
endc
else && если один из признаков ЦЕЛИКОМ пропущен
do case
case supervise=1
corr[j,i]=0
corr[i,j]=7
case supervise=3
corr[j,i]=9
case supervise=2
corr[j,i]=0
endc
endi
if j=1 .and. i=2
do elapsed with k1*(k1-1)/2,15,colpo1
endi
@ 13,3 say subs(done,1,arrcnt*100/(k1*(k1-1)))
next
a=inke()
if a=-9 .or. a=301
do way_out
endi
next
endi
Конец процедуры PAS_DIST
Процедура
PAS_W02
* Статистика
whiskers=0
vybros=0
prod=1
mm=0
helper=20
set cursor off
qu=adir("*.pdd")
if qu>0
sw_Esc(0)
k=0
fname=' '
definefile(@k,@fname)
if !empty(fname)
decl prnames[k], hibond[k], lobond[k], misscode[k], scale[k], prlabels[k]
go top
np=0
do whil !eof()
np++
prnames[np]=prnam
hibond[np]=high
lobond[np]=low
misscode[np]=missi
scale[np]=scala
prlabels[np]=prlab
skip
endd
use (fname+".pmf")
n=recc()
if n>1
sw_F1(0)
decl y[max(n,14)], y1[n], y71[20], y72[20], x1[n]
@ 2,3 clea to 22,54
in=fldat()
@ 2,3 say in
outputfile("STA")
set devi to prin
heading()
@ 3,1 say "Простые статистики и частоты. "+dtvr()
@ 5,1 say in
@ 6,1 say " "
set devi to scre
sw_Esc(0)
@ 4,3 say "С какого признака начать ( PgUp PgDn, затем
-+)"
g1=1
ye=1
ye=selevars(6,19,k,@g1,prnames)
sw_f10(0)
* -------- обработка признаков ---------
m=g1
do while m>0
sw_Esc(0)
set colo to W/
@ 3,3 clea to 22,54
set colo to W+*/
@ 10,8 say ''
set colo to W/
@ 10,10 say 'Подсчет статистик для признака '+prnames[m]
afill(y,0)
afill(y1,0)
e0=0
s=0
s1=0
s5=0
h7=0
s6=0
misscnt=0 && missings counter
use (fname+".pmf")
for j=1 to n
go j
ppp='p'+ltrim(str(m,3))
z=&ppp
if z#misscode[m] .and. z#-999.7
h7++
y1[h7]=z
s+=z
s1+=z*z
else
misscnt++
endi
next
if h7>1 && если есть нормальные наблюдения (и их
>1)
if scale[m]>1
s1=sqrt((s1*h7-s*s)/(h7-1)/h7)
s/=h7
s2=s1/sqrt(h7)
endi
show_time()
@ 10,10 say 'Сортировка значений признака '+prnames[m]+"
"
asort(y1,1,h7)
keyb ''
show_time()
priv kvn[9]
kvn[1]=.05
kvn[2]=.1
kvn[3]=.15
kvn[4]=.25
kvn[5]=.5
kvn[6]=.75
kvn[7]=.85
kvn[8]=.9
kvn[9]=.95
for ik=1 to 9
npl=(1+h7)*kvn[ik]
ka=int(npl)
de=npl-ka
if ka=h7
kd=ka
else
kd=ka+1
endi
if ka=0
kd=1
ka=1
endi
kvn[ik]=y1[ka]+de*(y1[kd]-y1[ka])
next
s8=kvn[4]
s9=kvn[6]
s3=kvn[5]
e=1
d=1
y[1]=1
x1[1]=y1[1]
for i=2 to h7
if y1[i]#y1[i-1]
d++
x1[d]=y1[i]
endi
y[d]++
if e<=y[d]
e=y[d]
s4=i
endi
next
s4=y1[s4]
m1=y1[h7]
m2=y1[1]
@ 3,3 clea to 22,54
helper=21
sw_F10(1)
sw_Esc(1)
set devi to prin
@ prow()+2,1 say "+----------------+ +"+repl('-',52)+"+"
@ prow()+1,1 say "¦ ПРИЗНАК "+prnames[m]+" ¦ ¦
"+prlabels[m]+" ¦"
@ prow()+1,1 say "+----------------+ +"+repl('-',52)+"+"
@ prow()+1,2 say "Шкала
"+if(scale[m]=1,"Номинальная",if(scale[m]=2,"Ранговая",if(scale[m]=3,"Абсолютная","не
задана")))
if scale[m]>1
@ prow(),34 say "Средняя "+str(s,11,4)
endi
@ prow(),59 say "Мода "+str(s4,11,4)
@ prow()+1,0 say ""
if scale[m]>1
@ prow(),23 say "Ср. кв. отклонение "+str(s1,11,4)
endi
@ prow(),60 say "Min "+str(m2,11,4)
@ prow()+1,0 say ""
if scale[m]>1
@ prow(),23 say "Ст. ошибка средней "+str(s2,11,4)
endi
@ prow(),60 say "Max "+str(m1,11,4)
@ prow()+1,34 say "Медиана "+str(s3,11,4)
@ prow(),58 say "Q3-Q1 "+str(s9-s8,11,4)
@ prow()+2,4 say "Различных значений "+alltrim(str(d,4))+"
, пропусков "+if(misscnt>0,alltrim(str(misscnt,4))+"
("+alltrim(str(100*misscnt/n,6,2))+"%)","нет")
@ prow()+1,4 say 'Метка Значение Частота Накопл. Доля %
Накопл. Общий %'
set devi to scre
exup=s9+3*(s9-s8)
exlo=s8-3*(s9-s8)
set colo to W/
@ 4,3 say ltrim(str(m,3))+"-й признак "
set colo to W+/
@ 4,col() say prnames[m]
set colo to W/
if scale[m]>2
@ 6,3 say " Средняя "+str(s,9,2)
@ 7,3 say "Ср. кв. отклонение "+str(s1,9,2)
@ 8,3 say "Ст. ошибка средней "+str(s2,9,2)
endi
@ 9,3 say " Медиана "+str(s3,9,2)
@ 5,35 say " Мода "+str(s4,9,2)
@ 6,35 say " Min "+str(m2,9,2)
@ 7,35 say " Max "+str(m1,9,2)
@ 8,35 say " Q3-Q1 "+str(s9-s8,9,2)
@ 9,35 say "Значений "+tran(d,'9999')
if misscnt>0
@ 4,35 say "Пропусков "+tran(misscnt,'9999')
endi
@ 11,3 say ' Значение Частота Накопл. Доля % Накопл.'
np=0
chna=0
if file (fname+'.pnv')
use (fname+'.pnv')
bor=.T.
else
bor=.F.
endi
vb=0
for i=1 to d
np++
if np>10
show_time()
np=1
if prod#3
set colo to +*/W
@ 24,3 say ""
set colo to /W
@ 24,5 say "Н"+ajmite_+"любую "+klavishu+" для
продолжения просмотра (ESC - выход)"
keyb ''
a=inke(10)
@ 24,0 say repl(' ',80)
if a=27
exit
endi
if a=-9 .or. a=301
way_out()
endi
endi
set colo to W/
@ 12,3 clea to 21,54
endi
chna+=y[i]
@ 11+np,3 say str(x1[i],11,4)+" "+tran(y[i],'9999')+"
"+tran(chna,'9999')+" "+str(100*y[i]/h7,7,2)+"
"+str(100*chna/h7,7,2)
set devi to prin
vlna=' '
if bor
loca for varnum=m .and. varval=x1[i]
if !eof()
vlna=valname
endi
endi
ev=x1[i]<exlo .or. x1[i]>exup
if ev
vb++
endi
@ prow()+1,4 say vlna+' '+str(x1[i],11,4)+" "+str(y[i],4)+"
"+str(chna,4)+" "+str(100*y[i]/h7,7,2)+"
"+str(100*chna/h7,7,2)+" "+str(100*y[i]/n,7,2)+"
"+if(ev,">>","")
set devi to scre
next
set devi to prin
if vb>0
vybros=1
endi
if vb>0 .and. prod#3 && если были выбросы и прогон не
пакетный
vb=0
@ prow()+2,4 say "Резкие выбросы: значения и номера
объектов"
@ prow()+2,0 say ""
use (fname+".pmf")
for i=1 to n
go i
z=&ppp
if z#misscode[m] .and. z#-999.7
if z<exlo .or. z>exup
vb++
if vb=4
@ prow()+1,0 say ""
vb=1
endi
@ prow(),pcol()+4 say str(z,11,4)+tran(i," (объект 99999)")
endi
endi
next
endi
if d>1
if scale[m]>1
@ prow()+2,4 say "Квантили распределения
(приблизительно)"
@ prow()+2,10 say ".05 .10 .15 .25 .50 .75 .85 .90 .95"
@ prow()+1,0 say ""
for i=1 to 9
@ prow(),i*10-6 say kvn[i] pict '999999.99'
next
* рисование ящика-с-усами
whiskers=1
@ prow()+2,1 say tran(m2,'99999.99')
@ prow(),19 say tran((m1+m2*3)/4,'99999.99')
@ prow(),37 say tran((m1+m2)/2,'99999.99')
@ prow(),55 say tran((m1*3+m2)/4,'99999.99')
@ prow(),73 say tran(m1,'99999.99')
@ prow()+1,6 say
"+-----------------+-----------------+-----------------+-----------------+"
nor=70/(m1-m2)
qq3=(y1[(h7*3+3)/4]+y1[(h7*3+4.9)/4])/2
qq1=(y1[(h7+3)/4]+y1[(h7+5)/4])/2
qq2=(y1[(h7+1)/2]+y1[(h7+2)/2])/2
sere=repl(" ",int((qq3-qq1)*nor))
te=int((qq2-qq1)*nor)
if len(sere)>0
sere=stuff(sere,te,1,":")
endi
ser1=repl("-",int((qq3-qq1)*nor))
@ prow()+1,6 say " "+repl("
",int((qq1-m2)*nor))+"+"+ser1+"+"
@ prow()+1,6 say
"+"+repl("-",int((qq1-m2)*nor))+"¦"+sere+"+"+repl("-",int((m1-qq3)*nor))+"¦"
@ prow()+1,6 say " "+repl("
",int((qq1-m2)*nor))+"+"+ser1+"+"
endi
else
if misscnt=0
@ prow()+2,4 say "*** Признак не содержит информации -
рекомендуется исключить его !"
@ prow()+1,0 say ""
endi
endi
set devi to scre
if prod#3
set colo to +*/W
@ 24,3 say ""
set colo to /W
if m1#m2
@ 24,5 say "Н"+ajmite_+"любую "+klavishu+" для
продолжения "
keyb ''
a=inke(0)
@ 24,5 say "Строится гистограмма (оценочная, для
выбора) "
set colo to W/
g1=int(5*log(h7)) && кол-во интервалов
r1=(m1-m2)/g1 && длина одного
hys=array(g1)
hmax=0
for u=1 to g1
np=0
ngr=m2+(u-1)*r1
vgr=m2+u*r1
if u=g1
vgr=vgr+1
endi
for t=1 to d
if ngr<=x1[t] .and. x1[t]<vgr
np=np+y[t]
endi
next
hys[u]=np
hmax=max(hmax,np)
next
each=if(hmax>14,hmax/14,1)
@ 11,3 clea to 21,54
for u=1 to g1
hm=hys[u]/each
@ 20,3+u say '-'
for v=1 to int(hm/2)
@ 20-v,3+u say '_'
next
if hm-int(hm/2)>0
@ 20-int(hm/2)-1,3+u say '_'
endi
next
@ 21,4 say 'В каждом значке (_) в среднем '+tran(each,'9999.9')+'
объекта.'
set colo to /W
if scale[m]>1 && гистограмма общего вида
each=if(hmax>8,hmax/8,1)
if each>1 && печатается, если отличается от
гисто с номерами
set devi to prin
@ prow()+1,0 say ""
for u=hmax/each to 1 step -1
@ prow()+1,0 say ""
for v=1 to g1
hm=hys[v]/each
if int(hm+.5)>=u
@ prow(),v+6 say "_"
else
if hm>0 .and. hm<1 .and. u=1
@ prow(),v+6 say "_"
endi
endi
next
next
@ prow()+1,4 say
"min"+repl("-",g1)+"max"+if(each>1," каждый
значок '_' содержит до "+ltrim(tran(each,'9999.9'))+'
объекта.',"")
@ prow()+1,0 say ''
set devi to scre
endi
endi
endi
helper=22
use (fname+".pmf")
@ 24,5 say "Н"+ajmite_+"любую "+klavishu+" "
if m1#m2
@ 24,27 say "(F2 - гистограмма с номерами объектов)"
endi
keyb ''
a=inke(0)
if a=-9 .or. a=301
way_out()
endi
@ 24,0 say repl(' ',80)
set colo to W/
@ 16,59 clea to 17,78
if a=-1 .and. m1#m2 && если нужна гисто и признак не
константа
@ 11,3 clea to 21,54
set colo to W+*/
@ 14,8 say ''
set colo to W/
@ 14,10 say 'Построение гистограммы'
done=repl('_',50)
@ 15,3 say repl('_',50)
ass="Числа на гистограмме представляют номера
объектов"
if scale[m]=3 .or. d>12 && если абсолютная шкала или
значений больше 12
g1=int(5*log(h7))
if g1>12
g1=12
endi
for ee=1 to g1
y[ee]=0
next
decl v[g1*7]
afill(v,0)
h=(g1-1)/(m1-m2)
ppp='p'+ltrim(str(m,3))
for i=1 to n
@ 15,3 say subs(done,1,i*50/n)
go i
z=&ppp
if z#misscode[m] .and. z#-999.7 && system missing
r=int((z-m2)*h)+1
y[r]=y[r]+1
if y[r]<=7
v[(y[r]-1)*g1+r]=i && (I-1)*k+J, j=1,k
endi
endi
next
* начало печати гистограммы
@ 14,3 clea to 22,54
set colo to /W
for i=1 to g1
for j=1 to min(y[i],7)
@ 20-j,i*4 say v[(j-1)*g1+i] pict '9999'
next
if y[i]>7
set colo to W/
@ 12,i*4 say " "
set colo to /W
endi
next
set colo to /W
@ 24,3 say "F2 - просмотр выбросов (крайних об'ектов)
на гистограмме"
set colo to W/
@ 11,4 say ass
helper=23
set devi to prin
@ prow()+2,0 say ""
mx=0
for i=1 to g1
if y[i]>7
@ prow(),i*4 say " ++"
endi
mx=max(y[i],mx)
next
for j=min(mx,7) to 1 step -1
@ prow()+1,0 say ""
for i=1 to g1
@ prow(),i*4 say iif(v[(j-1)*g1+i]>0,tran(v[(j-1)*g1+i],'9999')," ")
next
next
@ prow()+1,4 say repl('----',int((s3-m2)*h))+'medi'+repl('----',g1-int((s3-m2)*h)-1)
@ prow()+1,4 say 'min'
@ prow(),g1*4 say ' max'
@ prow()+2,0 say " "
set devi to scre
@ 20,4 say repl('----',g1)
@ 20,4*(int((s3-m2)*h)+1) say 'medi'
@ 21,4 say 'min'
@ 21,g1*4 say ' max'
set colo to +*/
@ 16,59 say ""
set colo to W/
@ 16,61 say "Н"+ajmite_+klavishu
@ 17,61 say "для продолжения"
keyb ''
a=inke(0)
if a=-9 .or. a=301
do way_out
endi
@ 16,59 clea to 17,78
set colo to /W
@ 24,0 say repl(' ',80)
set colo to W/
if a=-1
yc=1
decl ch[28]
helper=24
do whil yc#0
set colo to +*/W
@ 24,3 say ""
set colo to /W
@ 24,5 say "В"+yberite_+"об'ект ( , , -+ ) _ Esc - Выход из
гистограммы"
yc=0
set colo to /W,/W*
cnm=0
for i=1 to g1
if i<3 .or. i>g1-2
for j=min(y[i],7) to 1 step -1
cnm=cnm+1
ch[cnm]=v[(j-1)*g1+i]
@ 20-j,i*4 prom tran(ch[cnm],'9999')
next
endi
next
menu to yc
@ 24,0 say repl(' ',80)
set colo to W/
set colo to W/ , /W
if yc=0
exit
endi
go ch[yc]
save scre
@ 16,59 clea to 17,78
@ 4,13,8+min(k,11),43 box '________ '
set colo to W+/
@ 5,15 say "Об'ект No. "+ltrim(str(ch[yc],4))
set colo to W/
@ 7,15 say "No. Признак Значение"
np=0
helper=25
for t=1 to k
np++
if np=12
np=1
set colo to +*/W
@ 24,3 say ""
set colo to /W
@ 24,5 say "Н"+ajmite_+"любую "+klavishu+" для
продолжения просмотра"
keyb ''
a=inke(0)
@ 24,0 say repl(' ',80)
set colo to W/
@ 8,14 clea to 7+min(k,11),41
@ 16,59 clea to 17,78
if a=27
exit
endi
if a=-9 .or. a=301
do way_out
endi
endi
@ 7+np,15 say ltrim(str(t,3))
if t=m
set colo to W+/
endi
@ row(),20 say prnames[t]
set colo to W/
ppp='p'+ltrim(str(t,3))
z=&ppp
@ row(),col()+2 say "= "+tran(z,'999999.9999')
next
set colo to +*/W
@ 24,3 say ""
set colo to /W
helper=26
@ 24,5 say "Н"+ajmite_+"любую "+klavishu+" для
продолжения"
keyb ''
a=inke(0)
@ 24,0 say repl(' ',80)
set colo to W/
if a=27
exit
endi
rest scre
if a=-9 .or. a=301
do way_out
endi
endd
@ 16,59 clea to 17,78
endi
else && для номинальных и ранговых шкал
decl v[d*7]
afill(v,0)
e0=d
ppp='p'+ltrim(str(m,3))
afill(y,0)
for i=1 to n
@ 15,3 say subs(done,1,i*50/n)
go i
z=&ppp
if z#misscode[m] .and. z#-999.7 && system missing
r=ascan(x1,z)
y[r]=y[r]+1
if y[r]<=7
v[(y[r]-1)*e0+r]=i
endi
endi
next
@ 11,3 clea to 22,54
g1=e0
set colo to W/
for i=1 to g1
for j=1 to min(y[i],7)
@ 20-j,i*4 say v[(j-1)*g1+i] pict '9999'
next
if y[i]>7
@ 12,i*4 say " "
endi
next
@ 11,4 say ass
set devi to prin
@ prow()+2,0 say ""
mx=0
for i=1 to g1
if y[i]>7
@ prow(),i*4 say " ++"
endi
mx=max(y[i],mx)
next
for j=min(mx,7) to 1 step -1
@ prow()+1,0 say " "
for i=1 to g1
@ prow(),i*4 say iif(v[(j-1)*g1+i]>0,tran(v[(j-1)*g1+i],'9999')," ")
next
next
@ prow()+1,4 say repl('----',e0)
for i=1 to 3
prw=prow()+1
for j=i to e0 step 3
@ prw,j*4 say ltrim(tran(x1[j],'999999.9'))
next
next
@ prow()+1,0 say " "
set devi to scre
@ 20,4 say repl('----',g1)
for i=1 to g1
@ 21,i*4 say '#'+ltrim(str(i,2))
next
set colo to +*/
@ 16,59 say ""
set colo to W/
@ 16,61 say "Н"+ajmite_+klavishu
@ 17,61 say "для продолжения"
keyb ''
a=inke(0)
if a=-9 .or. a=301
do way_out
endi
@ 16,59 clea to 17,78
set colo to /W
@ 24,0 say repl(' ',80)
set colo to W/
endi
endi
endi
else
do show_time
@ 10,3 clea to 10,54
if h7=0
@ 10,4 say 'Все значения признака '+prnames[m]+' пропущены !'
else
@ 10,4 say 'Обнаружен только один нормальный объект !'
endi
a=inke(3)
endi
if psound
tone(800)
endi
@ 3,3 clea to 22,54
if prod=3 .and. m=mm
prod=1
endi
if prod#3
prod=1
if m=k
prod=4
endi
helper=27
do sw_F1 with 1
do show_time
@ 5,3 prom " Следующий признак
("+ltrim(str(iif(m<k,m+1,1),3))+"-й,
"+prnames[iif(m<k,m+1,1)]+") "
@ 6,3 prom " Признак с заданным номером "
@ 7,3 prom " Начиная с ... и кончая ... "
@ 8,3 prom " Закончить обработку "
menu to prod
if prod=1
if m<k
m++
else
m=1
endi
endi
do sw_F1 with 0
if prod=4 .or. prod=0
exit
endi
if prod=2
helper=28
@ 7,3 say spac(37)
@ 5,3 say spac(37)
@ 8,3 say "В"+yberite_+"нужный признак ( , затем -+
)"
g2=m+1
if g2>k
g2=1
endi
ye=1
set colo to /W
@ 10,19 say tran(g2,'999')+": "+prnames[g2]
do whil ye#13
ye=inke(0)
if ye=-9 .or. a=301
do way_out
endi
if ye=5
if g2<k
g2++
else
g2=1
endi
if g2=m
if g2<k
g2++
else
g2=1
endi
endi
endi
if ye=18
if g2<k-10
g2+=10
endi
if g2=m
if g2<k
g2++
else
g2=1
endi
endi
endi
if ye=3
if g2>10
g2-=10
endi
if g2=m
if g2>1
g2--
else
g2=k
endi
endi
endi
if ye=24
if g2>1
g2--
else
g2=k
endi
if g2=m
if g2>1
g2--
else
g2=k
endi
endi
endi
@ 10,19 say tran(g2,'999')+": "+prnames[g2]
endd
m=g2
endi
if prod=3 && Начиная ... и кончая ...
@ 5,3 clea to 8,50
@ 8,3 say "В"+yberite_+"начальный (от) признак ( , -+
)"
g2=m+1
if g2>k
g2=1
endi
set colo to /W
@ 10,19 say tran(g2,'999')+": "+prnames[g2]
ye=1
do whil ye#13
ye=inke(0)
if ye=-9 .or. a=301
do way_out
endi
if ye=5
if g2<k
g2++
else
g2=1
endi
if g2=m
if g2<k
g2++
else
g2=1
endi
endi
endi
if ye=18
if g2<k-10
g2+=10
endi
if g2=m
if g2<k
g2++
else
g2=1
endi
endi
endi
if ye=3
if g2>10
g2-=10
endi
if g2=m
if g2>1
g2--
else
g2=k
endi
endi
endi
if ye=24
if g2>1
g2--
else
g2=k
endi
if g2=m
if g2>1
g2--
else
g2=k
endi
endi
endi
@ 10,19 say tran(g2,'999')+": "+prnames[g2]
endd
m=g2
set colo to W/
@ 12,3 say "В"+yberite_+"конечный (до) признак ( , -+
)"
g2=m+1
if g2>k
g2=k
endi
set colo to /W
@ 14,19 say tran(g2,'999')+": "+prnames[g2]
ye=1
do whil ye#13
ye=inke(0)
if ye=-9 .or. a=301
do way_out
endi
if ye=5 .and. g2<k
g2++
endi
if ye=24 .and. g2>m
g2--
endi
if ye=18 .and. g2<k-10
g2+=10
endi
if ye=3 .and. g2>m+10
g2-=10
endi
@ 14,19 say tran(g2,'999')+": "+prnames[g2]
endd
mm=g2
endi
else
if m<mm
m++
endi
endi
a=inke()
if a=-9 .or. a=301
do way_out
endi
endd
do sw_F1 with 1
set devi to prin
if vybros=1
@ prow()+2,1 say "Выбросом считается значение
(отмечено >>), отстоящее от нижнего"
@ prow()+1,1 say "или верхнего квартиля более чем на
ТРИ межквартильных размаха."
endi
@ prow()+2,1 say "Коды пропусков при вычислениях
игнорируются."
if whiskers=1
@ prow()+2,1 say "Ящик-с-усами: ящик содержит 50%
объектов, усы - по 25%. Медиана - двоеточие (:)."
endi
@ prow()+2,1 say "(c) "+rel_year+" "+grp_+panda
@ prow()+1,0 say " "
set printer to
set devi to scre
set cursor on
endi
endi
else
set curs off
do filenotfound
endi
release all like y*
release all like x1
set printer to
retu
Конец процедуры PAS_W02
Процедура PAS_W03
* КРОСС-ТАБУЛЯЦИЯ
helper=30
set colo to W/
qu=adir("*.pdd")
set cursor off
if qu>0
save scre to crosst
sw_Esc(0)
sw_F1(0)
k=0
fname=' '
do definefile with k,fname
if !empty(fname)
decl prnames[k], misscode[k], scale[k], prlabels[k]
go top
np=0
do whil !eof()
np++
prnames[np]=prnam
misscode[np]=missi
scale[np]=scala
prlabels[np]=prlab
skip
endd
use (fname+".pmf")
n=recc()
@ 2,3 clea to 22,54
@ 2,3 say fldat()
set curs off
@ 4,3 prom " Трехмерная сопряженность "
@ 4,34 prom " двумерная "
upr=2
menu to upr
uup=0
if upr=1
uup=1
@ 4,3 say "Укажите управляющий признак ( , затем -+
)"
g1=1
set colo to /W
@ 6,19 say tran(g1,'999')+": "+prnames[g1]
ye=1
do whil ye#13
ye=inke(0)
if ye=-9 .or. ye=301
way_out()
endi
if ye=18
if g1<k-10
g1+=10
endi
endi
if ye=5
if g1<k
g1++
else
g1=1
endi
endi
if ye=3
if g1>10
g1-=10
endi
endi
if ye=24
if g1>1
g1--
else
g1=k
endi
endi
@ 6,19 say tran(g1,'999')+": "+prnames[g1]
endd
upr=g1
mx0=20
set colo to W/
@ 2,3 clea to 22,54
decl f[40], y[n], yy[40]
afill(yy,0)
@ 10,8 say '' color "W+*/"
srt='Сортировка значений признака '
@ 10,10 say srt+prnames[g1]
pp1='p'+ltrim(str(g1,3))
h71=0
afill(y,0)
for t=1 to n
go t
z1=&pp1
if z1#misscode[g1] .and. z1#-999.7
h71++
y[h71]=z1
endi
next
asort(y,1,h71)
d=1
yy[1]=1
f[1]=y[1]
do show_time
@ 10,10 clea to 10,54
pds='Подсчет количества различных значений'
@ 10,10 say pds
for i=2 to h71
if y[i]#y[i-1]
if d<mx0 && по управляющему признаку (вертикаль)
- до mx0 значений
d++
f[d]=y[i] && различные значения
endi
endi
yy[d]++ && сколько d-тых значений
next
do show_time
d0=d
contr:=aclone(f)
endi
mx1=20
mx2=18
@ 2,3 clea to 22,54
@ 4,4 say "Определение серии признаков ..."
@ 4,16 say "ПЕРВОЙ" color "W+/"
@ 6,3 prom " взять все признаки "
@ 6,26 prom " взять только некоторые "
nec=1
menu to nec
priv ch[k],realch1[k],realch2[k]
afill(ch,0)
np=0
for i=1 to k
if scale[i]<3
ch[++np]=i
endi
next
if np>0
acopy(ch,realch1)
acopy(ch,realch2)
chcnt=np
if nec=2
do markvars with np,0,np,8,ch,realch1,chcnt
endi
k1=chcnt
@ 5,3 clea to 22,54
@ 4,16 say "ВТОР" color "W+/"
@ 6,3 prom " взять все признаки "
@ 6,26 prom " взять только некоторые "
nec=1
menu to nec
np=0
for i=1 to k
if scale[i]<3
ch[++np]=i
endi
next
chcnt=np
if nec=2
do markvars with np,0,np,8,ch,realch2,chcnt
endi
k2=chcnt
release ch
if k1>1 .or. k2>1
h3=.5
else
h3=3
endi
@ 2,3 clea to 22,54
@ 10,4 say "Выбирайте стрелками, затем -+ :
печатать"
@ 12,27 say "или"
@ 12,3 prom " только частоту и долю "
@ 12,31 prom " все параметры "
menu to pat
do outputfile with "CRO"
beg=5
wid=11
ofs=4
wid1=12
ofs1=13
set devi to prin
do heading
@ 3,1 say "Таблица сопряженности. "+dtvr()
@ 5,1 say fldat()
set devi to scre
if uup=1
for gf0=1 to d0
pp0='p'+ltrim(str(upr,3))
pas_crss()
next
else
pas_crss()
endi
set devi to prin
@ prow()+2,1 say "АСО - асимптотическая стандартная
ошибка меры."
@ prow()+1,1 say "Уровень (Р) - вероятность случайного
получения данной величины"
@ prow()+1,1 say "статистики (меры), реальный уровень
значимости. Если выбранный"
@ prow()+1,1 say "заранее уровень значимости больше,
то гипотеза о независимости"
@ prow()+1,1 say "признаков отвергается."
@ prow()+2,1 say "См. А.Афифи,С.Эйзен. Статистический
анализ. - М.: Мир, 1982."
@ prow()+1,1 say " Г.Аптон. Анализ таблиц
сопряженности. - М.: ФиС, 1982."
@ prow()+2,1 say "(c) "+rel_year+" "+grp_+panda
@ prow()+1,0 say " "
set devi to scre
set printer to
@ 13,4 say sdelan_+rezu+"CRO"
else
@ 8,4 say "Все признаки - в абсолютной шкале !" color
"W+/"
endi
toto()
naj()
endi
else
do filenotfound
endi
release all like f*
release all like y
release all like b
set cursor on
set printer to
retu
Конец процедуры PAS_W03
Процедура PAS_W04
* CORRELATION
w1=1
do whil w1#6
helper=40
do show_time
do sw_F1 with 1
do sw_F10 with 1
do sw_Esc with 1
vybmenu(1)
@ 14,57 prom " Простая линейная "
@ 15,57 prom " Ранговая корреляция "
@ 16,57 prom " Ассоциация "
@ 17,57 prom " Нелинейная "
@ 18,57 prom " Проверка значимости "
@ 19,57 prom " Выход "
menu to w1
if w1=0
exit
endi
@ 24,0 say repl(' ',80) color "/W"
save scre to corrscre
qu=adir("*.pdd")
if w1=1 && Простая линейная
if qu>0
set cursor off
do sw_Esc with 0
do sw_F1 with 0
k=0
fname=' '
do definefile with k,fname
if empty(fname)
loop
endi
release prnames, misscode, scale
priv prnames[k], misscode[k]
go top
np=0
qabsscale=0
decl ch[k] && номера числовых признаков
do whil !eof()
np=np+1
prnames[np]=prnam
misscode[np]=missi
if scala=3
qabsscale=qabsscale+1
ch[qabsscale]=np
endi
skip
endd
do show_time
chcnt=qabsscale
use (fname+".pmf")
n=recc()
@ 2,3 clea to 22,54
@ 2,3 say "Файл данных: "+fname+". Об'ектов
"+ltrim(str(n,5))+", признаков "+ltrim(str(k,3))+","
if qabsscale>1 && если есть признаки с абсолютной
шкалой
priv realch[min(50,qabsscale)]
for t=1 to min(chcnt,50)
realch[t]=ch[t]
next
@ 3,3 say "из них "+ltrim(str(qabsscale,3))+" в абсолютной
шкале."
if qabsscale>10 && больше, чем нужно - выбрать до 50
@ 4,3 say "Можно выбрать не более 50 признаков !"
do markvars with qabsscale,1,50,6,ch,realch,chcnt
endi
k1=chcnt
release x1, mco, s1, corr
priv x1[k1], mco[k1], s1[k1], corr[k1,k1]
supervise=1
supex=0
do pas_corr
@ 12,3 clea to 21,54
@ 13,3 say rezu+"COR"+iif(k1<9," и на
экране","")
@ 15,3 say "Задайте уровень (%), тогда будут
выводиться только"
alfa=5.0
@ 16,3 say "значимые на этом уровне
коэффициенты" get alfa pict '99.9' range .1,90
helper=41
do sw_F1 with 1
set curs on
read
set curs off
do sw_F1 with 0
do show_time
do outputfile with "COR"
set devi to prin
do heading
@ 3,1 say "Линейная корреляция. "+dtvr()
@ 5,1 say fldat()
@ 6,1 say "Признаков в абсолютной шкале
"+ltrim(str(qabsscale,3))+", из них выбрано
"+ltrim(str(k1,3))+"."
@ 9,1 say "Матрица коэффициентов корреляции
Пирсона"
maxpri=10
tp=k1
if tp>maxpri
tp=maxpri
endi
yp=0
do whil .T.
if yp>0
do p66 with (yp),yp-maxpri
endi
do p44 with (yp)
do p55 with (yp),(yp)
yp=yp+maxpri
if k1<=yp
exit
endi
tp=k1-yp
if tp>maxpri
tp=maxpri
endi
endd
@ prow()+2,1 say "Печатаются только значимые
коэффициенты. Уровень alfa= "+ltrim(tran(alfa,'99.9'))+"
%"
@ prow()+1,1 say "Звездочки обозначают, что данный
признак - константа (или весь пропущен)."
@ prow()+1,1 say 'Отдельные пропущенные значения при
подсчете заменяются "чистыми" средними.'
@ prow()+2,1 say "См. Э.Ферстер,Б.Ренц. Методы
корреляции и регрессии. - М.:ФиС, 1983."
@ prow()+2,1 say "(c) "+rel_year+" "+grp_+panda
@ prow()+1,0 say " "
set printer to
set devi to scre
if k1<=8
set colo to +/
@ 1,0,23,79 box '________ '
set colo to W/
@ 2,3 say "Файл данных: "+fname+". Об'ектов
"+ltrim(str(n,5))+", признаков "+ltrim(str(k,3))+","
@ 3,3 say "из них "+ltrim(str(qabsscale,3))+" в абсолютной
шкале. Выбрано "+ltrim(str(k1,3))+"."
for d=1 to k1
@ 5,8*d+3 say prnames[realch[d]]
next
@ 6,3 say ""
for d=1 to k1
set colo to W/
@ row()+1,3 say prnames[realch[d]]
for f=1 to d
d7=corr[f,d]
d8=corr[d,f]
do case
case d7<alfa/200
set colo to W+/
case d7>=alfa/200
set colo to W/
endc
if d8=7
set colo to R/
endi
@ row(),8*f+3 say iif(d8=7,' *',tran(d8,'99.999'))
next
next
@ row()+2,3 say "Значимые на уровне
"+ltrim(tran(alfa,'99.9'))+" % коэффициенты "
set colo to W+/
@ row(),col() say "выделены"
set colo to W/
@ row(),col() say "."
set colo to R/
@ row()+1,3 say '* '
set colo to W/
@ row(),col() say "обозначает, что данный признак -
константа (или весь пропущен)."
release corr, realch, mco, x1, s1, yesch, ch
do naj
if a=-9 .or. a=301
do way_out
endi
endi
else
set colo to W+*/
@ 4,4 say ""
set colo to W/
@ 4,6 say "Мало признаков в абсолютной шкале !"
keyb ''
a=inke(2)
if a=-9 .or. a=301
do way_out
endi
endi
else
do filenotfound
endi
endi
if w1=2 && Ранговая корреляция ____________________________
if qu>0
set cursor off
do sw_Esc with 0
do sw_F1 with 0
k=0
fname=' '
do definefile with k,fname
if empty(fname)
loop
endi
release lobond, hibond, scale, prnames, misscode
priv prnames[k], misscode[k], ch[k]
go top
np=0
qrangscale=0
do whil !eof()
np=np+1
prnames[np]=prnam
misscode[np]=missi
if scala>1
qrangscale=qrangscale+1
ch[qrangscale]=np
endi
skip
endd
do show_time
chcnt=qrangscale
use (fname+".pmf")
n=recc()
@ 2,3 clea to 22,54
@ 2,3 say "Файл данных: "+fname+". Об'ектов
"+ltrim(str(n,5))+", признаков "+ltrim(str(k,3))+","
if qrangscale>1 && если есть признаки
неноминальные
release realch
priv realch[min(chcnt,20)]
for t=1 to min(chcnt,20)
realch[t]=ch[t]
next
@ 3,3 say "из них "+ltrim(str(qrangscale,3))+" в допустимых
шкалах."
if qrangscale>10 && больше, чем нужно - выбрать до 20
@ 4,3 say "Можно выбрать не более 20 признаков !"
do markvars with qrangscale,1,20,6,ch,realch,chcnt
endi
set colo to W*/
@ 12,3 say ""
set colo to W/
@ 12,5 say "Считывание данных"
@ 13,3 say repl('_',50)
@ 15,3 say vrzavoc_
colpo1=col()
k1=chcnt
release vi, mco
priv vi[k1,n], mco[k1]
* поддиагональные элементы corr содержат
корреляции, а
* наддиагональные - вероятности хвостов
afill(mco,0) && кол-вА нормальных значений
set colo to /W
@ 24,0 clea to 24,79
set colo to W/
zamer=seconds()
for i=1 to n
@ 13,3 say subs(done,1,i*50/n)
go i
for j=1 to k1
ppp="p"+ltrim(str(realch[j],3))
z=&ppp
if z#misscode[realch[j]] .and. z#-999.7
mco[j]=mco[j]+1
endi
vi[j,i]=z
next
do show_time
a=inke()
if a=-9 .or. a=301
do way_out
endi
if i=1
do elapsed with n,15,colpo1
endi
do show_time
next
* полученную матрицу выбранных данных
рассортировать
* (приписать ранги) на месте (по признакам).
coef=n*(n*n-1)/6
@ 12,5 say "Вычисление коэффициентов "
colp=col()
@ 13,3 say repl('_',50)
@ 15,colpo1 say " "
release corr, xi, vr, ac
priv corr[k1,k1], xi[n], vr[n], ac[n/2], v_[n], w_[n]
corr[k1,k1]=1
zamer=seconds()
for u=1 to k1-1
@ 12,colp say k1-u pict '999'
corr[u,u]=1
@ 13,3 say subs(done,1,u*50/(k1-1))
do show_time
ppp1="p"+ltrim(str(realch[u],3))
mcod1=misscode[realch[u]]
for w=1 to n
xi[w]=vi[u,w]
next
asort(xi)
vc=0
z1=1
xx=xi[1]
vr[1]=1
for i=2 to n
@ 12,34 say n-i pict '99999'
if xx=xi[i]
z1=z1+1
else
xx=xi[i]
vr[i]=i
if z1>1
vc=vc+1
ac[vc]=z1
sr=0
for t=1 to z1
sr=sr+i-t
next
sr=sr/z1
for t=1 to z1
vr[i-t]=sr
next
z1=1
endi
endi
next
if z1>1
vc=vc+1
ac[vc]=z1
endi
* получили : в VR - ранги
* в AC - связки
for j=1 to n
@ 12,41 say n-j pict '99999'
v_[j]=vr[ascan(xi,vi[u,j])]
next
a_=0
for t=1 to vc
a_=a_+ac[t]*(ac[t]^2-1)
next
a_=a_/12
for v=u+1 to k1
@ 12,41 say k1-v pict '99999'
do show_time
for w=1 to n
xi[w]=vi[v,w]
next
asort(xi)
vc=0
z1=1
xx=xi[1]
vr[1]=1
for i=2 to n
@ 12,48 say n-i pict '999'
if xx=xi[i]
z1=z1+1
else
xx=xi[i]
vr[i]=i
if z1>1
vc=vc+1
ac[vc]=z1
z1=1
sr=0
for t=1 to z1
sr=sr+i-t
next
sr=sr/z1
for t=1 to z1
vr[i-t]=sr
next
endi
endi
next
if z1>1
vc=vc+1
ac[vc]=z1
endi
for j=1 to n
w_[j]=vr[ascan(xi,vi[v,j])]
next
b_=0
for t=1 to vc
b_=b_+ac[t]*(ac[t]^2-1)
next
b_=b_/12
* теперь вычисление коэффициента
r_=0
for t=1 to n
r_=r_+(v_[t]-w_[t])^2
next
r_=(coef-r_-a_-b_)/sqrt((coef-a_-a_)*(coef-b_-b_))
corr[v,u]=r_
prb=0
if abs(r_)<1
tstat=r_*sqrt((n-2)/(1-r_*r_))
do student with (tstat),(n-2),prb
endi
corr[u,v]=prb && НАДдиагональ
if v=2
do elapsed with k1*(k1-1)/2,15,colpo1
endi
next
next
@ 12,3 clea to 21,54
if psound
tone(1046,1)
endi
@ 13,3 say rezu+"COR"+iif(k1<9," и на
экране","")
@ 15,3 say "Задайте уровень (%), тогда будут
выводиться только"
alfa=5.0
@ 16,3 say "значимые на этом уровне
коэффициенты" get alfa pict '99.9' range .1,90
helper=41
do sw_F1 with 1
set curs on
read
set curs off
do sw_F1 with 0
do show_time
do outputfile with "COR"
set devi to prin
do heading
@ 3,1 say "Ранговая корреляция. "+dtvr()
@ 5,1 say fldat()
@ 6,1 say "Признаков в допустимых шкалах
"+ltrim(str(qrangscale,3))+", из них выбрано
"+ltrim(str(k1,3))+"."
@ 9,1 say "Матрица коэффициентов ранговой
корреляции Спирмена"
maxpri=10
tp=k1
if tp>maxpri
tp=maxpri
endi
yp=0
do whil .T.
if yp>0
do p66 with (yp),yp-maxpri
endi
do p44 with (yp)
do p55 with (yp),(yp)
yp=yp+maxpri
if k1<=yp
exit
endi
tp=k1-yp
if tp>maxpri
tp=maxpri
endi
endd
@ prow()+2,1 say "Печатаются только значимые
коэффициенты. Уровень alfa= "+ltrim(tran(alfa,'99.9'))+"
%"
@ prow()+1,1 say "Звездочки обозначают, что данный
признак - константа (или весь пропущен)."
@ prow()+2,1 say "См. Э.Ферстер,Б.Ренц. Методы
корреляции и регрессии. - М.: ФиС, 1983."
@ prow()+2,1 say "(c) "+rel_year+" "+grp_+panda
@ prow()+1,0 say " "
set printer to
set devi to scre
if k1<=8
set colo to +/
@ 1,0,23,79 box '________ '
set colo to W/
@ 2,3 say "Файл данных: "+fname+". Об'ектов
"+ltrim(str(n,5))+", признаков "+ltrim(str(k,3))+","
@ 3,3 say "из них "+ltrim(str(qrangscale,3))+" в порядковой
шкале. Выбрано "+ltrim(str(k1,3))+"."
for d=1 to k1
@ 5,8*d+3 say prnames[realch[d]]
next
@ 6,3 say ""
for d=1 to k1
set colo to W/
@ row()+1,3 say prnames[realch[d]]
for f=1 to d
d7=corr[f,d]
d8=corr[d,f]
do case
case d7<alfa/200
set colo to W+/
case d7>=alfa/200
set colo to W/
endc
if d8=7
set colo to R/
endi
@ row(),8*f+3 say iif(d8=7,' *',tran(d8,'99.999'))
next
next
@ row()+2,3 say "Значимые на уровне
"+ltrim(tran(alfa,'99.9'))+" % коэффициенты "
set colo to W+/
@ row(),col() say "выделены"
set colo to W/
@ row(),col() say "."
endi
release xi, ac, realch, mco, yesch, w_, corr
release all like v*
do naj
if a=-9 .or. a=301
do way_out
endi
else
set colo to W+*/
@ 4,4 say ""
set colo to W/
@ 4,6 say "Мало неноминальных признаков !"
keyb ''
a=inke(2)
if a=-9 .or. a=301
do way_out
endi
endi
else
do filenotfound
endi
endi
if w1=3 && Качественная
endi
if w1=4 && Нелинейная
endi
if w1=5 && Проверка значимости
endi
rest scre from corrscre
set printer to
endd
retu
*********************************************************
Конец процедуры PAS_W04
Процедура PAS_W05
* ДИАГРАММЫ
w1=1
do whil w1#5
show_time()
sw_F1(1)
sw_F10(1)
sw_Esc(1)
helper=50
vybmenu(1)
@ 14,57 prom " Диаграмма рассеяния "
@ 15,57 prom " Графики "
@ 16,57 prom " Гистограммы "
@ 17,57 prom " Пробит-график "
@ 18,57 prom " Выход "
menu to w1
if w1=0 .or. w1=5
exit
endi
setcolor("/W")
@ 24,0 say repl(' ',80)
setcolor("W/")
qu=adir("*.pdd")
if qu>0
set cursor off
save scre to diascreen
sw_Esc(0)
sw_F10(0)
sw_F1(0)
k=0
fname=' '
definefile(@k,@fname)
if empty(fname)
loop
endi
decl prnames[k], hibond[k], lobond[k], misscode[k], scale[k]
go top
np=0
do whil !eof()
np=np+1
prnames[np]=prnam
hibond[np]=high
lobond[np]=low
misscode[np]=missi
scale[np]=scala
skip
endd
use (fname+".pmf")
n=recc()
@ 2,3 clea to 22,54
@ 2,3 say fldat()
set curs off
@ 4,3 say "По горизонтали откладывается"
@ 5,19 say "по вертикали"
g1=1
set colo to /W
@ 4,32 say tran(g1,'999')+": "+prnames[g1]
ye=1
do whil ye#13
ye=inke(0)
if ye=5 .and. g1<k
g1=g1+1
endi
if ye=24 .and. g1>1
g1=g1-1
endi
@ 4,32 say tran(g1,'999')+": "+prnames[g1]
endd
if w1=1
g2=g1+1
if g2>k
g2=1
endi
@ 5,32 say tran(g2,'999')+": "+prnames[g2]
ye=1
do whil ye#13
ye=inke(0)
if ye=5 .and. g2<k
g2=g2+1
if g2=g1
if g2<k
g2=g2+1
else
g2=1
endi
endi
endi
if ye=24 .and. g2>1
g2=g2-1
if g2=g1
if g2>1
g2=g2-1
else
g2=k
endi
endi
endi
@ 5,32 say tran(g2,'999')+": "+prnames[g2]
endd
else
do case
case w1=4
@ 5,32 say "аргумент СНФР"
endc
endi
setcolor("W/")
if w1=1 && Диаграммы рассеяния
sw_F1(1)
helper=51
@ 7,3 prom " Диаграмма с частотами "
@ 7,30 prom " с номерами объектов "
menu to diatype
qint1=20
if diatype#2
ciph=2
qint2=30
else
helper=511
ciph=6
qint2=15
decl dianum[qint1*qint2],varg3[n]
afill(dianum,0)
g3=1
@ 9,3 say "Вывести также значение признака #" get g3
pict '999' range 1,k
set cursor on
read
set cursor off
endi
codestring=" 123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ_"
do sw_F1 with 0
decl dia[qint1*qint2]
afill(dia,0)
go top
mx1=-999
mx2=-999
mn1=9999
mn2=9999
ppp1="p"+ltrim(str(g1,3))
ppp2="p"+ltrim(str(g2,3))
if diatype=2
ppp3="p"+ltrim(str(g3,3))
endi
mcod1=misscode[g1]
mcod2=misscode[g2]
@ 11,3 say "Определяется разброс признаков"
@ 12,3 say repl('_',50)
do whil !eof()
@ 12,3 say repl('_',recno()*50/n)
z1=&ppp1
z2=&ppp2
if z1#mcod1 .and. z1#-999.7 .and. z2#mcod2 .and. z2#-999.7
mx1=max(z1,mx1)
mn1=min(z1,mn1)
mx2=max(z2,mx2)
mn2=min(z2,mn2)
endi
skip
endd
if mx1>mn1
if mx2>mn2
r1=(qint2-1)/(mx1-mn1)
r2=(qint1-1)/(mx2-mn2)
miscnt=0
@ 11,3 say "Происходит размещение объектов на
диаграмме"
@ 12,3 say repl('_',50)
for i=1 to n
@ 12,3 say repl('_',i*50/n)
go i
z1=&ppp1
z2=&ppp2
if z1#mcod1 .and. z1#-999.7 .and. z2#mcod2 .and. z2#-999.7
j1=int((z1-mn1)*r1+1) && по горизонтали
j2=int((z2-mn2)*r2+1) && по вертикали
j3=(j2-1)*qint2+j1
dia[j3]=dia[j3]+1
if diatype=2
do case
case dia[j3]=1
dianum[j3]=i
case dia[j3]=2
dianum[j3]=varg3[dianum[j3]]+&ppp3
case dia[j3]>2
dianum[j3]=dianum[j3]+&ppp3
endc
endi
else
miscnt=miscnt+1
endi
if diatype=2
varg3[i]=&ppp3
endi
next
if miscnt>0
@ 14,3 say "Количество объектов, не включенных в
вывод = "+ltrim(tran(miscnt,'9999'))+"."
endi
@ 11,3 clea to 11,54
@ 11,3 say "Происходит вывод диаграммы в файл"
@ 12,3 say repl('_',50)
sec="----"+repl("+-----",15)
fir="--"+repl("+---",14)+"-+---"
outputfile("DIA")
set devi to prin
heading()
@ 3,1 say "Диаграмма рассеяния. "+dtvr()
@ 5,1 say fldat()
if miscnt>0
@ 6,1 say "Количество объектов, не включенных в
диаграмму из-за пропусков, равно
"+ltrim(tran(miscnt,'9999'))+"."
else
@ 6,1 say "Пропусков в заданных признаках не
обнаружено."
endi
leftmarg=10
prw=prow()+2
if diatype=2
for v=2 to 1 step -1
prw=prw+1
for u=v to qint2 step 2
@ prw,leftmarg+u*ciph-5 say tran(mn1+(u-1)*(mx1-mn1)/(qint2-1),'9999.9')
next
next
else
for v=3 to 1 step -2
prw=prw+1
for u=v to qint2-2 step 4
@ prw,leftmarg+u*ciph-3 say tran(mn1+(u-1)*(mx1-mn1)/(qint2-1),'9999.9')
next
if v=1
@ prw,leftmarg+qint2*ciph-4 say tran(mx1,'9999.9')
endi
next
endi
prw=prw+1
@ prw,leftmarg say iif(diatype=2,sec,fir)
sle=0
for i=qint1 to 1 step -1
set devi to scre
@ 12,3 say repl('_',(qint1-i+1)*50/qint1)
set devi to prin
prw=prw+1
if i<7
sle=sle+1
@ prw,1 say subs(prnames[g2],sle,1)
endi
@ prw,3 say tran(mn2+(i-1)*(mx2-mn2)/(qint1-1),'9999.9')
@ prw,leftmarg say "-"
for j=1 to qint2
j3=(i-1)*qint2+j
if diatype=2
qp=dia[j3]
if qp=1
np=dianum[j3]
@ prw,leftmarg+j*ciph-4 say "N "+tran(np,'9999')
else
if qp>1
@ prw,leftmarg+j*ciph-4 say "S "+tran(qp,'9999')
endi
endi
else
@ prw,leftmarg+j*ciph say subs(codestring,min(dia[j3]+1,37),1)
endi
next
@ prw,leftmarg+qint2*ciph+2+int(ciph/6) say "-"
@ prw,pcol()+1 say ltrim(tran(mn2+(i-1)*(mx2-mn2)/(qint1-1),'9999.9'))
if diatype=2
prw=prw+1
@ prw,leftmarg say "|"
for j=1 to qint2
j3=(i-1)*qint2+j
qp=dia[j3]
np=dianum[j3]
if qp=1
@ prw,leftmarg+j*ciph-4 say tran(varg3[np],'9999.9')
else
if qp>1
@ prw,leftmarg+j*ciph-4 say tran(np/qp,'9999.9')
endi
endi
next
@ prw,leftmarg+qint2*ciph+2+int(ciph/6) say "|"
endi
next
prw=prw+1
@ prw,leftmarg say iif(diatype=2,sec,fir)
if diatype=2
for v=1 to 2
prw=prw+1
for u=v to qint2 step 2
@ prw,leftmarg+u*ciph-5 say tran(mn1+(u-1)*(mx1-mn1)/(qint2-1),'9999.9')
next
next
else
for v=1 to 3 step 2
prw=prw+1
for u=v to qint2-2 step 4
@ prw,leftmarg+u*ciph-3 say tran(mn1+(u-1)*(mx1-mn1)/(qint2-1),'9999.9')
next
if v=1
@ prw,leftmarg+qint2*ciph-4 say tran(mx1,'9999.9')
endi
next
endi
prw=prw+2
for u=1 to 6
@ prw,leftmarg+(u-1)*(ciph/2+1) say subs(prnames[g1],u,1)
next
if diatype=2
@ prow()+2,1 say "Ячейку диаграммы составляют два
параметра. Если сюда попал только один объект, то
печатается его номер"
@ prow()+1,1 say "(после N) и ниже соответствующее
значение признака "+prnames[g3]+", иначе -
количество объектов (после S)"
@ prow()+1,1 say "и средняя признака "+prnames[g3]+" для
данных объектов."
else
@ prow()+2,1 say "На диаграмме частоты заменены
символами от 1 до Z, причем A обозначает 10,"
@ prow()+1,1 say "B - 11, C - 12, ... , Y - 34, Z - 35. Числа, большие 35
заменяются символом _."
endi
@ prow()+2,1 say "(c) "+rel_year+" "+grp_+panda
@ prow()+1,0 say ""
else
set colo to W+/
@ 14,3 say "Признак "+prnames[g2]+" - константа или
весь пропущен !"
set colo to W/
a=inke(3)
endi
set devi to scre
release all like dia*
else
set colo to W+/
@ 14,3 say "Признак "+prnames[g1]+" - константа или
весь пропущен !"
set colo to W/
a=inke(3)
endi
endi
if w1=2
endi
if w1=3
endi
if w1=4 && Пробит-график **************************************
setcolor("W*/")
@ 8,3 say ""
setcolor("W/")
@ 8,5 say "Сортировка значений признака "+prnames[g1]
h7=0
misscnt=0 && missings counter
priv y1[n],y[n]
afill(y,0)
ppp='p'+ltrim(str(g1,3))
for j=1 to n
go j
z=&ppp
if z#misscode[g1] .and. z#-999.7
h7++
y1[h7]=z
else
misscnt++
endi
next
if h7>1
asort(y1,1,h7)
show_time()
d=1
y[1]=1
for i=2 to h7
if y1[i]#y1[i-1]
d++ && кол-во pазличных значений
y1[d]=y1[i] && значения аргумента
endi
y[d]++ && частоты
next
priv pg[20,30]
for i=1 to 20
for j=1 to 30
pg[i,j]=0
next
next
sm=y[1]
for i=2 to d
sm+=y[i]
y[i]=sm && накопленные частоты
next
@ 8,5 clea to 8,54
@ 8,5 say "Вычисление обратной СНФР ..."
xmin=y1[1]
xmax=y1[d]
zmin=i_normal(1/h7)
zmax=i_normal(1-1/h7)
xr=29/(xmax-xmin)
zr=19/(zmax-zmin)
for t=1 to d
i1=(if(t<d,i_normal(y[t]/h7),zmax)-zmin)*zr+1
j1=(y1[t]-xmin)*xr+1
pg[i1,j1]++
next
@ 8,5 clea to 8,54
@ 8,5 say "Вывод пробит-графика в файл ..."
outputfile("PRP")
set devi to prin
heading()
@ 3,1 say "Пробит-график (проверка на
нормальность). "+dtvr()
@ 5,1 say fldat()
@ 7,1 say "По горизонтали откладываются значения
признака "+prnames[g1]+", а по вертикали"
@ 8,1 say "- соответствующие (по вероятностям)
значения аргумента стандартной"
@ 9,1 say "нормальной функции распределения
(СНФР)."
if misscnt>0
@ 11,1 say tran(misscnt,"Обнаpужено пpопущенных значений
99999")
endi
@ prow()+2,1 say "Различных значений "+ltrim(tran(d,'99999'))
@ prow()+3,2 say tran(xmin,'999999.99')
@ prow(),17 say tran((xmin*3+xmax)/4,'999999.99')
@ prow(),32 say tran((xmin+xmax)/2,'999999.99')
@ prow(),47 say tran((xmin+xmax*3)/4,'999999.99')
@ prow(),62 say tran(xmax,'999999.99')
diar="+--------------+--------------+--------------+--------------+"
@ prow()+1,8 say diar
for i=1 to 20
@ prow()+1,1 say tran(zmax-(i-1)/zr,'999.99 |')
for j=1 to 30
@ prow(),j*2+7 say if(pg[21-i,j]>0,"o","")
next
@ prow(),68 say tran(zmax-(i-1)/zr,'| 999.99')
next
@ prow()+1,8 say diar
@ prow()+1,2 say tran(xmin,'999999.99')
@ prow(),17 say tran((xmin*3+xmax)/4,'999999.99')
@ prow(),32 say tran((xmin+xmax)/2,'999999.99')
@ prow(),47 say tran((xmin+xmax*3)/4,'999999.99')
@ prow(),62 say tran(xmax,'999999.99')
@ prow()+2,1 say "Если вид графика отличается от
прямой, то можно попробовать преобразовать"
@ prow()+1,1 say "исходный признак, например,
прологарифмировать его."
@ prow()+2,1 say "(c) "+rel_year+" "+grp_+panda
@ prow()+2,1 say ""
set devi to scre
@ 8,3 say " "+sdelan_+rezu+"PRP"
else
@ 8,4 say "Признак - константа или весь пропущен
!"
endi
endi
toto()
naj()
release all like dia*
release pg,y,y1
else
filenotfound()
endi
set printer to
endd
retu
********************************************************
Конец процедуры PAS_W05
Процедура PAS_W61
* Кластерный анализ ****************
w2=1
soh="Сохранить классифицирующий признак ?"
opis="описывающий полученное разбиение. Его
можно"
ispo="использовать при выводе диаграмм, в ANOVA и
т.д."
do whil w2#6
helper=610
show_time()
sw_F1(1)
sw_F10(1)
sw_Esc(1)
vybmenu(2)
@ 14,57 prom " Таксономия (Форель) "
@ 15,57 prom " К-средних Мак-Кина "
@ 16,57 prom " Сгущения Уолтца "
@ 17,57 prom " КНП (Краб) "
@ 18,57 prom " ISODATA "
@ 19,57 prom " Выход "
menu to w2
if w2=0 .or. w2=6
exit
endi
set colo to /W
@ 24,0 say repl(' ',80)
set colo to W/
qu=adir("*.pdd")
if qu>0
sw_Esc(0)
sw_F1(0)
k=0
fname=' '
do definefile with k,fname
if !empty(fname)
decl prnames[k], hibond[k], lobond[k], misscode[k], scale[k]
go top
np=0
do whil !eof()
np=np+1
prnames[np]=prnam
hibond[np]=high
lobond[np]=low
misscode[np]=missi
scale[np]=scala
skip
endd
sw_Esc(1)
sw_F1(1)
@ 2,3 clea to 22,54
use (fname+".pmf")
n=recc()
@ 2,3 say fldat()
if w2=1 && ____________ Forel _____________________________
priv gp[n],c[k,2],op[n],t9[15,k],fp[n],f1p[n],f9[k],o1[k],owe[k]
afill(o1,0)
afill(f9,0)
afill(owe,1)
r=100
u=.5
p=5
@ 4,3 prom " Назначить веса признаков "
@ 4,31 prom " Признаки равноценны "
helper=620
wes=2
menu to wes
helper=619
set curs on
if wes=1
for w=1 to k
@ 6,5 say "Вес признака "+prnames[w]+" равен" get
owe[w] pict '9.9' range 0,1
read
next
endi
helper=611
@ 4,3 clea to 6,54
@ 4,3 say "Желательное число таксонов" get p pict '99'
range 1,min(15,n-1)
@ 5,5 say "Начальный радиус таксона" get r pict '999999.9'
@ 6,3 say "Коэффициент его уменьшения" get u pict '9.9'
range 0.1,1
read
set curs off
do sw_Esc with 0
set colo to W+*/
@ 10,8 say ''
set colo to W/
@ 10,10 say 'Подсчет дисперсий ... '
cl=col()
@ 11,3 say repl('_',50)
@ 13,3 say vrzavoc_
colpo1=col()
zamer=seconds()
for i=1 to n
@ 11,3 say subs(done,1,i*50/n)
@ 10,cl say n-i pict '9999'
go i
for j=1 to k
ppp="p"+ltrim(str(j,3))
z=&ppp
o1[j]+=z
f9[j]+=z*z
next
do show_time
a=inke()
if a=-9
do way_out
endi
if i=1
do elapsed with n,13,colpo1
endi
next
for j=1 to k && вычисление средних и ст. отклонений
f9[j]=sqrt((n*f9[j]-o1[j]^2)/(n*(n-1)))
o1[j]/=n
next
@ 10,3 clea to 22,54
set colo to W+*/
@ 11,8 say ''
set colo to W/
@ 11,10 say 'Идет классификация ... шаг '
clu=col()
t=0
@ 8,3 say "Текущее количество таксонов "
tt=col()
@ 9,16 say "Текущий радиус "
tr=col()
r=r/u
set colo to W+/
shag=0
do whil t<p
r=r*u
t=0
afill(fp,0)
afill(op,0)
do whil .T.
@ 11,clu say ltrim(tran(++shag,'9999'))+" "
@ 8,tt say ltrim(tran(t,'99'))+" "
@ 9,tr say ltrim(tran(r,'999999.9'))+" "
do show_time
a=inke() && проверка на запрос выхода
if a=-9
do way_out
endi
v1=0
for j=1 to n
if op[j]#1
v1++
endi
next
if v1=0
exit
endi
i=ascan(op,0)
if i=0
tone()
@ 10,60 say "Ошибка счета !!!"
a=inke(0)
endi
go i
for j=1 to k
ppp="p"+ltrim(str(j,3))
c[j,1]=&ppp
next
do whil .T.
afill(gp,0)
afill(f1p,0)
v=0
for i=1 to n
if op[i]#1
r1=0
go i
for j=1 to k
ppp="p"+ltrim(str(j,3))
r1+=owe[j]*(c[j,1]-&ppp)^2
next
r1=sqrt(r1)
if r1<=r
gp[i]=1
v++
endi
endi
next
r1=0
for j=1 to k
c[j,2]=0
for i=1 to n
go i
if gp[i]#0
ppp="p"+ltrim(str(j,3))
c[j,2]=c[j,2]+&ppp
endi
next
c[j,2]=c[j,2]/v
r1+=owe[j]*(c[j,1]-c[j,2])^2
next
r1=sqrt(r1)
if r1<=.00001
t++
for e=1 to n
op[e]=op[e]+gp[e]
f1p[e]=t*gp[e]
fp[e]=fp[e]+f1p[e]
next
exit && GOTO 6
endi
for j=1 to k
c[j,1]=c[j,2]
next
endd
endd
endd
set colo to W/
@ 11,8 say ' Классификация завершена. '
if psound
tone(1046*2,1)
endi
set colo to W+*/
@ 13,3 say ''
set colo to W/
@ 13,5 say 'Происходит вывод результатов в файл ...'
outputfile("TAX")
set devi to prin
heading()
@ 3,1 say "Таксономия (метод Форель). "+dtvr()
@ 5,1 say fldat()
@ prow()+1,1 say "Число таксонов "+tran(t,'99')+",
радиус каждого "+tran(r,'999999.9')
for i=1 to t
@ prow()+2,1 say "***** Таксон номер "+tran(i,'99')+" .
Номера объектов :"
l=0
@ prow()+2,0 say ""
for sp=1 to k
c[sp,1]:=c[sp,2]:=0
next
for j=1 to n
if fp[j]=i
l++
go j
for sp=1 to k
ppp="p"+ltrim(str(sp,3))
c[sp,1]+=&ppp
c[sp,2]+=(&ppp*&ppp)
next
if l%20=0
@ prow()+1,0 say ""
endi
@ prow(),pcol()+1 say alltrim(tran(j,'9999'))
endi
next
@ prow()+1,1 say "Всего объектов "+tran(l,'9999')+", или
"+tran(l/n*100,'999.9 %')
l0=max(l-1,1)
@ prow()+2,1 say "Признак Вес Средняя Ср.кв.от.
F-отношение"
for sp=1 to k
c0:=c[sp,2]:=(c[sp,2]*l-c[sp,1]^2)/l/l0
c[sp,2]:=sqrt(c[sp,2])
if abs(c0)<0.0000001
c0=-1
endi
o1[sp]=f9[sp]^2/c0
c[sp,1]/=l
if t<=15
t9[i,sp]=c[sp,1]
endi
if owe[sp]>0
@ prow()+1,1 say prnames[sp]+" "+tran(owe[sp],'9.9')+"
"+tran(c[sp,1],' 999999.9')+" "+tran(c[sp,2],'999999.99')+"
"+tran(o1[sp],'99999999.9')
endi
next
@ prow()+1,1 say repl('-',48)
next
if t<=15
@ prow()+2,1 say "Расстояния между центрами таксонов
(веса учитываются)"
for i=1 to k
s=0
for j=1 to t
s+=t9[j,i]
next
o1[i]=s/t
next
prw=prow()+2
for i=0 to t-1
@ prw,i*9+8 say tran(i,'99-й')
next
for i=1 to t
@ prow()+1,0 say tran(i,'99-й')
s=0
for l=1 to k
s+=owe[l]*(t9[i,l]-o1[l])^2
next
s=sqrt(s)
@ prow(),5 say tran(s,'9999999.9')
for j=1 to i-1
s=0
for l=1 to k
s+=owe[l]*(t9[i,l]-t9[j,l])^2
next
s=sqrt(s)
@ prow(),j*9+5 say tran(s,'9999999.9')
next
next
endi
@ prow()+3,1 say "'Нулевой' таксон содержит все
имеющиеся объекты (это просто обозначение)."
@ prow()+2,1 say "Средние для каждого таксона - это
координаты его геометрического центра."
@ prow()+1,1 say "F-критерий - отношение общей
дисперсии признака к дисперсии внутри
таксона,"
@ prow()+1,1 say "(отрицательное значение, если
внутритаксонная дисперсия равна нулю)."
@ prow()+2,1 say "Признаки, не участвовавшие в
классификации (с весом 0), не печатаются."
@ prow()+1,1 say "Данная процедура не предназначена
для обработки признаков с пропусками,"
@ prow()+1,1 say "поэтому предварительно
рекомендуется заполнить все пропуски средними
!"
@ prow()+3,1 say "См. Н.Загоруйко и др. ППП ОТЭКС (для
анализа данных). - М.: ФиС, 1986."
@ prow()+2,1 say "(c) "+rel_year+" "+grp_+panda
@ prow()+1,1 say ""
set devi to scre
@ 13,3 say sdelan_+rezu+"TAX "
if psound
tone(1046*2,1)
endi
@ 15,3 say soh
@ 15,41 prom " Да "
@ 15,47 prom " Нет "
helper=612
sp=1
menu to sp
if sp=1
setcolor("W+*/")
@ 17,3 say ""
setcolor("W/")
@ 17,5 say "Пpизнак сохpаняется ..."
* добавление признака ТАКСОН
use (fname+".pdd")
appe blan
repl prnam with "ТАКСОН", high with t, low with 1, missi with 0, scala
with 1
use (fname+".pmf")
copy stru extended to tmpstru
use tmpstru
appe blan
ppp="p"+ltrim(str(k+1,3))
repl Field_name with ppp, Field_type with "N", Field_len with 3, Field_dec with
0
use
crea tmpbody from tmpstru
erase tmpstru.dbf
use tmpbody
appe from (fname+".pmf")
for i=1 to n
go i
repl &ppp with fp[i]
next
use
erase (fname+".pmf")
rename tmpbody.dbf to (fname+".pmf")
setcolor("BG/")
@ 15,3 clea to 17,54
@ 15,3 say "В файл "+alltrim(fname)+" добавлен признак
ТАКСОН,"
@ 16,3 say opis
@ 17,3 say ispo
setcolor("W/")
endi
release all like g*
release all like o*
release all like t*
release all like f*
release all like c*
endi
if w2=2 && ____________ k-means ___________________________
priv e1[15,k],wp[15],x[k],y[15],e0[15,k],np[n],owe[k],z[n],zp[n]
@ 4,3 prom " Назначить веса признаков "
@ 4,31 prom " Признаки равноценны "
helper=620
afill(owe,1)
wes=2
menu to wes
helper=619
set curs on
if wes=1
for w=1 to k
@ 6,5 say "Вес признака "+prnames[w]+" равен" get
owe[w] pict '9.9' range 0,1
read
next
endi
helper=621
@ 4,3 clea to 6,54
r=3
@ 4,3 say "Сколько кластеров необходимо
получить" get r pict '99' range 1,min(15,n-1)
read
set curs off
setcolor("W+*/")
@ 7,8 say ''
setcolor("W/")
@ 7,10 say 'Идет классификация ... шаг '
clu=col()
afill(wp,1)
afill(np,0)
for i=1 to r
np[i]=i
go i
for j=1 to k
ppp="p"+ltrim(str(j,3))
e1[i,j]=&ppp
next
next
e0=aclone(e1)
* 4
o=0
setcolor("W+/")
do while .T.
o++
@ 7,clu say ltrim(tran(o,'999'))+" "
for i=1 to n
if np[i]=0
go i
for j=1 to k
ppp="p"+ltrim(str(j,3))
x[j]=&ppp
next
* sub 24 --------------------
s=0
for l=1 to k
s+=owe[l]*(e1[1,l]-x[l])^2
next
a=1
for m=2 to r
p=0
for l=1 to k
p+=owe[l]*(e1[m,l]-x[l])^2
next
if s>p
s=p
a=m
endi
next
* return from sub 24 --------
wp[a]++
for j=1 to k
e1[a,j]=(e1[a,j]*(wp[a]-1)+x[j])/wp[a]
next
np[i]=a
endi
next
s1=0
for j=1 to r
s=0
for i=1 to k
s+=(e1[j,i]-e0[j,i])^2
next
y[j]=sqrt(s)
s1+=y[j]
next
if s1<=r*.01 .or. o>100
exit
endi
e0=aclone(e1)
afill(wp,0)
afill(np,0)
endd
* 6
setcolor("W/")
@ 7,8 say ' Классификация завершена. '
if psound
tone(1046*2,1)
endi
setcolor("W+*/")
@ 10,3 say ''
setcolor("W/")
@ 10,5 say 'Происходит вывод результатов в файл ...'
do outputfile with "CLU"
set devi to prin
do heading
@ 3,1 say "Кластерный анализ (К-средних Мак-Кина).
"+dtvr()
@ 5,1 say fldat()
@ prow()+1,1 say "Число кластеров "+tran(r,'99')
for i=1 to r
@ prow()+2,1 say "***** КЛАСТЕР "+ltrim(tran(i,'99'))+":
объектов "+ltrim(tran(wp[i],'9999'))+". Мера сходимости =
"+ltrim(tran(y[i],'9999.999'))+". Номера объектов :"
@ prow()+2,0 say ""
if wp[i]>0
afill(x,0)
t=0
for j=1 to n
if np[j]=i
t++
zp[t]=j
s=0
go j
for l=1 to k
ppp="p"+ltrim(str(l,3))
x[l]+=&ppp*&ppp
s+=owe[l]*(e1[i,l]-&ppp)^2
next
z[j]=sqrt(s)
endi
next
if t#1
for j=1 to t-1
for l=j+1 to t
if z[zp[j]]>=z[zp[l]]
a=zp[j]
zp[j]=zp[l]
zp[l]=a
endi
next
next
endi
for j=1 to t
if j%20=0
@ prow()+1,0 say ""
endi
@ prow(),pcol()+1 say alltrim(tran(zp[j],'9999'))
next
@ prow()+2,1 say "Признак Вес Эталон Разброс"
s=max(t-1,1)
for l=1 to k
x[l]=sqrt((x[l]-t*e1[i,l]^2)/s)
if owe[l]>0
@ prow()+1,1 say prnames[l]+" "+tran(owe[l],'9.9')+"
"+tran(e1[i,l],'999999.9')+" "+tran(x[l],'99999.9')
endi
next
@ prow()+1,1 say repl('-',72)
endi
next
@ prow()+2,1 say "Объекты внутри кластера ранжированы
по удалению от эталона."
@ prow()+1,1 say "Мера сходимости характеризует
устойчивость кластера (чем меньше, тем лучше)."
@ prow()+2,1 say "Расстояния между эталонами
кластеров (веса учитываются) :"
t=r
for i=1 to k
s=0
for j=1 to t
s+=e1[j,i]
next
x[i]=s/t
next
prw=prow()+2
for i=0 to t-1
@ prw,i*10+8 say tran(i,'99-й')
next
for i=1 to t
@ prow()+1,0 say tran(i,'99-й')
s=0
for l=1 to k
s+=owe[l]*(e1[i,l]-x[l])^2
next
s=sqrt(s)
@ prow(),5 say tran(s,'9999999.9')
for j=1 to i-1
s=0
for l=1 to k
s+=owe[l]*(e1[i,l]-e1[j,l])^2
next
s=sqrt(s)
@ prow(),5+j*10 say tran(s,'9999999.9')
next
next
@ prow()+3,1 say "'Нулевой' кластер содержит все
имеющиеся объекты (это просто обозначение)."
@ prow()+1,1 say "Эталон для каждого таксона - это
координаты его геометрического центра."
@ prow()+2,1 say "Признаки, не участвовавшие в
классификации (с весом 0), не печатаются."
@ prow()+1,1 say "Данная процедура не предназначена
для обработки признаков с пропусками,"
@ prow()+1,1 say "поэтому предварительно
рекомендуется заполнить все пропуски средними
!"
@ prow()+3,1 say "См. И.Мандель. Кластерный анализ. - М.:
ФиС, 1988."
@ prow()+2,1 say "(c) "+rel_year+" "+grp_+panda
@ prow()+1,1 say ""
set devi to scre
@ 10,3 say sdelan_+rezu+"CLU "
if psound
tone(1046*2,1)
endi
@ 15,3 say soh
@ 15,41 prom " Да "
@ 15,47 prom " Нет "
helper=612
sp=1
menu to sp
if sp=1
setcolor("W+*/")
@ 17,3 say ""
setcolor("W/")
@ 17,5 say "Пpизнак сохpаняется ..."
* добавление признака КЛАСТР
use (fname+".pdd")
appe blan
repl prnam with "КЛАСТР", high with t, low with 1, missi with 0, scala
with 1
use (fname+".pmf")
copy stru extended to tmpstru
use tmpstru
appe blan
ppp="p"+ltrim(str(k+1,3))
repl Field_name with ppp, Field_type with "N", Field_len with 3, Field_dec with
0
use
crea tmpbody from tmpstru
erase tmpstru.dbf
use tmpbody
appe from (fname+".pmf")
for i=1 to n
go i
repl &ppp with np[i]
next
use
erase (fname+".pmf")
rename tmpbody.dbf to (fname+".pmf")
setcolor("BG/")
@ 15,3 clea to 17,54
@ 15,3 say "В файл "+alltrim(fname)+" добавлен признак
КЛАСТР,"
@ 16,3 say opis
@ 17,3 say ispo
setcolor("W/")
endi
release all like e*
release all like wp
release all like x*
release all like y*
release all like n*
release all like o*
release all like z*
endi
if w2=3 && ____________ сгущения
__________________________
priv npr[n],y[n],yp[n],e[15,k],owe[k],ep[15],x[k]
@ 4,3 prom " Назначить веса признаков "
@ 4,31 prom " Признаки равноценны "
helper=620
afill(owe,1)
wes=2
menu to wes
helper=619
set curs on
if wes=1
for w=1 to k
@ 6,5 say "Вес признака "+prnames[w]+" равен" get
owe[w] pict '9.9' range 0,1
read
next
endi
helper=642
@ 4,3 clea to 6,54
r=0
do whil r<=0
r=3
@ 4,3 say "Задайте радиус кластера (порог)" get r pict
'99999.9'
read
endd
set curs off
j0=int(n/2+1)
setcolor("W+*/")
@ 7,3 say ''
setcolor("W/")
@ 7,5 say 'Идет классификация ... шаг '
setcolor("W+/")
clu=col()
rad=.T.
cnt=0
do whil rad
rad=.F.
afill(ep,0)
afill(npr,0)
for i=1 to 15
for j=1 to k
e[i,j]=0
next
next
npr[j0]=1
epr:= ++ep[1]
go j0
for c=1 to k
ppp="p"+ltrim(str(c,3))
e[1,c]=(e[1,c]*(epr-1)+&ppp)/epr
next
g=1
for j=1 to n
cnt++
@ 7,clu say cnt pict '999999'
if npr[j]=0
s=0
go j
for l=1 to k
ppp="p"+ltrim(str(l,3))
x[l]=&ppp
s+=owe[l]*(e[1,l]-x[l])^2
next
a=1
for m=1 to g
p=0
for l=1 to k
p+=owe[l]*(e[m,l]-x[l])^2
next
if s>p
s=p
a=m
endi
next
if sqrt(s)>r
a=++g
if g>15
j0=int(n/2+1)
r+=r
setcolor("W/")
@ 8,3 say "Радиус мал. Пересчет с радиусом
"+ltrim(tran(r,'999999.9'))
setcolor("W+/")
rad=.T.
exit
endi
endi
epr:= ++ep[a]
go j
for c=1 to k
ppp="p"+ltrim(str(c,3))
e[a,c]=(e[a,c]*(epr-1)+&ppp)/epr
next
npr[j]=a
endi
next
if g=1
r=r/2
rad=.T.
setcolor("W/")
@ 8,3 clea to 8,54
@ 8,3 say "Радиус велик. Пересчет с радиусом
"+ltrim(tran(r,'999999.9'))
setcolor("W+/")
endi
endd
setcolor("W/")
@ 8,3 clea to 22,54
@ 7,5 say "Классификация завершена. Вывод
результатов ..."
outputfile("BUN")
set devi to prin
heading()
@ 3,1 say "Кластерный анализ (Области сгущения
Уолтца). "+dtvr()
@ 5,1 say fldat()
@ prow()+1,1 say "Областей сгущения "+tran(g,'99')+",
каждая радиусом до "+ltrim(tran(r,'999999.9'))
for i=1 to g
@ prow()+2,1 say "Область номер "+tran(i,'99')+"
включает объекты :"
@ prow()+1,0 say ""
t=0
for j=1 to n
if npr[j]=i
yp[++t]=j
s=0
go j
for l=1 to k
ppp="p"+ltrim(str(l,3))
s+=owe[l]*(e[i,l]-&ppp)^2
next
y[j]=sqrt(s)
endi
next
if t>1
for j=1 to t-1
for l=j+1 to t
if y[yp[j]]>y[yp[l]]
a=yp[j]
yp[j]=yp[l]
yp[l]=a
endi
next
next
endi
for j=1 to t
if j%20=0
@ prow()+1,0 say ""
endi
@ prow(),pcol()+1 say ltrim(tran(yp[j],'99999'))
next
@ prow()+1,1 say "Всего объектов
"+tran(t,'999')+tran(t*100/n," (999.99%)")
next
@ prow()+2,1 say "Объекты ранжированы по расстоянию
до локального центра."
t1:=t:=g
if t>10
t1=10
endi
@ prow()+2,1 say "Расстояния между центрами сгущений
(нулевой - глобальный центр)."
for i=1 to k
s=0
for j=1 to t
s+=e[j,i]
next
x[i]=s/t
next
@ prow()+2,0 say ""
for i=0 to t1-1
@ prow(),i*11+10 say tran(i,'99-й')
next
@ prow()+1,0 say ""
for i=1 to t1
@ prow()+1,1 say tran(i,'99-й')
s=0
for l=1 to k
s+=(e[i,l]-x[l])^2
next
s=sqrt(s)
@ prow(),6 say round(s,2) pict '9999999.99'
for j=1 to i-1
s=0
for l=1 to k
s+=(e[i,l]-e[j,l])^2
next
s=sqrt(s)
@ prow(),6+j*11 say round(s,2) pict '9999999.99'
next
next
if t>t1
@ prow()+2,11 say "0-й"
for i=1 to t-t1-1
@ prow(),10+i*11 say tran(i+t1,'99-й')
next
@ prow()+1,0 say ""
for i=t1+1 to t
@ prow()+1,1 say tran(i,'99-й')
s=0
for l=1 to k
s+=(e[i,l]-x[l])^2
next
s=sqrt(s)
@ prow(),6 say round(s,2) pict '9999999.99'
for j=t1+1 to i-1
s=0
for l=1 to k
s+=(e[i,l]-e[j,l])^2
next
s=sqrt(s)
@ prow(),6+(j-t1)*11 say round(s,2) pict '9999999.99'
next
next
@ prow()+2,0 say ""
for i=1 to t-t1
@ prow(),10+(i-1)*11 say tran(i+t1,'99-й')
next
@ prow()+1,0 say ""
for i=1 to t1
@ prow()+1,1 say tran(i,'99-й')
for j=t1+1 to t
s=0
for l=1 to k
s+=(e[i,l]-e[j,l])^2
next
s=sqrt(s)
@ prow(),6+(j-t1-1)*11 say round(s,2) pict '9999999.99'
next
next
endi
@ prow()+3,1 say "См. И.Мандель. Кластерный анализ. - М.:
ФиС, 1988."
@ prow()+2,1 say "(c) "+rel_year+" "+grp_+panda
@ prow()+1,1 say ""
set devi to scre
@ 7,3 clea to 22,54
@ 7,3 say sdelan_+rezu+"BUN"
if psound
tone(1046*2,1)
endi
@ 15,3 say soh
@ 15,41 prom " Да "
@ 15,47 prom " Нет "
helper=612
sp=1
menu to sp
if sp=1
setcolor("W+*/")
@ 17,3 say ""
setcolor("W/")
@ 17,5 say "Пpизнак сохpаняется ..."
* добавление признака СГУЩЕН
use (fname+".pdd")
appe blan
repl prnam with "СГУЩЕН", high with t, low with 1, missi with 0, scala
with 1
use (fname+".pmf")
copy stru extended to tmpstru
use tmpstru
appe blan
ppp="p"+ltrim(str(k+1,3))
repl Field_name with ppp, Field_type with "N", Field_len with 3, Field_dec with
0
use
crea tmpbody from tmpstru
erase tmpstru.dbf
use tmpbody
appe from (fname+".pmf")
for i=1 to n
go i
repl &ppp with npr[i]
next
use
erase (fname+".pmf")
rename tmpbody.dbf to (fname+".pmf")
setcolor("BG/")
@ 15,3 clea to 17,54
@ 15,3 say "В файл "+alltrim(fname)+" добавлен признак
СГУЩЕН,"
@ 16,3 say opis
@ 17,3 say ispo
setcolor("W/")
endi
release npr,y,yp,e,owe,ep,x
endi
if w2=4 && ____________ КРАБ _____________
set curs off
* [t,1] - No.перв объекта, 2]-второго, 3]-расст между ними
priv shway[n-1,3], x[k]
crea struct_
appe blan
repl Field_name with "fir", Field_type with "N", Field_len with 4,
Field_dec with 0
appe blan
repl Field_name with "sec", Field_type with "N", Field_len with 4,
Field_dec with 0
appe blan
repl Field_name with "dis", Field_type with "N", Field_len with 15,
Field_dec with 4
use
sele 2
crea dista from struct_
erase struct_.dbf
sele 1
use (fname+".pmf")
setcolor("W*/")
@ 12,3 say ""
setcolor("W/")
@ 12,5 say "Вычисляются расстояния ... "
colpo=col()
@ 13,3 say repl('_',50)
@ 15,3 say vrzavoc_
colpo1=col()
zamer=seconds()
for i=2 to n
show_time()
@ 13,3 say subs(done,1,i*50/n)
@ 12,colpo say tran(n-i+1,'99999')
go i
for t=1 to k
ppp="p"+ltrim(str(t,3))
x[t]=&ppp
next
for j=1 to i-1
go j
rasst=0
for t=1 to k
ppp="p"+ltrim(str(t,3))
rasst+=(&ppp-x[t])^2
next
sele 2
appe blan
repl fir with i, sec with j, dis with sqrt(rasst)
sele 1
next
if i=2
elapsed((n-1)*n/2+1,15,colpo1)
endi
next
* определение минимального расстояния
sele 2
go top
md=999999999
do whil !eof()
show_time()
if md>dis
md=dis
mf=fir
ms=sec
endi
skip
endd
@ 12,5 say "Определяется кратчайший путь (КНП) ...
"
colpo=col()
@ 13,3 say repl('_',50)
@ 15,colpo1 say " "
zamer=seconds()
np=0
do whil np<n-1
show_time()
dele for fir=mf .and. sec=ms
pack
np++
@ 13,3 say subs(done,1,np*50/(n-1))
@ 12,colpo say tran(n-1-np,'99999')
shway[np,1]=mf
shway[np,2]=ms
shway[np,3]=md
if np=n-1
exit
endi
md=999999999
for j=1 to np
sj1=shway[j,1]
sj2=shway[j,2]
loca for fir=sj1 .or. fir=sj2 .or. sec=sj1 .or. sec=sj2
do whil found()
if md>dis
if ascan(shway,{|x| x[1]=fir .or. x[2]=fir})=0 .or. ascan(shway,{|x| x[1]=sec .or.
x[2]=sec})=0
md=dis
mf=fir
ms=sec
endi
endi
cont
endd
next
if np=1
elapsed((n-1)*n/2-1,15,colpo1)
endi
endd
erase dista.dbf
@ 12,5 clea to 12,54
@ 13,3 clea to 22,54
if psound
tone(1046*2,1)
endi
hmc=int(n/5)+1
set curs on
helper=641
@ 12,5 say "Задайте количество кластеров" get hmc pict
'99' range 1,n-1
read
set curs off
@ 12,5 clea to 22,54
@ 12,5 say "Сортировка ..."
asort(shway,,, { |x,y| x[3] > y[3] })
for i=1 to hmc-1
shway[i,1]=-shway[i,1]
shway[i,2]=-shway[i,2]
next
* до этого места абсолютно верно
np=0
release x
priv defi[n]
afill(defi,0)
for i=1 to hmc-1
for g=1 to 2
a=inke()
if a=-9
way_out()
endi
np++
if np>hmc
exit
endi
fc=abs(shway[i,g])
if fc=0
np--
loop
endi
shway[i,g]=0
defi[fc]=np
* определение точек, непосредственно связанных с
данной
for j=1 to n-1
if shway[j,1]=fc
shway[j,1]=0
defi[shway[j,2]]=np
shway[j,2]=0
endi
if shway[j,2]=fc
shway[j,2]=0
defi[shway[j,1]]=np
shway[j,1]=0
endi
next
* определение всех остальных
qu=1
do whil qu>0
qu=0
for j=1 to n
if defi[j]>0
for l=1 to n-1
if abs(shway[l,1])=j
if shway[l,1]=j
defi[shway[l,2]]=defi[j]
shway[l,2]=0
qu++
endi
shway[l,1]=0
endi
if abs(shway[l,2])=j
if shway[l,2]=j
defi[shway[l,1]]=defi[j]
shway[l,1]=0
qu++
endi
shway[l,2]=0
endi
next
endi
next
endd
next
next
for i=1 to n
if defi[i]=0
defi[i]=hmc
endi
next
@ 12,5 say 'Происходит вывод результатов в файл ...'
outputfile("KRB")
set devi to prin
heading()
@ 3,1 say "Классификация методом кратчайшего
незамкнутого пути (Краб)."
@ 4,1 say dtvr()
@ 6,1 say fldat()
@ prow()+1,1 say "Число кластеров "+tran(hmc,'99')
for i=1 to hmc
@ prow()+2,1 say "+------------------+"
@ prow()+1,1 say tran(i,'¦ Кластер номер 99 ¦')
@ prow()+1,1 say "+------------------+"
l=0
@ prow()+1,0 say ""
for j=1 to n
if defi[j]=i
l++
if l%20=0
@ prow()+1,0 say ""
endi
@ prow(),pcol()+1 say alltrim(tran(j,'9999'))
endi
next
@ prow()+1,1 say "Всего объектов "+tran(l,'9999')+", или
"+tran(l/n*100,'999.9 %')
next
@ prow()+3,1 say "См. И.Мандель. Кластерный анализ. - М.:
ФиС, 1988."
@ prow()+1,1 say " Дж.Ту, Р.Гонсалес. Принципы
распознавания образов. - М.: Мир, 1978."
@ prow()+2,1 say "(c) "+rel_year+" "+grp_+panda
@ prow()+1,1 say ""
set devi to scre
@ 12,3 say sdelan_+rezu+"KRB "
if psound
tone(1046*2,1)
endi
@ 15,3 say soh
@ 15,41 prom " Да "
@ 15,47 prom " Нет "
helper=612
sp=1
menu to sp
if sp=1
setcolor("W+*/")
@ 17,3 say ""
setcolor("W/")
@ 17,5 say "Пpизнак сохpаняется ..."
* добавление признака КРАБ
use (fname+".pdd")
appe blan
repl prnam with "КРАБ", high with t, low with 1, missi with 0, scala with 1
use (fname+".pmf")
copy stru extended to tmpstru
use tmpstru
appe blan
ppp="p"+ltrim(str(k+1,3))
repl Field_name with ppp, Field_type with "N", Field_len with 3, Field_dec with
0
use
crea tmpbody from tmpstru
erase tmpstru.dbf
use tmpbody
appe from (fname+".pmf")
for i=1 to n
go i
repl &ppp with defi[i]
next
use
erase (fname+".pmf")
rename tmpbody.dbf to (fname+".pmf")
setcolor("BG/")
@ 15,3 clea to 17,54
@ 15,3 say "В файл "+alltrim(fname)+" добавлен признак
КРАБ,"
@ 16,3 say opis
@ 17,3 say ispo
setcolor("W/")
endi
release defi, shway
endi
if w2=5 && ____________ ISODATA
endi
toto()
naj()
if a=-9 .or. a=301
way_out()
endi
endi
else
set curs off
filenotfound()
endi
set printer to
endd
retu
******************************************
Конец процедуры PAS_W61
Процедура PAS_W62
* Регpессия ****************
w2=1
do whil w2#5
helper=6200
show_time()
sw_F1(1)
sw_F10(1)
sw_Esc(1)
vybmenu(2)
@ 14,57 prom " Линейная регрессия "
@ 15,57 prom " Пошаговая регрессия "
@ 16,57 prom " На главные факторы "
@ 17,57 prom " Нелинейная "
@ 18,57 prom " Выход "
menu to w2
if w2=0 .or. w2=5
exit
endi
set colo to /W
@ 24,0 say repl(' ',80)
set colo to W/
qu=adir("*.pdd")
if qu>0
do sw_Esc with 0
do sw_F1 with 0
k=0
fname=' '
definefile(@k,@fname)
if !empty(fname)
decl prnames[k], hibond[k], lobond[k], misscode[k], scale[k]
go top
np=0
do whil !eof()
np=np+1
prnames[np]=prnam
hibond[np]=high
lobond[np]=low
misscode[np]=missi
scale[np]=scala
skip
endd
@ 2,3 clea to 22,54
use (fname+".pmf")
n=recc()
@ 2,3 say fldat()
if w2=1 && ____________ linear _____________________________
set cursor off
@ 4,3 say "В"+yberite_+"зависимый признак ( , затем
-+ )"
dep=1
setcolor("/W")
@ 6,20 say tran(dep,'999')+": "+prnames[dep]
ye=1
do whil ye#13
ye=inke(0)
if ye=-9
do way_out
endi
if ye=5 .and. dep<k
dep=dep+1
endi
if ye=24 .and. dep>1
dep=dep-1
endi
@ 6,20 say tran(dep,'999')+": "+prnames[dep]
endd
priv realch[k-1],ch[k-1]
a=0
for y=1 to k
if y#dep
ch[++a]=y
endi
next
setcolor("W/")
@ 8,3 say "Отметьте независимые
признаки"+iif(k>50," (не более 50 !)","")
chcnt=0
markvars(k-1,0,50,10,@ch,@realch,@chcnt)
m=chcnt
priv indep[m],am[m+1],xtx[m+1,m+1],x[n],y[n]
for i=1 to m
indep[i]=realch[i]
next
release realch, ch
setcolor("W+*/")
@ 11,3 say ''
setcolor("W/")
@ 11,5 say 'Идет счет Методом Наименьших Квадратов ...'
@ 12,3 say repl('_',50)
@ 14,5 say "перемножение матриц, осталось "
colp=col()
ppp="p"+ltrim(str(dep,3))
mean_y=0.00000000
for l=1 to n
go l
y[l]=&ppp
mean_y+=&ppp
next
am[1]=mean_y
crm=0
* вычисляется только нижний треугольник и
диагональ
xtx[1,1]=n
for i=1 to m
ppp="p"+ltrim(str(indep[i],3))
z=0
zz1=0
for l=1 to n
go l
x[l]=&ppp
z+=x[l]*y[l]
zz1+=x[l]
next
am[i+1]=z
xtx[1,i+1]=zz1
xtx[i+1,1]=zz1
for j=1 to i
crm++
@ 14,colp say m*(m+1)/2-crm pict '9999'
ppp="p"+ltrim(str(indep[j],3))
s=0
for l=1 to n
go l
s+=x[l]*&ppp
next
xtx[i+1,j+1]=s
next
show_time()
a=inke()
if a=-9
way_out()
endi
next
@ 12,3 say subs(done,1,25)
@ 14,5 clea to 18,54
@ 14,5 say "обращение матрицы, осталось "
colp=col()
* обращение xtx (заполняется и верхний треугольник)
right=1
f_=m+1
priv h[f_]
for l=f_ to 1 step -1
p=xtx[1,1]
if p<=0
right=0
if psound
tone(100,10)
endi
@ 11,3 say ' '
setcolor("/W*")
@ 24,3 say ""
setcolor("/W")
@ 24,5 say "Нажмите любую клавишу (АНАЛИЗ
ПРЕКРАЩАЕТСЯ - НЕВЕРНЫЕ ДАННЫЕ)"
setcolor("W+/R")
sscr=savescreen(3,13,18,64)
dispbox(3,13,18,64,'+-+¦+-+¦ ')
@ 4,23 say "матрица X'X (модель регрессии Y=XA) не"
@ 5,23 say "удовлетворяет условиям применения
МНК."
@ 6,23 say "Если это не результат плохой
обусловлен-"
@ 7,23 say "ности матрицы вследствие погрешности
вы-"
@ 8,23 say "числений чисел с плавающей точкой,
то"
@ 9,23 say "значит столбцы матрицы X (признаки)
яв-"
@ 10,23 say "ляются линейно зависимыми."
@ 12,23 say "Проверьте ваши НЕЗАВИСИМЫЕ
(предсказы-"
@ 13,23 say "вающие) признаки и оставьте только
дей-"
@ 14,23 say "ствительно независимые."
@ 16,23 say "Подробнее о данной проблеме см. в
Руко-"
@ 17,23 say "водстве по использованию пакета."
setcolor("W+*/R")
@ 4,15 say "Ошибка:"
a=inke(2)
setcolor("W+/R")
@ 4,15 say "Ошибка:"
a=inke(0)
restscreen(3,13,18,64,sscr)
exit
endi
for i=2 to f_
q=xtx[i,1]
h[i]=iif(i>l,q/p,-q/p)
for j=2 to i
xtx[j-1,i-1]:=xtx[i-1,j-1]:=xtx[i,j]+q*h[j]
next
next
xtx[f_,f_]=1/p
for i=2 to f_
xtx[i-1,f_]:=xtx[f_,i-1]:=h[i]
next
@ 14,colp say l-1 pict '999'
show_time()
next
if right=1 && если не было ошибок счета
@ 12,3 say subs(done,1,40)
@ 14,5 clea to 18,54
@ 14,5 say "вычисление коэффициентов, осталось "
colp=col()
x=aclone(am)
for i=1 to m+1
z=0
for j=1 to m+1
z+=xtx[i,j]*x[j]
next
am[i]=z
@ 14,colp say f_-i pict '999'
next
@ 12,3 say subs(done,1,50)
@ 13,3 clea to 22,54
if psound
tone(1046*2,1)
endi
set colo to W/
@ 11,5 say 'Происходит вывод результатов в файл ... '
outputfile("REG")
set devi to prin
heading()
@ 3,1 say "Множественная линейная регрессия.
"+dtvr()
@ 4,1 say fldat()
@ 6,1 say "Уравнение регрессии имеет следующий вид
:"
@ prow()+2,5 say trim(prnames[dep])+" = "+alltrim(str(am[1],12,4))
for u=2 to m+1
if u%5=0
@ prow()+1,len(trim(prnames[dep]))+8 say ""
endi
sg=iif(am[u]<0,"-","+")
@ prow(),pcol() say
sg+alltrim(str(abs(am[u]),12,4))+"*("+trim(prnames[indep[u-1]])+")"
next
set devi to scre
@ 12,3 say repl('_',5)
x=array(n)
ssr=0.00000000
ssd=0.00000000
for i=1 to n
ss=am[1]
go i
for j=2 to f_
ppp="p"+ltrim(str(indep[j-1],3))
ss+=(am[j]*&ppp)
next
x[i]=ss && предсказанное значение (оценка по
регрессии)
ssr+=(y[i]-x[i])**2
ssd+=x[i]**2
next
@ 12,3 say repl("_",20)
set devi to prin
ssd-=(mean_y*mean_y/n)
sst=ssd+ssr
prb=0
fotn=ssd/ssr*(n-m-1)/m
f_distr(fotn,m,n-m-1,@prb)
@ prow()+3,15 say "Сумма квадратов Степ.свободы
Дисперсия F-отношение P-значение"
@ prow()+2,1 say "Регрессия "+str(ssd,15,3)+"
"+str(m,5)+" "+tran(ssd/m,'99999999.999')+"
"+iif(ssr>0.001,tran(fotn,'999999.999')," >>")+"
"+tran(prb,'999999.999')
@ prow()+1,1 say "Остатки "+str(ssr,15,3)+"
"+str(n-m-1,5)+" "+tran(ssr/(n-m-1),'99999999.999')
@ prow()+1,1 say "Исх.данные "+str(sst,15,3)+" "+str(n-1,5)
@ prow()+2,1 say "Доля дисперсии "+trim(prnames[dep])+",
объясненная регрессией = "+ltrim(tran(ssd/sst*100,'999.9 %'))
if n>m+1
if ssr#0
@ prow()+2,1 say ' Проверка нулевой гипотезы'
@ prow()+1,1 say ' Доверит. границы (95%-е) "признак не
улучшает предсказание"'
@ prow()+1,1 say "Признак Коэффициент Ст.ошибка
нижняя верхняя Т-статистика Р-значение"
@ prow()+1,1 say ""
for i=1 to f_
sai=sqrt(xtx[i,i]*ssr/(n-f_))
@ prow()+1,1 say if(i>1,prnames[indep[i-1]]+"
","Св.член")+" "+tran(am[i],'9999999.9999')+" "
@ prow(),pcol() say tran(sai,'999999.999')+"
"+tran(am[i]-sai*inv_st5(n-f_),'9999999.9999')+"
"+tran(am[i]+sai*inv_st5(max(n-f_,1)),'9999999.9999')
prb=0
if sai>0
qam=abs(am[i])/sai
student(qam,n-f_,@prb)
@ prow(),pcol()+7 say if(qam>9999999," >>>
",tran(qam,'9999999.9999'))+" "+tran(prb,'99.9999')
else
@ prow(),pcol()+7 say " >> 0.0000"
endi
next
endi
else
@ prow()+2,1 say "В модели регрессии число объектов
должно превосходить число независимых признаков
по кр.мере на 2."
endi
n1:=x1:=x[1] && расчетное значение
n2:=x2:=y[1]-x[1] && остаток
for i=1 to n
y[i]-=x[i] && в y[] - остатки
if n1>x[i] && n1=min x[]
n1=x[i]
endi
if x1<x[i] && x1=max x[]
x1=x[i]
endi
if n2>y[i] && n2=min y[]
n2=y[i]
endi
if x2<y[i] && x2=max y[]
x2=y[i]
endi
next
set devi to scre
@ 12,3 say repl("_",35)
set devi to prin
il=9/(x2-n2)
@ prow()+3,1 say "Гистограмма остатков. Длина
интервала "+ltrim(str(1/il,10,2))
@ prow()+2,1 say "Набл. Н.Гр."
release all like h
priv h[10]
afill(h,0)
for i=1 to n
h[(y[i]-n2)*il+1]++
next
for i=1 to 10
e=h[i]
@ prow()+1,1 say tran(e,'9999 ')+tran((n2*(10-i)+x2*(i-1))/9,'999999.9')+" ¦
"+repl("*",if(e>50,50,e))+iif(e>50," >>
","")
next
set devi to scre
@ 12,3 say repl("_",45)
set devi to prin
@ prow()+3,1 say "Зависимость остатков (отложены по
вертикали) от :"
@ prow()+2,17 say "Расчетных
значений"+spac(46)+"Положения (времени)"
h1=24/(x1-n1)
h2=14/(x2-n2)
h4=1/h2
release all like xtx
priv xt1[15,25], xt2[15,25]
for i=1 to 15
for j=1 to 25
xt1[i,j]:=xt2[i,j]:=0
next
next
for i=1 to n
xt1[(y[i]-n2)*h2+1,(x[i]-n1)*h1+1]++
xt2[(y[i]-n2)*h2+1,(i-1)*24/(n-1)+1]++
next
rpl="+------------+-----------+------------+------------+"
rpl=rpl+" "+rpl
@ prow()+1,1 say rpl
for i=15 to 1 step -1
@ prow()+1,1 say "|"
for j=1 to 25
if xt1[i,j]>0
@ prow(),j*2 say symprint(xt1[i,j])
endi
next
@ prow(),52 say "| "+tran(n2+(i-1)*h4,'999999.999')+" |"
for j=1 to 25
if xt2[i,j]>0
@ prow(),j*2+65 say symprint(xt2[i,j])
endi
next
@ prow(),116 say "|"
next
@ prow()+1,1 say rpl
@ prow()+1,0 say ltrim(tran(n1,'99999.99'))
@ prow(),9 say tran((x1+3*n1)/4,'99999.99')
@ prow(),21 say tran((x1+n1)/2,'99999.99')
@ prow(),34 say tran((x1*3+n1)/4,'99999.99')
@ prow(),47 say tran(x1,'99999.99')
@ prow(),66 say "1"
@ prow(),75 say tran((n+3)/4,'99999')
@ prow(),87 say tran((1+n)/2,'99999')
@ prow(),100 say tran((3*n+1)/4,'99999')
@ prow(),111 say n pict '99999'
@ prow()+2,1 say "Обозначения: . - 1 объект,"
@ prow()+1,1 say " : - 2 объекта,"
@ prow()+1,1 say " + - от 3 до 10,"
@ prow()+1,1 say " * - от 11 до 100,"
@ prow()+1,1 say " # - более 100."
@ prow()+3,1 say "См. А.Афифи, С.Эйзен. Статистический
анализ. - М.: Мир, 1982."
@ prow()+1,1 say " Н.Дрейпер, Г.Смит. Прикладной
регрессионный анализ. Кн.1,2. - М.: ФиС, 1987."
@ prow()+2,1 say "(c) "+rel_year+" "+grp_+panda
@ prow()+1,1 say ""
set devi to scre
@ 12,3 say repl("_",50)
@ 11,3 say sdelan_+rezu+"REG "
if psound
tone(1046*2,1)
endi
@ 14,3 say "Сохранить остатки ?"
@ 14,31 prom " Да "
@ 14,37 prom " Нет "
sw_F1(1)
helper=622
sp=1
menu to sp
sw_F1(0)
if sp=1
setcolor("W+*/")
@ 17,3 say ""
setcolor("W/")
@ 17,5 say "Остатки сохpаняются ..."
* добавление признака ОСТххх
use (fname+".pdd")
appe blan
ppt="ОСТ"+str(dep,3)
repl prnam with ppt, high with x2, low with n2, missi with 9999, scala with 3
use (fname+".pmf")
copy stru extended to tmpstru
use tmpstru
appe blan
ppp="p"+ltrim(str(k+1,3))
repl Field_name with ppp, Field_type with "N", Field_len with 8, Field_dec with
2
use
crea tmpbody from tmpstru
erase tmpstru.dbf
use tmpbody
appe from (fname+".pmf")
for i=1 to n
go i
repl &ppp with y[i]
next
use
erase (fname+".pmf")
rename tmpbody.dbf to (fname+".pmf")
setcolor("BG/")
@ 14,3 clea to 17,54
@ 15,3 say "В файл "+alltrim(fname)+" добавлен признак
"+ppt+","
@ 16,3 say "содержащий округленные до второго
знака после"
@ 17,3 say "точки остатки признака
"+prnames[dep]+"."
setcolor("W/")
endi
release all like in*
release all like am
release all like h*
release all like y*
release all like x*
toto()
naj()
if a=-9
way_out()
endi
endi
endi
if w2=2 && ____________ Stepwise ___________________________
endi
if w2=3 && ____________ на главные компоненты ___________
endi
if w2=4 && ____________ NONlinear ___________________________
endi
endi
else
set curs off
filenotfound()
endi
set printer to
endd
retu
******************************************
Конец процедуры PAS_W62
Процедура PAS_W63
* Дисперсионный ****************
w2=1
do whil w2#4
helper=630
do show_time
do sw_F1 with 1
do sw_F10 with 1
do sw_Esc with 1
do vybmenu with 2
@ 14,57 prom " One-Way ANOVA "
@ 15,57 prom " Two-Way ANOVA "
@ 16,57 prom " MANOVA "
@ 17,57 prom " Выход "
menu to w2
if w2=0 .or. w2=4
exit
endi
set colo to /W
@ 24,0 say repl(' ',80)
set colo to W/
qu=adir("*.pdd")
if qu>0
do sw_Esc with 0
do sw_F1 with 0
k=0
fname=' '
do definefile with k,fname
if !empty(fname)
decl prnames[k], hibond[k], lobond[k], misscode[k], scale[k]
go top
np=0
do whil !eof()
np=np+1
prnames[np]=prnam
hibond[np]=high
lobond[np]=low
misscode[np]=missi
scale[np]=scala
skip
endd
@ 2,3 clea to 22,54
sele 1
use (fname+".pmf")
n=recc()
@ 2,3 say fldat()
if w2=1 && ____________ oneway ANOVA _______________
set cursor off
@ 4,3 say "Укажите факторный признак ( , затем -+ )"
dep=1
setcolor("/W")
@ 6,20 say tran(dep,'999')+": "+prnames[dep]
ye=1
do whil ye#13
ye=inke(0)
if ye=-9
do way_out
endi
if ye=5 .and. dep<k
dep++
endi
if ye=24 .and. dep>1
dep--
endi
@ 6,20 say tran(dep,'999')+": "+prnames[dep]
endd
setcolor("W/")
@ 4,3 clea to 8,54
@ 4,3 say "Укажите исследуемый признак ( , и -+ )"
iss=if(dep#1,1,k)
setcolor("/W")
@ 6,20 say tran(iss,'999')+": "+prnames[iss]
ye=1
do whil ye#13
ye=inke(0)
if ye=-9
do way_out
endi
if ye=5 .and. iss<k
iss++
if iss=dep
if iss<k
iss++
else
iss=1
endi
endi
endi
if ye=24 .and. iss>1
iss--
if iss=dep
if iss>1
iss--
else
iss=k
endi
endi
endi
@ 6,20 say tran(iss,'999')+": "+prnames[iss]
endd
setcolor("W/")
priv yy[n], y[n], f[n]
afill(yy,0)
set colo to W+*/
@ 10,3 say ''
set colo to W/
@ 10,5 say 'Сортировка значений признака '+prnames[dep]+"
"
pp1='p'+ltrim(str(dep,3))
pp2='p'+ltrim(str(iss,3))
h71=0
afill(y,0)
fmea=0
fdis=0
gmin=99999
gmax=-99999
for t=1 to n
go t
z1=&pp1
z2=&pp2
if z1#misscode[dep] .and. z1#-999.7 .and. z2#misscode[iss] .and. z2#-999.7
h71++
y[h71]=z1
fmea+=z2
fdis+=z2*z2
if gmin>z2
gmin=z2
endi
if gmax<z2
gmax=z2
endi
endi
next
fdis=(h71*fdis-fmea^2)/(h71*(h71-1))
asort(y,1,h71)
d=1
yy[1]=1
f[1]=y[1]
do show_time
@ 10,5 clea to 10,54
@ 10,5 say 'Подсчет количества различных значений'
for i=2 to h71
if y[i]#y[i-1]
if d<50
d=d+1
f[d]=y[i] && различные значения
endi
endi
yy[d]++ && сколько d-тых значений
next
show_time()
d1=d
* основной счет по группам
@ 10,5 say 'Вычисление основных статистик для групп '
priv vol[d1], mea[d1], dis[d1]
afill(vol,0)
afill(mea,0)
afill(dis,0)
if file (fname+'.pnv')
sele 2
use (fname+'.pnv')
bor=.T.
else
bor=.F.
endi
outputfile("OWA")
set devi to prin
heading()
@ 3,1 say "Однофакторный дисперсионный анализ
(One-Way ANOVA). "+dtvr()
@ 4,1 say fldat()
@ prow()+2,1 say "Группирующий признак (фактор) -
"+prnames[dep]
@ prow()+1,1 say "Исследуемый (зависимый) признак -
"+prnames[iss]
@ prow()+2,1 say grp_+"Метка Объем Средняя Ст.ошибка
95%-ные границы средней"
@ prow()+1,1 say ""
pp2='p'+ltrim(str(iss,3))
priv q1q2q3[d1,5]
for i=1 to d1
sele 1
release tem
tem:={}
for j=1 to n
go j
z=&pp1
z2=&pp2
if z=f[i] .and. z2#misscode[iss] .and. z2#-999.7
vol[i]++
mea[i]+=z2
dis[i]+=z2*z2
aadd(tem,z2)
endi
next
hm=len(tem)
asort(tem)
q1q2q3[i,4]=tem[1]
q1q2q3[i,5]=tem[hm]
q1q2q3[i,1]=(tem[(hm+3)/4]+tem[(hm+5)/4])/2
q1q2q3[i,3]=(tem[(hm*3+3)/4]+tem[(hm*3+4.9)/4])/2
q1q2q3[i,2]=(tem[(hm+1)/2]+tem[(hm+2)/2])/2
if vol[i]>1
dis[i]=(vol[i]*dis[i]-mea[i]^2)/(vol[i]*(vol[i]-1))
else
dis[i]=0
endi
mea[i]/=vol[i]
* find label
vlna=tran(f[i],' 99999.99 ')
if bor
sele 2
loca for varnum=dep .and. varval=f[i]
if !eof()
vlna=valname
endi
endi
sto=sqrt(dis[i]/vol[i])
@ prow()+1,1 say tran(i,'9999')+" "+vlna+"
"+tran(vol[i],'99999')+" "+tran(mea[i],'99999.99')+"
"+tran(sto,'99999.999')
@ prow(),pcol()+6 say tran(mea[i]-sto*1.96,'99999.99')+"
"+tran(mea[i]+sto*1.96,'99999.99')
next
fsto=sqrt(fdis/h71)
@ prow()+2,1 say "В целом"+spac(13)+tran(h71,'99999')+"
"+tran(fmea/h71,'99999.99')+" "+tran(fsto,'99999.999')
@ prow(),pcol()+6 say tran(fmea/h71-fsto*1.96,'99999.99')+"
"+tran(fmea/h71+fsto*1.96,'99999.99')
ssb=0
ssw=0
cq=0
for t=1 to d1
ssb+=vol[t]*(mea[t]-fmea/h71)^2
ssw+=dis[t]*(vol[t]-1)
cq+=vol[t]^2
next
prb=0
if ssw>0
fst=ssb*(h71-d1)/(d1-1)/ssw
f_distr(fst,d1-1,h71-d1,@prb)
else
fst=999999.99
prb=0
endi
@ prow()+3,18 say "Сумма квадратов Ст.свободы Средний
квадрат F-отношение P-значение"
@ prow()+2,1 say "Межгрупповая "+str(ssb,13,2)+"
"+str(d1-1,3)+" "+str(ssb/(d1-1),11,2)+" "+str(fst,9,2)+"
"+tran(prb,'9.9999')
@ prow()+1,1 say "Внутригрупповая "+str(ssw,13,2)+"
"+str(h71-d1,5)+" "+str(ssw/(h71-d1),11,2)
@ prow()+1,1 say "Полная "+str(ssw+ssb,13,2)+" "+str(h71-1,5)
cq=(h71-cq/h71)/(d1-1)
cq=(ssb/(d1-1)-ssw/(h71-d1))/cq
cq=abs(cq)
@ prow()+2,1 say "Оценка дисперсии распределения
случайных эффектов = "+alltrim(str(cq,11,2))
release y,yy
if prb<.1 && когда мы, возможно, отвергнем
нулевую
@ prow()+3,1 say "Множественные сравнения Шеффе
(значимые на уpовне .05 pазличия отмечены
плюсом)"
@ prow()+2,2 say "Сpедние Гpуппы"
priv sheffe[d1], pointer[d1]
acopy(mea,sheffe)
for j=1 to d1
pointer[j]=j
next
for j=1 to d1-1
for l=j+1 to d1
if sheffe[pointer[j]]>sheffe[pointer[l]]
a0=pointer[j]
pointer[j]=pointer[l]
pointer[l]=a0
endi
next
next
for i=1 to d1
@ prow(),18+i*4 say pointer[i] pict '99'
next
@ prow()+1,1 say ""
fst=0
i_f_dist(5,d1-1,h71-d1,@fst)
fst5=fst
coef=(d1-1)*ssw/(h71-d1)
coef5=sqrt(coef*fst5)
for i=1 to d1
@ prow()+1,1 say tran(mea[pointer[i]],'99999.99')+tran(pointer[i],' 99')
for j=1 to i-1
dcs=mea[pointer[j]]-mea[pointer[i]]
sc=1/vol[pointer[j]]+1/vol[pointer[i]]
what=" "
if (dcs-coef5*sc)*(dcs+coef5*sc)>0
what="+"
endi
@ prow(),19+j*4 say what
next
next
endi
* ящики-с-усами для групп
@ prow()+3,1 say "Ящики-с-усами (для значений пpизнака
"+alltrim(prnames[iss])+")"
@ prow()+2,21 say tran(gmin,'99999.99')
@ prow(),39 say tran((gmax+gmin*3)/4,'99999.99')
@ prow(),57 say tran((gmax+gmin)/2,'99999.99')
@ prow(),75 say tran((gmax*3+gmin)/4,'99999.99')
@ prow(),93 say tran(gmax,'99999.99')
@ prow()+1,26 say
"+-----------------+-----------------+-----------------+-----------------+"
@ prow()+1,1 say ""
nor=70/(gmax-gmin)
for i=1 to d1
vlna=tran(f[i],' 99999.99 ')
if bor
sele 2
loca for varnum=dep .and. varval=f[i]
if !eof()
vlna=valname
endi
sele 1
endi
sere=repl(" ",int((q1q2q3[i,3]-q1q2q3[i,1])*nor))
te=int((q1q2q3[i,2]-q1q2q3[i,1])*nor)
if len(sere)>0
sere=stuff(sere,te,1,":")
endi
ser1=repl("-",int((q1q2q3[i,3]-q1q2q3[i,1])*nor))
@ prow()+1,26+int((q1q2q3[i,4]-gmin)*nor) say " "+repl("
",int((q1q2q3[i,1]-q1q2q3[i,4])*nor))+"+"+ser1+"+"
@ prow()+1,1 say tran(i,"Группа 99 ")+vlna
@ prow(), 26+int((q1q2q3[i,4]-gmin)*nor) say
"+"+repl("-",int((q1q2q3[i,1]-q1q2q3[i,4])*nor))+"¦"+sere+"+"+repl("-",int((q1q2q3[i,5]-q1q2q3[i,3])*nor))+"¦"
@ prow()+1,26+int((q1q2q3[i,4]-gmin)*nor) say " "+repl("
",int((q1q2q3[i,1]-q1q2q3[i,4])*nor))+"+"+ser1+"+"
next
@ prow()+2,26 say
"+-----------------+-----------------+-----------------+-----------------+"
@ prow()+1,21 say tran(gmin,'99999.99')
@ prow(),39 say tran((gmax+gmin*3)/4,'99999.99')
@ prow(),57 say tran((gmax+gmin)/2,'99999.99')
@ prow(),75 say tran((gmax*3+gmin)/4,'99999.99')
@ prow(),93 say tran(gmax,'99999.99')
@ prow()+3,1 say "Ящик содержит 50% наблюдений, медиана
обозначена двоеточием."
@ prow()+2,1 say "См. А.Афифи, С.Эйзен. Статистический
анализ. - М.: Мир, 1982."
@ prow()+2,1 say "(c) "+rel_year+" "+grp_+panda
@ prow()+1,1 say ""
set devi to scre
@ 10,3 clea to 10,54
@ 10,4 say sdelan_+rezu+"OWA"
endi
if w2=2 && ____________ twoway ANOVA _______________
endi
if w2=3 && ____________ MANOVA _______________
endi
release all like f*
release all like q*
release all like h*
release all like s*
release all like p*
do toto
do naj
if a=-9
do way_out
endi
endi
else
set curs off
do filenotfound
endi
set printer to
endd
retu
******************************************
Конец процедуры PAS_W63