Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
607 views
in Technique[技术] by (71.8m points)

python - Using R to fit data from a csv with a gamma function?

Using the following data:

time_stamp,secs,count
2013-04-30 23:58:55,1367366335,32
2013-04-30 23:58:56,1367366336,281
2013-04-30 23:58:57,1367366337,664
2013-04-30 23:58:58,1367366338,1255
2013-04-30 23:58:59,1367366339,2011
2013-04-30 23:59:00,1367366340,2819
2013-04-30 23:59:01,1367366341,3489
2013-04-30 23:59:02,1367366342,4230
2013-04-30 23:59:03,1367366343,4892
2013-04-30 23:59:04,1367366344,5185
2013-04-30 23:59:05,1367366345,5415
2013-04-30 23:59:06,1367366346,5428
2013-04-30 23:59:07,1367366347,5255
2013-04-30 23:59:08,1367366348,4897
2013-04-30 23:59:09,1367366349,4548
2013-04-30 23:59:10,1367366350,4154
2013-04-30 23:59:11,1367366351,3718
2013-04-30 23:59:12,1367366352,3521
2013-04-30 23:59:13,1367366353,3193
2013-04-30 23:59:14,1367366354,3015
2013-04-30 23:59:15,1367366355,2722
2013-04-30 23:59:16,1367366356,2545
2013-04-30 23:59:17,1367366357,2368
2013-04-30 23:59:18,1367366358,2265
2013-04-30 23:59:19,1367366359,2162
2013-04-30 23:59:20,1367366360,1967
2013-04-30 23:59:21,1367366361,1952
2013-04-30 23:59:22,1367366362,1750
2013-04-30 23:59:23,1367366363,1720
2013-04-30 23:59:24,1367366364,1624
2013-04-30 23:59:25,1367366365,1497
2013-04-30 23:59:26,1367366366,1452
2013-04-30 23:59:27,1367366367,1486
2013-04-30 23:59:28,1367366368,1328
2013-04-30 23:59:29,1367366369,1284
2013-04-30 23:59:30,1367366370,1266
2013-04-30 23:59:31,1367366371,1220
2013-04-30 23:59:32,1367366372,1152
2013-04-30 23:59:33,1367366373,1217
2013-04-30 23:59:34,1367366374,1144
2013-04-30 23:59:35,1367366375,1113
2013-04-30 23:59:36,1367366376,1051
2013-04-30 23:59:37,1367366377,1036
2013-04-30 23:59:38,1367366378,976
2013-04-30 23:59:39,1367366379,1018
2013-04-30 23:59:40,1367366380,974
2013-04-30 23:59:41,1367366381,946
2013-04-30 23:59:42,1367366382,910
2013-04-30 23:59:43,1367366383,912
2013-04-30 23:59:44,1367366384,853
2013-04-30 23:59:45,1367366385,824
2013-04-30 23:59:46,1367366386,908
2013-04-30 23:59:47,1367366387,800
2013-04-30 23:59:48,1367366388,841
2013-04-30 23:59:49,1367366389,828
2013-04-30 23:59:50,1367366390,775
2013-04-30 23:59:51,1367366391,761
2013-04-30 23:59:52,1367366392,827
2013-04-30 23:59:53,1367366393,725
2013-04-30 23:59:54,1367366394,752
2013-04-30 23:59:55,1367366395,730
2013-04-30 23:59:56,1367366396,701
2013-04-30 23:59:57,1367366397,685
2013-04-30 23:59:58,1367366398,680
2013-04-30 23:59:59,1367366399,704
2013-05-01 00:00:00,1367366400,642
2013-05-01 00:00:01,1367366401,688
2013-05-01 00:00:02,1367366402,682
2013-05-01 00:00:03,1367366403,660
2013-05-01 00:00:04,1367366404,602
2013-05-01 00:00:05,1367366405,672
2013-05-01 00:00:06,1367366406,663
2013-05-01 00:00:07,1367366407,649
2013-05-01 00:00:08,1367366408,609
2013-05-01 00:00:09,1367366409,594
2013-05-01 00:00:10,1367366410,554
2013-05-01 00:00:11,1367366411,554
2013-05-01 00:00:12,1367366412,594
2013-05-01 00:00:13,1367366413,552
2013-05-01 00:00:14,1367366414,565
2013-05-01 00:00:15,1367366415,553
2013-05-01 00:00:16,1367366416,550
2013-05-01 00:00:17,1367366417,513
2013-05-01 00:00:18,1367366418,525
2013-05-01 00:00:19,1367366419,518
2013-05-01 00:00:20,1367366420,485
2013-05-01 00:00:21,1367366421,554
2013-05-01 00:00:22,1367366422,527
2013-05-01 00:00:23,1367366423,490
2013-05-01 00:00:24,1367366424,509
2013-05-01 00:00:25,1367366425,460
2013-05-01 00:00:26,1367366426,473
2013-05-01 00:00:27,1367366427,491
2013-05-01 00:00:28,1367366428,511
2013-05-01 00:00:29,1367366429,490
2013-05-01 00:00:30,1367366430,438
2013-05-01 00:00:31,1367366431,487
2013-05-01 00:00:32,1367366432,449
2013-05-01 00:00:33,1367366433,426
2013-05-01 00:00:34,1367366434,482
2013-05-01 00:00:35,1367366435,484
2013-05-01 00:00:36,1367366436,433
2013-05-01 00:00:37,1367366437,455
2013-05-01 00:00:38,1367366438,441
2013-05-01 00:00:39,1367366439,450
2013-05-01 00:00:40,1367366440,440
2013-05-01 00:00:41,1367366441,446
2013-05-01 00:00:42,1367366442,470
2013-05-01 00:00:43,1367366443,421
2013-05-01 00:00:44,1367366444,456
2013-05-01 00:00:45,1367366445,442
2013-05-01 00:00:46,1367366446,424
2013-05-01 00:00:47,1367366447,416
2013-05-01 00:00:48,1367366448,455
2013-05-01 00:00:49,1367366449,366
2013-05-01 00:00:50,1367366450,407
2013-05-01 00:00:51,1367366451,413
2013-05-01 00:00:52,1367366452,420
2013-05-01 00:00:53,1367366453,396
2013-05-01 00:00:54,1367366454,428
2013-05-01 00:00:55,1367366455,344
2013-05-01 00:00:56,1367366456,372
2013-05-01 00:00:57,1367366457,388
2013-05-01 00:00:58,1367366458,394
2013-05-01 00:00:59,1367366459,383
2013-05-01 00:01:00,1367366460,391
2013-05-01 00:01:01,1367366461,385
2013-05-01 00:01:02,1367366462,410
2013-05-01 00:01:03,1367366463,359
2013-05-01 00:01:04,1367366464,328
2013-05-01 00:01:05,1367366465,396
2013-05-01 00:01:06,1367366466,388
2013-05-01 00:01:07,1367366467,402
2013-05-01 00:01:08,1367366468,349
2013-05-01 00:01:09,1367366469,355
2013-05-01 00:01:10,1367366470,315
2013-05-01 00:01:11,1367366471,360
2013-05-01 00:01:12,1367366472,357
2013-05-01 00:01:13,1367366473,369
2013-05-01 00:01:14,1367366474,337
2013-05-01 00:01:15,1367366475,358
2013-05-01 00:01:16,1367366476,367
2013-05-01 00:01:17,1367366477,357
2013-05-01 00:01:18,1367366478,362
2013-05-01 00:01:19,1367366479,302
2013-05-01 00:01:20,1367366480,346
2013-05-01 00:01:21,1367366481,338
2013-05-01 00:01:22,1367366482,366
2013-05-01 00:01:23,1367366483,341
2013-05-01 00:01:24,1367366484,332
2013-05-01 00:01:25,1367366485,317
2013-05-01 00:01:26,1367366486,310
2013-05-01 00:01:27,1367366487,327
2013-05-01 00:01:28,1367366488,333
2013-05-01 00:01:29,1367366489,322
2013-05-01 00:01:30,1367366490,344
2013-05-01 00:01:31,1367366491,330
2013-05-01 00:01:32,1367366492,312
2013-05-01 00:01:33,1367366493,292
2013-05-01 00:01:34,1367366494,315
2013-05-01 00:01:35,1367366495,326
2013-05-01 00:01:36,1367366496,313
2013-05-01 00:01:37,1367366497,307
2013-05-01 00:01:38,1367366498,325
2013-05-01 00:01:39,1367366499,308
2013-05-01 00:01:40,1367366500,299
2013-05-01 00:01:41,1367366501,293
2013-05-01 00:01:42,1367366502,285
2013-05-01 00:01:43,1367366503,295
2013-05-01 00:01:44,1367366504,310
2013-05-01 00:01:45,1367366505,309
2013-05-01 00:01:46,1367366506,289
2013-05-01 00:01:47,1367366507,307
2013-05-01 00:01:48,1367366508,311
2013-05-01 00:01:49,1367366509,294
2013-05-01 00:01:50,1367366510,279
2013-05-01 00:01:51,1367366511,280
2013-05-01 00:01:52,1367366512,247
2013-05-01 00:01:53,1367366513,285
2013-05-01 00:01:54,1367366514,288
2013-05-01 00:01:55,1367366515,250
2013-05-01 00:01:56,1367366516,284
2013-05-01 00:01:57,1367366517,281
2013-05-01 00:01:58,1367366518,250
2013-05-01 00:01:59,1367366519,266
2013-05-01 00:02:00,1367366520,305
2013-05-01 00:02:01,1367366521,257
2013-05-01 00:02:02,1367366522,252
2013-05-01 00:02:03,1367366523,261
2013-05-01 00:02:04,1367366524,253
2013-05-01 00:02:05,1367366525,273
2013-05-01 00:02:06,1367366526,290
2013-05-01 00:02:07,1367366527,267
2013-05-01 00:02:08,1367366528,289
2013-05-01 00:02:09,1367366529,259
2013-05-01 00:02:10,1367366530,281
2013-05-01 00:02:11,1367366531,242
2013-05-01 00:02:12,1367366532,280
2013-05-01 00:02:13,1367366533,300
2013-05-01 00:02:14,1367366534,284
2013-05-01 00:02:15,1367366535,244
2013-05-01 00:02:16,1367366536,236
2013-05-01 00:02:17,1367366537,269
2013-05-01 00:02:18,1367366538,292
2013-05-01 00:02:19,1367366539,279
2013-05-01 00:02:20,1367366540,234
2013-05-01 00:02:21,1367366541,262
2013-05-01 00:02:22,1367366542,245
2013-05-01 00:02:23,1367366543,258
2013-05-01 00:02:24,1367366544,230
2013-05-01 00:02:25,1367366545,252
2013-05-01 00:02:26,1367366546,273
2013-05-01 00:02:27,1367366547,216
2013-05-01 00:02:28,1367366548,245
2013-05-01 00:02:29,1367366549,254
2013-05-01 00:02:30,1367366550,293
2013-05-01 00:02:31,1367366551,331
2013-05-01 00:02:32,1367366552,320
2013-05-01 00:02:33,1367366553,303
2013-05-01 00:02:34,1367366554,378
2013-05-01 00:02:35,1367366555,348
2013-05-01 00:02:36,1367366556,347
2013-05-01 00:02:37,1367366557,316
2013-05-01 00:02:38,1367366558,333
2013-05-01 00:02:39,1367366559,299
2013-05-01 00:02:40,1367366560,284
2013-05-01 00:02:41,1367366561,279
2013-05-01 00:02:42,1367366562,261
2013-05-01 00:02:43,1367366563,259
2013-05-01 00:02:44,1367366564,223
2013-05-01 00:02:45,1367366565,229
2013-05-01 00:02:46,1367366566,243
2013-05-01 00:02:47,1367366567,217
2013-05-01 00:02:48,1367366568,225
2013-05-01 00:02:49,1367366569,215
2013-05-01 00:02:50,1367366570,265
2013-05-01 00:02:51,1367366571,240
2013-05-01 00:02:52,1367366572,248
2013-05-01 00:02:53,1367366573,226
2013-05-01 00:02:54,1367366574,188
2013-05-01 00:02:55,1367366575,213
2013-05-01 00:02:56,1367366576,216
2013-05-01 00:02:57,1367366577,222
2013-05-01 00:02:58,1367366578,282
2013-05-01 00:02:59,1367366579,234
2013-05-01 00:03:00,1367366580,216
2013-05-01 00:03:01,1367366581,195
2013-05-01 00:03:02,1367366582,224
2013-05-01 00:03:03,1367366583,227
2013-05-01 00:03:04,1367366584,256
2013-05-01 00:03:05,1367366585,213
2013-05-01 00:03:06,1367366586,225
2013-05-01 00:03:07,1367366587,229
2013-05-01 00:03:08,1367366588,217
2013-05-01 00:03:09,1367366589,220
2013-05-01 00:03:10,1367366590,209
2013-05-01 00:03:11,1367366591,231
2013-05-01 00:03:12,1367366592,195
2013-05-01 00:03:13,1367366593,231
2013-05-01 00:03:14,1367366594,185
2013-05-01 00:03:15,1367366595,173
2013-05-01 00:03:16,1367366596,190
2013-05-01 00:03:17,1367366597,196
2013-05-01 00:03:18,1367366598,190
2013-05-01 00:03:19,1367366599,209
2013-05-01 00:03:20,1367366600,206
2013-05-01 00:03:21,1367366601,227
2013-05-01 00:03:22,1367366602,199
2013-05-01 00:03:23,1367366603,182
2013-05-01 00:03:24,1367366604,202
2013-05-01 00:03:25,1367366605,208
2013-05-01 00:03:26,1367366606,190
2013-05-01 00:03:27,1367366607,220
2013-05-01 00:03:28,1367366608,190
2013-05-01 00:03:29,1367366609,205
2013-05-01 00:03:30,1367366610,196
2013-05-01 00:03:31,1367366611,182
2013-05-01 00:03:32,1367366612,187
2013-05-01 00:03:33,1367366613,186
2013-05-01 00:03:34,1367366614,195
2013-05-01 00:03:35,1367366615,224
2013-05-01 00:03:36,1367366616,210
2013-05-01 00:03:37,1367366617,188
2013-05-01 00:03:38,1367366618,167
2013-05-01 00:03:39,1367366619,194
2013-05-01 00:03:40,1367366620,187
2013-05-01 00:03:41,1367366621,213
2013-05-01 00:03:42,1367366622,193
2013-05-01 00:03:43,1367366623,187
2013-05-01 00:03:44,1367366624,209
2013-05-01 00:03:45,1367366625,204
2013-05-01 00:03:46,1367366626,199
2013-05-01 00:03:47,1367366627,186
2013-05-01 00:03:48,1367366628,197
2013-05-01 00:03:49,1367366629,196
2013-05-01 00:03:50,1367366630,182
2013-05-01 00:03:51,1367366631,169
2013-05-01 00:03:52,1367366632,196
2013-05-01 00:03:53,1367366633,198
2013-05-01 00:03:54,1367366634,187
2013-05-01 00:03:55,1367366635,178
2013-05-01 00:03:56,1367366636,197
2013-05-01 00:03:57,1367366637,172
2013-05-01 00:03:58,1367366638,186
2013-05-01 00:03:59,1367366639,181
2013-05-01 00:04:00,1367366640,167
2013-05-01 00:04:01,1367366641,201
2013-05-01 00:04:02,1367366642,172
2013-05-01 00:04:03,1367366643,190
2013-05-01 00:04:04,1367366644,167
2013-05-01 00:04:05,1367366645,152
2013-05-01 00:04:06,1367366646,186
2013-05-01 00:04:07,1367366647,193
2013-05-01 00:04:08,1367366648,207
2013-05-01 00:04:09,1367366649,182
2013-05-01 00:04:10,1367366650,182
2013-05-01 00:04:11,1367366651,183
2013-05-01 00:04:12,1367366652,204
2013-05-01 00:04:13,1367366653,163
2013-05-01 00:04:14,1367366654,164
2013-05-01 00:04:15,1367366655,156
2013-05-01 00:04:16,1367366656,177
2013-05-01 00:04:17,1367366657,162
2013-05-01 00:04:18,1367366658,168
2013-05-01 00:04:19,1367366659,178
2013-05-01 00:04:20,1367366660,201
2013-05-01 00:04:21,1367366661,165
2013-05-01 00:04:22,1367366662,181
2013-05-01 00:04:23,1367366663,152
2013-05-01 00:04:24,13

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)

-1 to the OP for not showing work, but +1 for posting an interesting problem. Here's a stab at it. The idea is to fit a curve to the (x, y) data, where the curve is a (scaled) gamma density. Since the data is a table of counts over a grid, a poisson error distribution seems to be a reasonable assumption. I've implicitly rescaled the problem so that the grid is 1...365 rather than 1367366395...1367366699 to make life easier; translating the results back to the original scale is an exercise for the reader.

negll <- function(par, x, y)
{
    shape <- par[1]
    rate <- par[2]
    mu <- dgamma(x, shape, rate) * sum(y)
    -2 * sum(dpois(y, mu, log=TRUE))
}


optim(c(1, 1), negll, x=seq_along(g$count), y=g$count, method="L-BFGS-B", lower=c(.001, .001))
$par
[1] 0.73034879 0.00698288

$value
[1] 62983.18

$counts
function gradient 
      32       32 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

This is actually a pretty poor fit (shape less than 1 indicates the fitted curve goes to infinity at the origin, in contrast to the observed data). Looking at the observed trend suggests that a gamma density is a poor choice for this data; the peak is too sharp and decays too slowly. A more heavy-tailed distribution like the (generalised) Pareto might be a better bet.

Couple more things:

  1. This is fitted to the entire table rather than just the first 30 points. Using only 30 points would make it even more difficult, since you get minimal information on the shape and rate from that part of the curve alone.

  2. There's a noticeable bump in the curve 2/3rds of the way down. Whether this is materially significant is for the OP to decide, but I hope they actually noticed this rather than skipping directly to the curve-fitting stage.


与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...