Форум программистов, компьютерный форум CyberForum.ru

Оптимизация алгоритма - C++

Восстановить пароль Регистрация
 
vlad_light
4 / 4 / 0
Регистрация: 24.09.2012
Сообщений: 178
10.03.2013, 18:35     Оптимизация алгоритма #1
Условие:
Дана выборка http://www.cyberforum.ru/cgi-bin/latex.cgi?(X_i, Y_i)_{i=1}^N. Предполагается, что она была построена по следующему закону:
http://www.cyberforum.ru/cgi-bin/latex.cgi?\begin{cases}<br />
Y=\beta \xi ^2+\varepsilon \\<br />
X=\xi +\delta \\<br />
\end{cases}<br />
где http://www.cyberforum.ru/cgi-bin/latex.cgi?\beta\in\mathbb R - константа, http://www.cyberforum.ru/cgi-bin/latex.cgi?\xi\sim N(\mu ,\sigma ^2), \varepsilon\sim N(0,\sigma _\varepsilon ^2), \delta\sim N(0,\sigma _\delta ^2), причем все http://www.cyberforum.ru/cgi-bin/latex.cgi?\xi ,\varepsilon ,\delta - независимы.
Из всего этого нам известны: http://www.cyberforum.ru/cgi-bin/latex.cgi?\sigma _\varepsilon ,\sigma _\delta и нужно найти (оценить) параметр http://www.cyberforum.ru/cgi-bin/latex.cgi?\beta.
Для этого существуют (как минимум) 2 метода: предварительное оценивание и совместное оценивание. Рассмотрим каждый из них по-отдельности.

1. Предварительное оценивание
Сперва оценим http://www.cyberforum.ru/cgi-bin/latex.cgi?\mu и http://www.cyberforum.ru/cgi-bin/latex.cgi?\sigma:
http://www.cyberforum.ru/cgi-bin/latex.cgi?\mu = \hat\mu =\frac 1N\sum\limits _{i=1}^NX_i; http://www.cyberforum.ru/cgi-bin/latex.cgi?\sigma = \hat\sigma =\sqrt{\frac{1}{N-1}\sum\limits _{i=1}^N((X_i-\mu)^2)-\sigma _\delta ^2}
Для оценки http://www.cyberforum.ru/cgi-bin/latex.cgi?\beta мы используем функцию http://www.cyberforum.ru/cgi-bin/latex.cgi?q(x,y,\beta )=\frac{y-m(x,\beta )}{v(x,\beta )}\cdot \frac{dm}{d\beta}(x,\beta )
где http://www.cyberforum.ru/cgi-bin/latex.cgi?m(x,\beta )=E(Y|X)=\beta ((\frac{\mu\sigma _\delta ^2+x\sigma ^2}{\sigma _\delta ^2+\sigma ^2})^2+\frac{\sigma ^2\sigma _\delta ^2}{\sigma ^2+\sigma _\delta ^2}) -условное матожидание
http://www.cyberforum.ru/cgi-bin/latex.cgi?v(x,\beta )=Var(Y|X)=\beta ^2 (4(\frac{\mu\sigma _\delta ^2+x\sigma ^2}{\sigma _\delta ^2+\sigma ^2})^2\frac{\sigma ^2\sigma _\delta ^2}{\sigma ^2+\sigma _\delta ^2}+2(\frac{\sigma ^2\sigma _\delta ^2}{\sigma ^2+\sigma _\delta ^2})^2)+\sigma _\varepsilon ^2 -условная дисперсия
Тогда решение уравнения http://www.cyberforum.ru/cgi-bin/latex.cgi?\sum\limits _{i=1}^N q(X_i,Y_i,\hat\beta)=0 относительно http://www.cyberforum.ru/cgi-bin/latex.cgi?\hat\beta и будет нашей оценкой для http://www.cyberforum.ru/cgi-bin/latex.cgi?\beta, т.е. http://www.cyberforum.ru/cgi-bin/latex.cgi?\hat\beta\to\beta ,N\to\infty почти наверное.

2. Совместное оценивание
Пусть http://www.cyberforum.ru/cgi-bin/latex.cgi?\theta = (\mu ,\sigma ,\beta )^T. Тогда функцию http://www.cyberforum.ru/cgi-bin/latex.cgi?q перепишем в виде:
http://www.cyberforum.ru/cgi-bin/latex.cgi?q(x,y;\theta )=\frac{y-m(x;\theta )}{v(x;\theta )}\cdot \frac{\partial m}{\partial\theta}(x;\theta )+\frac{\partial \ln f_X}{\partial\theta}(x;\theta)
где http://www.cyberforum.ru/cgi-bin/latex.cgi?f_X (x;\theta )=\frac{1}{\sqrt{2\pi (\sigma ^2+\sigma _\delta ^2)}}\exp (-\frac 12\cdot\frac{(x-\mu )^2}{\sigma ^2+\sigma _\delta ^2}) - плотность http://www.cyberforum.ru/cgi-bin/latex.cgi?X.
Тогда для оценки параметра http://www.cyberforum.ru/cgi-bin/latex.cgi?\theta нужно решить систему из 3 уравнений с 3 неизвестными: http://www.cyberforum.ru/cgi-bin/latex.cgi?\mu ,\sigma ,\theta:
http://www.cyberforum.ru/cgi-bin/latex.cgi?\begin{cases}<br />
\sum\limits _{i=1}^N \frac{y-m(x;\theta )}{v(x;\theta )}\cdot \frac{\partial m}{\partial\mu}(x;\theta )+\frac{\partial \ln f_X}{\partial\mu}(x;\theta) = 0 \\<br />
\sum\limits _{i=1}^N \frac{y-m(x;\theta )}{v(x;\theta )}\cdot \frac{\partial m}{\partial\sigma}(x;\theta )+\frac{\partial \ln f_X}{\partial\sigma}(x;\theta) =0 \\<br />
\sum\limits _{i=1}^N \frac{y-m(x;\theta )}{v(x;\theta )}\cdot \frac{\partial m}{\partial\beta}(x;\theta )=0 \\<br />
\end{cases}<br />
Выражения http://www.cyberforum.ru/cgi-bin/latex.cgi?\frac{\partial m}{\partial\mu}(x;\theta ),\frac{\partial m}{\partial\sigma}(x;\theta ),\frac{\partial m}{\partial\beta}(x;\theta ),\frac{\partial \ln f_X}{\partial\mu}(x;\theta),\frac{\partial \ln f_X}{\partial\sigma}(x;\theta) в явном виде я здесь выписывать не буду, поскольку их достаточно просто посчитать.
При решении данной системы мы получим оценку на параметры http://www.cyberforum.ru/cgi-bin/latex.cgi?\mu ,\sigma ,\beta, которая также будет стремиться к настоящим параметрам почти наверное.

Задание:
Моя задача -- реализовать эти два метода, чтоб эмпирически преверить гипотезу о неравенстве дисперсий:
http://www.cyberforum.ru/cgi-bin/latex.cgi?\Sigma _1>\Sigma _2 (http://www.cyberforum.ru/cgi-bin/latex.cgi?\Leftrightarrow \Sigma _1 - \Sigma _2 - положительно определенная матрица)
где http://www.cyberforum.ru/cgi-bin/latex.cgi?\Sigma _1 и http://www.cyberforum.ru/cgi-bin/latex.cgi?\Sigma _2 - ковариационные матрицы оценок для метода 1 и 2 соответственно, т.е. показать, что дисперсия оценки 2 будет меньше, чем дисперсия оценки 1 (http://www.cyberforum.ru/cgi-bin/latex.cgi?\Leftrightarrowоценка 2 -- точнее)
Либо опровергнуть её, показав, что
http://www.cyberforum.ru/cgi-bin/latex.cgi?\Sigma _1=\Sigma _2

Я написал следующий код для 1-ого метода. Проблема в том, что работает он очень медленно. Как его можно оптимизировать? Заранее огромное спасибо!
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
// main.cpp
 
#include "Sample.h"
#include "LS.h"
#include <ctime>
 
int main ()
{
    srand (time (0));
 
    const int N = 20;
 
    const int SIZE = 100;
    const double beta = 10;
    const double mu = 1975;
    const double sigma = 10;
    const double sigma_epsilon = 10;
    const double sigma_delta = 10;
 
    const double left = 1;
    const double right = 20;
    const double precision = 0.01;
 
    std::vector<double> solutions;
    solutions.reserve (N);
 
    for (int i = 0; i < N; ++i)
    {
        Sample A (SIZE, beta, mu, sigma, sigma_epsilon, sigma_delta);
        A.initSample ();
 
        double solution = solverLS (left, right, precision, A);
        solutions.push_back (solution);
 
        std::cout << solution << std::endl;
    }
 
    std::cout << "Mean: " << sampleMean (solutions) << std::endl;
    std::cout << "Deviation: " << std::sqrt ( sampleDispersion (solutions) ) << std::endl;
 
    std::cout << "Done";
    std::cin.get ();
 
    return 0;
}
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
// Sample.h
 
#pragma once
 
#include <vector>
#include <cstdlib>
#include <iostream>
 
class Sample
{
private:
    std::vector<double> Y;
    std::vector<double> X;
 
    int N;
 
    double beta;
    
    double mu;
    double sigma;
        
    double sigma_epsilon;
    double sigma_delta;
 
    double box_muller (const double M, const double D) const;
 
public:
    Sample ();
    Sample (const Sample& copy);
    Sample& operator= (const Sample& copy);
    Sample (const int N, const double beta, const double mu, const double sigma, 
            const double sigma_epsilon, const double sigma_delta);
 
    ~Sample ();
 
    void initSample ();
 
    const std::vector<double> returnY () const;
    const std::vector<double> returnX () const;
 
    const double returnSigma_epsilon () const;
    const double returnSigma_delta () const;
 
    // #debug_begin
 
    void displayY () const;
    void displayX () const;
 
    const double returnBeta () const;
    const double returnMu () const;
    const double returnSigma () const;
 
    // #debug_end
};
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
//Sample.cpp
 
#include "Sample.h"
 
double Sample::box_muller (const double M, const double D) const
{
    double x;
    double y;
 
    double s = 0;
 
    while (s == 0 || s > 1)
    {
        x = ( ( (double) rand () ) / RAND_MAX - 0.5 ) * 2;
        y = ( ( (double) rand () ) / RAND_MAX - 0.5 ) * 2;
 
        s = x * x + y * y;
    }
 
    return ( x * std::sqrt (-2 * std::log (s) / s) ) * D + M;
}
 
Sample::Sample ()
{
 
}
 
Sample::Sample (const Sample& copy)
{
 
}
 
Sample& Sample::operator= (const Sample& copy)
{
    return *this;
}
 
Sample::Sample (const int init_N, const double init_beta, const double init_mu, 
                const double init_sigma, const double init_sigma_epsilon, 
                const double init_sigma_delta)
{
    N = init_N;
    beta = init_beta;
    mu = init_mu;
    sigma = init_sigma;
    sigma_epsilon = init_sigma_epsilon;
    sigma_delta = init_sigma_delta;
}
 
Sample::~Sample ()
{
 
}
 
void Sample::initSample ()
{
    Y.reserve (N);
    X.reserve (N);
    for (int i = 0; i < N; ++i)
    {
        double xi = box_muller (mu, sigma);
 
        Y.push_back ( beta * xi * xi + box_muller (0, sigma_epsilon) );
        X.push_back ( xi + box_muller (0, sigma_delta) );
    }
}
 
const std::vector<double> Sample::returnY () const
{
    return Y;
}
 
const std::vector<double> Sample::returnX () const
{
    return X;
}
 
const double Sample::returnSigma_epsilon () const
{
    return sigma_epsilon;
}
 
const double Sample::returnSigma_delta () const
{
    return sigma_delta;
}
 
void Sample::displayY () const
{
    std::cout << "Y - sample:" << std::endl;
 
    for (int i = 0; i < N; ++i)
        std::cout << Y.at (i) << std::endl;
}
 
void Sample::displayX () const
{
    std::cout << "X - sample:" << std::endl;
 
    for (int i = 0; i < N; ++i)
        std::cout << X.at (i) << std::endl;
}
 
const double Sample::returnBeta () const
{
    return beta;
}
 
const double Sample::returnMu () const
{
    return mu;
}
 
const double Sample::returnSigma () const
{
    return sigma;
}
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
//LS.h
 
#pragma once
 
#include "Sample.h"
 
const double sampleMean (const std::vector<double>& sample);
const double sampleDispersion (const std::vector<double>& sample);
 
const double LS (const Sample& sample, const double beta);
const double fastLS (const Sample& sample, const double beta);
 
const double solverLS (const double a, const double b, const double precision, const Sample& sample);
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
// LS.cpp
 
#include "LS.h"
 
const double sampleMean (const std::vector<double>& sample)
{
    double mean = 0;
 
    for (auto iter = sample.cbegin (); iter != sample.cend (); ++iter)
        mean += *iter;
 
    return mean / ( sample.size () );
}
 
const double sampleDispersion (const std::vector<double>& sample)
{
    double mean = sampleMean (sample);
    double deviation = 0;
 
    for (auto iter = sample.cbegin (); iter != sample.cend (); ++iter)
        deviation += ( *iter - mean) * ( *iter - mean);
 
    return std::sqrt ( deviation / ( sample.size () - 1 ) );
}
 
const double LS (const Sample& sample, const double beta)
{
    const double mu = sampleMean ( sample.returnX () );
 
    const double sigma_deltaSquared = sample.returnSigma_delta () * sample.returnSigma_delta ();
    const double sigma_epsilonSquared = sample.returnSigma_epsilon () * sample.returnSigma_epsilon ();
 
    const double sigmaSquared = sampleDispersion ( sample.returnX () ) - sigma_deltaSquared;
 
    // c2 - right summand
 
    const double c22 = sigma_deltaSquared + sigmaSquared;
    const double c21 = sigma_deltaSquared * sigmaSquared;
 
    const double c2 = c21 / c22;
 
    // c1 - left summand
 
    const double c12 = c22 * c22;
 
    std::vector<double> c1;
    c1.reserve ( sample.returnX ().size () );
 
    for (int i = 0; i < sample.returnX ().size (); ++i)
    {
        const double c = mu * sigma_deltaSquared + sample.returnX ().at (i) * sigmaSquared;
        const double c11 = c * c;
 
        c1.push_back (c11 / c12);
    }
 
    // dm - derivative of m with respect to beta
 
    std::vector<double> dm;
    dm.reserve ( sample.returnX ().size () );
 
    for (int i = 0; i < sample.returnX ().size (); ++i)
        dm.push_back ( c1.at (i) + c2 );
 
    // m - expected value E(Y|X)
 
    std::vector<double> m;
    m.reserve ( sample.returnX ().size () );
 
    for (int i = 0; i < sample.returnX ().size (); ++i)
        m.push_back ( dm.at (i) * beta );
 
    // v - dispersion D(Y|X)
 
    std::vector<double> v;
    v.reserve ( sample.returnX ().size () );
 
    for (int i = 0; i < sample.returnX ().size (); ++i)
        v.push_back ( ( 4 * c1.at (i) * c2 + 2 * c2 * c2 ) * beta * beta + sigma_epsilonSquared);
 
    // q - LinearScore function
 
    std::vector<double> q;
    q.reserve ( sample.returnX ().size () );
 
    for (int i = 0; i < sample.returnX ().size (); ++i)
        q.push_back ( dm.at (i) * ( sample.returnY ().at (i) - m.at (i) ) / v.at (i) );
 
    // Q - sum of q's
 
    double Q = 0;
 
    for (int i = 0; i < q.size (); ++i)
        Q += q.at (i);
 
    return Q;
}
 
const double fastLS (const Sample& sample, const double beta)
{
    const double mu = sampleMean ( sample.returnX () );
    const double sigma_deltaSquared = sample.returnSigma_delta () * sample.returnSigma_delta ();
    const double sigma_epsilonSquared = sample.returnSigma_epsilon () * sample.returnSigma_epsilon ();
    const double sigmaSquared = sampleDispersion ( sample.returnX () ) - sigma_deltaSquared;
    const double c22 = sigma_deltaSquared + sigmaSquared;
    const double c2 = ( sigma_deltaSquared * sigmaSquared ) / ( c22 );
    const double c12 = c22 * c22;
    double Q = 0;
    for (int i = 0; i < sample.returnX ().size (); ++i)
    {
        const double c = mu * sigma_deltaSquared + sample.returnX ().at (i) * sigmaSquared;
        const double c1 = c * c / c12;
        const double dm = c1 + c2;
        const double m = dm * beta;
        const double v = ( 4 * c1 * c2 + 2 * c2 * c2 ) * beta * beta + sigma_epsilonSquared;
        const double q = dm * (sample.returnY ().at (i) - m ) / v;
        Q += q;
    }   
    return Q;
}
 
const double solverLS (const double a, const double b, const double precision, const Sample& sample)
{
    double left = a;
    double right = b;
 
    double LSLeft = fastLS (sample, left);  // = LS (...)
    double LSRight = fastLS (sample, right);    // = LS (...)
 
    double distance = std::abs ( LSLeft - LSRight ); 
 
    while (distance > precision)
    {
        double mid = ( left + right ) / 2;
        double LSMid = fastLS (sample, mid);    // = LS (...)
 
        if ( ( LSLeft > 0 && LSMid < 0 ) || ( LSLeft < 0 && LSMid > 0 ) )
        {
            right = mid;
            LSRight = LSMid;
        }
        else
        {
            left = mid;
            LSLeft = LSMid;
        }
 
        distance = std::abs ( LSRight - LSLeft );
    }
 
    return ( left + right ) / 2;
}
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
10.03.2013, 18:35     Оптимизация алгоритма
Посмотрите здесь:

Оптимизация C++
Оптимизация алгоритма вычисления определителя матрицы C++
Оптимизация алгоритма C++
C++ Оптимизация программы на С++
Оптимизация функции C++
После регистрации реклама в сообщениях будет скрыта и будут доступны все возможности форума.
Avazart
 Аватар для Avazart
6901 / 5141 / 252
Регистрация: 10.12.2010
Сообщений: 22,602
Записей в блоге: 17
10.03.2013, 18:45     Оптимизация алгоритма #2
То что угляделось:
C++
1
 c1.at (i)
Используй [] либо итераторы они быстрее...
C++
1
2
for (int i = 0; i < sample.returnX ().size (); ++i)
        m.push_back ( dm.at (i) * beta );
Можно вынести
C++
1
2
3
const int size= sample.returnX ().size ();
for (int i = 0; i < size; ++i)
        m.push_back ( dm.at (i) * beta );
Добавлено через 1 минуту
Все остальное возможно дело математики.
Да и может оказаться что использовать массивы намного быстрее чем вектора.
vlad_light
4 / 4 / 0
Регистрация: 24.09.2012
Сообщений: 178
10.03.2013, 18:56  [ТС]     Оптимизация алгоритма #3
Я вот тоже думал, что обычные массивы быстрее будут. Просто я никогда их не использовал, поэтому побоялся изначально с их помощью писать Сейчас попробую всё заменить и протестировать.
Всю математику я описал выше -- для первого метода всё в явном виде. Если тебя не затруднит -- можешь запустить этот код у себя и проверить быстродействие: оно действительно должно примерно столько времени считать или должно быть в раз 50-100 быстрее? Спасибо!
diagon
Higher
 Аватар для diagon
1920 / 1186 / 49
Регистрация: 02.05.2010
Сообщений: 2,925
Записей в блоге: 2
10.03.2013, 19:27     Оптимизация алгоритма #4
Какой у вас компилятор? Если какой-нибудь древний, то мучайтесь(основные тормоза могут быть как раз из-за глупого компилятора). Если что-то современное, то просто врубите профилирование и посмотрите, какая часть кода жрет больше всего времени.
vlad_light
4 / 4 / 0
Регистрация: 24.09.2012
Сообщений: 178
10.03.2013, 19:41  [ТС]     Оптимизация алгоритма #5
у меня мс вс 2010 экспресс.
врубите профилирование
Это что такое?

Добавлено через 8 минут
нашел одну ошибку:
в файле LS.cpp функция sampleDispersion должна возвращать результат БЕЗ корня. (строка 23)
diagon
Higher
 Аватар для diagon
1920 / 1186 / 49
Регистрация: 02.05.2010
Сообщений: 2,925
Записей в блоге: 2
10.03.2013, 20:02     Оптимизация алгоритма #6
Цитата Сообщение от vlad_light Посмотреть сообщение
у меня мс вс 2010 экспресс.
В экспрэсс профайлера вроде нету... Ультимейт студия есть на торрентах.
Там достаточно удобный профайлер, для каждой строчки кода показывает время, которое заняло ее выполнение.
Yandex
Объявления
10.03.2013, 20:02     Оптимизация алгоритма
Ответ Создать тему
Опции темы

Текущее время: 20:20. Часовой пояс GMT +3.
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2016, vBulletin Solutions, Inc.
Рейтинг@Mail.ru