книги / Статистика и анализ геологических данных

		READ		( 5 ,1 0 0 0 )				IORD
	C	IO R D I= IO R D				+	1
		READ			X - Y	DATA		AND	P R IN T	I T OUT
	C'
	C	CALL		R E A D M (A ,N ,M ,1 0 0 ,2 )
		CALL		P R IN T M (A ,N ,M ,1 0 0 ,2 )
	C	W RITE			( 6 ,1 9 8 5 )
		CALCULATE				SUMS		FOR	LEAST	SQUARES SOLUTION
	C . . .
	C	DO 10 0 1 = 1 , IO R D I
		C ( I ) = 0 . 0
		DO	101		J - l , I O R D I
	101	B ( I , J ) = 0 . 0
		CONTINUE
	1 0 0	CONTINUE
		DO	102		I * I , N
		X P ( l ) = I . O
		DO	103		J = 2 ,IO R D I
	103	X P ( J ) = X P ( J - l ) * A ( I , 1 )
		CONTINUE
		DO	104		J - l , I O R D I
		DO	105		K = l,IO R D I
	105	B ( J ,K ) = B (J ,K ) + X P ( J ) * X P ( K >
		CONTINUE
	1 0 4	C ( J ) » C ( J ) + X P ( J ) * A ( . I ,2 )
		CONTINUE
	102	CONTINUE
	C	SOLVE			SLE
	C . . .
	C	CALL		P R IN T M (B ,IO R D I,IO R D I,2 0 ,2 0 )
		W RITE			( 6 ,1 0 0 1 )
		CALL		P R IN T M (C ,1 . IO R D I, I , 2 0 )
		W RITE			( 6 ,1 0 0 2 )
		CALL		S L E ( B ,C ,IO R D I, 2 0 , I . 0 E - 0 6 )
		CALL		PRINTM			( С , I , I O R D I , 1 ,2 0 )
	C	W RITE			( 6 ,1 0 0 3 )
		CALCULATE				PREDICTED			VALUES
	C . . .
	C	DO	106		1 = 1 ,N
		D ( I , I ) « A ( I , I )
		D ( I , 2 ) = A ( I , 2 )
		X X P = I. 0
		Y Y P = 0 .0
		DO	107		J = l,IO R D I

Y Y P = Y Y P + X X P *C (J )

X X P = X X P *D (I, I )

107CONTINUE

D ( I , 3 )= Y Y P

D ( I , 4 ) * D ( I , 2 ) - YYP

106 CONTINUE

CALL P R IN T M (D ,N .4 ,1 0 0 ,4 ) W RITE ( 6 , I 0 0 4 )

C

C . . . CALCULATE ERROR MEASURES C

S Y - 0 . 0

S Y 2 - 0 .0

S Y C - 0 .0

S Y C 2 = 0 .0

DO 108 I = I ,N

S Y = S Y + D (1 , 2 )

S Y 2 = S Y 2 + D (1 , 2 ) * D ( 1 , 2 ) S Y C = S Y C + D (1 , 3 )

S Y C 2 = S Y C 2 + D (1 , 3 ) * D ( 1 , 3 ) 108 CONTINUE

О О О

О О О

. . .	CALCULATE			P O IN TS	EOUALLY SPACED	ALONG REGRESSION L IN E
	X M IN = A (I, 1 )
	XM AX=XM IN
	DO	109	1 = 1,N
	IF ( A ( I , 1 ) . L T . X M IN ) X M IN = A (1 , 1 )
109	I F	( A ( I	, I )	.C T .	XMAX) ХЛ1АХ=А( I ,4	)
	CONTINUE
	D X = (X M A X -X M IN )/I0 0 . 0
	X = X M IN		1 = 1 ,1 0 1
	DO	110
	X X P = I. 0
	Y Y P = 0 .0		J = l,IO R D I
	DO	I I 1

Y Y P = Y Y P + X X P *C (J )

XXP=XXP*X

ИI CONTINUE E ( I , 1 )= X

E ( I , 2 )=Y Y P X=X+DX

110 CONTINUE

. . . PLOT	O R IG IN A L DATA	AND REGRESSION L IN E
CALL	P L O T E R (A ,N ,I 0	0 , E . 1 0 1 , 1 U I )

О О О

. . .

S T A T IS T IC A L A N A LYSIS

S S T = S Y 2 - S Y * S Y / F L 0 A T ( N )

S S R = S Y C 2 - S Y C * S Y C / F L 0 A T ( N )

S S D = S S T - S S R

R 2 = S S R / S S T

R = S 0 R T ( R 2 )

W R I T E ( 6 , 1 9 9 0 )

W R I T E ( 6 , 1 9 9 5 ) IO R D

W R I T E ( 6 , 2 0 0 0 ) N

W R I T E ( 6 , 2 0 0 I )

SST

W R I T E ( 6 , 2 0 0 2 )

SSR

W R I T E ( 6 , 2 0 0 3 )

SSD

W R I T E ( 6 , 2 0 0 4 ) R2

W R I Т Е ( 6 , 2 0 0 5 ) R

C A L L

E X I T

I 0 0 0

FORM AT

( 1 2 )

1 00 1

F O R M AT( / '

C O E F .

M A T R I X

OF UNKNOWN

P AR A M ET ER S I N

NORMAL E O N S ' )

1 0 0 2

F O R M A T ( / '

VECTO R

OF

C R OSSPR ODUCTS

OF

X

AND Y ' )

1 0 0 3 F O R M A T t / ' TH E P AR A M ET ER S O F T H E R E G R E S S IO N E Q U A T I O N ' )

1 0 0 4 F O R M A T ( / ' CO L I = X V A R I A B L E ' / ' C O L 2 = Y V A R I A B L E ' /

I ' CO L 3 = Y V A L U E B AS E D ON R E G R E S S IO N E Q U A T I O N ' /

2 ' COL

4

=

C O L

2

-

CO L

3 ' )

1 9 8 5

F O R M A T C / '

O R I G I N A L

X

AND

Y D A T A ' )

1 9 9 0 FO RM AT ( I H I )

OF

E Q U A T IO N

=

, 1 5 )

1 9 9 5

FORMAT

( 2 1 H 00R D E R

2 0 0 0

FO RM AT

( 2 1H0NUMBER

OF

S A M P L E S =

, 1 5 )

2 0 0 1

FO RM AT

( 2 5 H 0 T 0 T A L

SUMS

O F SOUARES

=

, F I 5 . 4 )

2 0 0 2

FO RM AT

( 3 7 H 0 S U M S

OF

SQUA RES

DUE

TO

R E G R E S S IO N

=

, F I 5 . 4 )

2 0 0 3

FORMAT

( 3 6 H 0 S U M S

OF

SQUARES

DUE

TO

D E V I A T I O N

=

. F I 5 . 4 )

2 0 0 4

FORMAT

( I9 H 0 G 0 0 D N E S S

OF

F I T

=

, F I 5 . 6 )

2 0 0 5

FORM AT

( 2 7 H O C O R R E L A T I O N

C O E F F I C I E N T

=

, F I 5 . 6 )

END

Программа 5.5. POLYD

Программа 5.5, называемая POLYD, предназначена для построения кривой полиномиальной аппроксимации по методу наименьших квадратов. В ней используется ряд программ, рас­ смотренных в гл. 4, а именно READM (программа 4.1), PRINTM (программа 4.2) и SLE (программа 4.9). Она включает также уже описанное построение графика, (программа 5.4, PLOTER). Главная программа предназначена для ввода числа М, являю­ щегося степенью полиномиального аппроксимирующего уравне­ ния, и использует для ввода программу 4.1 (READM), поэтому за ней должны следовать описания формата, предназначенные для этой программы. Программа вычисляет матрицу смешан­ ных произведений, а соответствующее матричное уравнение ре­ шается по программе SLE. Коэффициенты печатаются и исполь­ зуются для получения -массива предсказанных значений и остатков. Последние в свою очередь печатаются. Основная про­ грамма выдает также различные статистические характери­ стики, необходимые для проведения регрессионного анализа. Наконец, печатаются графики исходных данных и построенной линии аппроксимации. Программа составлена таким образом, что в нее легко ввести изменения с целью получения других статистик или вывода данных в другой форме, если это нужно для конкретной задачи. За подробным изложением вопроса нелинейной регрессии мы отсылаем читателя к книге Ли [15]. Подпрограммы нелинейной регрессии содержатся также в неко­

торых источниках вычислительных программ,	перечисленных
в приложении. Обратите особое внимание на BMD05R в системе
BMD и KWIKR8 в Computer Contribution, 28.	Как правило,
программы вычерчивания графиков, такие, как	PLOTER (про­

грамма 5.4), имеются в каждом вычислительном центре. Для этих целей используется либо АЦПУ (автоматическое цифровое печатающее устройство), либо графопостроитель.

Мы можем теперь использовать эту программу для построе­ ния нелинейной регрессии для совокупности данных, указанных в табл. 5.7 и 5.10. Построим квадратичную аппроксимацию, что позволит нам убедиться в том, что повышение степени приводит к значительному улучшению качества аппроксимации. Полино­ миальная кривая второй степени, подобранная к этим данным, изображена на фиг. 5.12. В данном случае уравнение регрессии имеет следующий вид:

YI = PO+ P JX1+ P 2XI = 122,9- 7 .9 Х .+ 0 , IX?.

В этом примере необходимые для выполнения дисперсион­ ного анализа статистики принимают следующие значения:

SST = 2 1 363,0	SSR	=20673,2
SSD = 689,8	SSPE= 126,0
R2 = 0 ,9 7	R	= 0,98

Фиг. 5.12. Кривая полиномиальной регрессии второй степени, полученная по замерам влажности, взятым из табл. 5.7 и 5.10.

Легко	убедиться,	что значения S S T и
в случае	линейной	аппроксимации, так

S S P E такие же, как и как они не содержат

оценок величин Y. Как можно было ожидать, более гибкая квадратичная кривая ближе подходит к наблюдаемым данным, чем прямая линия. Сумма квадратов отклонений относительно линии регрессии уменьшилась с 5177,8 до 689,8. Это большое уменьшение,' но такая разница не всегда столь значительна. Таблицу дисперсионного анализа можно снова расширить для проверки этой характеристики' (табл. 5.12).

Как можно увидеть из этой таблицы, сумма квадратов по­ лучается в результате вычитания суммы квадратов для линей­ ной функции S S R I и з аналогичной суммы квадратов для квад­ ратичной функции SS R2. Эта новая сумма квадратов является мерой улучшения качества аппроксимации в результате введе­ ния дополнительного члена в уравнение регрессии. В крите­ рии6 табл. 5.12 эта величина, которая обозначена через SS2-i, используется для оценки дополнительной дисперсии в регрессии. Ее значимость проверяется точно таким же образом, как и для самого уравнения регрессии. Если окончательное значение F -критерия попадает в критическую область, то добавочный член дает существенный вклад в регрессию и его следует сохра-

								Т а б л и ц а 5.12
		Дисперсионный анализ для определения значимости дополнительных
				членов в нелинейной регрессии
	Источник изменчивости				С умма	Число степе­	Средние	F -критерий
				квадратов		ней свободы	квадраты
	Линейная регрессия				s s RI	1	M S RI	M SR2/MSji,2
	Квадратичная регрессия				s s R2	2	M S R2
	Квадратичные добавки				S S 2__J	1	M S 2_ !	M S 2_ ,/ M S £ 2
	Квадратичное отклонение				s s D2	n — 3	M S D2
	Общая дисперсия				s s T	n — 1
	а	Критерии	зн ач и м ости	д л я	к вад р ати ч н ого п риближ ения.
	б	К ритерии	зн ач и м ости	д л я	проверки	к ач е ств а квад р ати чн ой аппрокси м аци и по с р а в ­
	нению	с линейной .

нить. Если же полученное значение несущественно, то дополни­ тельный член не дает значимого вклада в регрессию. Необхо­ димо отметить, что критерий а может быть существенным, в то время как критерий 6 таким не является. Это происходит потому, что критерий а предназначен для проверки гипотезы о равенстве нулю комбинации линейных и квадратичных чле­ нов. При этом линейная часть может оказаться высоко значи­ мой, в то время как вклад квадратичного члена будет очень низким. Тогда значение критерия а будет существенным благо­ даря значительному влиянию одного линейного члена. Иногда может оказаться, что значимыми являются либо один из двух, либо оба, либо ни один из двух членов.

Читатель уже мог заметить, что корреляционная зависи­ мость всегда увеличивается при добавлении новых членов полинома. Если число членов полинома достигает (п — 1), то коэффициент корреляции становится равным 1,00 независимо от степени разброса данных точек. Однако приведенные выше критерии показывают, что увеличение коэффициента корреля­ ции не имеет статистических оснований. При этом, если средние значения квадратов отклонений увеличиваются, то F-отношение, характеризующее значимость аппроксимации, уменьшается. Оценка дисперсии частично зависит от числа наблюдений,

использованных для ее вычисления, или от числа степеней сво­ боды. Последнее постоянно уменьшается по мере увеличения числа коэффициентов в уравнении регрессии. Напомним (см. гл. 3), что мы теряем в точности на каждом оцениваемом пара­ метре одну степень свободы и что коэффициенты bi полиноми­ ального уравнения являются оценками коэффициентов регрес­ сии Pi.

При выполнении определенных статистических требований процедуру проверки гипотезы о незначимости добавляемых чле­ нов можно распространить на члены более высоких степеней в уравнении полиномиальной регрессии. Если получение допол­ нительных наблюдений реально, то эти критерии можно исполь­ зовать в комбинации с критерием проверки отсутствия прибли­ жения или критерием чисто случайной компоненты. В табл. 5.13 приведена дополненная схема дисперсионного анализа для квад­ ратичной регрессии и объединенной выборки значений влаж­ ности осадка.

Т а б л и ц а 5.13

Результаты дисперсионного анализа для определения значимости квадратичной регрессии содержания воды в осадке на глубине

	Источник изменчивости	С ум м а	Число степ е­	Средние	F -критерий
		квадратов	ней свободы	квадраты

	Линейная регрессия		16185,19	1	16185,19	M SR9
						M SD,	277,83“
	Квадратичная регрессия-		20 673,24	2	10 336,62
	Квадратичные добавки		4488,05	1	4488,05	M S,	,
						M S d;	= 9* . э б
	Квадратичное отклонение		689,76	13	45,37
	Общая дисперсия	121 363,00		15
	а К в а д р ати ч н а я	регр есси я	вы соко зн ач и м а .
	<5 К в а д р ат и ч н а я	регр есси я	п риводит к	зн ач и т ел ьн о м у	улучш ен и ю	к а ч е с т в а	ап п р ок ­

си м ац ии по сравн ен и ю с линейной регресси ей .

Следуя изложенной выше схеме и используя программу 5.5, проведите дисперсионный анализ значений влажности осадка для полиномиальной регрессии третьего порядка. Является ли член третьей степени значимым? Улучшает ли регрессию добав­ ление членов более высоких степеней?

В некоторых случаях нас может интересовать не только по­ лучение оценок в заданных точках или отклонения от линии регрессии, но и ее наклон и значение, при котором этот наклон изменяется. В качестве примера задачи такого типа мы рассмот­ рим стратиграфическую последовательность по скважине, пробу­ ренной в породах нижнего палеозоя в восточной Оклахоме. Целью геологического исследования являлось восстановление условий, существовавших во время формирования отложений. Изученная толща состоит из слоев алевролита и песчаника, относительно которых допускалось, что они имеют морское про­ исхождение. Геологами было сделано предположение, что оса­ дочный бассейн постепенно наполнялся, и, по мере того как бе­ реговая линия продвигалась по направлению к местоположению скважины, мощность последовательно образовывавшихся слоев песчаника увеличивалась. Толща насчитывала тысячи слоев и было бы крайне обременительно измерять каждый из них. Вме­ сто этого была измерена мощность каждого слоя песчаника через интервал в 10 футов. Эти измерения приведены в табл. 5.14. Гео­ логу интересно знать, существует ли зависимость между мощ­ ностью отдельного слоя и общей мощностью накопленного осадка. Отметим, что суммарная мощность X измеряется в фик­ сированных точках, а мощность индивидуальных слоев рассмат­ ривается как случайная величина, что позволяет использовать регрессионную модель. Иными словами, геолог должен прове­ рить гипотезу, заключающуюся в том, что коэффициент регрес­ сии bi значимо отличается от нуля.

После того как описанным выше методом получено уравне­ ние регрессии Y i=b0 + biXi, можно оценить дисперсию относи­ тельно линии регрессии, используя величину M S D . Последнюю в свою очередь можно использовать для вычисления t-ста-

ТИСТИКИ:	* -	, Ь'		(5.29)
	yM SD/SSx
Среднее значение	квадратов,	связанное	с отклонением ( M S D ) ,
равно S S D , деленному на п — 2 степеней свободы,				как это ука­
зано в табл. 5.8. Суммы квадратов SSx			находятся	по формуле

п( £ х , ) 2

ss = 2 x >— Vl=1n	•	(5.зо)
i= i

Эта величина используется для проверки одной из гипотез:

1.Н0: Pi= 0 ,

2.H0:P i< 0 ,

3.Н0: Pi ^ 0.

в ней можно сделать незначительные изменения, которые поз­ волят вычислить величину SSx. Некоторые из результатов, ис­ пользуемых в этом критерии, указаны ниже.

Уравнение регрессии: Y= 4,25 + 0,020 X

SST =730,94		SSR =	425,18
SSD=	305,76	SSX =	1 142 450,00
R2 =	0 ,5 8	R = 0 ,7 6

Исходные данные содержат 50 наблюдений, поэтому вели­ чина SS D соответствует 48 степеням свободы. При пятипроцент­ ном уровне значимости (а = 0,05) критическое значение t-стати­ стики при v = 46 степенях свободы равно 1,68. Значение t-стати­ стики равно

У 6,37/1 142 450.00

Полученное число лежит в критической области, и поэтому мы должны отклонить гипотезу о том, что наклон прямой регрес­ сии отрицателен или равен нулю. Последовательности присуща небольшая, но вполне определенная тенденция увеличения мощ­ ности отдельных слоев.

Приведенный только что критерий является частным слу­

чаем критерия
t	bi — Pi	(5.31)

V M S D/ S S X

для проверки гипотезы, заключающейся в том, что наклон ли­

нии	регрессии имеет	некоторое заранее заданное значение Рь
В	рассмотренном	случае (критерий 5.29) значение Pi равно

нулю. Этот критерий, точнее некоторая его разновидность, имеет важное применение в анализе временных рядов. Методы ана­ лиза временных рядов основаны на предположении об отсут­ ствии тренда в изучаемых данных, т. е. что наклон линии регрес­ сии по отношению к временной оси (или оси расстояний) равен нулю. Если тренд существует, то его нужно устранить, иначе анализ временных рядов теряет свою силу. Ряды, не имеющие значительного линейного тренда, называются стационарными. Если в данных имеется устойчивый или направленный тренд, то он называется эволюционным или нестационарным.

Одним из предположений в теории линейной регрессии явля­ ется предположение о том, что дисперсия относительно линии регрессии постоянна. Это предположение можно проверить с помощью исследования остатков между данными и их оцен­ ками. Если дисперсия постоянна, то остатки образуют более или менее равномерную полосу около линии регрессии. Если же

16 З а к а з № 455

имеется постепенное изменение ширины полосы отклонений, то дисперсия не может быть постоянной. Эти два условия даже получили свои несколько устрашающие названия: гомоседастичности для постоянной дисперсии и гетероседастичности для изменяющейся дисперсии. Быстрый и не вполне точный способ определения типа изменения дисперсии относительно линии рег­ рессии состоит в построении линейной регрессии для абсолют­ ных значений отклонений. Изменение дисперсии в последова­ тельности будет проявляться как значительный наклон.

Другим допущением регрессионного анализа является сле­ дующее: отклонения от линии регрессии не коррелируются между собой. Под автокорреляцией в данном случае подразу­ мевается стремление остатков к образованию групп близких отклонений в одну и ту же сторону по отношению к линии рег­ рессии. Присутствие ряда последовательностей автокоррелированных остатков может указывать на то, что регрессионная модель не соответствует исходным данным. Может также слу­ читься, что автокоррелированные отклонения свидетельствуют о существовании явлений, представляющих геологический ин­ терес. Эти вопросы будут рассмотрены подробнее в главе о тренд-анализе, где автокоррелированные положительные остат­ ки выбираются в качестве показателей экономического потен­ циала запасов нефти и других полезных ископаемых. Проверка наличия автокорреляции производится с помощью критерия скачков, применяемого к последовательности знаков отклоне­ ний от линии регрессии, или одного из методов, рассмотренных в разделе об автокорреляции.

Интересна зависимость дисперсии и автокорреляции остат­ ков, выявленная при изучении горных разработок в северном Квебеке. На месторождении золота бульдозером была проде­ лана длинная траншея. Вдоль нее с некоторым интервалом были отобраны пробы, в которых определялось содержание золота. Тренд полученных значений был очевидным, а на одном конце траншеи были отмечены значительные отклонения от линии регрессии. Обычно это предвещает богатое месторождение зо­ лота. Именно такие месторождения часто обладают крайне низкими значениями содержаний золота в большей части мине­ рализованной зоны, но наряду с этим попадаются и богатые жилы. Кроме того, в них группами или скачками встречаются большие положительные отклонения, которые также указывают на то, что траншея пересекает зону минерализованных жил. По данным табл. 5.15 проверьте наличие тренда в значениях содержаний золота и исследуйте поведение дисперсии вдоль профиля.

Используя критерий знаков, ответьте на вопросы: можно ли считать, что отклонения распределены случайно относительно