Помогите написать линейную регрессию - страница 5

Новый комментарий
 
Rosh писал (а) >>

А вот что выдает Excel 2007



Так что, проверять, возможно, нужно Маткад.

Если я правильно понял, то эксель выдал 3-й результат, отличный от первых двух ). А истина то где ? В совпал с маткадом, а вот коэф А другой

 
Prival писал (а) >>

Если я правильно понял, то эксель выдал 3-й результат, отличный от первых двух ). А истина то где ? В совпал с маткадом, а вот коэф А другой

Вообще при таком количестве значащих цифр и таком диапазоне вычисления идут где-то на задворках мантиссы. То есть в ответ могут привнестись даже элементы случайности. Я думаю на корректность можно рассчитывать только для какого-нибудь специального алгоритма повышенной точности. А в данном случае лучше таки сдвинуть начало координат ближе к диапазону X.



P.S. Особенно когла сумма X*X считается, просто прямо в унитаз информация идёт :)

 
lna01 писал (а) >>

Вообще при таком количестве значащих цифр и таком диапазоне вычисления идут где-то на задворках мантиссы. То есть в ответ могут привнестись даже элементы случайности. Я думаю на корректность можно рассчитывать только для какого-нибудь специального алгоритма повышенной точности. А в данном случае лучше таки сдвинуть начало координат ближе к диапазону X.

Дело в том, что начал подготовку к чемпионату. И свои наработки в маткаде начал переводить на MQL. Если помнишь мы строили АКФ (автокорреляционную функцию), я начал с неё, решил рассчитывать по прямым формулам, т.к. через преобразования Фурье сильно нагружаю процессор.

Первая попытка и сразу грабли АКФ в маткаде не совпадает с АКФ в MQL, вот и начал разбираться искать откуда грабли растут (

Попробую сместить X (Time) в точку 0. Но это опять все по новой перепроверять. И так уже от 50% идей приходиться отказываться

 
Prival писал (а) >>

Дело в том, что начал подготовку к чемпионату. И свои наработки в маткаде начал переводить на MQL. Если помнишь мы строили АКФ (автокорреляционную функцию), я начал с неё, решил рассчитывать по прямым формулам, т.к. через преобразования Фурье сильно нагружаю процессор.

Первая попытка и сразу грабли АКФ в маткаде не совпадает с АКФ в MQL, вот и начал разбираться искать откуда грабли растут (

Попробую сместить X (Time) в точку 0. Но это опять все по новой перепроверять. И так уже от 50% идей приходиться отказываться

МТ держит в мантиссе 15 цифр. Извлекаем корень, получаем 10^7. То есть возводить в квадрат, да ещё после этого суммировать числа больше 10000000 - см. постскрипт к предыдущему посту :). К счастью этот предел как раз соответствует числу минутных баров в реальной истории, так что если речь о нём, то ещё должно работать. А если это время, то сдвиг начала координат просто напрашивается.


P.S. Кстати, если ты свою функцию на чемпионате собираешься использовать, добавь защиту от деления на ноль. Иначе есть риск что твой индикатор просто встанет посредине чемпионата. Или в начале. Помнишь, была такая фишка с Фурье.

 

Тот же самый алгоритм в Java

import java.util.ArrayList;

public class Prival {
public static void main(String arg[]){
int N = 6;
double Y[];
double X[];
ArrayList<Double> Parameters = new ArrayList<Double>();
Parameters.add(0.0);
Parameters.add(0.0);
X = new double[6];
Y = new double[6];
for ( int i = 0; i < N; i ++ )
{
// массивы Y и X для проверки работоспособности
// intercept = -3.33333333 slope = 5.00000000

X[i]=i;
Y[i]=i*i;
}

LinearRegr(X, Y, N,Parameters);
System.out.println("intercept = "+Parameters.get(0)+" slope = "+ Parameters.get(1));

// вторая проверка
X[0]=1216640160;
X[1]=1216640100;
X[2]=1216640040;
X[3]=1216639980;
X[4]=1216639920;
X[5]=1216639860;

Y[0]=1.9971;
Y[1]=1.9970;
Y[2]=1.9967;
Y[3]=1.9969;
Y[4]=1.9968;
Y[5]=1.9968;


LinearRegr(X, Y, N, Parameters);
System.out.println("intercept = "+Parameters.get(0)+" slope = "+ Parameters.get(1));

}
public static void LinearRegr(double X[], double Y[], int N, ArrayList<Double> Parameters){
double sumY = 0.0, sumX = 0.0, sumXY = 0.0, sumX2 = 0.0;
double A=0,B=0;
for ( int i = 0; i < N; i ++ ){
sumY +=Y[i];
sumXY +=X[i]*Y[i];
sumX +=X[i];
sumX2 +=X[i]*X[i];
}
B=(sumXY*N-sumX*sumY)/(sumX2*N-sumX*sumX);
A=(sumY-sumX*B)/N;
Parameters.set(0, A);
Parameters.set(1, B);
}
}


Результат:


intercept = -3.3333333333333335 slope = 5.0
intercept = -1102.169141076954 slope = 9.075536028198574E-7

Process finished with exit code 0

 

Rosh

Согласен, что эти формулы дадут один и тот же результат вот маткад

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

Я сделал проверку произвел сортировку по X и снова рассчитал коэффициенты

РЕЗУЛЬТАТ поменялся !!!, а такого не должно быть. Скорее всего ошибка связана с накоплением ошибки из-за возведения в квадрат больших чисел (Candid прав). Я поковырялся в литературе и нашел более простую формулу, нет квадратов и вроде меньше вычислений.

Результат совпадает с маткадом, и от сортировки не зависит.

Рекомендую использовать эту формулу для расчетов коэффициентов линейной регрессии.

//+------------------------------------------------------------------+
//|                                                       LinReg.mq4 |
//|                                                    Привалов С.В. |
//|                                             Skype -> privalov-sv |
//+------------------------------------------------------------------+
#property copyright "Привалов С.В."
#property link      "Skype -> privalov-sv"

//+------------------------------------------------------------------+
//| script program start function                                    |
//+------------------------------------------------------------------+
int start()
  {
//----
   int      N=6;                 // Размер массива
   double   Y[],X[],A=0.0, B=0.0;
   
  ArrayResize(X,N);
  ArrayResize(Y,N);
      
// проверка 
    X[0]=1216640160;
    X[1]=1216640100;
    X[2]=1216640040;
    X[3]=1216639980;
    X[4]=1216639920;
    X[5]=1216639860;
    
    Y[0]=1.9971;
    Y[1]=1.9970;    
    Y[2]=1.9967;
    Y[3]=1.9969;    
    Y[4]=1.9968;    
    Y[5]=1.9968;
    
    
  LinearRegr(X, Y, N, A, B);
  
  Print("A = ", DoubleToStr(A,8)," B = ",DoubleToStr(B,8));
           
//----
   return(0);
  }
//+------------------------------------------------------------------+
//| Рассчет коэффициентов A и B в уравнении                          |
//| y(x)=A*x+B                                                       |
//| используються формулы https://forum.mql4.com/ru/10780/page5&nbsp;      |
//+------------------------------------------------------------------+

void LinearRegr(double X[], double Y[], int N, double& A, double& B)
{
      double mo_X = 0.0, mo_Y = 0.0, var_0 = 0.0, var_1 = 0.0;
      
    for ( int i = 0; i < N; i ++ )
      {
        mo_X +=X[i];
        mo_Y +=Y[i];
      }
    mo_X /=N;
    mo_Y /=N;
        
    for ( i = 0; i < N; i ++ )
      {
        var_0 +=(X[i]-mo_X)*(Y[i]-mo_Y);
        var_1 +=(X[i]-mo_X)*(X[i]-mo_X);
      }
        A = var_0 / var_1;
        B = mo_Y - A * mo_X;
}

Скрипт прилагаю если кто-то причешет LinearRegr (на предмет ошибок при работе на реале и увеличит быстродействие) то будет вообще хорошо. Я поменял местами A и B, т.к.

запись y(x)=a*x+b более привычна (для меня по книгам).

Файлы:
linreg_1.mq4  2 kb
 

Я не вижу каким образом результат может зависеть от сортировки. Сортировка в формулах нигде явным образом не используется.


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


Если требуются массовые вычисления линейной регрессии на ряде ценовых последовательностей, то лучше выделить отдельные буфера в индикаторе и считать методом накопления итогов. Это позволяет ускорить расчеты на порядки. Пример - Оптимизированная АМА Кауфмана : Perry Kaufman AMA optimized

 
Rosh писал (а) >>

Я не вижу каким образом результат может зависеть от сортировки. Сортировка в формулах нигде явным образом не используется.


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

Если требуются массовые вычисления линейной регрессии на ряде ценовых последовательностей, то лучше выделить отдельные буфера в индикаторе и считать методом накопления итогов. Это позволяет ускорить расчеты на порядки. Пример - Оптимизированная АМА Кауфмана : Perry Kaufman AMA optimized

1. В том то и дело что результат не должен зависеть от сортировки, а в том алгоритме зависит. Вот проверьте.

//+------------------------------------------------------------------+
//| script program start function                                    |
//+------------------------------------------------------------------+
int start()
  {
//----
   int      N=6;                 // Размер массива
   double   Y[],X[],Y1[],X1[],A=0.0, B=0.0;
   
  ArrayResize(X,N);
  ArrayResize(Y,N);
  ArrayResize(X1,N);
  ArrayResize(Y1,N);
      
// проверка 
    X[0]=1216640160;
    X[1]=1216640100;
    X[2]=1216640040;
    X[3]=1216639980;
    X[4]=1216639920;
    X[5]=1216639860;
    
    Y[0]=1.9971;
    Y[1]=1.9970;    
    Y[2]=1.9967;
    Y[3]=1.9969;    
    Y[4]=1.9968;    
    Y[5]=1.9968;
    

// отсортируем массив по возрастанию X (исходный массив был по убыванию)
  for (int i = 0; i < N; i++)
   {
   X1[i]=X[N-i-1];
   Y1[i]=Y[N-i-1];
//   Print(X[i], " ", X1[i], " ", Y[i], " ", Y1[i]);
   }            
//----
// 
  LinearRegr(X, Y, N, A, B);
  Print("A = ", DoubleToStr(A,8)," B = ",DoubleToStr(B,8));
  LinearRegr(X1, Y1, N, A, B);
  Print(" A = ", DoubleToStr(A,8)," B = ",DoubleToStr(B,8));

  LinearRegr1(X, Y, N, A, B);
  Print("A = ", DoubleToStr(A,8)," B = ",DoubleToStr(B,8));
  LinearRegr1(X1, Y1, N, A, B);
  Print(" A = ", DoubleToStr(A,8)," B = ",DoubleToStr(B,8));

   return(0);
  }

//-------------------------------------------------------------------------------
// использование этой формулы приводит к ошибкам если X=Time
// формула предложена вот тут https://forum.mql4.com/ru/10780/page4
//| y(x)=A+B*x  

void LinearRegr(double X[], double Y[], int N, double& A, double& B)
{
      double sumY = 0.0, sumX = 0.0, sumXY = 0.0, sumX2 = 0.0;
      
    for ( int i = 0; i < N; i ++ )
    {
        sumY   +=Y[i];
        sumXY  +=X[i]*Y[i];
        sumX   +=X[i];
        sumX2  +=X[i]*X[i];
    }
   B=(sumXY*N-sumX*sumY)/(sumX2*N-sumX*sumX);
   A=(sumY-sumX*B)/N;
}

//+------------------------------------------------------------------+
//| Формула предлагаемая мной                                        |
//| Рассчет коэффициентов A и B в уравнении                          |
//| y(x)=A*x+B                                                       |
//| используються формулы https://forum.mql4.com/ru/10780/page5&nbsp;      |
//+------------------------------------------------------------------+

void LinearRegr1(double X[], double Y[], int N, double& A, double& B)
{
      double mo_X = 0.0, mo_Y = 0.0, var_0 = 0.0, var_1 = 0.0;
      
    for ( int i = 0; i < N; i ++ )
      {
        mo_X +=X[i];
        mo_Y +=Y[i];
      }
    mo_X /=N;
    mo_Y /=N;
        
    for ( i = 0; i < N; i ++ )
      {
        var_0 +=(X[i]-mo_X)*(Y[i]-mo_Y);
        var_1 +=(X[i]-mo_X)*(X[i]-mo_X);
      }
        A = var_0 / var_1;
        B = mo_Y - A * mo_X;
}

Результат

2008.07.30 13:51:08 LinReg EURUSD,M1: A = 0.00000090 B = -1098.77264952

2008.07.30 13:51:08 LinReg EURUSD,M1: A = 0.00000090 B = -1098.77264952

2008.07.30 13:51:08 LinReg EURUSD,M1: A = -1078.77267965 B = 0.00000089

2008.07.30 13:51:08 LinReg EURUSD,M1: A = -1102.16914108 B = 0.00000091

А такого быть не должно.

2. То что появляются два цикла я вижу, поэтому и просил как-то увеличить быстродействие. Быстрее возможно будет вот этот алгоритм 'Регрессия: что это такое?', но его надо тоже оптимизировать ( думаю Vinin это уже сделал)

3. За Кауфмана спасибо хороший индикатор. Если Вы не забыли перед вторым чемпионатом я отлавливал в нем неточности работы. Спасибо что вы их поправили.

З.Ы. Просьба у кого есть Маtlab. Введите эти массивы и рассчитайте встроенными формулами (насколько я помню, там есть), и выложите результат сюда. Что бы прийти к единому мнению. Спасибо. Помогайте )). Rosh-а довольно трудно убедить, но и у меня лоб военный ))

 
Prival писал (а) >>

2. То что появляются два цикла я вижу, поэтому и просил как-то увеличить быстродействие. Быстрее возможно будет вот этот алгоритм 'Регрессия: что это такое?', но его надо тоже оптимизировать ( думаю Vinin это уже сделал)

LWMA действительно безопаснее чем X*X, так что ваша с Mathemat'ом работа обретает новый смысл :). Но мою первую рекомендацию (сдвиг начала координат) я продолжаю считать лучшим вариантом. Неужели формальная замена повсюду Time[pos] на Time[pos]-Time[Bars-1] несёт такой риск внести ошибку?

 
Prival писал (а) >>

1. В том то и дело что результат не должен зависеть от сортировки, а в том алгоритме зависит. Вот проверьте.

Результат

2008.07.30 13:51:08 LinReg EURUSD,M1: A = 0.00000090 B = -1098.77264952

2008.07.30 13:51:08 LinReg EURUSD,M1: A = 0.00000090 B = -1098.77264952

2008.07.30 13:51:08 LinReg EURUSD,M1: A = -1078.77267965 B = 0.00000089

2008.07.30 13:51:08 LinReg EURUSD,M1: A = -1102.16914108 B = 0.00000091

А такого быть не должно.



Вставьте в ваш код распринтокву:

//-------------------------------------------------------------------------------
// использование этой формулы приводит к ошибкам если X=Time
// формула предложена вот тут https://forum.mql4.com/ru/10780/page4
//| y(x)=A+B*x  
 
void LinearRegr(double X[], double Y[], int N, double& A, double& B)
{
      double sumY = 0.0, sumX = 0.0, sumXY = 0.0, sumX2 = 0.0;
      
    for ( int i = 0; i < N; i ++ )
    {
        sumY   +=Y[i];
        sumXY  +=X[i]*Y[i];
        sumX   +=X[i];
        sumX2  +=X[i]*X[i];
    }
   Print("sumY = ",DoubleToStr(sumY,8)," sumX = ",DoubleToStr(sumX,8)," sumXY = ",DoubleToStr(sumXY,8)," sumX2 = ",DoubleToStr(sumX2,8));
   Print("sumXY*dN-sumX*sumY = ",DoubleToStr(sumXY*dN-sumX*sumY,8));    
   Print("sumX2*dN-sumX*sumX = ",DoubleToStr(sumX2*dN-sumX*sumX,8));    
 
   B=(sumXY*N-sumX*sumY)/(sumX2*N-sumX*sumX);
   A=(sumY-sumX*B)/N;
}

Получите что-то вроде этого:

первый вызов
sumY = 11.98130000 sumX = 7299840060.00000000 sumXY = 14576928951.87000100 sumX2 = 8881277483596863500.00000000
sumXY*dN-sumX*sumY = 0.34199524
sumX2*dN-sumX*sumX = 376832.00000000
A = -1102.16914108 B = 0.00000091
второй вызов
sumY = 11.98130000 sumX = 7299840060.00000000 sumXY = 14576928951.87000300 sumX2 = 8881277483596864500.00000000
sumXY*dN-sumX*sumY = 0.34202576
sumX2*dN-sumX*sumX = 385024.00000000
A = -1078.77267965 B = 0.00000089

Это еще один подводный камень компьютерных вычислений и округлений. С одной стороны я сам не ожидал таких граблей, но с другой стороны такая разница вполне объяснима, когда два ряда значений (X и Y) имеют слишком большую разницу в порядке значений.

1234567
Новый комментарий