LinuxForum.Fans
Тема

Octave и мат.статистика

Автор: Osanve 30 сообщений 9988 просмотров
21.05.2011 07:00 #163916

Здравствуйте.
Требуется сделать курсовую работу по теории вероятности и математической статистике.
Поставил Octave (на лекциях давали MathCad (под Linux вроде тоже есть, но если не ошибаюсь - платный)) и прочитал этот курс статей: http://xgu.ru/wiki/GNU_Octave.
Несмотря на это, застопорился на первом же задании, в котором для выборок надо найти оценки. В MathCad выборка осуществляется с помощью функции rnorm. В Octave её нет.
Что можно прочитать, чтобы все-таки разобраться и сделать эту курсовую? :)
Заранее спасибо.

23.05.2011 11:23 · отредактировано #164302

Что-то я не вижу в матлабе функции rnorm. Возможно тебе нужна normrnd, так она есть и в octave.

23.05.2011 19:57 #164460
Vascom написал:
в матлабе

Не в матлабе функция, а в маткаде, который нам дают на лекции :)

23.05.2011 20:05 #164462

Да, не заметил.
В любом случае Octave - замена именно Matlab, а не Mathcad.

23.05.2011 20:13 #164465

Насчет этого я в курсе, но на OpenNet было указано, что заменой MathCad является GAP. Но там опять же было несовпадение функций с MathCad, после чего мне посоветовали Octave.

23.05.2011 20:37 #164467

Хм...
Для формирования выборки использовал normrnd, далее, для нахождения выборочного среднего - mean. По идее должно было получиться одно число, т.к., грубо говоря, находим среднее арифмитическое (по формуле).
А вот то, что получилось:

Показать

При этом через MathCad получается все нормально (работа девушки):
http://dl.dropbox.com/u/20415838/%D1%82%D0%B2%D0%B8%D0%BC%D1%81.pdf
http://dl.dropbox.com/u/20415838/%E2%84%966.pdf

24.05.2011 11:01 · отредактировано #164541

Дело в том, что mean считает среднее для вектора, а ты подсунул ему матрицу 100х100, вот он для каждого столбца и подсчитал среднее.
Читай help mean

Чтобы создать вектор, команда должна выглядеть так:
normrnd(m,s,ns,nc)
где
ns - число строк
nc - число стобцов.
Если nc не указывать, то будет создана квадратная матрица ns x ns.

24.05.2011 13:54 #164583
Vascom написал:
Читай help mean

Читал прилагающуюся к Octave pdf-ку - там про вектор не заметил

Vascom написал:
Чтобы создать вектор, команда должна выглядеть так:
normrnd(m,s,ns,nc)
где
ns - число строк
nc - число стобцов.
Если nc не указывать, то будет создана квадратная матрица ns x ns.

Была такая догадка, но также была уверенность в том, что я перевел описание неправильно.

Сейчас вызвал таким образом: normrnd(m, s, n, 1); - вроде все нормально.
Благодарю.

26.05.2011 14:11 #165091

Кстати, для мат-статистики вроде люди рекомендуют R. Сам его пока особо не распробовал - не до того было.

26.05.2011 14:13 #165093

R - это язык программирования. Матлаб и Octave - несколько другое. А уж тем более Mathcad

26.05.2011 23:41 #165297
Battle Coder написал:
Кстати, для мат-статистики вроде люди рекомендуют R. Сам его пока особо не распробовал - не до того было.

В обозримом будущем мне предстоит сделать курсовик по мат.статистике и в следующем семестре у нас вычислительная математика и вроде как все, поэтому пока не вижу смысла сильно заморачиваться по этому поводу.

26.05.2011 23:54 #165298

Я сейчас не могу понять одного: стало интересно и я посмотрел какие числа заносятся в выборку.
Выборки получаются следующим образом:

text
X = normrnd(m, s, n, 1);
Y = normrnd(m, s, n, 1);

Где m = 0, s = 1, n = 180.
При этом в выборке X вроде как все нормально (на таком объеме проверить сложновато), а выборка Y состоит из 180-ти одинаковых чисел. Я думаю, что все-таки так не должно быть.
В чем ошибка?

27.05.2011 00:12 #165300

Извиняюсь. Была моя ошибка. Случайно нажал не ту кнопку на клаве, в результате чего в цикле выполнялось не вычитание, а присваивание :)

28.05.2011 16:00 #165748

И снова прошу вашей помощи.

Имеется код:

text
# -------------- [Четвертое задание] --------------
# Для X построить эмпирическую кривую плотности
# распределения, разбив интервал (X(1), X(n)) на 10
# интервалов. На этом же графике изобразить
# теоретическую кривую

function tmp = xb(u, n)
  if (u<0)
    tmp = 0;
  else
    tmp = 1/n;
  endif
endfunction


function tmp = fb(u, n, X)
  temp = 0;
  for i = 1:n
    temp += xb(u - X(i), n);
  endfor
  tmp = temp;
endfunction


# Функция для построения эмпирической кривой
# плотности распределения:
function y = fEmp(x, Sx2, Xsr, n, X)
  tmp = 0;
  for i = 1:9
    tmp += (fb(sqrt(Sx2)*0.8*(i-4)+Xsr, n, X) - fb(sqrt(Sx2)*0.8*(i-5)+Xsr, n, X))*(xb(x-(sqrt(Sx2)*0.8*(i-5)+Xsr), n)-xb(x-(sqrt(Sx2)*0.8*(i-4)+Xsr), n))*(n/(0.8*sqrt(Sx2)));
  endfor
  y = tmp;
endfunction


# Строим эмпирическую кривую плотности распределения:
Xe = 0:1:5;
Ye = fEmp(Xe, Sx2, Xsr, n, X);


# Строим теоретическую кривую плотности распределения:
Yt = normpdf(Xe, Xsr, sqrt(Sx2));


# Прорисовываем график:
plot(Xe, Yt, Xe, Ye)
grid("on")
# -------------------------------------------------

Формулы взяты из учебника. Некоторые переменные были определены выше.

Собственно проблема в том, что график, строящийся функцией normpdf, представляет собой прямую линию, а график, строящийся с помощью самописной функции fEmp вообще не выводится.

К сожалению никак не могу понять в чем проблема.

Заранее спасибо.

P.S. На всякий случай выкладываю полностью скрипт, который есть на данный момент:

Показать
text
# -------------- [Первое задание] --------------
# Для выборок X, Y найти оценки: выборочное среднее
# выборочная дисперсия, оценка центрального момента
# порядка 3 и 4, коэффициент эксцесса, коэффициент
# ассимитрии, оценка корреляционного момента, оценка
# коэффициента корреляции, размах выборки X и Y

n = 180; # Объем выборки
m = 0; # Коэффициент сдвига
s = 1; # Коэффициент масштаба
a = 0.05; # Уровень значимости


# Формируем выборки с нормальным распределением:
X = normrnd(m, s, n, 1);
Y = normrnd(m, s, n, 1);


# Находим выборочное среднее:
Xsr = mean(X)
Ysr = mean(Y)


# Находим выборочную дисперсию:
Sx2 = 0;
Sy2 = 0;
for i = 1:n
  Sx2 += X(i) - Xsr;
  Sy2 += Y(i) - Ysr;
endfor
Sx2 /= n - 1
Sy2 /= n - 1


# Находим оценку центрального момента порядка 3 и 4:
M3X = 0;
M3Y = 0;
M4X = 0;
M4Y = 0;
for i = 1:n
  M3X += (X(i) - Xsr)^3;
  M3Y += (Y(i) - Ysr)^3;
  M4X += (X(i) - Xsr)^4;
  M4Y += (Y(i) - Ysr)^4;
endfor
M3X /= n - 1
M3Y /= n - 1
M4X /= n - 1
M4Y /= n - 1


# Находим коэффициент эксцесса:
Ex = (M4X / Sx2^2) - 3
Ey = (M4Y / Sy2^2) - 3


# Находим коэффициент ассимитрии:
Sx = M3X / Sx2^(3/2)
Sy = M3Y / Sy2^(3/2)


# Находим оценку корреляционного момента:
Rxy = 0;
for i = 1:n
  Rxy += (X(i) - Xsr) * (Y(i) - Ysr);
endfor
Rxy /= n - 1


# Находим оценку коэффициента корреляции:
ro = Rxy / (Sx2^(1/2) * Sy2^(1/2))


# Находим размах выборки:
wX = X(n) - X(1)
wY = Y(n) - Y(1)
# ----------------------------------------------


# -------------- [Второе задание] --------------
# Для оценок по выбоке X построить доверительные
# доверительные интервалы для математического
# ожидания (считая дисперсию известной и
# неизвестной) и дисперсии


# Находим доверительный интервал для математического ожидания
# в случае с известной дисперсией для выборки X:
dims1 = Xsr - ((Sx2^(1/2) / n^(1/2)) * norminv(1-a/2, m, s))
dims2 = Xsr + ((Sx2^(1/2) / n^(1/2)) * norminv(1-a/2, m, s))


# Находим доверительный интервал для математического ожидания
# в случае с неизвестной дисперсией для выборки X:
dimns1 = Xsr - ((Sx2^(1/2) / n^(1/2)) * tinv(1-a/2, n-1))
dimns2 = Xsr + ((Sx2^(1/2) / n^(1/2)) * tinv(1-a/2, n-1))


# Находим доверительный интервал для дисперсии
# для выборки X:
dis1 = (n-1)*Sx2/chi2inv(1-a/2, n-1)
dis2 = (n-1)*Sx2/chi2inv(a/2, n-1)
# ----------------------------------------------


# -------------- [Третье задание] --------------
# По выборкам X, Y построить доверительный интервал
# для коэффициента корреляции


# Определяем доверительный интервал по коэффициенту
# корреляции для выборок X, Y:
z = (1/2)*log((1+Rxy)/(1-Rxy));
z1 = z - tinv(1-a/2, n) * (1/(n-3)^(1/2));
z2 = z + tinv(1-a/2, n) * (1/(n-3)^(1/2));
ro1 = tanh(z1)
ro2 = tanh(z2)
# ----------------------------------------------


# -------------- [Четвертое задание] --------------
# Для X построить эмпирическую кривую плотности
# распределения, разбив интервал (X(1), X(n)) на 10
# интервалов. На этом же графике изобразить
# теоретическую кривую


function tmp = xb(u, n)
  if (u<0)
    tmp = 0;
  else
    tmp = 1/n;
  endif
endfunction


function tmp = fb(u, n, X)
  temp = 0;
  for i = 1:n
    temp += xb(u - X(i), n);
  endfor
  tmp = temp;
endfunction


# Функция для построения эмпирической кривой
# плотности распределения:
function y = fEmp(x, Sx2, Xsr, n, X)
  tmp = 0;
  for i = 1:9
    tmp += (fb(sqrt(Sx2)*0.8*(i-4)+Xsr, n, X) - fb(sqrt(Sx2)*0.8*(i-5)+Xsr, n, X))*(xb(x-(sqrt(Sx2)*0.8*(i-5)+Xsr), n)-xb(x-(sqrt(Sx2)*0.8*(i-4)+Xsr), n))*(n/(0.8*sqrt(Sx2)));
  endfor
  y = tmp;
endfunction


# Строим эмпирическую кривую плотности распределения:
Xe = 0:1:5;
Ye = fEmp(Xe, Sx2, Xsr, n, X);


# Строим теоретическую кривую плотности распределения:
Yt = normpdf(Xe, Xsr, sqrt(Sx2));


# Прорисовываем график:
plot(Xe, Yt, Xe, Ye)
grid("on")
# -------------------------------------------------


# -------------- [Пятое задание] --------------
# ---------------------------------------------


# -------------- [Шестое задание] --------------
# ----------------------------------------------

30.05.2011 14:03 #166222

У тебя Xe = 0:1:5;
это что вообще такое? Если интервалы то их должно быть 10. А если данные для которых надо построить кривую плотности распределения, то конечно у тебя будет горизонтальная прямая с такими исходными данными.

По твоей рукописной функции ни чего не строится, потому что Ye у тебя представляет собой одну точку. Вот она и строится, да и то она нулю равна.

В общем как-то всё плохо у тебя в четвёртом задании.

31.05.2011 19:43 #166650
Vascom написал:
R - это язык программирования. Матлаб и Octave - несколько другое. А уж тем более Mathcad

Я бы не назвал R именно языком программирования. То есть, конечно, это язык программирования - но по той же логике и matlab - язык программирования.
R - среда для статистических расчётов. Логично, что у неё свой синтаксис и свой язык программирования... Есть и много оболочек для R.

01.06.2011 10:56 #166749
Vascom написал:
У тебя Xe = 0:1:5;
это что вообще такое? Если интервалы то их должно быть 10. А если данные для которых надо построить кривую плотности распределения, то конечно у тебя будет горизонтальная прямая с такими исходными данными.

А как правильно надо задать? Это данные по координате x для построения графика.

Vascom написал:
По твоей рукописной функции ни чего не строится, потому что Ye у тебя представляет собой одну точку.

Т.е. я так понимаю определение координаты y надо поставить в цикл? Или можно использовать оператор "." перед вызовом функции?

01.06.2011 11:02 #166754

Если Xe это для построения графика, то зачем ты их используешь для вычисления Ye и Yt? Я так понимаю, что для этого тебе надо использовать изначальный массив случайных величин X или какую-то функцию от него (разбитие на участки или что-то около того).

Я не вникал глубоко в твою функцию, какая-то она монстроузная, так что не могу ответить на второй вопрос.

Войдите или зарегистрируйтесь, чтобы ответить.