Matlab  |  Mathcad  |  Maple  |  Mathematica  |  Statistica  |  Другие пакеты
Internet-класс  |  Примеры  |  Методики  |  Банк задач  |  Консультации & Форум  |  Download  |  Ссылки  |  Конкурсы
 
Пакет MATHEMATICA: Первые уроки
Казань-2001
Э.Ю.Лернер, О.А.Кашина

Вернуться на страницу <Методические разработки>
archive.gif (75 bytes) Архив пособия (397 Кб, формат RTF)

Рецензенты: доц. кафедры теории относительности и гравитации КГУ, к.ф.-м.н. Р.А. Даишев,
                    доц. кафедры геометрии КГУ, к.ф.-м.н. М.А. Малахальцев

Цель издания – выработка начальных навыков работы с пакетом MATHEMATICA.

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

Готовится к печати:

Учебное пособие Э.Ю.Лернер, О.А. Кашина. Экономическое моделирование и прогнозирование на компьютере (в рамках совместного проекта Казанского и Гиссенского университетов). Рекомендуем его всем, кто заинтересован в совместном применении пакетов MATHEMATICA и STATISTICA в социально-экономической сфере. В основе изложения – реальные и игровые задачи, в том числе и те, что рассмотрены в настоящей брошюре.

Введение

"Лиха беда – начало", – гласит народная мудрость. Действительно, (как не раз подтверждала история), от правильно сделанных первых шагов порой зависит успех всего мероприятия. Что же касается освоения пакета MATHEMATICA, то здесь, на наш взгляд, наилучшим стартом будет ознакомление с примерами применения пакета к решению некоторых практических задач. Это поможет читателю накопить ту критическую массу знаний, которая затем позволит ему более глубоко изучить пакет.

Мы не стремимся максимально охватить здесь возможности пакета и не объясняем его внутреннюю логику – для этого рекомендуем нашему читателю русскоязычные источники: [1] (а также [2] и [3]). Однако, по нашему убеждению, читать эти книги целесообразно уже тогда, когда Вы сталкиваетесь с какой-то конкретной задачей и уже имеете некоторые начальные навыки работы с пакетом. Начать же мы предлагаем с пары наглядных примеров, хорошо иллюстрирующих возможности применения пакета, и, кроме того, представляющих самостоятельный интерес.

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

В конце брошюры приведен список использованных горячих клавиш и команд пакета. Подробнее они описаны в Help’е.

Использованные шрифты несут определенную смысловую нагрузку:

"Times New Roman, курсив, разрядка" -  используется для математических терминов, знание которых необходимо для понимания материала.

"Courier New, жирный" используются для команд пакета MATHEMATICA.

"Times New Roman, курсив" – для важных понятий.

"Times New Roman, курсив, жирный" – для особо важных понятий.

"Times New Roman, жирный, но не курсив" – для обозначения опций пакета MATHEMATICA.

Мы благодарны научному сотруднику Казанского Физико-Технического Института Юрию Кандрашкину за помощь, оказанную нам в освоении пакета MATHEMATICA, а также за ценные замечания и профессору Гиссенского университета Армину Бонету (Armin Bohnet) – за приобретение пакета.

Эдуард Лернер (, Ольга Кашина ()

1. Считаем собственные доходы…

Представьте себе, что Вы положили 1 млн. рублей в банк "Золотые горы" под 100% годовых. Для простоты будем считать, что это было 1 января 2000 года. Придя через год, Вы получите в 2 раза большую сумму, – то есть на 1 января 2001 года Вы будете иметь 2 млн. рублей. Если же Вы придете не через год, а через полгода, то есть 1 июля 2000 года, то получите сумму в 1,5 раза большую первоначальной, то есть 1,5 млн. рублей. Еще через полгода она увеличится в 1,5 раза и составит 2,25 млн. рублей. А если приходить в банк каждый месяц? 1 февраля 2001 года сумма Вашего вклада составит 1 млн., увеличенный в раз, то есть млн. рублей. Нетрудно видеть, что каждый месяц мы будем увеличивать сумму вклада в раз. Таким образом, 1 января 2001 года сумма вклада составит , то есть, как нетрудно посчитать на калькуляторе, примерно 2 млн. 613 305 рублей. Вы видите, что это значительно больше, чем сумма, полученная в первом или во втором случае. Не ошибка ли это? Ведь при работе с дробями калькулятор производит округления. Проверим наши расчеты в пакете MATHEMATICA.

Мы считаем, что пакет MATHEMATICA версии 3.0 или выше уже установлен на Вашем компьютере. Если ее фирменного значка еще нет на Вашем рабочем столе, найдите файл Mathematica.exe с соответствующей пиктограммой и перетащите ее на рабочий стол. Двойной щелчок по пиктограмме – и пакет запущен. Перед Вами ненадолго появляется заставка, а затем – чистый лист (Notebook, записная книжка) и меню с панелью инструментов наверху.

Просто пишем на нем (13/12)^12и нажимаем Shift+Enter.

Замечание. В пакете MATHEMATICA существуют множество приемов, облегчающих ввод формул. Так, например, для возведения в степень некоторого выражения его нужно выделить и нажать Ctrl+^. Скобки рисуются сами собой и курсор переходит в верхний индекс, куда следует ввести показатель степени. Ctrl+пробел возвращает курсор в строку. Вообще, многие клавиши, будучи нажаты вместе с Ctrl, создают шаблон формулы – поэкспериментируйте.

Вы увидите следующую картину, изображенную на рис.1.1. Ну как тут не вспомнить анекдот про математика, который на вопрос “Где мы?”, заданный ему с воздушного шара заблудившимися путешественниками, немного подумав, дал абсолютно точный, но совершенно бесполезный ответ. (Догадайтесь, какой).

Рисунок 1.1. Фрагмент рабочего документа пакета MATHEMATICA.

Не удивляйтесь – перед Вами – дробь, которая является абсолютно точным значением суммы Вашего вклада.

Но не спешите обвинять пакет MATHEMATICA – Вы задали вопрос в терминах целых чисел – в них же и выражен ответ. Если бы после числа 13 Вы поставили точку, ответ был бы в 10-ичных дробях: 2.61304, то есть примерно то же, что и было получено с помощью калькулятора. Но мы помним, что это число получено в результате расчетов с округленными значениями.

Как же представить точное значение (дробь) в удобоваримом виде? В пакете MATHEMATICA есть возможность записать результаты вычислений со сколь угодно высокой степенью точности (она ограничивается лишь возможностями компьютера). Здесь нужно иметь в виду следующее: в пакете MATHEMATICA хранится протокол всех сделанных в процессе данного сеанса работы вычислений. (Как говорил один известный персонаж Ильфа и Петрова, "У меня все ходы записаны"). In[номер] обозначает ввод, Out[номер] – вывод при обращении с указанным номером. Так как мы обращались к системе лишь однажды, Out[1] содержит полученную нами дробь. Обращение: N[Out[1]] будет означать, что результат первого вычисления переводится в десятичную дробь. Точность при этом выбирается автоматически. Обращение: N[Out[1],10] будет означать, что мы хотим получить в результате 10 знаков (всего). Такой же результат был получен нами с помощью калькулятора. Как и при работе с калькулятором, мы можем использовать результаты предыдущего вычисления. Обращение к предыдущей выходной ячейке осуществляется с помощью знака процента: %. Так, ввод N[%,10] на втором шаге эквивалентен N[Out[1],10].

Замечание. Чтобы заставить несколько команд восприниматься вместе (и тем самым избежать вывода результатов выполнения (Out’ов) каждой из них), нужно разделять их символом ";" .

Мы не увидели пока преимуществ пакета MATHEMATICA перед калькулятором, кроме возможности задавать число цифр для вывода результата. Но не будем спешить с выводами.

Вернемся в наш банк "Золотые горы". Что если пересчитывать вклад каждый день, а еще лучше – каждый чаc, каждую минуту, каждую секунду? Cоздается впечатление, что если пересчет будет производиться очень часто, то через год сумма возрастет в очень много раз. Хватит ли принесенных нами мешков, чтобы унести все деньги? Прикинем: при ежесекундном пересчете вклада итоговая сумма составит 2 млн. 718 281 рубль. Хм, мы думали, будет больше. Возникает смутное подозрение – не существует ли предела этой суммы при сколь угодно частом посещении банка? Разберемся: пусть n – количество наших временных интервалов в году (дней, минут, секунд,…). Соответствующие формулы имеют вид: . Существует ли при ? Спросим об этом пакет MATHEMATICA. Сначала выясним, знает ли MATHEMATICA пределы. Это можно сделать с помощью электронного учебника, встроенного в пакет и вызываемого при помощи опции Help .

Нужно сказать, что MATHEMATICA – очень дружественный пакет. В нем существует несколько способов поиска вспомогательной информации. Так, достаточно набрать ??Limit, (или, например, только ??Li* или ??*imit), и перед нами появится список всех ключевых слов, содержащих данную подстроку. Выделим мышкой нужное нам слово и нажмем F1 – откроется соответствующая страница Help. Можно было сделать и так: HelpHelp – ввести слово Limit в предлагаемое поле и нажать Enter.

Help содержит не только описание синтаксиса обращения к функции, но и, что очень удобно, примеры ее использования. Так, мы увидим пример вычисления предела, общий вид обращения к функции вычисления предела: Limit[expr, ], описание действия оператора и (в рубрике Futher Examples) соответствующие примеры.

Замечание. В пакете MATHEMATICA стремление аргумента обозначается знаком "" (дефис и "строго больше"), бесконечность – словом Infinity. Стрелочку и бесконечность можно изобразить покрасивее, если до и после "" и "inf" нажать Esc. Если одни и те же операции вызываются неоднократно, удобно использовать входящие панели инструментов пакета (Palettes). Их вызов осуществляется через меню: FilesPalettes. (См. рис. 1.2.) В палитре 2 Basic CalculationsCalculusCommon Operations существует шаблон для задания предела, в 5 Complete Characters – знаки бесконечности и стремления.
Очень важное замечание. Для ввода с помощью шаблона лучше предварительно задать стиль входной ячейки как Standard: Cell – Convert To – StandartForm.

На рисунке 1.2 показан результат применения двух способов написания выражения для вычисления предела: ручного, требующего по словам известного физика Р. Йоста [4], "рудиментарного владения латинским и греческим алфавитами", и шаблонного. Там же изображены две инструментальные палитры – отыщите на них шаблончик для показателя степени и символ бесконечности. (Еще более наглядный пример использования шаблонов Вы увидите на странице 18.)

Рисунок 1.2. Фрагмент рабочего документа и две палитры.

Итак, сколь бы часто мы ни пересчитывали сумму нашего (извините, конечно же, Вашего) вклада, существует предел ее увеличения. Он равен числу е, свойства которого были установлены в 18 веке великим математиком Леонардом Эйлером. (Подробности – см., например, в [5]).

Здесь речь шла о 100% годовых. К сожалению, такая ставка нереальна даже для банка "Золотые горы". Предположим, что процентная ставка равна x. (Ранее x=1.) Будет ли существовать предел суммы вклада в этом случае? Опять обратимся к пакету MATHEMATICA. Результатом вычисления: Limit[(1+x/n)^n,n->Infinity] является . Таким образом, независимо от процентной ставки сумма вклада будет стабилизироваться при сколь угодно большом значении n. При каждом значении x мы будем получать свое значение предела. Это и будет та предельная сумма, которую Вы получили бы через год после оформления вклада, если бы вообще не покидали стен банка и непрерывно осуществляли пересчет вклада.

Что же делать тем вкладчикам, у кого нет под рукой компьютера с пакетом MATHEMATICA и даже нет калькулятора? Давайте сделаем для них следующее: распечатаем в пакете MATHEMATICA график предельной суммы как функции от процентной ставки x. С этим графиком они для любого значения x легко узнают предел своих финансовых возможностей. Построение графиков (функция Plot) – одна из наиболее используемых возможностей пакета MATHEMATICA. (Ее параметры описаны на стр.8.) Набираем: Plot[E^x,{x, 0,1.2}]> и получаем картинку, изображенную на рисунке 1.3.

Рисунок 1.3. График зависимости предельного годового дохода (т.е. дохода, полученного в конце года при бесконечно частом пересчете суммы вклада) от процентной ставки.

Процентная ставка изменяется здесь от 0 до 120%.

График хорош, но можно сделать его еще более понятным и легко читаемым. Для этого существует множество параметров функции Plot. Им можно присваивать значения с помощью значка "->" (rule, правило): Plot[присвоения через запятую]. В таблице 1.1 описаны некоторых из параметров функции Plot.

Таблица 1.1. Параметры функции Plot.

PlotLabel->"INCOME DEPENDING ON X RATE"
Заголовок графика
PlotRange->{0,3.5}
Диапазон изменения переменной y на графике. (По умолчанию он выбирается так, чтобы показать характер поведения функции на всем заданном интервале изменения аргумента.) Если точка y=0 не попадает в этот диапазон, ось абсцисс переносится в нижнюю или верхнюю (в зависимости от функции) его границу.
AxesLabel->{"x","E^x"}
Надписи на осях абсцисс и ординат.
Ticks->{{0.0,0.2,0.4,0.6,0.8,1.0,1.2,1.4},

{0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5}}

Разметка осей абсцисс (первый список) и ординат (второй список). Точки, не попадающие в диапазон изменения переменных x, y, не отображаются.
GridLines –>

{0.0,0.4,0.8,1.2},{0.0,0.5,1.0,1.5,2.0,2.5}}

Линии вспомогательной решетки – абсциссы ее узлов заданы компонентами первого списка, ординаты – второго. Отображаются лишь те участки линий решетки, которые попадают в диапазон изменений переменных xy.

Заметим, что значения свойств Ticks и GridLines не связаны друг с другом.

PlotStyle->{Thickness[0.009],RGBColor[1,0,0]}
Стиль линии графика. Здесь заданы параметры: толщина линии как 0.009 ширины рисунка и цвет графика (красный).
AxesStyle->{Thickness[0.005],RGBColor[0,1,0]}
Стиль осей. Здесь задан один только список параметров, он используется для обеих осей. Можно было задать разные стили для оси абсцисс и оси ординат, (создав для этого два списка параметров), но мы поленились.

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

Сохраним записную книжку с помощью опции Save меню File. При этом создается файл с расширением .nb. Пакет MATHEMATICA позволяет работать одновременно с несколькими файлами (записными книжками). При этом нумерация всех “ходов” является сквозной. Не пугайтесь, если Вы ненароком закрыли все записные книжки и видите меню MATHEMATICи парящим прямо на рабочем столе или на фоне окна какого-нибудь другого приложения наподобие улыбки Чеширского кота. Откройте новую записную книжку опцией: FileNew.

Получим наконец материализованный итог нашей работы – распечатаем график. Печать в пакете осуществляется при выборе опции Print меню File, но прежде установим стиль печати (основные параметры страницы) опцией Printing Settings меню File. Не будем перечислять параметры печати – они традиционны для Windows. В пакете MATHEMATICA существуют также так называемые стилевые таблицы (stylesheets), задающие режим печати: рабочий, черновой, режим презентации и пр. Режимы различаются размером шрифта и плотностью печати. Выбор режима осуществляется с помощью опции Printing Style Environment меню Format.

Теперь Вы можете получить протокол всего Вашего диалога с пакетом. Однако, вся "портянка", как правило, бывает не нужна. Нужно убрать из диалогового окна лишние фрагменты, написать заголовок, комментарии, привести соответствующие формулы. Обратите внимание: справа есть вертикальные линии, выделяющие на экране отдельные ячейки и их группы. Вы можете делать с ячейками все, что хотите – удалять их, свертывать, объединять, управлять их стилем – тем самым содержимое экрана может быть оформлено так, как Вы того пожелаете. Пока мы имели дело с ячейками типа Input, Output, Standard . Существуют и другие типы: заголовки (Title, Subtitle,...), текст, формулы и пр. Они выполняют функции соответствующих опций текстовых редакторов, позволяя красиво оформить содержимое записной книжки. Основные возможности работы с ячейками в пакете MATHEMATICA мы перечислили в таблице 1.2.

Таблица 1.2. Работа с ячейками

Что сделать с ячейкой Как это сделать
Выделить Подвести курсор (он изменит свою форму на тонкую стрелку с вертикальной чертой впереди) к скобке ячейки и щелкнуть мышью.
Удалить Нажать Delete. Выделенная ячейка будет удалена.
Сгруппировать Существует два режима группировки: автоматический и ручной. Выбор режима осуществляется опцией CellCell Grouping. При автоматическом режиме (Automatic Grouping) группируются входная и выходная ячейки. При ручном режиме (Manual Grouping) можно сгруппировать любое количество смежных ячеек. Для этого их следует выделить и выбрать опцию CellCell GroupingGroup Cells.
Свернуть/Развернуть

(для группы ячеек)

Двойной щелчок на соответствующей скобке.
Задать стиль ячейки FormatStyle – соответствующая опция.

Итак, мы рассмотрели один пример применения пакета MATHEMATICA. Отметим некоторые свойства пакета, на которые Вы, должно быть, уже обратили внимание:

а) Все выражения (имена функций, параметров, обозначения известных констант) пишутся с большой буквы. Например: Е (наш злополучный предел), Infinity (бесконечность), Pi (число ), GoldenRatio (отношение золотого сечения: кстати, именно таким является принятое по умолчанию соотношение длины и ширины прямоугольной области, на которой размещается график).

б) Все аргументы функций заключены в квадратные скобки. Строки заключены в кавычки. Некоторые параметры функций определяются совокупностью свойств (см., например, PlotStyle в таблице 1.1), значения которых могут задаваться в произвольном порядке.

в) Ввести одну и ту же формулу можно разными способами: набором слов с клавиатуры или с помощью шаблона соответствующей панели (Pallete). (См., например, рисунок 1.2).

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

2. Считаем чужие доходы (что еще занятнее).

Вы слышите по местному радио, как мэр Вашего городка утверждает, что средний доход его жителей составляет 1 000 условных единиц. Вы возмущены! Самым "популярным" среди Ваших многочисленных родственников, знакомых, знакомых ваших родственников и родственников Ваших знакомых является доход в 500  у.е. Хотя, чуть поразмыслив, Вы понимаете, что такое возможно. Если представители высшего света (т.е. 1% населения) получают, скажем, по 100 000 у.е., то даже при нулевом доходе остального населения (99%), средний доход будет составлять: 0.01100 000=1 000 (у.е.). Говоря языком математической статистики, наиболее распространенное значение дохода  – мода, – не совпадает со средним его значением, или математическим ожиданием.

Через некоторое время Вы получаете повышение по службе, и Ваш доход вырастает аж до 800 у.е. Вы замечаете, что половина Ваших сослуживцев Вам завидует, а половина смотрит свысока. Значит, Вы – типичный представитель среднего класса! Но доход-то у Вас гораздо ниже среднего! Это Вам не дает покоя, и Вы достаете учебники, где излагаются основные сведения из теории вероятностей и математической статистике. (Хорошо если у Вас есть книги по прикладной статистике, например, [6], – прочитайте для начала первые две главы. Если же Вы одолеете главы 9 и 10 книги [7], то восприятие последующего материала этой брошюры не составит для Вас труда.) Теперь мы уверены в том, что Вы владеете основными понятиями из теории вероятностей (непрерывные и дискретные случайные величины (с.в.), плотность и функция распределения, числовые характеристики с.в.: среднее, дисперсия, среднее квадратическое отклонение, квантили), а также из математической статистики (выборка, генеральная совокупность, выборочные оценки параметров распределения).

Итак, 800 у.е. – это медиана распределения дохода, или 50% квантиль. Перед нами три (!) различные характеристики с.в., причем все они в некотором смысле характеризуют средний ее уровень. (Аналогичная ситуация описана в [8, с.146]). Это все понятно. Но среднее значение дохода в обществе – это все равно, что средняя температура пациентов в больнице. Чтобы оценить состояние общества, нужно знать закон распределения дохода, а не только его отдельные параметры.

Опишем закон распределения дохода для Вашего городка и проанализируем его взаимосвязь с модой, медианой и математическим ожиданием.

Одним из наиболее распространенных распределений вероятностей является так называемое нормальное, или Гауссово распределение. Плотность нормального распределения с параметрами (a, s) имеет вид:

.

По нормальному закону распределены многие с.в., совершенно разные по своей физической сути. Не подойдет ли оно для описания распределения доходов?

Как известно, мат. ожидание, мода и медиана нормального распределения совпадают с параметром a. (Вы можете легко убедиться в этом, используя очевидное равенство .) Для нашего же случая мат. ожидание, мода и медиана дохода различны. Значит, доход в нашем городе явно распределен ненормально! Не будем призывать к перераспределению доходов, а займемся лучше исследованием свойств других распределений.

Как известно, нормальным является распределение суммы большого числа независимых факторов, – в этом случае говорят об аддитивном характере их воздействия на результирующую с.в. В экономике же характер совместного действия многочисленных случайных факторов нередко является не аддитивным, а мультипликативным: значение результирующего признака h, достигнутое за счет действия случайного фактора x, пропорционально текущему значению h с коэффициентом (1+x). (Так, в примере из раздела 1, сумма вклада в банке “Золотые горы” на 1 февраля 2000 года пропорциональна сумме вклада на 1 января 2000 года с коэффициентом (1+x), где x равна 1/12 ставки процента на 2000 год.) Если число случайных факторов достаточно велико, воздействие каждого из них незначительно и имеет мультипликативный характер, то результат будет иметь не нормальное, а так называемое логарифмически-нормальное (логнормальное) распределение

Плотность логнормального распределения имеет вид:  
(2.1)
Функция распределения:  
(2.2)
 
Упражнение 2.1. Докажите, что логарифм с.в., распределенной по логнормальному закону с параметрами (a, s), имеет нормальное распределение с параметрами (ln a, s). Кроме того, покажите с помощью пакета MATHEMATICA (см. указание к данному упражнению), что плотность логнормального распределения есть производная функции (2.2) и начертите ее график при значениях параметров a=800 и s =0.685568 и сравните его с графиком на рисунке 2.2.
Указание (причем очень важное): MATHEMATICA – высоко интеллектуальный пакет, умеющий дифференцировать и интегрировать аналитически заданные функции. Соответствующие команды:

D[дифференцируемая функция, аргумент],

Integrate[подинтегральная функция, {аргумент, нижний предел, верхний предел}].

Замечание. Если Вам лень набирать команды на клавиатуре, используйте шаблоны из панели: 2 Basic Calculations (см. рис. 1.2). Вписывать формулы в квадратики шаблона – работа почти ювелирная. Облегчить ее можно за счет “внедрения” подынтегральной функции в placeholder. Для этого нужно выделить формулу и нажать шаблон интеграла. Перемещение курсора между нижним и верхним пределами интегрирования осуществляется с помощью клавиш: Ctrl+% (Ctrl+5).
Другие характеристики логнормального распределения:
Среднее значение (математическое ожидание):  

(2.3)

Мода:
Медиана:
     

Из перечисленных выше формул видно, что куб медианы логнормального распределения равен произведению квадрата среднего на моду. В нашем примере это равенство выполняется приблизительно: .

Проведем статистический анализ распределения дохода.

В пакете MATHEMATICA есть отдельные модули (packages) для выполнения различных специальных вычислений. Чтобы выполнить статистические расчеты, подключим "мастера", который будет сам осуществлять поиск нужного модуля: <<Statistics`Maste`. (Обратите внимание, что здесь используются обратные апострофы – они расположены на клавиатуре под знаком "тильда".)

Замечание. Чтобы избежать синтаксических ошибок при вводе названий законов распределений, лучше использовать Help и копировать оттуда целые выражения. Удобно также не набирать длинные ключевые слова целиком – после нескольких введенных символов нажмите Ctrl+K и выберете нужный пункт из возникающего при этом списка. Однако это не срабатывает для ключевых слов отдельных модулей.

Подготовим "почву" для последующего анализа распределения дохода. Так как мы будем неоднократно обращаться к логнормальному закону, его удобно было бы как-то запомнить. Пакет MATHEMATICA предоставляет нам такую возможность – переменными в нем могут быть не только числа или строки, но и другие объекты. Найдем нужное нам распределение в общем списке понятных MATHEMATICе законов: HelpHelpAdd-onsStandard packagesStatisticsContiniousDistributions: LogNormalDistribution[mu,sigma]. Найдем значения параметров. Будем считать, что значения параметров совпадают с их точечными оценками. Для нашего примера число 800 является оценкой для медианы, то есть параметра а в формулах (2.3). Величина lna является значением параметра mu логнормального распределения в обозначениях пакета MATHEMATICA (а вовсе не значением математического ожидания логнормально распределенной с.в., как написано в Help’е 3-ей версии пакета MATHEMATICA – в 4-ой версии эта ошибка уже исправлена). Значение параметра sigma найдем из условия: . Получаем: sigma=0,685568. Запомним выбранный закон с указанными параметрами как значение некоторой переменной (назовем ее mylog):

mylog=LogNormalDistribution[Log[800], 0.685568]>

Замечание. Имя переменной в пакете MATHEMATICA может быть любой последовательностью латинских строчных или прописных букв, цифр и знака $. Строчные и прописные буквы различаются.

Теперь мы можем ответить на многие вопросы. Например, нас интересует, какой процент населения имеет доход выше среднего. Иначе говоря, нас интересует вероятность превышения случайной величиной значения х, равного среднему доходу. Из определения функции распределения с.в. () эта вероятность есть . Так как средний доход равен 1 000 у.е., запишем: 1-CDF[mylog,1000]. Здесь CDF[distribution,x] – значение функции выбранного распределения (distribution) в точке x.) Получим: 0.372406. Следовательно, доход выше среднего имеют 37,24% жителей города.

А сколько человек имеют очень маленький доход, скажем, не превышающий 400 у.е.? CDF[mylog,400]: почти 16%.

Какой процент населения имеет доход свыше 10 000? Получается, лишь чуть больше 0,1%. Это противоречит Вашим представлениям о высшем свете, – в него, по-вашему, входит каждый сотый житель, и половина из них имеет доход свыше 10 000.

Какой же доход является "проходным баллом" в высшее общество? Эта задача обратна предыдущей. Мы ищем такое значение дохода x, что вероятность иметь меньший доход составляет 99%. Иначе говоря, мы ищем корень уравнения, или 99% квантиль:

Quantile[mylog,0.99] дает 3 942 "с хвостиком". Вам же кажется, что "проходной балл" должен быть существенно выше.

Попробуем разобраться. Мы видим, что распределение дохода в Вашем городке, возможно, подчинено логнормальному закону. По крайней мере, если рассматривать лишь ту категорию населения, доход которой расположен левее медианы (условимся называть их бедными), то все полученные результаты неплохо согласуются с Вашими представлениями. Что же касается людей, чей доход превышает медиану (их будем называть богатыми), то здесь дело обстоит несколько иначе, чем Вам представляется. Богатым, похоже, логнормальный "закон не писан". Отметим, что согласно нашему определению бедные (как и богатые) составляют половину населения.

Считать чужие деньги – занятие интересное не только для нас с Вами. Еще в конце прошлого века итальянский экономист Вильфредо Парето рассматривал распределение дохода налогоплательщиков (см. [9]). Поскольку часть населения налогов не платит (ввиду того что их доход не превышает установленного значения), Парето имел дело лишь с теми значениями, которые превышают этот порог. (Позднее такие распределения были названы усеченными.) Парето обнаружил, что закон распределения высоких доходов описывается формулой:

(2.4)

Здесь параметр с обозначает минимальное значение с.в., а параметр a характеризует кривизну функции.

Итак, мы хотим описать распределение дохода среди богатого населения. Значение параметра с плотности распределения Парето, таким образом, равно медиане (800 у.е.), а значение параметра a нам нужно оценить. Получить необходимые для этого промежуточные результаты Вы можете, сделав следующие упражнения:

Упражнение 2.2. Найдите вручную формулы мат. ожидания и функции распределения Парето и проверьте результат с помощью пакета MATHEMATICA.
Указание: для нахождения математического ожидания удобно использовать функцию Mean статистического модуля пакета MATHEMATICA.
Упражнение 2.3. Для данных нашего примера убедитесь с помощью пакета MATHEMATICA в том, что средний доход бедной части населения примерно равен моде.
Указание: плотность логнормального распределения нужно умножить на 2х (т.к. бедное население составляет лишь половину общей численности) и взять от нее интеграл в пределах: 0 – 800.

Оценим значение параметра a, исходя из среднего дохода богатой части населения. Поскольку бедная часть населения имеет средний доход, примерно равный 500 у.е. (см. упражнение 2.3), общий средний доход равен 1 000 у.е., численность богатого и бедного населения одинакова, то средний доход богатых равен 1 500 у.е. Формула для мат. ожидания распределения Парето (см. упражнение 2.2) имеет вид:

  (2.5)

Отсюда получаем: .

Запомним закон Парето с нашими параметрами:

myPareto=ParetoDistribution[800,2.1486].

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

Важное замечание. Как известно, плотность любого распределения ограничивает на плоскости фигуру единичной площади. Мы используем "половину" логнормальной кривой для описания распределения доходов бедного населения. Плотность же распределения Парето перед "стыковкой" необходимо уполовинить, то есть домножить его на долю богатого населения. Тогда результирующая кривая будет ограничивать фигуру площадью Ѕ+Ѕ=1.

Вычислим значение плотности распределения Парето в медиане (x=800): PDF[myPareto,800]/2. Получаем число 0.00134288. Вычислим при x=800 значение плотности логнормального распределения:

PDF[mylog,800] есть 0.000727394, что существенно меньше значения плотности распределения Парето. Таким образом, ожидаемой "стыковки" не происходит.

Чем же это объяснить? Напомним, что у нас нет реальных данных, и все выводы о параметрах распределения дохода мы делаем на основании: а) известного математического ожидания (мы не смеем подвергать сомнению правильность утверждений мэра); б) приблизительного значения моды (полученного по Вашим оценкам); в) еще гораздо более приблизительной оценки медианы (основанной лишь на эмоциональной окраске Ваших контактов с сослуживцами). Похоже, Вы не совсем точны в своих оценках  – медиана, видимо, расположена левее, (то есть по нашей двоичной классификации Вы уже относитесь к категории богатых).

Уточним медиану распределения дохода.

Напомним, что бедняки подчиняются логнормальному закону, а богачи живут по закону Парето. В качестве исходных данных будем использовать моду и математическое ожидание. Из формул (2.3) получим выражение для медианы логнормального распределения: .

Для обозначения моды, математического ожидания, медианы и параметра s будем использовать соответственно переменные mod, mean, med, sig. Для чистоты эксперимента очистим их от возможных прежних значений: Clear[mod, mean, med, sig].

Замечание. Имена очищаемых переменных могут содержать знак *, при этом их следует заключать в кавычки. Например, Clear["m*"] очистит переменные mod, mean и med, Clear["*"] – сразу все переменные. Еще более радикальный способ начать расчеты с "чистого листа" – закрыть ядро (Kernel) пакета MATHEMATICA.

Выполним присвоения:

mod=500; mean=1000; med=mod E^sig^2.

Напомним, что в точке med графики плотностей, описывающих доход бедного (логнормальный закон) и богатого населения (половина от плотности распределения Парето) должны сомкнуться. Мы отказались от предположения, что медиана составляет 800 у.е. и пытаемся найти ее численное значение. Очевидно, что для этого достаточно найти значение переменной sig. Но сначала нужно исключить переменную . Приравняем значения плотности (2.1) и уполовиненной плотности (2.4) при x, равном значению медианы. (x=a для (2.1) и x=c для (2.4)). Получим уравнение: , которое решим относительно :

sol=Solve[1/(Sqrt[2 Pi] sig)==al/2,al].

(Синтаксис пакета требует в записи уравнений двойного равенства.)

Замечание. Знак умножения (*) в пакете MATHEMATICA можно заменить пробелом, что мы только что сделали. Вообще, что выгодно отличает язык MATHEMATICи от универсальных языков программирования – это правильное восприятие естественных нотаций, например, 2x.
 
Еще одно замечание. Мы могли избежать выписывания формул для плотностей следующим образом:

mylog=LogNormalDistribution[Log[med],sig]
myPar=ParetoDistribution[med,al]
sol=Solve[PDF[mylog,med]==0.5 PDF[myPar,med],al]
.

 
И еще одно замечание. MATHEMATICA позволяет использовать традиционное обозначение для греческих букв: (набираем: Esc–a–Esc), (Esc–b–Esc),…, (Esc-p-Esc),…. Символ обозначает известную константу (3.14….), остальные буквы могут использоваться для обозначения переменных.

Nota bene! Решение получается не численное (так как al неизвестно), а в виде подстановки, точнее – списка, хотя и состоящего из единственного элемента. (Дело в том, что MATHEMATICA не знала заранее, сколько корней будет получено.) Для упрощения изложения запишем этот первый и единственный элемент списка sol в переменную s:

s= sol[[1]]. Напомним, что al – величина переменная, поэтому при дальнейших обращениях к ней найденное только что значение будет утрачено. Чтобы зафиксировать его, введем новую переменную (alpha), в которую сделаем подстановку s. Синтаксис этой операции имеет вид:

alpha=al/.s>.

Вычислим средний доход бедной части населения (как функцию от sig). Введем переменную:

mylog=LogNormalDistribution[Log[med],sig].

Согласно указанию к упражнению 2.3 получаем:

Средний доход богатого населения отличается от общего среднего дохода на разность между последним и доходом бедного населения: meanrich=mean+(mean-meanpoor).

Таким образом, мы получили для среднего дохода богатого населения сложное выражение от sig. Вместе с тем, этот же средний доход может быть найден по формуле мат. ожидания распределения Парето (2.5), где параметр с совпадает со значением переменной med: cPareto=med; meanPareto=alpha/(alpha-1)cPareto.

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

Введем функцию f=meanPareto–meanrich. Прежде чем строить ее график, найдем область ее определения. Построим систему ограничений на значения переменной sig, исходя из экономического смысла параметров. Подключим модуль для решения алгебраических неравенств:

<<Algebra`InequalitySolve`

и запишем функцию для решения неравенств:

InequalitySolve[mod<=med&&med<=mean&&sig>0&&alpha>1,sig].

Здесь значок && означает "логическое и", то есть одновременное выполнение условий. Решением этого неравенства является интервал: . Построим на нем (исключая границы) график функции f Plot[f,{sig,0.1,Sqrt[2/Pi]-0.1}].

 

Рисунок 2.1. График функции f.

Очевидно, что уравнение f=0 имеет единственный корень.

На графике видно, что единственный ноль функции f близок к 0.5 – возьмем эту точку в качестве начального приближения к решению: solution=FindRoot[f==0,{sig,0.5}]. В отличие от функции Solve, действующей по принципу “огласите весь список, пожалуйста”, то есть выдающей весь набор найденных решений, функция FindRoot находит только один корень, к которому сходится итерационный процесс, начатый в указанной точке. Получаем значение sig, равное 0.476 011. Вычислим теперь значения параметра alpha и медианы (точки сопряжения плотностей) с помощью найденного значения sig

{alpha,med}={alpha,med}/.solution. Они соответственно равны 1.676 19 и 627.155.

Замечание. Запомните полученное значение a =1.676 19. Позднее мы еще вернемся к нему.

Построим графики обоих распределений с найденными параметрами: mylog=mylog/.solution; g1=Plot[PDF[mylog,x],{x,0,med}].

  Рисунок 2.2. График плотности логнормального распределения. Описывает распределение дохода бедной части населения.

График будет храниться в переменной g1.

myPareto=ParetoDistribution[med, alpha],

g2=Plot[0.5 PDF[myPareto,x],{x,med,1500}]

Рисунок 2.3. График плотности распределения. Парето. Описывает распределение дохода богатой части населения.

 График будет храниться в переменной g2.

Теперь – самый ответственный момент. Изобразим оба графика на одной координатной плоскости: Show[g1,g2].

Рисунок 2.4. График плотности распределения дохода всего населения Вашего города.

Ура! Графики сомкнулись (рисунок 2.4), причем линия получилась достаточно гладкой, – это означает, что производные обеих функций мало различаются в точке стыковки.

Последнее и самое приятное замечание. Построение закона распределения дохода жителей Вашего города опирается на нахождение корня уравнения f=0. Без пакета MATHEMATICA доказательства существования и единственности решения было бы делом весьма непростым. Мы же поступили подобно древним египтянам, которые вместо математического доказательства рисовали на своих папирусах чертеж (рис. 2.1) и писали "смотри!".

Итак, мы описали распределение дохода. Это позволит нам ответить на многие вопросы – причем не только математические, но и социально-экономические.

Напомним, что в начале нашего небольшого исследования мы уже пытались оценить долю населения с фиксированным диапазоном дохода. При этом мы опирались на предположения о логнормальном распределении дохода для всего населения и о том, что медиана распределения составляет 800 у.е. Поскольку полученные результаты нас не вполне удовлетворили (они противоречили нашим наблюдениям в отношении богатого населения), мы уточнили математическую модель – построили плотность распределения дохода как функцию, состоящую из двух "кусков": до медианы она совпадает с плотностью логнормального распределения, а после – с плотностью распределения Парето. Посмотрим, как согласуется с реальностью наша уточненная модель.

Для начала выясним, сколько жителей городка беднее Вас. Понятно, что это все, кто относится к категории бедных (50%) плюс некоторая часть богатого населения. Имеем:

      (2.6)

Здесь – плотность распределения Парето, – функция этого распределения (с только что вычисленными значениями параметров).

Вычислим последнее слагаемое суммы (2.6) в пакете MATHEMATICA: CDF[myPareto,800]/2. Получаем 0.167 515. Таким образом, беднее Вас 50%+16.7%, то бишь две трети населения Вашего городка!

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

(2.7)

Для x=1 000 вычислим: 0.5–0.5 CDF[myPareto, 1000]. Получаем: 0.228 735. Таким образом, доход выше среднего имеют 22,87% горожан (но никак не 37,24%, как предполагалось ранее).

Сколько же человек имеют очень маленький доход, не превышающий 400 у.е.? CDF[mylog,400] составляет 17,24%, что мало отличается от результата, полученного в начале нашего исследования (16%).

Какой процент населения имеет доход свыше 10 000? Согласно (2.7) получаем: 0.5–0.5 CDF[myPareto,10000], что составляет примерно 0,48%, то есть почти в 5 раз выше предыдущего значения. Вы были абсолютно правы в своих оценках высшего света!

Ну и наконец, выясним, какой же доход является "проходным баллом" в высшее общество? Найдем его (как значение переменной x) из условия: . Получаем:

.

Корень этого уравнения является 98%-ым квантилем распределения Парето. Quantile[myPareto,0.98] есть 6 470.91, что чуть ли не вдвое выше предыдущего значения, в истинности которого Вы сильно сомневались.

Итак, 1% населения (составляющий высшее общество Вашего городка) имеет доход не ниже 6 470 у.е.; подняв планку до 10 000 у.е., мы “отсеем” 0,52% всего населения, то есть 52% высшего света (составляющего, как уже было сказано, 1% общества). Проверим высший свет на прочность. Посмотрим, что от него останется, если поднять планку еще во столько же раз (то есть примерно в 1.55).

Выведем формулу "отсева" при изменении минимального дохода с величины x на y. Поскольку "отсев" составляет , относительная величина “отсева” равна . Используя формулу (2.7), проделаем вычисления в пакете MATHEMATICA:

z=0.5-0.5 CDF[myPareto,10000];

w=0.5-0.5 CDF[myPareto,15500]; (z-w)/z.

Получаем 0.520 32, то есть те же самые 52%.

Упражнение 2.4. Убедитесь, с помощью пакета MATHEMATICA, что относительная величина "отсева" при увеличении суммы дохода на 1% постоянна и равна параметру распределения Парето.
Указание. Положите переменную x равной некоторому значению дохода, например, 1 000. Тогда y составит 1,01 x, то есть 1 010. Проведите эксперименты для разных значений x.

Пусть N(x) – количество людей с доходом, не меньшим x. Ясно, что =. Таким образом, параметр характеризует эластичность N(x) по доходу x (подробнее см. [9]).

Нетрудно показать, что

, (2.8)

где n – общее число жителей.

Упражнение 2.5. С помощью формулы для функции распределения Парето (см. упражнение 2.2) докажите равенство (2.8).

Итак, абсолютная величина "отсева" пропорциональна N(x), а значит, обратно пропорциональна . Вспомним, что x – это доход, а =1.676 19. Значит, при поднятии планки на 1% “отсеяно” будет тем меньше людей, чем больше x. Такая "непотопляемость" богатых подтверждает известную народную мудрость "деньги к деньгам", причем при больших значениях эта тенденция проявляется ярче.

По мнению последователей Парето, значение , близкое к 1,5, характеризует состояние социального равновесия; увеличение свидетельствует о социальной нестабильности (вследствие большой поляризации общества). Таким образом, как поется в популярной некогда песне, "городок наш ничего", хотя некоторое превышение над идеальным значением все же заставляет задуматься… Заметьте, все началось с праздного, казалось бы, любопытства относительно чужих доходов, а в итоге мы получили инструмент для экономического и социологического анализа.

Литература

1. Воробьев Е.М. Введение в систему “МАТЕМАТИКА”. // М.: "Финансы и статистика", 1998.
2. Дьяконов В.П. Mathematica 4: учебный курс. // Издательский дом "Питер", 2001.
3. Аладьев В.З., Шишаков М.Л. Введение в среду пакета MATHEMATICA 2.2. // М.: Изд-во "Филин", 1997.
4. Рид М., Саймон Б. Методы современной математической физики. Т.1. Функциональный анализ.// М.: Мир, 1977.
5. Гарднер М. А ну-ка, догадайся! М.: Мир, 1984.
6. Тюрин Ю.Н., Макаров А.А. Статистический анализ данных на компьютере. // М.: "ИНФРА–М", 1998.
7. Феллер В. Введение в теорию вероятностей и ее приложения. Т.1. // М. Мир, 1984.
8. Гарднер М. Математические досуги. М.: Оникс, 1995.
9. Ланге О. Введение в эконометрику. М.: "Прогресс", 1964

Приложение. Список использованных возможностей пакета MATHEMATICA

Ключ Действие Стр.
Shift+Enter Завершение ввода строки; отправка на выполнение. 4
F1 Вызов HELP. 6
Cell Обращение к меню для работы с ячейками. 6,10
Palettes Обращение к меню для работы с палитрами. 6,7
Ctrl+^ Возведение (выделенного выражения) в степень. 4
Ctrl+пробел Возвращение курсора в основную строку. 4
Ctrl+% (Ctrl+5) Перемещение курсора между нижним и верхним пределами интегрирования. 13
Ctrl+K+подстрока Вызов списка ключевых слов, содержащих данную подстроку. 14
??подстрока Вызов списка ключевых слов, содержащих данную подстроку. В подстроке для замещения отсутствующих символов допустимо использование символа *. 6
(Esc–a–Esc) Символ . 18
(Esc–b–Esc) Символ . 18
. . . . . . . 18
(Esc–p–Esc) Символ . 18
     
. . . . . . 18
Команда Действие Стр.
Out[номер] Обращение к выходной ячейке с указанным номером. 5
% Обращение к предыдущей выходной ячейке 5
N[выражение, целое число] Преобразование выражения в 10-ичную дробь. Необязательный параметр – целое число – общее количество знаков результата. 5
Clear Очистка переменных. 17
/. Операция подстановки. 18
Limit Нахождение предела. 6
Plot Построение графика. 7,8
Show Отображает графики, являющиеся аргументами, на одной координатной плоскости. 20
D Дифференцирование. 13
Integrate Интегрирование. 13
FindRoot Нахождение корня уравнения. 19
Solve Нахождение всех корней уравнения. 17
<<Statistics`Master` Вызов модуля статистических расчетов. 13
PDF Вычисление плотности распределения. 16
CDF Вычисление функции распределения. 14
Mean Вычисление математического ожидания. 16
Quantile Вычисление квантиля распределения. 15
<<Algebra

`InequalitySolve`

Вызов модуля решения алгебраических неравенств. 19

Вернуться на страницу <Методические разработки>

Карта сайта | На первую страницу | Поиск |О проекте |Сотрудничество |
Exponenta Pro | Matlab.ru

Наши баннеры


Copyright © 2000-2003. Компания SoftLine. Все права защищены.

Дата последнего обновления информации на сайте: 15.04.03
Сайт начал работу 1.09.00

www.softline.ru

Призы для подписчиков научно-практического журнала: Exponenta Pro. Математика в приложениях