2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему
 
 Помогите решить уравнение
Сообщение14.12.2012, 21:03 


07/03/11
690
Даны все $\{X_i, Y_i, \sigma _i| i=\overline {1,n}\}$, и имеется явный вид для формул $m(X), v(X)$. Как найти параметр $\beta$ в следующем уравнении?$$\sum\limits _{i=1}^n\frac{Y_i-\beta m(X_i)}{\beta ^2v(X_i)+\sigma _i^2}m(X_i)=0$$Явного вида скорее всего не получится, хотя бы алгоритм подскажите. ($n$ - очень большое)

 Профиль  
                  
 
 Re: Помогите решить уравнение
Сообщение14.12.2012, 23:15 
Заслуженный участник


12/09/10
1547
Что в вашем понимании очень большое? Тысяча? Миллион? Квадриллион?
Порядок скажите...

 Профиль  
                  
 
 Re: Помогите решить уравнение
Сообщение14.12.2012, 23:45 


07/03/11
690
до миллиона, где-то 10.000-100.000

 Профиль  
                  
 
 Re: Помогите решить уравнение
Сообщение14.12.2012, 23:53 
Заслуженный участник


12/09/10
1547
Это немного.
Можете использовать любой алгоритм поиска корней нелинейного уравнения.
Даже методом половинного деления быстро посчитается.

 Профиль  
                  
 
 Re: Помогите решить уравнение
Сообщение15.12.2012, 23:43 


07/03/11
690
Написал методом половинного деления -- работает очень медленно. Начиная с 2000 считает около 10 секунд. Очень долго каждый раз вручную пересчитывать такую сумму. Кто ещё подскажет что-то?

 Профиль  
                  
 
 Re: Помогите решить уравнение
Сообщение16.12.2012, 00:17 
Заслуженный участник


09/05/12
25179
vlad_light в сообщении #658908 писал(а):
Написал методом половинного деления -- работает очень медленно. Начиная с 2000 считает около 10 секунд. Очень долго каждый раз вручную пересчитывать такую сумму. Кто ещё подскажет что-то?

По-видимому, Вы как-то очень неудачно написали дихотомию. Сейчас интереса ради набросал эту программу (для случайных данных) - для $n=10^6$ с итоговой погрешностью $10^{-8}$ счет занимает 0.1 секунды (на одном ядре приличного, но вполне обычного процессора).

 Профиль  
                  
 
 Re: Помогите решить уравнение
Сообщение16.12.2012, 01:11 


07/03/11
690
Сколько времени считает одну итерацию?
код: [ скачать ] [ спрятать ]
  1. #include <iostream> 
  2. #include <random> 
  3. #include <cstdlib> 
  4. #include <ctime> 
  5.  
  6. const int SAMPLE_SIZE = 1000; 
  7. const double PRECISION = 0.01; 
  8.  
  9. const double MU_TEST = -240.3; 
  10. const double SIGMA_TEST_SQUARED = 591.7; 
  11. const double BETA_TEST = 45.6; 
  12.  
  13. const double SIGMA_EPSILON_SQUARED = 7.6; 
  14. const double SIGMA_DELTA_SQUARED = 1.8; 
  15.  
  16. double Box_Muller (); 
  17. void init_sample (double* sample_x, double* sample_y); 
  18.  
  19. double sample_mean (const double* sample_x); 
  20. double sample_squared_variance (const double* sample_x); 
  21.  
  22. double sum_Q (const double* sample_x, const double* sample_y, const double beta); 
  23. double dichotomy_method (const double epsilon, const double a, const double b, const double* sample_x, const double* sample_y); // or any other solver 
  24.  
  25. int main () 
  26. srand (time (NULL)); 
  27.  
  28. double sample_X[SAMPLE_SIZE], sample_Y[SAMPLE_SIZE]; 
  29.  
  30. init_sample (sample_X, sample_Y); 
  31.  
  32. std::cout << dichotomy_method (PRECISION, -1000, 1000, sample_X, sample_Y); 
  33.  
  34. std::cin.get (); 
  35. return 0; 
  36.  
  37. double Box_Muller () 
  38. do 
  39. double r_1 = rand (); 
  40. r_1 = (r_1 / RAND_MAX - 0.5) * 2; 
  41. double r_2 = rand (); 
  42. r_2 = (r_2 / RAND_MAX - 0.5) * 2; 
  43.  
  44. double s = r_1 * r_1 + r_2 * r_2; 
  45.  
  46. if (s > 0 && s <= 1) 
  47. return r_1 * std::sqrt (-2 * std::log (s) / s); 
  48. while (true); 
  49.  
  50. void init_sample (double* sample_x, double* sample_y) 
  51. for (int i = 0; i < SAMPLE_SIZE; ++i) 
  52. const double EPSILON = std::sqrt (SIGMA_EPSILON_SQUARED) * Box_Muller (); 
  53. const double DELTA = std::sqrt (SIGMA_DELTA_SQUARED) * Box_Muller (); 
  54.  
  55. const double X = MU_TEST + std::sqrt (SIGMA_TEST_SQUARED) * Box_Muller () + DELTA; 
  56. const double Y = BETA_TEST * (MU_TEST + std::sqrt (SIGMA_TEST_SQUARED) * Box_Muller ()) *  
  57. (MU_TEST + std::sqrt (SIGMA_TEST_SQUARED) * Box_Muller ()) + EPSILON; 
  58.  
  59. sample_x[i] = X; 
  60. sample_y[i] = Y; 
  61.  
  62. double sample_mean (const double* sample_x) 
  63. double mean = 0; 
  64.  
  65. for (int i = 0; i < SAMPLE_SIZE; ++i) 
  66. mean += sample_x[i]; 
  67.  
  68. return mean / SAMPLE_SIZE; 
  69.  
  70. double sample_squared_variance (const double* sample_x) 
  71. const double mean = sample_mean (sample_x); 
  72.  
  73. double variance = 0; 
  74.  
  75. for (int i = 0; i < SAMPLE_SIZE; ++i) 
  76. variance += (sample_x[i] - mean) * (sample_x[i] - mean); 
  77.  
  78. return variance / (SAMPLE_SIZE - 1) - SIGMA_DELTA_SQUARED; 
  79.  
  80. double sum_Q (const double* sample_x, const double* sample_y, const double beta) 
  81. double sum = 0; 
  82.  
  83. for (int i = 0; i < SAMPLE_SIZE; ++i) 
  84. const double left_summand = (sample_mean (sample_x) * SIGMA_DELTA_SQUARED + sample_x[i] * sample_squared_variance (sample_x)) /  
  85. (SIGMA_DELTA_SQUARED + sample_squared_variance (sample_x)) *  
  86. (sample_mean (sample_x) * SIGMA_DELTA_SQUARED + sample_x[i] * sample_squared_variance (sample_x)) /  
  87. (SIGMA_DELTA_SQUARED + sample_squared_variance (sample_x)); 
  88. const double right_summand = SIGMA_DELTA_SQUARED * sample_squared_variance (sample_x) / (SIGMA_DELTA_SQUARED + sample_squared_variance (sample_x)); 
  89.  
  90. const double dm = left_summand + right_summand; 
  91.  
  92. const double v = beta * beta * (4 * left_summand * right_summand + 2 * right_summand * right_summand) + SIGMA_EPSILON_SQUARED; 
  93.  
  94. sum += (sample_y[i] - beta * dm) * dm / v; 
  95.  
  96. return sum; 
  97.  
  98. double dichotomy_method (const double epsilon, const double a, const double b, const double* sample_x, const double* sample_y) 
  99. double new_a = a, new_b = b; 
  100.  
  101. while (std::abs (new_a - new_b) > epsilon) 
  102. const double c = (new_a + new_b) / 2; 
  103.  
  104. if (sum_Q (sample_x, sample_y, new_a) * sum_Q (sample_x, sample_y, c) < 0) 
  105. new_b = c; 
  106. else 
  107. new_a = c; 
  108.  
  109. std::cout << c << std::endl; 
  110.  
  111. return (new_a + new_b) / 2; 

 Профиль  
                  
 
 Re: Помогите решить уравнение
Сообщение16.12.2012, 01:54 
Заслуженный участник


09/05/12
25179
vlad_light в сообщении #658943 писал(а):
Сколько времени считает одну итерацию?

Поскольку начальный интервал имел размер $10^1$, а в итоге сократился до $10^{-8}$, то это 30 итераций. Дихотомия же.

Так... послушайте, а зачем Вы при каждом вычислении суммы заново считаете все $m(X_i)$ и $v(X_i)$? И зачем при выборе правого или левого отрезка заново пересчитываете значение суммы в new_a? Да, Вы "идеологически правильно" разнесли все по куче функций без глобальных переменных, но именно это и привело к такой безумной потере производительности.

 Профиль  
                  
 
 Re: Помогите решить уравнение
Сообщение16.12.2012, 08:36 
Заслуженный участник


12/09/10
1547
Вместо функций sample_mean и sample_squared_variance сделайте 2 массива. Этого будет достаточно. Сложность сразу сократится с $O(n^2)$ до $O(n)$. Все остальное уже мелочи.

 Профиль  
                  
 
 Re: Помогите решить уравнение
Сообщение16.12.2012, 14:31 


07/03/11
690
Всем большое спасибо, учёл все комментарии -- работает за 0.1 сек :D Ещё раз спасибо!

(Оффтоп)

Не будет из меня програмиста :-(

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 10 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group