Browse files

Draft of tobit example

  • Loading branch information...
1 parent f5ff605 commit ad1e8716607d50f0eee62411bb69183dd93e5f99 @johnmyleswhite committed Mar 23, 2012
View
14 bugs/tobit/tobit.bugs
@@ -0,0 +1,14 @@
+model
+{
+ for (i in 1:N)
+ {
+ mu[i] <- a * x[i] + b
+ y[i] ~ dnorm(mu[i], tau) T(0, )
+ }
+
+ a ~ dnorm(0, 0.0001)
+ b ~ dnorm(0, 0.0001)
+
+ tau <- pow(sigma, -2)
+ sigma ~ dunif(0, 10000)
+}
View
202 data/tobit/tobit_regression.csv
@@ -0,0 +1,202 @@
+"X","YStar","Y"
+-10,-19.6264538107423,0
+-9.9,-18.6163566757779,0
+-9.8,-19.43562861241,0
+-9.7,-16.8047191978622,0
+-9.6,-17.8704922281846,0
+-9.5,-18.820468384118,0
+-9.4,-17.3125709475715,0
+-9.3,-16.8616752948708,0
+-9.2,-16.8242186483465,0
+-9.1,-17.5053883871564,0
+-9,-15.4882188315492,0
+-8.9,-16.4101567635886,0
+-8.8,-17.2212405805418,0
+-8.7,-18.6146998871775,0
+-8.6,-15.0750690818569,0
+-8.5,-16.0449336090152,0
+-8.4,-15.8161902630989,0
+-8.3,-14.6561637893147,0
+-8.2,-14.5787788049019,0
+-8.1,-14.6060986787825,0
+-8,-14.0810226283918,0
+-7.9,-14.0178636992689,0
+-7.8,-14.5254350166348,0
+-7.7,-16.3893516958634,0
+-7.6,-13.5801742521053,0
+-7.5,-14.056128739529,0
+-7.4,-13.9557955067053,0
+-7.3,-15.0707523838993,0
+-7.2,-13.8781500551086,0
+-7.1,-12.7820584398003,0
+-7,-11.641320448471,0
+-6.9,-12.902787727343,0
+-6.8,-12.2123283884406,0
+-6.7,-12.4538050405829,0
+-6.6,-13.5770595568286,0
+-6.5,-12.4149945632997,0
+-6.4,-12.1942899537103,0
+-6.3,-11.6593133967112,0
+-6.2,-10.2999746280161,0
+-6.1,-10.4368242515425,0
+-6,-11.1645235962536,0
+-5.9,-11.0533616801365,0
+-5.8,-9.90303662459526,0
+-5.7,-9.84333680132634,0
+-5.6,-10.8887556945495,0
+-5.5,-10.7074951569621,0
+-5.4,-9.43541803786317,0
+-5.3,-8.83146707548458,0
+-5.2,-9.51234621215023,0
+-5.1,-8.31889227354579,0
+-5,-8.60189411963293,0
+-4.9,-9.41202639325077,0
+-4.8,-8.25888030857558,0
+-4.7,-9.52936309608079,0
+-4.6,-6.76697629829896,0
+-4.5,-6.01960010149414,0
+-4.4,-8.16722147646651,0
+-4.3,-8.64413462631653,0
+-4.2,-6.83028037255759,0
+-4.1,-7.33505460388082,0
+-4,-4.59838223949522,0
+-3.9,-6.83924000273317,0
+-3.8,-5.91026063754922,0
+-3.7,-6.37199784121933,0
+-3.6,-6.9432732088824,0
+-3.5,-5.81120770048566,0
+-3.4,-7.60495862889104,0
+-3.3,-4.13444513843711,0
+-3.2,-5.2467466617881,0
+-3.1,-3.02738832963785,0
+-3,-4.52449047110034,0
+-2.9,-5.50994643092181,0
+-2.8,-3.98927364651094,0
+-2.7,-5.33409763164425,0
+-2.6,-5.4536334002391,0
+-2.5,-3.70855376448254,0
+-2.4,-4.24329187321843,0
+-2.3,-3.59889464836838,0
+-2.2,-3.32565867584833,0
+-2.1,-3.78952094618807,0
+-2,-3.5686687328185,0
+-1.9,-2.93517861512383,0
+-1.8,-1.42191300342679,0
+-1.7,-3.92356680042976,0
+-1.6,-1.60605381237158,0
+-1.5,-1.66704962878648,0
+-1.4,-0.736900162723638,0
+-1.3,-1.9041839236343,0
+-1.2,-1.02998119008371,0
+-1.1,-0.932901209227768,0
+-1,-1.54252003099165,0
+-0.9,0.407867805983172,0.407867805983172
+-0.799999999999999,0.560402615694954,0.560402615694954
+-0.699999999999999,0.300213649515,0.300213649515
+-0.6,1.38683345454085,1.38683345454085
+-0.5,0.558486425565304,0.558486425565304
+-0.399999999999999,-1.07659220845803,0
+-0.299999999999999,-0.173265414236884,0
+-0.199999999999999,-0.624612614898354,0
+-0.0999999999999996,0.326599363560689,0.326599363560689
+0,0.379633322775876,0.379633322775876
+0.100000000000001,1.24211587314424,1.24211587314424
+0.200000000000001,0.489078351447557,0.489078351447557
+0.300000000000001,1.75802877240408,1.75802877240408
+0.4,1.14541535608118,1.14541535608118
+0.5,3.76728726937265,3.76728726937265
+0.600000000000001,2.91670747601721,2.91670747601721
+0.700000000000001,3.31017422949523,3.31017422949523
+0.800000000000001,2.98418535782635,2.98418535782635
+0.9,4.48217608051942,4.48217608051942
+1,2.36426354605102,2.36426354605102
+1.1,2.73835526963944,2.73835526963944
+1.2,4.83228223854166,4.83228223854166
+1.3,2.94930364668963,2.94930364668963
+1.4,3.59261925639804,3.59261925639804
+1.5,3.60719207055802,3.60719207055802
+1.6,3.8800071314515,3.8800071314515
+1.7,4.12088669702344,4.12088669702344
+1.8,5.09418833126783,5.09418833126783
+1.9,4.62266951773039,4.62266951773039
+2,4.49404253788574,4.49404253788574
+2.1,6.54303882517041,6.54303882517041
+2.2,5.18542059145313,5.18542059145313
+2.3,5.42044346995661,5.42044346995661
+2.4,5.69980925878644,5.69980925878644
+2.5,6.71266630705141,6.71266630705141
+2.6,6.12643559587368,6.12643559587368
+2.7,6.36236582853295,6.36236582853295
+2.8,5.91833952124434,5.91833952124434
+2.9,6.47572972775368,6.47572972775368
+3,7.06016044043452,7.06016044043452
+3.1,6.61110551374034,6.61110551374034
+3.2,7.93149619263257,7.93149619263257
+3.3,6.08160591821321,6.08160591821321
+3.4,8.10655786078977,8.10655786078977
+3.5,6.46355017646241,6.46355017646241
+3.6,7.89902387316339,7.89902387316339
+3.7,7.871720095555,7.871720095555
+3.8,7.947905219319,7.947905219319
+3.9,8.74310322215261,8.74310322215261
+4,7.08564057431999,7.08564057431999
+4.1,10.3765833120186,10.3765833120186
+4.2,7.735027563788,7.735027563788
+4.3,9.13646959852761,9.13646959852761
+4.4,8.68407989495716,8.68407989495716
+4.5,9.24918099880655,9.24918099880655
+4.6,12.2871665456283,12.2871665456283
+4.7,10.4173956196933,10.4173956196933
+4.8,9.31369946956568,9.31369946956568
+4.9,9.15939446558142,9.15939446558142
+5,11.4501871012727,11.4501871012727
+5.1,11.1814401672854,11.1814401672854
+5.2,11.0819316254562,11.0819316254562
+5.3,10.6706378525463,10.6706378525463
+5.4,10.3125396898585,10.3125396898585
+5.5,10.9248077033843,10.9248077033843
+5.6,13.2000288037139,13.2000288037139
+5.7,11.7787333052032,11.7787333052032
+5.8,11.2155731526155,11.2155731526155
+5.9,14.6692906224236,14.6692906224236
+6,13.4251003773724,13.4251003773724
+6.1,12.961352899087,12.961352899087
+6.2,14.458483048709,14.458483048709
+6.3,14.4864226513749,14.4864226513749
+6.4,13.1807569517689,13.1807569517689
+6.5,16.2061024645405,16.2061024645405
+6.6,13.944972969859,13.944972969859
+6.7,12.9755053497872,12.9755053497872
+6.8,14.4556003980458,14.4556003980458
+6.9,15.0075383392323,15.0075383392323
+7,17.3079783990594,17.3079783990594
+7.1,15.3058023678937,15.3058023678937
+7.2,15.8569988054234,15.8569988054234
+7.3,15.5228470646435,15.5228470646435
+7.4,15.4659991576335,15.4659991576335
+7.5,15.9652739716887,15.9652739716887
+7.6,16.9876396056302,16.9876396056302
+7.7,18.4752450086523,18.4752450086523
+7.8,17.6273924387638,17.6273924387638
+7.9,18.0079083983867,18.0079083983867
+8,15.768676578442,15.768676578442
+8.1,18.1838955700534,18.1838955700534
+8.2,17.6199248036606,17.6199248036606
+8.3,16.1327499709078,16.1327499709078
+8.4,18.3210227426481,18.3210227426481
+8.5,17.841245395284,17.841245395284
+8.6,19.6645873119698,19.6645873119698
+8.7,17.6339180003953,17.6339180003953
+8.8,18.1697882460715,18.1697882460715
+8.9,17.8738905026226,17.8738905026226
+9,18.8228960385635,18.8228960385635
+9.1,19.6020117794863,19.6020117794863
+9.2,18.6682518268804,18.6682518268804
+9.3,20.4303731679817,20.4303731679817
+9.4,18.5919172136955,18.5919172136955
+9.5,18.9520155871923,18.9520155871923
+9.6,21.6411577068443,21.6411577068443
+9.7,19.3841525346954,19.3841525346954
+9.8,21.0119747123175,21.0119747123175
+9.9,20.4189239488911,20.4189239488911
+10,21.4094018396509,21.4094018396509
View
13 generators/tobit/tobit.R
@@ -0,0 +1,13 @@
+set.seed(1)
+
+a <- 2
+b <- 1
+
+x <- seq(-10, 10, by = 0.1)
+y.star <- a * x + b
+y.star <- y.star + rnorm(length(y.star), 0, 1)
+y <- ifelse(y.star >= 0, y.star, 0)
+
+write.csv(data.frame(X = x, YStar = y.star, Y = y),
+ file = file.path('data', 'tobit', 'tobit_regression.csv'),
+ row.names = FALSE)
View
BIN graphs/tobit/data_plot.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
BIN graphs/tobit/hidden_data_plot.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
40 scripts/tobit/tobit.R
@@ -0,0 +1,40 @@
+# Basic inference.
+library('rjags')
+
+df <- read.csv(file.path('data',
+ 'tobit',
+ 'tobit_regression.csv'))
+
+png(file.path('graphs', 'tobit', 'data_plot.png'))
+with(df, plot(X, Y))
+dev.off()
+
+png(file.path('graphs', 'tobit', 'hidden_data_plot.png'))
+with(df, plot(X, YStar))
+dev.off()
+
+jags <- jags.model(file.path('bugs',
+ 'tobit',
+ 'tobit.bugs'),
+ data = list('x' = with(df, X),
+ 'y' = with(df, Y),
+ 'N' = nrow(df)),
+ n.chains = 4,
+ n.adapt = 1000)
+
+mcmc.samples <- coda.samples(jags,
+ c('a', 'b', 'sigma'),
+ 1000)
+
+png(file.path('graphs', 'tobit', 'plot1.png'))
+plot(mcmc.samples)
+dev.off()
+
+summary(mcmc.samples)
+
+# Compare with MLE OLS solution.
+lm(Y ~ X, data = df)
+
+# Compare with MLE tobit solution.
+library('VGAM')
+fit <- vglm(Y ~ X, family = tobit(Lower = 0), df, trace = TRUE, crit = 'coeff')

0 comments on commit ad1e871

Please sign in to comment.