Как подогнать кривую в форме колокола к двумерным данным рассеяния?

У меня есть данные разброса x-y, которые демонстрируют поведение в форме колокола (то есть в форме нормального распределения) в течение года. Это первичные производственные данные из высоких широт (подробнее здесь, статья платная, но надеюсь цифры видны).

Вопрос

Как подогнать кривую формы нормального распределения к данным разброса в ggplot2?

Пример данных

x <- structure(list(yday = c(238, 238, 238, 242, 242, 250, 250, 253, 
254, 169, 199, 208, 230, 21, 37, 88, 94, 102, 125, 125, 95, 98, 
100, 101, 103, 93, 96, 97, 97, 99, 291, 300, 316, 332, 363, 12, 
27, 49, 68, 256, 263, 263, 264, 265, 266, 127, 127, 127, 127, 
133, 123, 127, 133, 127, 133, 141, 148, 155, 112, 120, 127, 134, 
169, 169, 169, 169, 169, 124, 124, 124, 124, 126, 126, 127, 130, 
132, 134, 134, 136, 138, 140, 142, 144, 146, 149, 152, 154, 127, 
130, 132, 134, 134, 136, 138, 140, 142, 144, 146, 149, 152, 154, 
127, 130, 132, 134, 134, 136, 138, 140, 142, 144, 146, 149, 152, 
154, 127, 130, 132, 134, 134, 136, 138, 140, 142, 144, 146, 149, 
152, 154, 127, 130, 132, 134, 134, 136, 138, 140, 142, 144, 146, 
149, 152, 154, 146, 154, 154, 154, 143, 156, 134, 162, 164, 168, 
166, 71, 38, 39), value = c(48802.2, 28869.36, 46370.31, 67936, 
91442, 89559.4, 46862.2, 7895.3, 19660.72, 273540.7, 254615, 
268930, 264310.56, 72561, 114520, 42950, 149151.15, 530610, 53289.2, 
45, 1776, 20504, 1768, 5740, 29497, 340762, 259837, 11576, 847238, 
1773275, 1555.92, 3108.48, 8579.1, 25677.8, 17697.32, 2887.56, 
5311.2, 13127.98, 38006.4, 135128, 71003, 10454.75, 41389.6, 
15266.5, 58601.7, 206984.918282083, 265165.058077198, 90485.4790849673, 
5705618.16993847, 1616527.31316346, 1610059.4788107, 5689427.93092749, 
3261840.85863376, 1911057.16943202, 1023812.55301328, 1579191.31813709, 
2241683.51873045, 4398531.75676259, 933143.183151504, 1771596.51236257, 
3000366.86522064, 1219826.2208944, 247538.595548984, 353927.523573691, 
323722.062546854, 278081.544235635, 601042.642308546, 2317070.57555887, 
963348.671707912, 8125168.04401668, 1860334.91955526, 18673.3716353477, 
682901.426071428, 348046.291238703, 1387947.38534056, 112176.06673827, 
203778.898538342, 304593.428222028, 1015454.26894711, 384172.102208766, 
1211065.9345086, 580449.092224899, 556147.163209095, 707840.652723421, 
2919016.89462558, 878518.35266303, 760837.632557093, 437441.609086177, 
3761984.12905246, 1008524.7172583, 153914.10863321, 209919.739543153, 
5165174.16501832, 592152.070785338, 754878.057858348, 548774.607567716, 
784679.488265372, 1191547.05905346, 867977.806474748, 1601076.47417622, 
1059665.24406883, 509654.672768973, 1878007.77720015, 217773.469093887, 
282571.399726361, 98438.9397685662, 889753.057501427, 564416.438455766, 
1843608.78521975, 727213.52622083, 689307.464580901, 1018069.45500141, 
1188687.56383149, 1352651.53225745, 726363.223839249, 420446.302222222, 
856363.289847527, 18015.1056535911, 229628.041636759, 165657.605285714, 
164394.955219614, 510449.457665504, 1778093.57209278, 610564.603533888, 
889187.420481627, 2762856.39975472, 863978.618292937, 1697540.61924469, 
859284.971006319, 197822.92972973, 389791.010663156, 144870.825015753, 
196128.64471631, 133023.96688172, 70787.1033258229, 518223.208383732, 
4051834.77590046, 628912.36192445, 378818.831615793, 413839.100579421, 
1091509.82410159, 1187325.98867099, 331226.406610866, 1729022.69104484, 
4663215.78870189, 2843159.21140248, 708207.236363227, 1436498.03405122, 
1158173.94324553, 448469.915666212, 901903.855484778, 1599625.0472896, 
1141633.00553421, 728670.952878351, 123982.148723477, 112304.540084388, 
4011.30312056738)), .Names = c("yday", "value"), row.names = c(NA, 
-157L), class = "data.frame")

Пример

ggplot(x, aes(x = yday, y = value)) + geom_point() + xlab("Day of year")

"Кривая" добавлена ​​вручную.

Что я пробовал

В предыдущих работах использовалось нормальное распределение. , распределение Вейбулла и логарифмически нормальное распространение. Я, очевидно, не могу напрямую подогнать распределение плотности к этим данным, так как данные двумерные. Вот где я борюсь. Я мог бы настроить высоту распределения плотности, но это не будет соответствовать модели, и у R должен быть лучший способ сделать это. Мое первоначальное ощущение состоит в том, что мне нужно сформулировать кривую из вышеупомянутых распределений и использовать nls, чтобы подогнать эти кривые. Я, однако, немного запутался в том, как формулировать уравнения и вызовы функций. Как только вызов nls сформулирован, его можно передать слоям stat_function или geom_smooth. Любая помощь будет оценена по достоинству.